最近一直在學HTML5和CSS3，Numpy的東西都有些生疏，那本書是已經看完了的，緊跟着相關的代碼也都敲了一遍，還是發現了一些問題，因為這樣的學習方式，總感覺太被動，緊緊跟着示例代碼，缺少了整體觀，即使你現在問我Numpy可以處理什么問題，我還是回答不出。所以，有必要回頭重來一遍，再一次審視代碼背后的意義，寫博客真的是一個很不錯的方式，畢竟，如果你不懂，寫出來的文字必然也是混亂的。
那，下面記錄一下Numpy學習筆記（二）
Example1
文件讀寫：數據不應該僅僅存在內存里，應該及時保存在硬盤上，以文件的形式
# -*- coding:utf-8 -*- # 導入numpy庫 import numpy as np
# 文件讀寫 # 先創建一個單位矩陣 a = np.eye(2) # numpy自帶的eye方法 print a # 使用savetxt函數保存 np.savetxt("D:\Learn\Code\python\exercise\eye.txt", a) |
結果如下：
查看相關文件夾，已經出現eye.txt，打開
Example2
CSV文件：逗號分隔值文件，很少聽見，但經常遇見，用Excel打開的效果幾乎與.xls文件一模一樣
Numpy是用來處理數據的，而CSV是用來存儲數據的，看起來淵源很深呢。
先來看一個CSV文件
里面記錄的是蘋果公司的股票，第一列是股票代碼，第二列是日期，第三列為空，下面依次是開盤價、最高價、最低價和收盤價和成交量
接下來我們將收盤價和成交量分別載入數組
# CSV文件 # delimiter定義分隔符，默認是空格；usecols定義選取的；unpack默認False，解包 c, v = np.loadtxt('D:\Learn\Code\python\exercise\data.csv', \ delimiter=",", usecols=(6, 7), unpack=True) print u"收盤價:", c print u"成交量:", v |
看結果：
並不困難
Example3
各種平均值的求法：書里面列舉了三個，成交量加權平均價格、算術平均值、時間加權平均價格
成交量平均價格（VWAP）表示成績量越大，該價格所占的權重就越大
算術平均值就是權重為1
時間加權平均價格（TWAP）表示日期越近的價格所占比重越大
# 成交量平均價格 c, v = np.loadtxt("D:\Learn\Code\python\exercise\data.csv",\ delimiter=",", usecols=(6, 7), unpack=True) # 計算以成交量加權的價格 vwap = np.average(c, weights=v) print "VWAP:", vwap
# 算術平均價格 mean = np.mean(c) print "Mean:", mean
# 時間加權平均價格 # 先構造一個時間序列 t = np.arange(len(c)) twap = np.average(c, weights=t) print "TWAP:", twap |
結果如下：
Example4
最大值與最小值
# 載入每日最高價和最低價的數據 h, l = np.loadtxt("D:\Learn\Code\python\exercise\data.csv",\ delimiter=",", usecols=(4, 5), unpack=True) # 直接調用max()和min()函數 print "highest:", np.max(h) print "lowest:", np.min(l)
# ptp函數可以計算取值范圍 print "Spred high price:", np.ptp(h) print "Spred low price:", np.ptp(l) |
結果如下：
無論怎么講都是相當便捷的
Example5
簡單統計分析，這里主要是中位數、和方差
# 簡單統計分析 c = np.loadtxt("D:\Learn\Code\python\exercise\data.csv",\ delimiter=",", usecols=(6, ), unpack=True) # 尋找中位數 print "median:", np.median(c) # 利用排序來檢驗正確與否 sorted_c = np.msort(c) print "sorted_c:", sorted_c N = len(c) if N%2 != 0: middle = sorted_c[N/2] else: middle = (sorted_c[N/2] + sorted_c[N/2-1])/2 print "middle:", middle
# 計算方差 variance = np.var(c) print "variance:", variance |
結果如下：
所需要的無非就只是調用個函數那么簡單！
Example6
計算股票收益率，分成這么幾個（都是經濟學的東西，看看就好，重點還是看numpy如何使用上）：簡單收益率、對數收益率和波動率
# 股票收益率 c = np.loadtxt("D:\Learn\Code\python\exercise\data.csv",\ delimiter=",", usecols=(6, ), unpack=True) # 簡單收益率 # diff函數返回相鄰數組元素差值組成的數組 diff = np.diff(c) print "diff:", diff # 收益率等於這一天與前一天的差值除以前一天的值 returns = diff/c[:-1] print "returns:", returns # 計算收益率的標注差 print "Standard deviation:", np.std(returns)
# 對數收益率 logreturns = np.diff(np.log(c)) print "logreturns:", logreturns print "Standard deviation:", np.std(logreturns)
# 計算收益率為正值的情況 posretindices = np.where(returns > 0) print "Indices with positive returns:", posretindices
# 波動率（對數收益率的標准差除以其均值，再除以交易日倒數的平方根，交易日通常取252天） annual_volatility = np.std(logreturns)/np.mean(logreturns)/np.sqrt(1.0/252) print "Annual volatility:", annual_volatility print "Month volatility:", annual_volatility * np.sqrt(1.0/12) |
結果如下：
Example7
日期分析：分析CSV文件，會看到有一列日期數據，格式為dd-mm-yyyy，Python自帶的時間處理模塊，可以很自如地處理，而Numpy本身是無法處理這樣的數據的，所以在解包的時候要先編寫一個轉換函數，然后將這個函數傳遞到loadtxt里的converters參數里
# 時間處理 # 先導入時間模塊 import datetime # 編寫轉換函數 def datetostr(s): truetime = datetime.datetime.strptime(s, '%d-%m-%Y').weekday() return truetime # 開始讀取數據，並添加converters參數 dates, close = np.loadtxt("D:\Learn\Code\python\exercise\data.csv",\ delimiter=",", usecols=(1, 6), converters={1: datetostr}, \ unpack=True) print "Dates:", dates # 創建一個包含五個元素的數組 averages = np.zeros(5) # where會根據條件返回滿足條件的元素索引，take可以從索引中取出數據 for i in range(5): indices = np.where(dates == i) prices = np.take(close, indices) avg = np.mean(prices) print "Day", i, "prices", prices, "Average", avg averages[i] = avg |
結果是這樣的：
這個練習有些困難，不僅使用到了Python里的datetime模塊，還涉及到了loadtxt轉換參數converters的使用，更不必說where和take了
Example8
Data.csv文件里，每一行代表每一天，如果數據量很大，我們可以考慮將它們進行壓縮，按周進行匯總
# 周匯總 import datetime def datetostr(s): truetime = datetime.datetime.strptime(s, '%d-%m-%Y').weekday() return truetime dates, open, high, low, close = np.loadtxt("D:\Learn\Code\python\exercise\data.csv",\ delimiter=",", usecols=(1, 3, 4, 5, 6), converters={1: datetostr}, \ unpack=True) # 為了方便計算，僅取前三周數據 dates = dates[:16]
# 尋找第一個星期一 first_monday = np.ravel(np.where(dates == 0))[0] # where返回的是個多維數組，需要展平 print "First Monday", first_monday # 尋找最后一個星期五 last_friday = np.ravel(np.where(dates == 4))[-1] print "Last Friday"
weeks_indices = np.arange(first_monday, last_friday+1) print "Weeks indices initial", weeks_indices
weeks_indices = np.split(weeks_indices, 3) print "Weeks indices after split", weeks_indices
# 為了后面的apply_along_axis def summarize(a, o, h, l, c): monday_open = o[a[0]] week_high = np.max(np.take(h, a)) week_low = np.min(np.take(l, a)) friday_close = c[a[-1]] return ("APPL", monday_open, week_high, week_low, friday_close) # apply_along_axis內涵很豐富 weeksummary = np.apply_along_axis(summarize, 1, weeks_indices, open, high, low, close) print "Week summary", weeksummary # savetxt參數其實有很多 np.savetxt("D:\Learn\Code\python\exercise\weeksummary.csv", weeksummary, delimiter=",", fmt="%s") |
看結果：
打開文件夾，查看weekssummay.csv文件
Example9
真實波動幅度均值（ATR），這里取20天的數據
# 真實波動幅度均值 h, l, c = np.loadtxt('D:\Learn\Code\python\exercise\data.csv', delimiter=',', usecols=(4, 5, 6), unpack=True)
N = 20 # 切片 h = h[-N:] l = l[-N:]
print "len(h)", len(h), "len(l)", len(l) print "Close", c # 計算前一日的收盤價 previousclose = c[-N-1:-1]
print "len(previousclose)", len(previousclose) print "Previous close", previousclose # maximum函數可以選擇出每個元素的最大值 truerange = np.maximum(h-l, h-previousclose, previousclose-l)
print "True range", truerange # zeros函數初始化數組為0 atr = np.zeros(N)
atr[0] = np.mean(truerange)
for i in range(1, N): atr[i] = (N-1)*atr[i-1] + truerange[i] atr[i] /= N
print "ATR", atr |
結果是這樣的
Example10
簡單移動平均線
這里會涉及到一個重要的函數：consolve函數，即卷積函數。卷積的概念百度百科上是這樣解釋的：卷積是兩個變量在某范圍內相乘后求和的結果。我在知乎上看到的一個答案，說的更簡潔，即加權求和。
# 簡單移動平均線 # N是移動窗口的大小 N = 5 # 權重是個平均值 weights = np.ones(N) / N print "Weights", weights
c = np.loadtxt('D:\Learn\Code\python\exercise\data.csv', delimiter=',', usecols=(6, ), unpack=True) # 要從卷積運算中取出與原數組重疊的區域 sma = np.convolve(weights, c)[N-1:-N+1] # 生成一個時間序列 t = np.arange(N-1, len(c)) # 用matplotlib繪圖 plt.plot(t, c[N-1:],'b-', lw=1.0) plt.plot(t, sma, 'r-', lw=2.0) plt.show() |
結果是這樣的
Matplotlib的繪圖效果相當不錯，這里還可以添加一些參數
plt.xlabel("Date") plt.ylabel("Close Price") plt.title(u"Simple Moving Average") plt.annotate('before convolve', xy=(12.8, 363), xytext=(15, 363), arrowprops=dict(facecolor='black',shrink=0.005)) plt.annotate('after convolve', xy=(15, 358), xytext=(17, 358), arrowprops=dict(facecolor='black',shrink=0.005)) |
Example11
指數移動平均線
指數移動平均線使用的權重是指數衰減的，其他的與Example10一樣
# 指數移動平均線 x = np.arange(5) # 對x求指數，exp函數 print "Exp", np.exp(x) # linspace函數實現等距分隔 print "Linspace", np.linspace(-1, 0 ,5) N = 5 # 上面是兩個示范，下面才是真的 weights = np.exp(np.linspace(-1, 0, N)) weights /= weights.sum() print "Weights", weights
c = np.loadtxt('D:\Learn\Code\python\exercise\data.csv', delimiter=',', usecols=(6, ), unpack=True) ema = np.convolve(weights, c)[N-1:-N+1] t = np.arange(N-1, len(c)) plt.plot(t, c[N-1:],'b-', lw=1.0) plt.plot(t, ema, 'r-', lw=2.0) plt.show() |
結果如下
與上一幅圖略有不同
Example12
布林帶
股票市場的一種常用指標，基本形態是由三條軌道線組成
中軌：簡單移動平均線
上軌：比簡單移動平均線高兩倍標准差的距離，標准差是簡單移動平均線的標准差
下軌：比簡單移動平均線低兩倍標准差的距離，標准差是簡單移動平均線的標准差
# 繪制布林帶 N = 5 # 計算權重 weights = np.ones(N)/N
c = np.loadtxt('D:\Learn\Code\python\exercise\data.csv', delimiter=',', usecols=(6, ), unpack=True) # 簡單移動平均線，注意切片 sma = np.convolve(weights, c)[N-1:-N+1] deviation = []
# 標准差為計算簡單移動平均線所用數據的標准差 for i in range(0, len(c)-N+1): dev = c[i:i+N]
averages = np.zeros(N) # fill函數可以將數組元素賦為單一值，平均值恰好為sma數組里的元素 averages.fill(sma[i]) dev = dev-averages dev = dev ** 2 dev = np.sqrt(np.mean(dev)) deviation.append(dev)
# 書上的代碼 ''' for i in range(N-1, len(c)): if i+N<len(c): dev = c[i:i+N] else: dev = c[-N:]
averages = np.zeros(N) averages.fill(sma[i-N-1]) dev = dev-averages dev = dev ** 2 dev = np.sqrt(np.mean(dev)) deviation.append(dev) '''
deviation = 2 * np.array(deviation) # 每個sma的元素應對應一個標注差 print len(deviation), len(sma) upperBB = sma + deviation lowerBB = sma - deviation
c_slice = c[N-1:] # 檢驗數據是否全都落入上軌和下軌內 between_bands = np.where((c_slice<upperBB)&(c_slice>lowerBB))
print lowerBB[between_bands] print c[between_bands] print upperBB[between_bands] between_bands = len(np.ravel(between_bands)) print "Ratio between bands", float(between_bands)/len(c_slice)
# 繪圖，這個就比較簡單了 t = np.arange(N-1, len(c)) plt.plot(t, c_slice, lw=1.0) plt.plot(t, sma, lw=2.0) plt.plot(t, upperBB, lw=3.0) plt.plot(t, lowerBB, lw=3.0) plt.show() |
結果如下：
書上的代碼略有不同，但我並沒有看懂，而是照着自己的理解，根據布林帶的規則寫的，如果大家發現有什么問題，或者我的寫法是錯誤的，希望能及時的提醒我，也好及時更改(*^_^*)
Example13
線性模型，Numpy里的linalg包是專門用來處理此類問題，今后還會接觸
# 用線性模型預測價格
# 用於預測所取的樣本量 N = 5
c = np.loadtxt("D:\Learn\Code\python\exercise\data.csv", delimiter=',', usecols = (6, ), unpack=True) # 取后N個數 b = c[-N:] # 倒序 b = b[::-1] print "b", b
# 初始化一個N*N的二維數組A A = np.zeros((N,N), float) print "Zeros N by N", A
# A[i]與b對應 for i in range(N): A[i, ] = c[-N-1-i:-1-i]
print "A", A
# lstsq函數擬合數據，返回值包括系數向量、殘差數組、A的秩以及A的奇異值 (x, residuals, rank, s) = np.linalg.lstsq(A, b)
print x, residuals, rank, s
# x提供系數，dot點積即可預測下一次股價 print np.dot(b, x) |
結果如下：
實際查得下一個交易日收盤價為353.56，說明股票市場還是很不規律的
Example14
數組的修剪和壓縮
# 數組的修剪和壓縮
# clip返回一個修剪過的數組，小於等於給定最小值的設為給定最小值，反之亦然 a = np.arange(5) print "a=", a print "Clipped", a.clip(1, 2)
# compress返回一個給定篩選條件后的數組 print "Compressed", a.compress(a>2) |
結果也是比較明顯：
Example15
階乘
# 階乘 b = np.arange(1, 9) print "b=", b # 一個prod()函數即可，省略了循環 print "Factorial", b.prod()
# cumprod函數可以計算累計乘積 print "Factorials", b.cumprod() |
看結果：
總結：從這次的筆記中可以看到，Numpy的確是很大程度上拓展了Python的統計功能，很多時候你所需要做的無非只是寫個函數名，直接調用即可，Matplotlib同時提供了個性化很強的繪圖功能，兩者集成就可以實現很好的數據可視化。當然，這里的幾個例子無非只是拋磚引玉，接下來的學習之路還很漫長。
源代碼在這里：https://github.com/Lucifer25/Learn-Python/blob/master/Numpy/exercise2.py
參考資料：http://docs.scipy.org/doc/numpy/reference/
http://matplotlib.org/