Numpy快速上手指南 --- 進階篇
文章目錄
歡迎點選原文連結，Fork到自己的工作區進行實現和除錯。
閱讀過本文的讀者可以點選
1. 廣播法則(rule)
廣播法則能使通用函式有意義地處理不具有相同形狀的輸入。
-
廣播第一法則：如果所有的輸入陣列維度不都相同，一個“1”將被重複地新增在維度較小的陣列上直至所有的陣列擁有一樣的維度。
-
廣播第二法則：確定長度為1的陣列沿著特殊的方向表現地好像它有沿著那個方向最大形狀的大小。對陣列來說，沿著那個維度的陣列元素的值理應相同。
應用廣播法則之後，所有陣列的大小必須匹配。更多細節可以從這個[文件找到]。
2. 花哨的索引和索引技巧
Numpy比普通Python序列提供更多的索引功能。除了索引整數和切片，正如我們之前看到的，陣列可以被整數陣列和布林陣列索引。
通過陣列索引
from numpy import *
a = arange(12)**2 # the first 12 square numbers
i = array( [ 1,1,3,8,5 ] ) # an array of indices
a[i] # the elements of a at the positions i
>>>
array([ 1, 1, 9, 64, 25])
j = array( [ [ 3, 4], [ 9, 7 ] ] ) # a bidimensional array of indices a[j] # the same shape as j >>> array([[ 9, 16], [81, 49]])
當被索引陣列a是多維的時，每一個唯一的索引數列指向a的第一維5。以下示例通過將圖片標籤用調色版轉換成色彩影象展示了這種行為。
palette = array( [ [0,0,0], # 黑色
[255,0,0], # 紅色
[0,255,0], # 綠色
[0,0,255], # 藍色
[255,255,255] ] ) # 白色
image = array( [ [ 0, 1, 2, 0 ], # each value corresponds to a color in the palette
[ 0, 3, 4, 0 ] ] )
palette[image] # the (2,4,3) color image
>>>
array([[[ 0, 0, 0],
[255, 0, 0],
[ 0, 255, 0],
[ 0, 0, 0]],
[[ 0, 0, 0],
[ 0, 0, 255],
[255, 255, 255],
[ 0, 0, 0]]])
我們也可以給出不不止一維的索引，每一維的索引陣列必須有相同的形狀。
a = arange(12).reshape(3,4)
a
>>>
array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
i = array( [ [0,1], # indices for the first dim of a
[1,2] ] )
j = array( [ [2,1], # indices for the second dim
[3,3] ] )
a[i,j] # i and j must have equal shape
>>>
array([[ 2, 5],
[ 7, 11]])
a[i,2]
>>>
array([[ 2, 6],
[ 6, 10]])
a[:,j] # i.e., a[ : , j]
>>>
array([[[ 2, 1],
[ 3, 3]],
[[ 6, 5],
[ 7, 7]],
[[10, 9],
[11, 11]]])
自然，我們可以把i和j放到序列中(比如說列表)然後通過list索引。
l = [i,j]
a[l] # 與 a[i,j] 相等
>>>
array([[ 2, 5],
[ 7, 11]])
然而，我們不能把i和j放在一個數組中，因為這個陣列將被解釋成索引a的第一維。
s = array( [i,j] )
a[s] # not what we want
# ---------------------------------------------------------------------------
# IndexError Traceback (most recent call last)
# <ipython-input-100-b912f631cc75> in <module>()
# ----> 1 a[s]
# IndexError: index (3) out of range (0<=index<2) in dimension 0
a[tuple(s)] # same as a[i,j]
另一個常用的陣列索引用法是搜尋時間序列最大值6。
time = linspace(20, 145, 5) # time scale
data = sin(arange(20)).reshape(5,4) # 4 time-dependent series
time
>>>
array([ 20. , 51.25, 82.5 , 113.75, 145. ])
data
>>>
array([[ 0. , 0.84147098, 0.90929743, 0.14112001],
[-0.7568025 , -0.95892427, -0.2794155 , 0.6569866 ],
[ 0.98935825, 0.41211849, -0.54402111, -0.99999021],
[-0.53657292, 0.42016704, 0.99060736, 0.65028784],
[-0.28790332, -0.96139749, -0.75098725, 0.14987721]])
ind = data.argmax(axis=0) # index of the maxima for each series
ind
>>>
array([2, 0, 3, 1])
time_max = time[ ind] # times corresponding to the maxima
data_max = data[ind, range(data.shape[1])] # => data[ind[0],0], data[ind[1],1]...
time_max
>>>
array([ 82.5 , 20. , 113.75, 51.25])
data_max
>>>
array([ 0.98935825, 0.84147098, 0.99060736, 0.6569866 ])
all(data_max == data.max(axis=0)) # True
>>>
True
你也可以使用陣列索引作為目標來賦值：
a = arange(5)
a
>>>
array([0, 1, 2, 3, 4])
a[[1,3,4]] = 0
a
>>>
array([0, 0, 2, 0, 0])
然而，當一個索引列表包含重複時，賦值被多次完成，保留最後的值：
a = arange(5)
a[[0,0,2]]=[1,2,3]
a
>>>
array([2, 1, 3, 3, 4])
這足夠合理，但是小心如果你想用Python的+=結構，可能結果並非你所期望：
a = arange(5)
a[[0,0,2]]+=1
a
>>>
array([1, 1, 3, 3, 4])
即使0在索引列表中出現兩次，索引為0的元素僅僅增加一次。這是因為Python要求a+=1和a=a+1等同。
通過布林陣列索引
當我們使用整數陣列索引陣列時，我們提供一個索引列表去選擇。通過布林陣列索引的方法是不同的我們顯式地選擇陣列中我們想要和不想要的元素。
我們能想到的使用布林陣列的索引最自然方式就是使用和原陣列一樣形狀的布林陣列。
a = arange(12).reshape(3,4)
b = a > 4
b # b is a boolean with a's shape
>>>
array([[False, False, False, False],
[False, True, True, True],
[ True, True, True, True]], dtype=bool)
a[b] # 1d array with the selected elements
>>>
array([ 5, 6, 7, 8, 9, 10, 11])
這個屬性在賦值時非常有用：
a[b] = 0 # All elements of 'a' higher than 4 become 0
a
>>>
array([[0, 1, 2, 3],
[4, 0, 0, 0],
[0, 0, 0, 0]])
你可以參考曼德博集合示例看看如何使用布林索引來生成曼德博集合的影象。
第二種通過布林來索引的方法更近似於整數索引；對陣列的每個維度我們給一個一維布林陣列來選擇我們想要的切片。
a = arange(12).reshape(3,4)
b1 = array([False,True,True]) # first dim selection
b2 = array([True,False,True,False]) # second dim selection
a[b1,:] # selecting rows
>>>
array([[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
a[b1] # same thing
>>>
array([[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
a[:,b2] # selecting columns
>>>
array([[ 0, 2],
[ 4, 6],
[ 8, 10]])
a[b1,b2] # a weird thing to do
>>>
array([ 4, 10])
注意一維陣列的長度必須和你想要切片的維度或軸的長度一致，在之前的例子中，b1是一個秩為1長度為三的陣列(a的行數)，b2(長度為4)與a的第二秩(列)相一致。
ix_()函式
ix_函式可以為了獲得多元組的結果而用來結合不同向量。例如，如果你想要用所有向量a、b和c元素組成的三元組來計算a+b*c：
a = array([2,3,4,5])
b = array([8,5,4])
c = array([5,4,6,8,3])
ax,bx,cx = ix_(a,b,c)
ax
array([[[2]],
[[3]],
[[4]],
[[5]]])
bx
array([[[8],
[5],
[4]]])
cx
array([[[5, 4, 6, 8, 3]]])
ax.shape, bx.shape, cx.shape
((4, 1, 1), (1, 3, 1), (1, 1, 5))
result = ax+bx*cx
result
array([[[42, 34, 50, 66, 26],
[27, 22, 32, 42, 17],
[22, 18, 26, 34, 14]],
[[43, 35, 51, 67, 27],
[28, 23, 33, 43, 18],
[23, 19, 27, 35, 15]],
[[44, 36, 52, 68, 28],
[29, 24, 34, 44, 19],
[24, 20, 28, 36, 16]],
[[45, 37, 53, 69, 29],
[30, 25, 35, 45, 20],
[25, 21, 29, 37, 17]]])
result[3,2,4]
17
a[3]+b[2]*c[4]
17
你也可以實行如下簡化：
def ufunc_reduce(ufct, *vectors):
vs = ix_(*vectors)
r = ufct.identity
for v in vs:
r = ufct(r,v)
return r
然後這樣使用它：
ufunc_reduce(add,a,b,c)
array([[[15, 14, 16, 18, 13],
[12, 11, 13, 15, 10],
[11, 10, 12, 14, 9]],
[[16, 15, 17, 19, 14],
[13, 12, 14, 16, 11],
[12, 11, 13, 15, 10]],
[[17, 16, 18, 20, 15],
[14, 13, 15, 17, 12],
[13, 12, 14, 16, 11]],
[[18, 17, 19, 21, 16],
[15, 14, 16, 18, 13],
[14, 13, 15, 17, 12]]])
這個reduce與ufunc.reduce(比如說add.reduce)相比的優勢在於它利用了廣播法則，避免了建立一個輸出大小乘以向量個數的引數陣列。
用字串索引
參見 RecordArray
3. 線性代數
簡單陣列運算
from numpy import *
from numpy.linalg import *
a = array([[1.0, 2.0], [3.0, 4.0]])
print (a)
>>>
[[ 1. 2.]
[ 3. 4.]]
a.transpose()
>>>
array([[ 1., 3.],
[ 2., 4.]])
inv(a)
>>>
array([[-2. , 1. ],
[ 1.5, -0.5]])
u = eye(2) # unit 2x2 matrix; "eye" represents "I"
u
>>>
array([[ 1., 0.],
[ 0., 1.]])
j = array([[0.0, -1.0], [1.0, 0.0]])
dot (j, j) # matrix product
>>>
array([[-1., 0.],
[ 0., -1.]])
trace(u) # trace
>>>2.0
y = array([[5.], [7.]])
solve(a, y)
>>>
array([[-3.],
[ 4.]])
eig(j)
# Parameters:
# square matrix
# Returns
# The eigenvalues, each repeated according to its multiplicity.
# The normalized (unit "length") eigenvectors, such that the
# column ``v[:,i]`` is the eigenvector corresponding to the
# eigenvalue ``w[i]``
(array([ 0.+1.j, 0.-1.j]),
>>>
array([[ 0.70710678+0.j , 0.70710678-0.j ],
[ 0.00000000-0.70710678j, 0.00000000+0.70710678j]]))
矩陣類
這是一個關於矩陣類的簡短介紹。
A = matrix('1.0 2.0; 3.0 4.0')
A
>>>
matrix([[ 1., 2.],
[ 3., 4.]])
type(A) # file where class is defined
>>>
numpy.matrixlib.defmatrix.matrix
A.T # transpose
>>>
matrix([[ 1., 3.],
[ 2., 4.]])
X = matrix('5.0 7.0')
Y = X.T
Y
>>>
matrix([[ 5.],
[ 7.]])
print (A*Y) # matrix multiplication
>>>
[[ 19.]
[ 43.]]
print (A.I) # inverse
>>>
[[-2. 1. ]
[ 1.5 -0.5]]
solve(A, Y) # solving linear equation
>>>
matrix([[-3.],
[ 4.]])
索引：比較矩陣和二維陣列
注意Numpy中陣列和矩陣有些重要的區別。Numpy提供了兩個基本的物件：一個N維陣列物件和一個通用函式物件。其它物件都是建構在它們之上 的。特別的，矩陣是繼承自Numpy陣列物件的二維陣列物件。對陣列和矩陣，索引都必須包含合適的一個或多個這些組合：整數標量、省略號 (ellipses)、整數列表;布林值，整數或布林值構成的元組，和一個一維整數或布林值陣列。矩陣可以被用作矩陣的索引，但是通常需要陣列、列表或者 其它形式來完成這個任務。
像平常在Python中一樣，索引是從0開始的。傳統上我們用矩形的行和列表示一個二維陣列或矩陣，其中沿著0軸的方向被穿過的稱作行，沿著1軸的方向被穿過的是列。9
讓我們建立陣列和矩陣用來切片：
A = arange(12)
A
>>>
array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11])
A.shape = (3,4)
M = mat(A.copy())
print (type(A)," ",type(M))
>>>
<class 'numpy.ndarray'> <class 'numpy.matrixlib.defmatrix.matrix'>
print (A)
>>>
[[ 0 1 2 3]
[ 4 5 6 7]
[ 8 9 10 11]]
print (M)
>>>
[[ 0 1 2 3]
[ 4 5 6 7]
[ 8 9 10 11]]
現在，讓我們簡單的切幾片。基本的切片使用切片物件或整數。例如，A[:]和M[:]的求值將表現得和Python索引很相似。然而要注意很重要的一點就是Numpy切片陣列不建立資料的副本;切片提供統一資料的檢視。
print (A[:])
print (A[:].shape)
>>>
[[ 0 1 2 3]
[ 4 5 6 7]
[ 8 9 10 11]]
(3, 4)
print (M[:])
print (M[:].shape)
>>>
[[ 0 1 2 3]
[ 4 5 6 7]
[ 8 9 10 11]]
(3, 4)
現在有些和Python索引不同的了：你可以同時使用逗號分割索引來沿著多個軸索引。
print (A[:,1])
print (A[:,1].shape)
>>>
[1 5 9]
(3,)
print (M[:,1]);
print (M[:,1].shape)
>>>
[[1]
[5]
[9]]
(3, 1)
注意最後兩個結果的不同。對二維陣列使用一個冒號產生一個一維陣列，然而矩陣產生了一個二維矩陣。
例如，一個M[2,:]切片產生了一個形狀為(1,4)的矩陣，相比之下，一個數組的切片總是產生一個最低可能維度11的陣列。例如，如果C是一個三維陣列，C[…,1]產生一個二維的陣列而C[1,:,1]產生一個一維陣列。從這時開始，如果相應的矩陣切片結果是相同的話，我們將只展示陣列切片的結果。
假如我們想要一個數組的第一列和第三列，一種方法是使用列表切片：
A[:,[1,3]]
>>>
array([[ 1, 3],
[ 5, 7],
[ 9, 11]])
稍微複雜點的方法是使用 take() 方法 method:
A[:,].take([1,3],axis=1)
>>>
array([[ 1, 3],
[ 5, 7],
[ 9, 11]])
如果我們想跳過第一行，我們可以這樣：
A[1:,].take([1,3],axis=1)
>>>
array([[ 5, 7],
[ 9, 11]])
或者我們僅僅使用A[1:,[1,3]]。還有一種方法是通過矩陣向量積(叉積)。
A[ix_((1,2),(1,3))]
>>>
array([[ 5, 7],
[ 9, 11]])
為了讀者的方便，在次寫下之前的矩陣：
A[ix_((1,2),(1,3))]
>>>
array([[ 5, 7],
[ 9, 11]])
現在讓我們做些更復雜的。比如說我們想要保留第一行大於1的列。一種方法是建立布林索引：
A[0,:]>1
>>>
array([False, False, True, True], dtype=bool)
A[:,A[0,:]>1]
>>>
array([[ 2, 3],
[ 6, 7],
[10, 11]])
就是我們想要的！但是索引矩陣沒這麼方便。
M[0,:]>1
>>>
matrix([[False, False, True, True]], dtype=bool)
indice太多造成無法索引
M[:,M[0,:]>1]
---------------------------------------------------------------------------
IndexError Traceback (most recent call last)
<ipython-input-65-79db943ff8d4> in <module>()
----> 1 M[:,M[0,:]>1]
~/anaconda/lib/python3.6/site-packages/numpy/matrixlib/defmatrix.py in __getitem__(self, index)
316
317 try:
--> 318 out = N.ndarray.__getitem__(self, index)
319 finally:
320 self._getitem = False
IndexError: too many indices for array
這個過程的問題是用“矩陣切片”來切片產生一個矩陣12，但是矩陣有個方便的A屬性，它的值是陣列呈現的。所以我們僅僅做以下替代：
M[:,M.A[0,:]>1]
>>>
matrix([[ 2, 3],
[ 6, 7],
[10, 11]])
如果我們想要在矩陣兩個方向有條件地切片，我們必須稍微調整策略，代之以：
A[A[:,0]>2,A[0,:]>1]
>>>
array([ 6, 11])
M[M.A[:,0]>2,M.A[0,:]>1]
>>>
matrix([[ 6, 11]])
我們需要使用向量積ix_:
A[ix_(A[:,0]>2,A[0,:]>1)]
>>>
array([[ 6, 7],
[10, 11]])
M[ix_(M.A[:,0]>2,M.A[0,:]>1)]
>>>
matrix([[ 6, 7],
[10, 11]])
4. 技巧和提示
下面我們給出簡短和有用的提示。
“自動” 改變形狀
更改陣列的維度，你可以省略一個尺寸，它將被自動推匯出來。
a = arange(30)
a.shape = 2,-1,3 # -1 means "whatever is needed"
a.shape
>>>
(2, 5, 3)
a
>>>
array([[[ 0, 1, 2],
[ 3, 4, 5],
[ 6, 7, 8],
[ 9, 10, 11],
[12, 13, 14]],
[[15, 16, 17],
[18, 19, 20],
[21, 22, 23],
[24, 25, 26],
[27, 28, 29]]])
向量組合(stacking)
我們如何用兩個相同尺寸的行向量列表構建一個二維陣列？在 MATLAB中這非常簡單：如果 x和 y是兩個相同長度的向量，你僅僅需要做 m=[x;y]。在 Numpy中這個過程通過函式 column_stack、dstack、hstack和 vstack來完成，取決於你想要在那個維度上組合。例如：
x = arange(0,10,2)
x
>>>
array([0, 2, 4, 6, 8])
y = arange(5)
y
>>>
array([0, 1, 2, 3, 4])
m = vstack([x,y])
m
>>>
array([[0, 2, 4, 6, 8],
[0, 1, 2, 3, 4]])
xy = hstack([x,y])
xy
array([0, 2, 4, 6, 8, 0, 1, 2, 3, 4])
直方圖(histogram)
Numpy中histogram函式應用到一個數組返回一對變數：直方圖陣列和箱式向量。注意：matplotlib也有一個用來建立直方圖的函式(叫作hist,正如matlab中一樣)與Numpy中的不同。
主要的差別是pylab.hist自動繪製直方圖，而numpy.histogram僅僅產生資料。
import numpy
import pylab
# Build a vector of 10000 normal deviates with variance 0.5^2 and mean 2
mu, sigma = 2, 0.5
v = numpy.random.normal(mu,sigma,10000)
# Plot a normalized histogram with 50 bins
pylab.hist(v, bins=50, normed=1)
pylab.title('Matplotlib Version')# matplotlib version (plot)
pylab.show()
# Compute the histogram with numpy and then plot it
(n, bins) = numpy.histogram(v, bins=50, normed=True) # NumPy version (no plot)
pylab.plot(.5*(bins[1:]+bins[:-1]), n)
pylab.title('Numpy Version')
pylab.show()
上述文中所有程式碼部分都可以使用線上資料分析協作工具K-Lab復現，點選連結一鍵直達~
K-Lab提供基於Jupyter Notebook的線上資料分析服務，涵蓋Python&R等主流程式語言，可滿足資料科學家、人工智慧工程師、商業分析師等資料工作者線上完成資料處理、模型搭建、程式碼除錯、撰寫分析報告等資料分析全過程。