標籤：

【翻譯】《利用Python進行數據分析·第2版》第4章（下）NumPy基礎：數組和矢量計算

04-27

作者：SeanCheney Python愛好者社區專欄作者

簡書專欄：jianshu.com/u/130f76596

前文傳送門：

【翻譯】《利用Python進行數據分析·第2版》第1章 準備工作

【翻譯】《利用Python進行數據分析·第2版》第2章（上）Python語法基礎，IPython和Jupyter

【翻譯】《利用Python進行數據分析·第2版》第2章（中）Python語法基礎，IPython和Jupyter

【翻譯】《利用Python進行數據分析·第2版》第2章（下）Python語法基礎，IPython和Jupyter

【翻譯】《利用Python進行數據分析·第2版》第3章（上）Python的數據結構、函數和文件

【翻譯】《利用Python進行數據分析·第2版》第3章（中）Python的數據結構、函數和文件

【翻譯】《利用Python進行數據分析·第2版》第3章（下）Python的數據結構、函數和文件

4.3 利用數組進行數據處理

NumPy數組使你可以將許多種數據處理任務表述為簡潔的數組表達式（否則需要編寫循環）。用數組表達式代替循環的做法，通常被稱為矢量化。一般來說，矢量化數組運算要比等價的純Python方式快上一兩個數量級（甚至更多），尤其是各種數值計算。在後面內容中（見附錄A）我將介紹廣播，這是一種針對矢量化計算的強大手段。

作為簡單的例子，假設我們想要在一組值（網格型）上計算函數sqrt(x^2+y^2)。np.meshgrid函數接受兩個一維數組，併產生兩個二維矩陣（對應於兩個數組中所有的(x,y)對）：

In [155]: points = np.arange(-5, 5, 0.01) # 1000 equally spaced points

In [156]: xs, ys = np.meshgrid(points, points)

In [157]: ys

Out[157]:

array([[-5. , -5. , -5. , ..., -5. , -5. , -5. ],

[-4.99, -4.99, -4.99, ..., -4.99, -4.99, -4.99],

[-4.98, -4.98, -4.98, ..., -4.98, -4.98, -4.98],

...,

[ 4.97, 4.97, 4.97, ..., 4.97, 4.97, 4.97],

[ 4.98, 4.98, 4.98, ..., 4.98, 4.98, 4.98],

[ 4.99, 4.99, 4.99, ..., 4.99, 4.99, 4.99]])

現在，對該函數的求值運算就好辦了，把這兩個數組當做兩個浮點數那樣編寫表達式即可：

In [158]: z = np.sqrt(xs ** 2 + ys ** 2)

In [159]: z

Out[159]:

array([[ 7.0711, 7.064 , 7.0569, ..., 7.0499, 7.0569, 7.064 ],

[ 7.064 , 7.0569, 7.0499, ..., 7.0428, 7.0499, 7.0569],

[ 7.0569, 7.0499, 7.0428, ..., 7.0357, 7.0428, 7.0499],

...,

[ 7.0499, 7.0428, 7.0357, ..., 7.0286, 7.0357, 7.0428],

[ 7.0569, 7.0499, 7.0428, ..., 7.0357, 7.0428, 7.0499],

[ 7.064 , 7.0569, 7.0499, ..., 7.0428, 7.0499, 7.0569]])

作為第9章的先導，我用matplotlib創建了這個二維數組的可視化：

In [160]: import matplotlib.pyplot as plt

In [161]: plt.imshow(z, cmap=plt.cm.gray); plt.colorbar()

Out[161]: <matplotlib.colorbar.Colorbar at 0x7f715e3fa630>

In [162]: plt.title("Image plot of $sqrt{x^2 + y^2}$ for a grid of values")

Out[162]: <matplotlib.text.Text at 0x7f715d2de748>

見圖4-3。這張圖是用matplotlib的imshow函數創建的。

圖4-3 根據網格對函數求值的結果

將條件邏輯表述為數組運算

numpy.where函數是三元表達式x if condition else y的矢量化版本。假設我們有一個布爾數組和兩個值數組：

In [165]: xarr = np.array([1.1, 1.2, 1.3, 1.4, 1.5])

In [166]: yarr = np.array([2.1, 2.2, 2.3, 2.4, 2.5])

In [167]: cond = np.array([True, False, True, True, False])

假設我們想要根據cond中的值選取xarr和yarr的值：當cond中的值為True時，選取xarr的值，否則從yarr中選取。列表推導式的寫法應該如下所示：

In [168]: result = [(x if c else y)

.....: for x, y, c in zip(xarr, yarr, cond)]

In [169]: result

Out[169]: [1.1000000000000001, 2.2000000000000002, 1.3, 1.3999999999999999, 2.5]

這有幾個問題。第一，它對大數組的處理速度不是很快（因為所有工作都是由純Python完成的）。第二，無法用於多維數組。若使用np.where，則可以將該功能寫得非常簡潔：

In [170]: result = np.where(cond, xarr, yarr)

In [171]: result

Out[171]: array([ 1.1, 2.2, 1.3, 1.4, 2.5])

np.where的第二個和第三個參數不必是數組，它們都可以是標量值。在數據分析工作中，where通常用於根據另一個數組而產生一個新的數組。假設有一個由隨機數據組成的矩陣，你希望將所有正值替換為2，將所有負值替換為－2。若利用np.where，則會非常簡單：

In [172]: arr = np.random.randn(4, 4)

In [173]: arr

Out[173]:

array([[-0.5031, -0.6223, -0.9212, -0.7262],

[ 0.2229, 0.0513, -1.1577, 0.8167],

[ 0.4336, 1.0107, 1.8249, -0.9975],

[ 0.8506, -0.1316, 0.9124, 0.1882]])

In [174]: arr > 0

Out[174]:

array([[False, False, False, False],

[ True, True, False, True],

[ True, True, True, False],

[ True, False, True, True]], dtype=bool)

In [175]: np.where(arr > 0, 2, -2)

Out[175]:

array([[-2, -2, -2, -2],

[ 2, 2, -2, 2],

[ 2, 2, 2, -2],

[ 2, -2, 2, 2]])

使用np.where，可以將標量和數組結合起來。例如，我可用常數2替換arr中所有正的值：

In [176]: np.where(arr > 0, 2, arr) # set only positive values to 2

Out[176]:

array([[-0.5031, -0.6223, -0.9212, -0.7262],

[ 2. , 2. , -1.1577, 2. ],

[ 2. , 2. , 2. , -0.9975],

[ 2. , -0.1316, 2. , 2. ]])

傳遞給where的數組大小可以不相等，甚至可以是標量值。

數學和統計方法

可以通過數組上的一組數學函數對整個數組或某個軸向的數據進行統計計算。sum、mean以及標準差std等聚合計算（aggregation，通常叫做約簡（reduction））既可以當做數組的實例方法調用，也可以當做頂級NumPy函數使用。

這裡，我生成了一些正態分布隨機數據，然後做了聚類統計：

In [177]: arr = np.random.randn(5, 4)

In [178]: arr

Out[178]:

array([[ 2.1695, -0.1149, 2.0037, 0.0296],

[ 0.7953, 0.1181, -0.7485, 0.585 ],

[ 0.1527, -1.5657, -0.5625, -0.0327],

[-0.929 , -0.4826, -0.0363, 1.0954],

[ 0.9809, -0.5895, 1.5817, -0.5287]])

In [179]: arr.mean()

Out[179]: 0.19607051119998253

In [180]: np.mean(arr)

Out[180]: 0.19607051119998253

In [181]: arr.sum()

Out[181]: 3.9214102239996507

mean和sum這類的函數可以接受一個axis選項參數，用於計算該軸向上的統計值，最終結果是一個少一維的數組：

In [182]: arr.mean(axis=1)

Out[182]: array([ 1.022 , 0.1875, -0.502 , -0.0881, 0.3611])

In [183]: arr.sum(axis=0)

Out[183]: array([ 3.1693, -2.6345, 2.2381, 1.1486])

這裡，arr.mean(1)是「計算行的平均值」，arr.sum(0)是「計算每列的和」。

其他如cumsum和cumprod之類的方法則不聚合，而是產生一個由中間結果組成的數組：

In [184]: arr = np.array([0, 1, 2, 3, 4, 5, 6, 7])

In [185]: arr.cumsum()

Out[185]: array([ 0, 1, 3, 6, 10, 15, 21, 28])

在多維數組中，累加函數（如cumsum）返回的是同樣大小的數組，但是會根據每個低維的切片沿著標記軸計算部分聚類：

In [186]: arr = np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])

In [187]: arr

Out[187]:

array([[0, 1, 2],

[3, 4, 5],

[6, 7, 8]])

In [188]: arr.cumsum(axis=0)

Out[188]:

array([[ 0, 1, 2],

[ 3, 5, 7],

[ 9, 12, 15]])

In [189]: arr.cumprod(axis=1)

Out[189]:

array([[ 0, 0, 0],

[ 3, 12, 60],

[ 6, 42, 336]])

表4-5列出了全部的基本數組統計方法。後續章節中有很多例子都會用到這些方法。

用於布爾型數組的方法

在上面這些方法中，布爾值會被強制轉換為1（True）和0（False）。因此，sum經常被用來對布爾型數組中的True值計數：

In [190]: arr = np.random.randn(100)In [191]: (arr > 0).sum() # Number of positive valuesOut[191]: 42

另外還有兩個方法any和all，它們對布爾型數組非常有用。any用於測試數組中是否存在一個或多個True，而all則檢查數組中所有值是否都是True：

In [192]: bools = np.array([False, False, True, False])

In [193]: bools.any()

Out[193]: True

In [194]: bools.all()

Out[194]: False

這兩個方法也能用於非布爾型數組，所有非0元素將會被當做True。

排序

跟Python內置的列表類型一樣，NumPy數組也可以通過sort方法就地排序：

In [195]: arr = np.random.randn(6)

In [196]: arr

Out[196]: array([ 0.6095, -0.4938, 1.24 , -0.1357, 1.43 , -0.8469])

In [197]: arr.sort()

In [198]: arr

Out[198]: array([-0.8469, -0.4938, -0.1357, 0.6095, 1.24 , 1.43 ])

多維數組可以在任何一個軸向上進行排序，只需將軸編號傳給sort即可：

In [199]: arr = np.random.randn(5, 3)

In [200]: arr

Out[200]:

array([[ 0.6033, 1.2636, -0.2555],

[-0.4457, 0.4684, -0.9616],

[-1.8245, 0.6254, 1.0229],

[ 1.1074, 0.0909, -0.3501],

[ 0.218 , -0.8948, -1.7415]])

In [201]: arr.sort(1)

In [202]: arr

Out[202]:

array([[-0.2555, 0.6033, 1.2636],

[-0.9616, -0.4457, 0.4684],

[-1.8245, 0.6254, 1.0229],

[-0.3501, 0.0909, 1.1074],

[-1.7415, -0.8948, 0.218 ]])

頂級方法np.sort返回的是數組的已排序副本，而就地排序則會修改數組本身。計算數組分位數最簡單的辦法是對其進行排序，然後選取特定位置的值：

In [203]: large_arr = np.random.randn(1000)

In [204]: large_arr.sort()

In [205]: large_arr[int(0.05 * len(large_arr))] # 5% quantile

Out[205]: -1.5311513550102103

更多關於NumPy排序方法以及諸如間接排序之類的高級技術，請參閱附錄A。在pandas中還可以找到一些其他跟排序有關的數據操作（比如根據一列或多列對表格型數據進行排序）。

唯一化以及其它的集合邏輯

NumPy提供了一些針對一維ndarray的基本集合運算。最常用的可能要數np.unique了，它用於找出數組中的唯一值並返回已排序的結果：

In [206]: names = np.array([Bob, Joe, Will, Bob, Will, Joe, Joe])

In [207]: np.unique(names)

Out[207]:

array([Bob, Joe, Will],

dtype=<U4)

In [208]: ints = np.array([3, 3, 3, 2, 2, 1, 1, 4, 4])

In [209]: np.unique(ints)

Out[209]: array([1, 2, 3, 4])

拿跟np.unique等價的純Python代碼來對比一下：

In [210]: sorted(set(names))

Out[210]: [Bob, Joe, Will]

另一個函數np.in1d用於測試一個數組中的值在另一個數組中的成員資格，返回一個布爾型數組：

In [211]: values = np.array([6, 0, 0, 3, 2, 5, 6])

In [212]: np.in1d(values, [2, 3, 6])

Out[212]: array([ True, False, False, True, True, False, True], dtype=bool)

NumPy中的集合函數請參見表4-6。

4.4 用於數組的文件輸入輸出

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)

4.5 線性代數

線性代數（如矩陣乘法、矩陣分解、行列式以及其他方陣數學等）是任何數組庫的重要組成部分。不像某些語言（如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.])

@符（類似Python 3.5）也可以用作中綴運算符，進行矩陣乘法：

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中列出了一些最常用的線性代數函數。

4.6 偽隨機數生成

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中的部分函數。在下一節中，我將給出一些利用這些函數一次性生成大量樣本值的範例。

4.7 示例：隨機漫步

我們通過模擬隨機漫步來說明如何運用數組運算。先來看一個簡單的隨機漫步的例子：從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 簡單的隨機漫步

不難看出，這其實就是隨機漫步中各步的累計和，可以用一個數組運算來實現。因此，我用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))

4.8 結論

雖然本書剩下的章節大部分是用pandas規整數據，我們還是會用到相似的基於數組的計算。在附錄A中，我們會深入挖掘NumPy的特點，進一步學習數組的技巧。


推薦閱讀：

※Python 2.7安裝常用的科學計算分析包
※python numpy的樣本標準差怎麼寫?
※ImagePy教程 —— 打開圖像
※ImagePy開發文檔 —— 圖像封裝類
※為什麼用 Numpy 還是慢, 你用對了嗎?

TAG:Python | 數據分析 | numpy |