Numpy能夠讀寫磁盤上的文本數據或二進制數據。這一小節只討論Numpy的內置二進制格式,因為更多的用戶會使用pandas或其它工具加載文本或表格數據(見第6章)。
np.save和np.load是讀寫磁盤數組數據的兩個主要函數。默認情況下,數組是以未壓縮的原始二進制格式保存在擴展名為.npy的文件中的:
In [213]: arr = np.arange(10)
In [214]: np.save('some_array', arr)
如果文件路徑末尾沒有擴展名.npy,則該擴展名會被自動加上。然后就可以通過np.load讀取磁盤上的數組:
In [215]: np.load('some_array.npy')
Out[215]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
通過np.savez可以將多個數組保存到一個未壓縮文件中,將數組以關鍵字參數的形式傳入即可:
In [216]: np.savez('array_archive.npz', a=arr, b=arr)
加載.npz文件時,你會得到一個類似字典的對象,該對象會對各個數組進行延遲加載:
In [217]: arch = np.load('array_archive.npz')
In [218]: arch['b']
Out[218]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
如果要將數據壓縮,可以使用numpy.savez_compressed:
In [219]: np.savez_compressed('arrays_compressed.npz', a=arr, b=arr)
線性代數
線性代數(如矩陣乘法、矩陣分解、行列以及其他方陣數學等)是任何數組庫的重要組成部分。不像某些語言(如matlab),通過*對兩個二維數組相乘得到的是一個元素級的積,而不是一個矩陣點積。因此,Numpy提供了一個用于矩陣乘法的dot函數(既是一個數組方法也是numpy命名空間的一個函數):
In [223]: x = np.array([[1., 2., 3.], [4., 5., 6.]])
In [224]: y = np.array([[6., 23.], [-1, 7], [8, 9]])
In [225]: x
Out[225]:
array([[ 1., 2., 3.],
[ 4., 5., 6.]])
In [226]: y
Out[226]:
array([[ 6., 23.],
[ -1., 7.],
[ 8., 9.]])
In [227]: x.dot(y)
Out[227]:
array([[ 28., 64.],
[ 67., 181.]])
x.dot(y)等價于np.dot(x, y)
In [228]: np.dot(x, y)
Out[228]:
array([[ 28., 64.],
[ 67., 181.]])
一個二維數組跟一個大小合適的一維數組的矩陣點積運算之后將會得到一個一維數組:
In [229]: np.dot(x, np.ones(3))
Out[229]: array([ 6., 15.])
@符也可以用來用作運算符,進行矩陣乘法:
In [230]: x @ np.ones(3)
Out[230]: array([ 6., 15.])
numpy.linalg中有一組標準的矩陣分解運算以及諸如求逆和行列式之類的東西。它們跟matlab和R等語言所使用的是相同的行業標準線性代數庫,如BLAS、LAPACK、Intel MKL(Math Kernel Library,可能有,取決于你的NumPy版本)等:
In [231]: from numpy.linalg import inv, qr
In [232]: X = np.random.randn(5, 5)
In [233]: mat = X.T.dot(X)
In [234]: inv(mat)
Out[234]:
array([[ 933.1189, 871.8258, -1417.6902, -1460.4005, 1782.1391],
[ 871.8258, 815.3929, -1325.9965, -1365.9242, 1666.9347],
[-1417.6902, -1325.9965, 2158.4424, 2222.0191, -2711.6822],
[-1460.4005, -1365.9242, 2222.0191, 2289.0575, -2793.422 ],
[ 1782.1391, 1666.9347, -2711.6822, -2793.422 , 3409.5128]])
In [235]: mat.dot(inv(mat))
Out[235]:
array([[ 1., 0., -0., -0., -0.],
[-0., 1., 0., 0., 0.],
[ 0., 0., 1., 0., 0.],
[-0., 0., 0., 1., -0.],
[-0., 0., 0., 0., 1.]])
In [236]: q, r = qr(mat)
In [237]: r
Out[237]:
array([[-1.6914, 4.38 , 0.1757, 0.4075, -0.7838],
[ 0. , -2.6436, 0.1939, -3.072 , -1.0702],
[ 0. , 0. , -0.8138, 1.5414, 0.6155],
[ 0. , 0. , 0. , -2.6445, -2.1669],
[ 0. , 0. , 0. , 0. , 0.0002]])
表達式X.T.dot(X)計算X和它的轉置X.T的點積。
表4-7列出了一些最常用的線性代數函數:
偽隨機數生成
numpy.random模塊對Python內置的random進行了補充,增加了一些用于高效生成多種概率分布的樣本值的函數。例如,你可以用normal來得到一個標準正態分布的4×4樣本數組:
In [238]: samples = np.random.normal(size=(4, 4))
In [239]: samples
Out[239]:
array([[ 0.5732, 0.1933, 0.4429, 1.2796],
[ 0.575 , 0.4339, -0.7658, -1.237 ],
[-0.5367, 1.8545, -0.92 , -0.1082],
[ 0.1525, 0.9435, -1.0953, -0.144 ]])
而Python內置的random模塊則只能一次生成一個樣本值。從下面的測試結果可以看出,如果需要產生大量樣本值,numpy.random快了不止一個數量級:
In [240]: from random import normalvariate
In [241]: N = 1000000
In [242]: %timeit samples = [normalvariate(0, 1) for _ in range(N)]
1.77 s +- 126 ms per loop (mean +- std. dev. of 7 runs, 1 loop each)
In [243]: %timeit np.random.normal(size=N)
61.7 ms +- 1.32 ms per loop (mean +- std. dev. of 7 runs, 10 loops each)
我們說這些都是偽隨機數,是因為它們都是通過算法基于隨機數生成器種子,在確定性的條件下生成的。你可以用Numpy的np.random.seed更改隨機數生成種子:
In [244]: np.random.seed(1234)
numpy.random的數據生成函數使用了全局的隨機種子。要避免全局狀態,你可以使用numpy.random.RandomState,創建一個與其它隔離的隨機數生成器:
In [245]: rng = np.random.RandomState(1234)
In [246]: rng.randn(10)
Out[246]:
array([ 0.4714, -1.191 , 1.4327, -0.3127, -0.7206, 0.8872, 0.8596,
-0.6365, 0.0157, -2.2427])
表4-8列出了numpy.random中的部分函數。在下一節中,我將給出一些利用這些函數一次性生成大量樣本值的范例。
示例:隨機漫步
我們通過模擬隨機漫步來說明如何運用數組運算,先來看一個簡單的隨機漫步的例子:從0開始,步長1和-1出現的概率相等。
下面是一個通過內置的random模塊以純Python的方式實現1000步的隨機漫步:
In [247]: import random
.....: position = 0
.....: walk = [position]
.....: steps = 1000
.....: for i in range(steps):
.....: step = 1 if random.randint(0, 1) else -1
.....: position += step
.....: walk.append(position)
.....:
圖4-4是根據前100個隨機漫步值生成的折線圖:
In [249]: plt.plot(walk[:100])
圖4-4 簡單的隨機漫步.png
不難看出,這其實就是隨機漫步中各步的累計和,可以用一個數組運算來實現。因此,我用np.random模塊一次性隨機產生1000個“擲硬幣”結果(即兩個數中任選一個),將其分別設置成1或-1,然后計算累計和:
In [251]: nsteps = 1000
In [252]: draws = np.random.randint(0, 2, size=nsteps)
In [253]: steps = np.where(draws > 0, 1, -1)
In [254]: walk = steps.cumsum()
有了這些數據之后,我們就可以沿著漫步路徑做一些統計工作了,比如求取最大值和最小值:
In [255]: walk.min()
Out[255]: -3
In [256]: walk.max()
Out[256]: 31
現在來看一個復雜點的統計任務——首次穿越時間,即隨機漫步過程中第一次到達某個特定值的時間。假設我們想要知道本次隨機漫步需要多久才能舉例初始0點至少10步遠(任一方向均可)。np.abs(walk)>=10可以得到一個布爾型數組,它表示的是距離是否達到或超過10,而我們想要知道的是第一個10或-10的索引。可以用argmax來解決這個問題,它返回的是該布爾型數組第一個最大值的索引(True就是最大值):
In [257]: (np.abs(walk) >= 10).argmax()
Out[257]: 37
注意,這里使用argmax并不是很高效,因為它無論如何都會對數組進行完全掃描。在本例中,只要發現了一個True,那我們就知道它是個最大值了。
一次模擬多個隨機漫步
如果你希望模擬多個隨機漫步過程(比如5000個),只需對上面的代碼做一點點修改即可生成所有的隨機漫步過程。只要給numpy.random的函數傳入一個二元元組就可以產生一個二維數組,然后我們就可以一次性計算5000個隨機漫步過程(一行一個)的累計和了:
In [258]: nwalks = 5000
In [259]: nsteps = 1000
In [260]: draws = np.random.randint(0, 2, size=(nwalks, nsteps)) # 0 or 1
In [261]: steps = np.where(draws > 0, 1, -1)
In [262]: walks = steps.cumsum(1)
In [263]: walks
Out[263]:
array([[ 1, 0, 1, ..., 8, 7, 8],
[ 1, 0, -1, ..., 34, 33, 32],
[ 1, 0, -1, ..., 4, 5, 4],
...,
[ 1, 2, 1, ..., 24, 25, 26],
[ 1, 2, 3, ..., 14, 13, 14],
[ -1, -2, -3, ..., -24, -23, -22]])
現在,我們來計算所有隨機漫步過程的最大值和最小值:
In [264]: walks.max()
Out[264]: 138
In [265]: walks.min()
Out[265]: -133
得到這些數據之后,我們來計算30或-30的最小穿越時間。這里稍微復雜些,因為不是5000個過程都到達了30。我們可以用any方法來對此進行檢查:
In [266]: hits30 = (np.abs(walks) >= 30).any(1)
In [267]: hits30
Out[267]: array([False, True, False, ..., False, True, False], dtype=bool)
In [268]: hits30.sum() # Number that hit 30 or -30
Out[268]: 3410
然后我們利用這個布爾型數組選出那些穿越了30(絕對值)的隨機漫步(行),并調用argmax在軸1上獲取穿越時間:
In [269]: crossing_times = (np.abs(walks[hits30]) >= 30).argmax(1)
In [270]: crossing_times.mean()
Out[270]: 498.88973607038122
請嘗試用其他分布方式得到漫步數據。只需使用不同的隨機數生成函數即可,如normal用于生成指定均值和標準差的正態分布數據:
In [271]: steps = np.random.normal(loc=0, scale=0.25,
.....: size=(nwalks, nsteps))