第4章 NumPy基礎:陣列和向量計算
參考
目錄
NumPy(Numerical Python的簡稱)是Python數值計算最重要的基礎包。大多數提供科學計算的包都是用NumPy的陣列作為構建基礎。
NumPy的部分功能如下:
- ndarray,一個具有向量算術運算和複雜廣播能力的快速且節省空間的多維陣列。
- 用於對整組資料進行快速運算的標準數學函式(無需編寫迴圈)。
- 用於讀寫磁碟資料的工具以及用於操作記憶體對映檔案的工具。
- 線性代數、隨機數生成以及傅立葉變換功能。
- 用於整合由C、C++、Fortran等語言編寫的程式碼的A C API。
由於NumPy提供了一個簡單易用的C API,因此很容易將資料傳遞給由低階語言編寫的外部庫,外部庫也能以NumPy陣列的形式將資料返回給Python。NumPy本身並沒有提供多麼高階的資料分析功能,理解NumPy陣列以及面向陣列的計算將有助於你更加高效地使用諸如pandas之類的工具。
對於大部分資料分析應用而言,大家應該關注的功能主要為以下幾點:
- 用於資料整理和清理、子集構造和過濾、轉換等快速的向量化陣列運算。
- 常用的陣列演算法,如排序、唯一化、集合運算等。
- 高效的描述統計和資料聚合/摘要運算。
- 用於異構資料集的合併/連線運算的資料對齊和關係型資料運算。
- 將條件邏輯表述為陣列表示式(而不是帶有if-elif-else分支的迴圈)。
- 資料的分組運算(聚合、轉換、函式應用等)。。
處理表格資料時,通常將pandas作為統計和分析工作的基礎。
NumPy之於數值計算特別重要的原因之一,是因為它可以高效處理大陣列的資料。
- NumPy是在一個連續的記憶體塊中儲存資料,獨立於其他Python內建物件。NumPy的C語言編寫的演算法庫可以操作記憶體,而不必進行型別檢查或其它前期工作。比起Python的內建序列,NumPy陣列使用的記憶體更少。
- NumPy可以在整個陣列上執行復雜的計算,而不需要Python的for迴圈。
正如下面的例子:
import numpy as np hust_arr = np.arange(1000000) hust_list = list(range(1000000)) %time for _ in range(10): hust_arr2 = hust_arr * 2
out put:
Wall time: 19 ms
%time for _ in range(10): hust_list2 = [x * 2 for x in hust_list]
out put:
Wall time: 750 ms
基於NumPy的演算法要比純Python快10到100倍(甚至更快),並且使用的記憶體更少。
4.1 NumPy的ndarray:一種多維陣列物件
NumPy最重要的一個特點就是其N維陣列物件(即ndarray),該物件是一個快速而靈活的大資料集容器。利用陣列對整塊資料執行數學運算,其語法跟標量元素之間的運算一樣。
要明白Python是如何利用與標量值類似的語法進行批次計算,首先引入NumPy,然後生成一個包含隨機資料的小陣列:
import numpy as np
# Generate some random data
data = np.random.randn(2, 4)
data
array([[-0.10339226, 1.00605839, -1.86735608, -0.40424423],
[ 0.06615611, -1.51357686, -0.32021485, -0.25574128]])
然後進行數學運算:
data * 10
array([[ -1.03392263, 10.06058393, -18.6735608 , -4.04244231],
[ 0.66156112, -15.1357686 , -3.20214852, -2.55741283]])
data + data
array([[-0.20678453, 2.01211679, -3.73471216, -0.80848846],
[ 0.13231222, -3.02715372, -0.6404297 , -0.51148257]])
第一個例子中,所有的元素都乘以10。第二個例子中,每個元素都與自身相加。
ndarray是一個通用的同構資料多維容器,其中的所有元素必須是相同型別的。每個陣列都有一個shape(一個表示各維度大小的元組)和一個dtype(一個用於說明陣列資料型別的物件):
In [7]: data.shape
Out[7]: (2, 4)
In [8]: data.dtype
Out[8]: dtype('float64')
本章將會介紹NumPy陣列的基本用法,這對於本書後面各章的理解基本夠用。雖然大多數資料分析工作不需要深入理解NumPy,但是精通面向陣列的程式設計和思維方式是成為Python科學計算牛人的一大關鍵步驟。
4.1.1 建立ndarray
建立陣列最簡單的辦法就是使用array函式。它接受一切序列型的物件(包括其他陣列),然後產生一個新的含有傳入資料的NumPy陣列。下面以一個列表的轉換為例:
In [13]: data1 = [5, 3.4, 6, 2, 7]
arr1 = np.array(data1)
arr1
Out[13]: array([5. , 3.4, 6. , 2. , 7. ])
巢狀序列(比如由一組等長列表組成的列表)將會被轉換為一個多維陣列:
In [14]: data2 = [[4, 3, 2, 1], [10, 9, 8, 7]]
In [15]: arr2 = np.array(data2)
In [16]: arr2
Out[16]: array([[ 4, 3, 2, 1],
[10, 9, 8, 7]])
因為data2是列表的列表,NumPy陣列arr2的兩個維度的shape是從data2引入的。可以用屬性ndim和shape驗證:
In [17]: arr2.ndim
Out[17]: 2
In [18]: arr2.shape
Out[18]: (2, 4)
除非特別說明(稍後將會詳細介紹),np.array會嘗試為新建的這個陣列推斷出一個較為合適的資料型別。資料型別儲存在一個特殊的dtype物件中。比如說,在上面的兩個例子中,我們有:
In [19]: arr1.dtype
Out[19]: dtype('float64')
In [20]: arr2.dtype
Out[20]: dtype('int32')
除np.array之外,還有一些函式也可以新建陣列。比如,zeros和ones分別可以建立指定長度或形狀的全0或全1陣列。empty可以建立一個沒有任何具體值的陣列。要用這些方法建立多維陣列,只需傳入一個表示形狀的元組即可:
In [21]: np.zeros(10)
Out[21]: array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
In [22]: # 定義4行6列的陣列
np.zeros((4, 6))
Out[22]:
array([[0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0.]])
In [23]: # 定義一個三維陣列
np.empty((2, 3, 2))
Out[23]:
array([[[9.39944447e-312, 2.86558075e-322],
[0.00000000e+000, 0.00000000e+000],
[0.00000000e+000, 3.76231868e+174]],
[[4.96125191e-090, 1.57211215e-076],
[3.61112012e+174, 1.31505911e-047],
[6.48224660e+170, 5.82471487e+257]]])
需要對空陣列進行初始化,否則元素的值不一定全為0
arange是Python內建函式range的陣列版:
In [26]: np.arange(10)
Out[26]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
表4-1列出了一些陣列建立函式。由於NumPy關注的是數值計算,因此,如果沒有特別指定,資料型別基本都是float64(浮點數)。
4.1.2 ndarray的資料型別
dtype(資料型別)是一個特殊的物件,它含有ndarray將一塊記憶體解釋為特定資料型別所需的資訊:
In [27]: arr1 = np.array([0, 1, 2], dtype = np.float64)
In [28]: arr2 = np.array([4, 5, 6], dtype = np.int32)
In [29]: arr1.dtype
Out[29]: dtype('float64')
In [30]: arr2.dtype
Out[30]: dtype('int32')
標準的雙精度浮點值(即Python中的float物件)需要佔用8位元組(即64位)。因此,該型別在NumPy中就記作float64。表4-2列出了NumPy所支援的全部資料型別。
可以通過ndarray的astype方法明確地將一個數組從一個dtype轉換成另一個dtype:
In [31]: arr = np.array([0, 1, 2, 3, 4, 5, 6])
In [32]: arr.dtype
Out[32]: dtype('int32')
In [33]: float_arr = arr.astype(np.float)
In [34]: float_arr.dtype
Out[34]: dtype('float64')
在下個例子中,整數被轉換成了浮點數。如果將浮點數轉換成整數,則小數部分將會被擷取刪除:
In [35]: arr = np.array([2.5, 3.9, 1.4])
In [36]: arr
Out[36]: array([2.5, 3.9, 1.4])
In [37]: arr.astype(np.int32)
Out[37]: array([2, 3, 1])
如果某字串陣列表示的全是數字,也可以用astype將其轉換為數值形式:
In [39]: numeric_string = np.array(['3.12', '-8.23', '1932'], dtype = np.string_)
In [40]: numeric_string.astype(float)
Out[40]: array([ 3.12, -8.23, 1932. ])
注意:使用numpy.string_型別時,一定要小心,因為NumPy的字串資料是大小固定的,發生擷取時,不會發出警告。pandas提供了更多非數值資料的便利的處理方法。
陣列的dtype還有另一個屬性:
In [48]: int_array = np.arange(9)
In [49]: calibers = np.array([.23, .45, 89.34], dtype = np.float64)
In [50]: int_array.astype(calibers.dtype)
Out[50]: array([0., 1., 2., 3., 4., 5., 6., 7., 8.])
可以用簡潔的型別程式碼來表示dtype:
In [51]: empty_uint32 = np.empty(8, dtype = 'u4')
In [52]: empty_uint32
Out[52]:
array([4054629320, 442, 0, 0, 4072783887,
32843595, 3782017648, 2147491258], dtype=uint32)
4.1.3 NumPy陣列的運算
In [53]: arr = np.array([[1., 2., 3.], [6., 5., 4.]])
In [54]: arr
Out[54]:
array([[1., 2., 3.],
[6., 5., 4.]])
In [55]: arr *arr
Out[55]:
array([[ 1., 4., 9.],
[36., 25., 16.]])
In [56]: arr - arr
Out[56]:
array([[0., 0., 0.],
[0., 0., 0.]])
陣列與標量的算術運算會將標量值傳播到各個元素:
In [57]: 1 / arr
Out[57]:
array([[1. , 0.5 , 0.33333333],
[0.16666667, 0.2 , 0.25 ]])
In [58]: arr ** 2
Out[58]:
array([[ 1., 4., 9.],
[36., 25., 16.]])
大小相同的陣列之間的比較會生成布林值陣列:
In [59]: arr2 = np.array([[0., 4., 1.], [7., 2., 10.]])
In [60]: arr2
Out[60]:
array([[ 0., 4., 1.],
[ 7., 2., 10.]])
In [61]: arr2 > arr
Out[61]:
array([[False, True, False],
[ True, False, True]])
不同大小的陣列之間的運算叫做廣播(broadcasting),下面對廣播進行簡單介紹:
4.1.4 陣列的廣播
如果兩個陣列的維數不相同,則元素到元素的操作是不可能的。 然而,在 NumPy 中仍然可以對形狀不相似的陣列進行操作,因為它擁有廣播功能。 較小的陣列會廣播到較大陣列的大小,以便使它們的形狀可相容。
如果滿足以下規則,可以進行廣播:
- ndim較小的陣列會在前面追加一個長度為 1 的維度。
- 輸出陣列的每個維度的大小是輸入陣列該維度大小的最大值。
- 如果輸入在每個維度中的大小與輸出大小匹配,或其值正好為 1,則可以在計算中使用該輸入。
- 如果輸入的某個維度大小為 1,則該維度中的第一個資料元素將用於該維度的所有計算。
如果上述規則產生有效結果,並且滿足以下條件之一,那麼陣列被稱為可廣播的。
- 陣列擁有相同形狀。
- 陣列擁有相同的維數,每個維度擁有相同長度,或者長度為 1。
- 陣列擁有極少的維度,可以在其前面追加長度為 1 的維度,使上述條件成立。
請看下面例項:
import numpy as np
a = np.array([[0.0,0.0,0.0],[10.0,10.0,10.0],[20.0,20.0,20.0],[30.0,30.0,30.0]])
b = np.array([1.0,2.0,3.0])
print ('第一個陣列:')
print (a)
print ('\n第二個陣列:')
print (b)
print ('\n第一個陣列加第二個陣列:')
print (a + b)
output
第一個陣列:
[[ 0. 0. 0.]
[ 10. 10. 10.]
[ 20. 20. 20.]
[ 30. 30. 30.]]
第二個陣列:
[ 1. 2. 3.]
第一個陣列加第二個陣列:
[[ 1. 2. 3.]
[ 11. 12. 13.]
[ 21. 22. 23.]
[ 31. 32. 33.]]
下面的圖片展示了陣列b如何通過廣播來與陣列a相容。
4.1.5 基本的索引和切片
NumPy陣列的索引是一個內容豐富的主題,因為選取資料子集或單個元素的方式有很多。一維陣列很簡單。從表面上看,它們跟Python列表的功能差不多:
In [1]: import numpy as np
arr = np.arange(10)
In [2]: arr
Out[2]:
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
In [3]: arr[5]
Out[3]:
5
In [4]: arr[5:8] #第6~8個元素,左閉右開
Out[4]:
array([5, 6, 7])
In [5]: arr[5:8] = 12
In [6]: arr
Out[6]:
array([ 0, 1, 2, 3, 4, 12, 12, 12, 8, 9])
將一個標量值賦值給一個切片時(如arr[5:8]=12),該值會自動傳播到整個選區。跟列表最重要的區別在於,陣列切片是原始陣列的檢視。這意味著資料不會被複制,檢視上的任何修改都會直接反映到源陣列上。
作為例子,先建立一個arr的切片:
In [7]: arr_slice = arr[5:8]
In [8]: arr_slice
Out[8]: array([12, 12, 12])
修改arr_slice中的值,變動也會體現在原始陣列arr中:
In [9]: arr_slice[1] = 12345
In [10]: arr
Out[10]:
array([ 0, 1, 2, 3, 4, 12, 12345, 12, 8,
9])
切片[ : ]會給陣列中的所有值賦值:
In [11]: arr_slice[:] = 64
In [13]: arr
Out[13]:
array([ 0, 1, 2, 3, 4, 64, 64, 64, 8, 9])
注意:由於NumPy的設計目的是處理大資料,假如NumPy堅持要將資料複製來複制去的話會產生何等的效能和記憶體問題。所以ndarray就是直接在原始資料上進行操作,降低記憶體。
對於高維度陣列,能做的事情更多。在一個二維陣列中,各索引位置上的元素不再是標量而是一維陣列:
In [14]: arr2d = np.array([[1, 1, 1], [2, 3, 5], [2, 4, 5]])
In [15]: arr2d[2]
Out[15]: array([2, 4, 5])
因此,可以對各個元素進行遞迴訪問,但這樣需要做的事情有點多。你可以傳入一個以逗號隔開的索引列表來選取單個元素。也就是說,下面兩種方式是等價的:
In [16]: arr2d[0][2]
Out[16]: 1
In [17]: arr2d[0, 2]
Out[17]: 1
下圖說明了二維陣列的索引方式。軸0作為行,軸1作為列。
在多維陣列中,如果省略了後面的索引,則返回物件會是一個維度低一點的ndarray(它含有高一級維度上的所有資料)。因此,在2×2×3陣列arr3d中:
In [21]: arr3d = np.array([[[1, 2, 3],[4, 5, 6]],[[7, 8, 9],[10, 11, 12]]])
In [22]: arr3d
Out[22]:
array([[[ 1, 2, 3],
[ 4, 5, 6]],
[[ 7, 8, 9],
[10, 11, 12]]])
arr3d[0]是一個2×3陣列
In [23]: arr3d[0]
Out[23]:
array([[1, 2, 3],
[4, 5, 6]])
標量值和陣列都可以被賦值給arr3d[0]:
In [24]: old_values = arr3d[0].copy()
In [25]: arr3d[0] = 42
In [26]: arr3d
Out[26]:
array([[[42, 42, 42],
[42, 42, 42]],
[[ 7, 8, 9],
[10, 11, 12]]])
In [27]: arr3d[0] = old_values
In [28]: arr3d
Out[28]:
array([[[ 1, 2, 3],
[ 4, 5, 6]],
[[ 7, 8, 9],
[10, 11, 12]]])
相似的,arr3d[1,0]可以訪問索引以(1,0)開頭的那些值(以一維陣列的形式返回):
In [29]: arr3d[1,0]
Out[29]:
array([7, 8, 9])
可以分解成兩步進行,但結果是一致的。
In [30]: x = arr3d[1]
In [31]: x
Out[31]: array([[ 7, 8, 9], [10, 11, 12]])
In [32]: x[0]
Out[32]: array([7, 8, 9])
注意,在上面所有這些選取陣列子集的例子中,返回的陣列都是檢視。
4.1.6 切片索引
ndarray的切片語法跟Python列表這樣的一維物件差不多:
In [33]: arr
Out[33]: array([ 0, 1, 2, 3, 4, 64, 64, 64, 8, 9])
In [34]: arr[1:6]
Out[34]: array([ 1, 2, 3, 4, 64])
對於之前的二維陣列arr2d,其切片方式稍顯不同:
In [36]: arr2d[:2]
Out[36]: array([[1, 1, 1], [2, 3, 5]])
可以看出,它是沿著第0軸(即第一個軸)切片的。也就是說,切片是沿著一個軸向選取元素的。表示式arr2d[:2]可以被認為是“選取arr2d的前兩行”。
也可以一次傳入多個切片,就像傳入多個索引那樣:
In [37]: arr2d[:2, 1:]
Out[37]:
array([[1, 1],
[3, 5]])
像這樣進行切片時,只能得到相同維數的陣列檢視。通過將整數索引和切片混合,可以得到低維度的切片。
例如,我可以選取第三行的前兩列:
In [38]: arr2d[2, :2]
Out[38]: array([2, 4])
下圖對此進行了說明。注意,“只有冒號”表示選取整個軸,因此可以像下面這樣只對高維軸進行切片:
In [39]: arr2d[:, :1]
Out[39]:
array([[1],
[2],
[2]])
自然,對切片表示式的賦值操作也會被擴散到整個選區:
In [40]: arr2d[:2,1:] = 0
In [41]: arr2d
Out[41]:
array([[1, 0, 0],
[2, 0, 0],
[2, 4, 5]])
4.1.7 布林型索引
假設我們有一個用於儲存資料的陣列以及一個儲存姓名的陣列(含有重複項)。在這裡,我將使用numpy.random中的randn函式生成一些正態分佈的隨機資料:
In [42]: names = np.array(['Bob', 'Joe', 'Will', 'Bob', 'Will', 'Joe', 'Joe'])
data = np.random.randn(7,4)
In [43]: names
Out[43]: array(['Bob', 'Joe', 'Will', 'Bob', 'Will', 'Joe', 'Joe'], dtype='<U4')
In [44]: data
Out[44]:
array([[ 0.73603062, -0.12614702, -0.54407036, 1.21298923],
[-0.59730696, 1.2702293 , 0.71270184, -0.74228686],
[-0.12988137, 0.23263983, -0.08371573, 0.42319046],
[ 1.78970064, 0.20477398, -0.12707496, 1.96050673],
[ 0.71197402, 1.13861239, 1.58047673, 1.48136366],
[-0.70338155, -1.15144029, -0.50996925, 0.79713611],
[ 1.41090564, 0.40505494, -0.04298217, 0.51167567]])
假設每個名字都對應data陣列中的一行,而我們想要選出對應於名字"Bob"的所有行。跟算術運算一樣,陣列的比較運算(如==)也是向量化的。因此,對names和字串"Bob"的比較運算將會產生一個布林型陣列:
In [45]: names == 'Bob'
Out[45]: array([ True, False, False, True, False, False, False])
這個布林型陣列可用於陣列索引:
In [46]: data[names == 'Bob']
Out[46]:
array([[ 0.73603062, -0.12614702, -0.54407036, 1.21298923],
[ 1.78970064, 0.20477398, -0.12707496, 1.96050673]])
布林型陣列的長度必須跟被索引的軸長度一致。此外,還可以將布林型陣列跟切片、整數(或整數序列,稍後將對此進行詳細講解)混合使用:
下面的例子,我選取了names == 'Bob'
的行,並索引了列:
In [47]: data[names == 'Bob', 2:]
Out[47]:
array([[-0.54407036, 1.21298923],
[-0.12707496, 1.96050673]])
In [48]: data[names == 'Bob', 3]
Out[48]: array([1.21298923, 1.96050673])
要選擇除"Bob"以外的其他值,既可以使用不等於符號(!=),也可以通過~對條件進行否定:
In [49]: names != 'Bob'
Out[49]: array([False, True, True, False, True, True, True])
In [50]: data[~(names == 'Bob')]
Out[50]:
array([[-0.59730696, 1.2702293 , 0.71270184, -0.74228686],
[-0.12988137, 0.23263983, -0.08371573, 0.42319046],
[ 0.71197402, 1.13861239, 1.58047673, 1.48136366],
[-0.70338155, -1.15144029, -0.50996925, 0.79713611],
[ 1.41090564, 0.40505494, -0.04298217, 0.51167567]])
~操作符用來反轉條件很好用:
In [51]: cond = names == 'Bob'
In [52]: data[~cond]
Out[52]:
array([[-0.59730696, 1.2702293 , 0.71270184, -0.74228686],
[-0.12988137, 0.23263983, -0.08371573, 0.42319046],
[ 0.71197402, 1.13861239, 1.58047673, 1.48136366],
[-0.70338155, -1.15144029, -0.50996925, 0.79713611],
[ 1.41090564, 0.40505494, -0.04298217, 0.51167567]])
選取這三個名字中的兩個需要組合應用多個布林條件,使用&(和)、|(或)之類的布林算術運算子即可:
In [56]: mask = (names == 'Bob') | (names == 'Will')
In [57]: mask
Out[57]: array([ True, False, True, True, True, False, False])
In [58]: data[mask]
Out[58]:
array([[ 0.73603062, -0.12614702, -0.54407036, 1.21298923],
[-0.12988137, 0.23263983, -0.08371573, 0.42319046],
[ 1.78970064, 0.20477398, -0.12707496, 1.96050673],
[ 0.71197402, 1.13861239, 1.58047673, 1.48136366]])
通過布林型索引選取陣列中的資料,將總是建立資料的副本,即使返回一模一樣的陣列也是如此。
注意:Python關鍵字and和or在布林型陣列中無效。要使用&與|。
通過布林型陣列設定值是一種經常用到的手段。為了將data中的所有負值都設定為0,我們只需:
In [59]: data[data < 0] = 0
In [60]: data
Out[60]:
array([[0.73603062, 0. , 0. , 1.21298923],
[0. , 1.2702293 , 0.71270184, 0. ],
[0. , 0.23263983, 0. , 0.42319046],
[1.78970064, 0.20477398, 0. , 1.96050673],
[0.71197402, 1.13861239, 1.58047673, 1.48136366],
[0. , 0. , 0. , 0.79713611],
[1.41090564, 0.40505494, 0. , 0.51167567]])
通過一維布林陣列設定整行或列的值也很簡單:
In [61]: data[names == 'Joe'] = 7
In [62]: data
Out[62]:
array([[0.73603062, 0. , 0. , 1.21298923],
[7. , 7. , 7. , 7. ],
[0. , 0.23263983, 0. , 0.42319046],
[1.78970064, 0.20477398, 0. , 1.96050673],
[0.71197402, 1.13861239, 1.58047673, 1.48136366],
[7. , 7. , 7. , 7. ],
[7. , 7. , 7. , 7. ]])
當然,這些操作也可以用Pandas來做。
4.1.8 花式索引
花式索引(Fancy indexing)是一個NumPy術語,它指的是利用整數陣列進行索引。假設我們有一個8×4陣列:
In [2]: import numpy as np
arr = np.empty((8,4))
In [3]: for i in range(8):
arr[i] = i
In [4]: arr
Out[4]:
array([[0., 0., 0., 0.],
[1., 1., 1., 1.],
[2., 2., 2., 2.],
[3., 3., 3., 3.],
[4., 4., 4., 4.],
[5., 5., 5., 5.],
[6., 6., 6., 6.],
[7., 7., 7., 7.]])
為了以特定順序選取行子集,只需傳入一個用於指定順序的整數列表或ndarray即可:
In [5]: arr[[4, 3, 0, 6]]
Out[5]:
array([[4., 4., 4., 4.],
[3., 3., 3., 3.],
[0., 0., 0., 0.],
[6., 6., 6., 6.]])
這段程式碼確實達到我們的要求了!使用負數索引將會從末尾開始選取行:
In [6]:barr[[-3, -5, -7]]
Out[6]:
array([[5., 5., 5., 5.],
[3., 3., 3., 3.],
[1., 1., 1., 1.]])
一次傳入多個索引陣列會有一點特別。它返回的是一個一維陣列,其中的元素對應各個索引元組:
In [7]: arr = np.arange(32).reshape((8,4))
In [8]: arr
Out[8]:
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, 30, 31]])
In [9]: arr[[1, 5, 7, 2], [0, 3, 1, 2]]
Out[9]: array([ 4, 23, 29, 10])
最終選出的是元素(1,0)、(5,3)、(7,1)和(2,2)。無論陣列是多少維的,花式索引總是一維的。
記住,花式索引跟切片不一樣,它總是將資料複製到新陣列中。
4.1.9 陣列轉置和軸對換
轉置是重塑的一種特殊形式,它返回的是源資料的檢視(不會進行任何複製操作)。陣列不僅有transpose方法,還有一個特殊的T屬性:
In [10]: arr = np.arange(15).reshape(3, 5)
In [11]: arr
Out[11]:
array([[ 0, 1, 2, 3, 4],
[ 5, 6, 7, 8, 9],
[10, 11, 12, 13, 14]])
In [12]: arr.T
Out[12]:
array([[ 0, 5, 10],
[ 1, 6, 11],
[ 2, 7, 12],
[ 3, 8, 13],
[ 4, 9, 14]])
在進行矩陣計算時,經常需要用到該操作,比如利用np.dot計算矩陣內積:
In [13]: arr = np.random.randn(6, 3)
In [14]: arr
Out[14]:
array([[ 0.76416581, -0.59371434, 0.68181681],
[ 0.30954994, 0.73782933, 1.17631282],
[ 0.06555348, 0.33505602, 0.22756087],
[-0.05010979, 2.00622115, -1.17668862],
[ 0.66066926, 1.20410624, -0.27191449],
[ 1.48285532, 0.1020095 , -0.74993385]])
In [15]: np.dot(arr.T, arr)
Out[15]:
array([[ 3.32192256, 0.64291292, -0.33265926],
[ 0.64291292, 6.49435244, -2.22525257],
[-0.33265926, -2.22525257, 3.92130431]])
對於高維陣列,transpose需要得到一個由軸編號組成的元組才能對這些軸進行轉置(比較費腦子):
In [16]: arr = np.arange(16).reshape((2, 2, 4))
In [17]: arr
Out[17]:
array([[[ 0, 1, 2, 3],
[ 4, 5, 6, 7]],
[[ 8, 9, 10, 11],
[12, 13, 14, 15]]])
In [18]: arr.transpose((1, 0, 2))
Out[18]:
array([[[ 0, 1, 2, 3],
[ 8, 9, 10, 11]],
[[ 4, 5, 6, 7],
[12, 13, 14, 15]]])
這裡,第一個軸被換成了第二個,第二個軸被換成了第一個,最後一個軸不變。
簡單的轉置可以使用.T,它其實就是進行軸對換而已。ndarray還有一個swapaxes方法,它需要接受一對軸編號:
In [19]: arr
Out[19]:
array([[[ 0, 1, 2, 3],
[ 4, 5, 6, 7]],
[[ 8, 9, 10, 11],
[12, 13, 14, 15]]])
In [20]: arr.swapaxes(1, 2)
Out[20]:
array([[[ 0, 4],
[ 1, 5],
[ 2, 6],
[ 3, 7]],
[[ 8, 12],
[ 9, 13],
[10, 14],
[11, 15]]])
swapaxes也是返回源資料的檢視(不會進行任何複製操作)。
4.2 通用函式:快速的元素級陣列函式
通用函式(即ufunc)是一種對ndarray中的資料執行元素級運算的函式。你可以將其看做簡單函式(接受一個或多個標量值,併產生一個或多個標量值)的向量化包裝器。
許多ufunc都是簡單的元素級變體,如sqrt和exp:
In [21]: arr = np.arange(10)
In [22]: arr
Out[22]:
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
In [23]: np.sqrt(arr)
Out[23]:
array([0. , 1. , 1.41421356, 1.73205081, 2. ,
2.23606798, 2.44948974, 2.64575131, 2.82842712, 3. ])
In [24]: np.exp(arr)
Out[24]:
array([1.00000000e+00, 2.71828183e+00, 7.38905610e+00, 2.00855369e+01,
5.45981500e+01, 1.48413159e+02, 4.03428793e+02, 1.09663316e+03,
2.98095799e+03, 8.10308393e+03])
這些都是一元(unary)ufunc。另外一些(如add或maximum)接受2個數組(因此也叫二元(binary)ufunc),並返回一個結果陣列:
In [25]: x = np.random.randn(8)
In [26]: y = np.random.randn(8)
In [27]: x
Out[27]:
array([-0.95361393, -1.61004574, 0.15676379, -0.21776008, 0.69062465,
-1.52815122, -0.04371557, 2.01002837])
In [28]: y
Out[28]:
array([ 0.41073896, -0.67278923, -0.15690505, -0.71940875, -1.68803873,
1.20356913, 0.13458085, -1.6300652 ])
In [30]: np.maximum(x, y)
Out[30]:
array([ 0.41073896, -0.67278923, 0.15676379, -0.21776008, 0.69062465,
1.20356913, 0.13458085, 2.01002837])
這裡,numpy.maximum計算了x和y中元素級別最大的元素。
雖然並不常見,但有些ufunc的確可以返回多個數組。modf就是一個例子,它是Python內建函式divmod的向量化版本,它會返回浮點數陣列的小數和整數部分:
In [31]: arr = np.random.randn(7) * 5
In [32]: arr
Out[32]:
array([-4.4759812 , -0.69628583, -9.28954211, 0.97787069, -9.80958373,
3.50748262, 1.95993912])
In [33]: remainder, whole_part = np.modf(arr)
In [34]: remainder
Out[34]:
array([-0.4759812 , -0.69628583, -0.28954211, 0.97787069, -0.80958373,
0.50748262, 0.95993912])
In [35]: whole_part
Out[35]:
array([-4., -0., -9., 0., -9., 3., 1.])
Ufuncs可以接受一個out可選引數,這樣就能在陣列原地進行操作:
print(arr)
print(np.sqrt(arr))
print(np.sqrt(arr, arr))
print(arr)
output
[ -2.88452855 -1.14582772 3.00420715 1.22281642 -0.38185437
-11.36241841 6.52291062]
[ nan nan 1.73326488 1.1058103 nan nan
2.55399895]
[ nan nan 1.73326488 1.1058103 nan nan
2.55399895]
[ nan nan 1.73326488 1.1058103 nan nan
2.55399895]
表4-3和表4-4分別列出了一些一元和二元ufunc。
4.3 利用陣列進行資料處理
NumPy陣列使你可以將許多種資料處理任務表述為簡潔的陣列表示式(否則需要編寫迴圈)。用陣列表示式代替迴圈的做法,通常被稱為向量化。
作為簡單的例子,假設我們想要在一組值(網格型)上計算函式sqrt(x^2+y^2)
。np.meshgrid函式接受兩個一維陣列,併產生兩個二維矩陣(對應於兩個陣列中所有的(x,y)對):
In [40]: points = np.arange(-5, 5, 0.01) #步長0.01 ,總計1000個離散點
In [41]: xs, ys = np.meshgrid(points, points)
In [42]: ys
Out[42]:
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 [45]: z = np.sqrt(xs ** 2 + ys ** 2)
In [46]: z
Out[46]:
array([[7.07106781, 7.06400028, 7.05693985, ..., 7.04988652, 7.05693985,
7.06400028],
[7.06400028, 7.05692568, 7.04985815, ..., 7.04279774, 7.04985815,
7.05692568],
[7.05693985, 7.04985815, 7.04278354, ..., 7.03571603, 7.04278354,
7.04985815],
...,
[7.04988652, 7.04279774, 7.03571603, ..., 7.0286414 , 7.03571603,
7.04279774],
[7.05693985, 7.04985815, 7.04278354, ..., 7.03571603, 7.04278354,
7.04985815],
[7.06400028, 7.05692568, 7.04985815, ..., 7.04279774, 7.04985815,
7.05692568]])
採用matplotlib建立了這個二維陣列的視覺化:
In [47]: import matplotlib.pyplot as plt
In [48]: plt.imshow(z, cmap = plt.cm.gray); plt.colorbar()
Out[48]: <matplotlib.colorbar.Colorbar at 0x16e18e7dc88>
In [49]: plt.title('Image plot of $\sqrt{x^2 + y^2}$ for a grid of values')
Out[49]: Text(0.5,1,'Image plot of $\\sqrt{x^2 + y^2}$ for a grid of values')
這張圖是用matplotlib的imshow函式建立的。
4.3.1 將條件邏輯表述為陣列運算
numpy.where函式是三元表示式x if condition else y的向量化版本。假設我們有一個布林陣列和兩個值陣列:
In [57]: xarr = np.array([1.1, 1.2, 1.3, 1.4, 1.5])
In [58]: yarr = np.array([2.1, 2.2, 2.3, 2.4, 2.5])
In [59]: cond = np.array([True, False, True, True, False])
假設我們想要根據cond中的值選取xarr和yarr的值:當cond中的值為True時,選取xarr的值,否則從yarr中選取。列表推導式的寫法應該如下所示:
In [63]: result = [(x if c else y)
for x, y, c in zip(xarr, yarr, cond)]
In [64]: result
Out[64]:
[1.1, 2.2, 1.3, 1.4, 2.5]
這有幾個問題。第一,它對大陣列的處理速度不是很快(因為所有工作都是由純Python完成的)。第二,無法用於多維陣列。若使用np.where,則可以將該功能寫得非常簡潔:
In [65]: result = np.where(cond, xarr, yarr)
In [66]: result
Out[66]:
array([1.1, 2.2, 1.3, 1.4, 2.5])
np.where的第二個和第三個引數不必是陣列,它們都可以是標量值。在資料分析工作中,where通常用於根據另一個數組而產生一個新的陣列。假設有一個由隨機資料組成的矩陣,你希望將所有正值替換為2,將所有負值替換為-2。若利用np.where,則會非常簡單:
In [67]: arr = np.random.randn(4, 4)
In [68]: arr
Out[68]:
array([[-0.64496779, 0.52414139, 1.18606729, 1.24215752],
[ 0.14011134, 0.68157075, -1.17944087, 0.77633086],
[-1.52888488, 0.74012192, 0.67633297, 1.07999609],
[-1.85749656, -0.98081351, 0.56145716, -0.19420933]])
In [69]: arr > 0
Out[69]:
array([[False, True, True, True],
[ True, True, False, True],
[False, True, True, True],
[False, False, True, False]])
In [70]: np.where(arr > 0, 2, -2 )
Out[70]:
array([[-2, 2, 2, 2],
[ 2, 2, -2, 2],
[-2, 2, 2, 2],
[-2, -2, 2, -2]])
使用np.where,可以將標量和陣列結合起來。例如,我可用常數2替換arr中所有正的值:
In [71]: np.where(arr > 0, 2, arr)
Out[71]:
array([[-0.64496779, 2. , 2. , 2. ],
[ 2. , 2. , -1.17944087, 2. ],
[-1.52888488, 2. , 2. , 2. ],
[-1.85749656, -0.98081351, 2. , -0.19420933]])
傳遞給where的陣列大小可以不相等,甚至可以是標量值。
4.3.2 數學和統計方法
可以通過陣列上的一組數學函式對整個陣列或某個軸向的資料進行統計計算。sum、mean以及標準差std等聚合計算(aggregation,通常叫做約簡(reduction))既可以當做陣列的例項方法呼叫,也可以當做頂級NumPy函式使用。
這裡,可以生成一些正態分佈隨機資料,然後做了聚類統計:
In [72]: arr = np.random.randn(5, 4)
In [73]: arr
Out[73]:
array([[-0.18969458, 2.42217634, -0.11111327, -0.57700033],
[-1.01329696, 0.00691129, -1.50832842, 0.24767229],
[-1.41086114, -1.17894482, -0.71458589, 0.18638283],
[-0.02377788, 0.83628221, -0.98033713, 1.26823669],
[ 1.18835246, 0.37427623, 0.57923965, -1.60212904]])
In [74]: arr.mean()
Out[74]: -0.11002697365467654
In [75]: np.mean(arr)
Out[75]: -0.11002697365467654
In [76]: arr.sum()
Out[76]: -2.2005394730935306
mean和sum這類的函式可以接受一個axis選項引數,用於計算該軸向上的統計值,最終結果是一個少一維的陣列:
In [77]: arr.mean(axis = 1)
Out[77]: array([ 0.38609204, -0.56676045, -0.77950226, 0.27510097, 0.13493483])
In [78]: arr.sum(axis = 0)
Out[78]: array([-1.4492781 , 2.46070125, -2.73512506, -0.47683756])
這裡,arr.mean(1)是“計算行的平均值”,arr.sum(0)是“計算每列的和”。
其他如cumsum和cumprod之類的方法則不聚合,而是產生一個由中間結果組成的陣列:
In [79]: arr = np.array([1, 2, 3, 4, 5, 6, 7, 8])
In [80]: arr.cumsum()
Out[80]: array([ 1, 3, 6, 10, 15, 21, 28, 36], dtype=int32)
在多維陣列中,累加函式(如cumsum)返回的是同樣大小的陣列,但是會根據每個低維的切片沿著標記軸計算部分聚類:
In [81]: arr = np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])
In [82]: arr
Out[82]:
array([[0, 1, 2],
[3, 4, 5],
[6, 7, 8]])
In [83]: arr.cumsum(axis = 0)
Out[83]:
array([[ 0, 1, 2],
[ 3, 5, 7],
[ 9, 12, 15]], dtype=int32)
In [84]: arr.cumprod(axis = 1)
Out[84]:
array([[ 0, 0, 0],
[ 3, 12, 60],
[ 6, 42, 336]], dtype=int32)
表4-5列出了全部的基本陣列統計方法。後續章節中有很多例子都會用到這些方法。
4.3.3 用於布林型陣列的方法
在上面這些方法中,布林值會被強制轉換為1(True)和0(False)。因此,sum經常被用來對布林型陣列中的True值計數:
In [85]: arr = np.random.randn(100)
In [86]: (arr > 0).sum()
Out[86]: 40
另外還有兩個方法any和all,它們對布林型陣列非常有用。any用於測試陣列中是否存在一個或多個True,而all則檢查陣列中所有值是否都是True:
In [87]: bools = np.array([False, False, True, False])
In [88]: bools.any()
Out[88]: True
In [89]: bools.all()
Out[89]: False
這兩個方法也能用於非布林型陣列,所有非0元素將會被當做True。
4.3.4 排序
跟Python內建的列表型別一樣,NumPy陣列也可以通過sort方法就地排序:
In [99]: arr = np.random.randn(6)
In [100]: arr
Out[100]:
array([ 0.89934768, 0.49889372, -1.80138268, 1.35137542, -0.29677817,
-1.04741995])
In [101]: arr.sort()
In [102]: arr
Out[102]:
array([-1.80138268, -1.04741995, -0.29677817, 0.49889372, 0.89934768,
1.35137542])
多維陣列可以在任何一個軸向上進行排序,只需將軸編號傳給sort即可:
In [95]: arr = np.random.randn(5, 3)
In [96]: arr
Out[96]:
array([[-0.44073589, -0.84379604, -1.74775132],
[-1.17503284, 0.54021733, -0.18219635],
[-0.22771804, 0.98648339, -0.71349615],
[ 0.5721064 , 0.38012503, 0.67521682],
[ 0.49641397, -0.50484124, -2.19343809]])
In [97]: arr.sort(1)
In [98]: print(arr)
[[-1.74775132 -0.84379604 -0.44073589]
[-1.17503284 -0.18219635 0.54021733]
[-0.71349615 -0.22771804 0.98648339]
[ 0.38012503 0.5721064 0.67521682]
[-2.19343809 -0.50484124 0.49641397]]
頂級方法np.sort返回的是陣列的已排序副本,而就地排序則會修改陣列本身。計算陣列分位數最簡單的辦法是對其進行排序,然後選取特定位置的值:
In [103]: large_arr = np.random.randn(1000)
In [104]: large_arr.sort()
In [105]: large_arr[int(0.05 * len(large_arr))]
Out[105]: -1.6334705078901595
NumPy中提供了各種排序相關功能。 這些排序函式實現不同的排序演算法,每個排序演算法的特徵在於執行速度,最壞情況效能,所需的工作空間和演算法的穩定性。 下表顯示了三種排序演算法的比較。
種類 | 速度 | 最壞情況 | 工作空間 | 穩定性 |
---|---|---|---|---|
'quicksort'(快速排序) | 1 | O(n^2) | 0 | 否 |
'mergesort'(歸併排序) | 2 | O(n*log(n)) | ~n/2 | 是 |
'heapsort'(堆排序) | 3 | O(n*log(n)) | 0 | 否 |
numpy.sort()
sort()函式返回輸入陣列的排序副本。 它有以下引數:
numpy.sort(a, axis, kind, order)
其中:
序號 | 引數及描述 |
---|---|
1 | a 要排序的陣列 |
2 | axis 沿著它排序陣列的軸,如果沒有陣列會被展開,沿著最後的軸排序 |
3 | kind 預設為'quicksort'(快速排序) |
4 | order 如果陣列包含欄位,則是要排序的欄位 |
示例:
In [90]: import numpy as np
a = np.array([[3, 7], [9, 2]])
In [94]:
print('我們的陣列:')
print(a)
print('\n')
print('呼叫sort()函式:')
print(np.sort(a))
print('\n')
print('沿軸 0 排序:')
print(np.sort(a, axis = 0))
print('\n')
# 對欄位排序
dt = np.dtype([('name', 'S10'), ('age', int)])
a = np.array([('raju',21),('anil', 25),('ravi', 17), ('amar', 27)], dtype = dt)
print('我們的陣列:')
print(a)
print('\n')
print('按 name 排序')
print(np.sort(a, order = 'name'))
output:
我們的陣列:
[[3 7]
[9 2]]
呼叫sort()函式:
[[3 7]
[2 9]]
沿軸 0 排序:
[[3 2]
[9 7]]
我們的陣列:
[(b'raju', 21) (b'anil', 25) (b'ravi', 17) (b'amar', 27)]
按 name 排序
[(b'amar', 27) (b'anil', 25) (b'raju', 21) (b'ravi', 17)]
numpy.argsort()
numpy.argsort()函式對輸入陣列沿給定軸執行間接排序,並使用指定排序型別返回資料的索引陣列。 這個索引陣列用於構造排序後的陣列。
In [6]: import numpy as np
In [7]: x = np.array([3, 1, 2])
In [9]:
print('我們的陣列是:')
print(x)
print('\n')
print('對 x 呼叫 argsort()函式:')
y = np.argsort(x)
print(y)
print('\n')
print('以排序後的順序重構原陣列')
print(x[y])
print('\n')
print('使用迴圈重構原陣列')
for i in y:
print (x[i])
putput:
我們的陣列是:
[3 1 2]
對 x 呼叫 argsort()函式:
[1 2 0]
以排序後的順序重構原陣列
[1 2 3]
使用迴圈重構原陣列
1
2
3
numpy.lexsort()
函式使用鍵序列執行間接排序。 鍵可以看作是電子表格中的一列。 該函式返回一個索引陣列,使用它可以獲得排序資料。 注意,最後一個鍵恰好是 sort 的主鍵。
import numpy as np
nm = ('raju','anil','ravi','amar')
dv = ('f.y.', 's.y.', 's.y.', 'f.y.')
ind = np.lexsort((dv,nm))
print ('呼叫 lexsort() 函式:' )
print (ind)
print ('\n')
print ('使用這個索引來獲取排序後的資料:')
print ([nm[i] + ", " + dv[i] for i in ind])
output
呼叫 lexsort() 函式:
[3 1 0 2]
使用這個索引來獲取排序後的資料:
['amar, f.y.', 'anil, s.y.', 'raju, f.y.', 'ravi, s.y.']
4.3.5 唯一化以及其它的集合邏輯
NumPy提供了一些針對一維ndarray的基本集合運算。最常用的可能要數np.unique了,它用於找出陣列中的唯一值並返回已排序的結果:
In [11]: names = np.array(['Bob', 'Joe', 'Will', 'Bob', "Will", 'Joe', "Joe"])
In [12]: np.unique(names)
Out[12]: array(['Bob', 'Joe', 'Will'], dtype='<U4')
In [13]: ints = np.array([3, 3, 3, 2, 2, 1, 1, 4, 4])
In [14]: np.unique(ints)
Out[14]: array([1, 2, 3, 4])
拿跟np.unique等價的純Python程式碼來對比一下:
In [15]: sorted(set(names))
Out[15]: ['Bob', 'Joe', 'Will']
另一個函式np.in1d用於測試一個數組中的值在另一個數組中的成員資格,返回一個布林型陣列:
In [16]: values = np.array([6, 0, 0, 3, 2, 5, 6])
In [18]: np.in1d(values, [2, 3, 6])
Out[18]: array([ True, False, False, True, True, False, True])
NumPy中的集合函式如下表所示:
4.4 用於陣列的檔案輸入輸出
NumPy能夠讀寫磁碟上的文字資料或二進位制資料。
np.save和np.load是讀寫磁碟陣列資料的兩個主要函式。預設情況下,陣列是以未壓縮的原始二進位制格式儲存在副檔名為.npy的檔案中的:
In [19]: arr = np.arange(10)
In [20]: np.save('some_array', arr)
In [21]: print(np.load('some_array.npy'))
out[22]: [0 1 2 3 4 5 6 7 8 9]
如果檔案路徑末尾沒有副檔名.npy,則該副檔名會被自動加上。然後就可以通過np.load讀取磁碟上的陣列。
通過np.savez可以將多個數組儲存到一個未壓縮檔案中,將陣列以關鍵字引數的形式傳入即可:
In [21]: np.savez('array_archive.npz', a = arr, b = arr)
載入.npz檔案時,你會得到一個類似字典的物件,該物件會對各個陣列進行延遲載入:
In [22]: arch = np.load('array_archive.npz')
In [23]: arch['b']
Out[23]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
如果要將資料壓縮,可以使用numpy.savez_compressed:
In [26]: np.savez_compressed('arrays_compressed.npz', a = arr, b = arr)
4.5 線性代數
線性代數(如矩陣乘法、矩陣分解、行列式以及其他方陣數學等)是任何陣列庫的重要組成部分。不像某些語言(如MATLAB),通過*對兩個二維陣列相乘得到的是一個元素級的積,而不是一個矩陣點積。因此,NumPy提供了一個用於矩陣乘法的dot函式(既是一個數組方法也是numpy名稱空間中的一個函式):
In [27]: x = np.array([[1., 2., 3.],[4., 5., 6]])
In [28]: y = np.array([[6., 23.],[-1, 7],[8, 9]])
In [29]: x
Out[29]:
array([[1., 2., 3.],
[4., 5., 6.]])
In [30]: y
Out[30]:
array([[ 6., 23.],
[-1., 7.],
[ 8., 9.]])
In [31]: x.dot(y)
Out[31]:
array([[ 28., 64.],
[ 67., 181.]])
x.dot(y)等價於np.dot(x, y):
In [32]: np.dot(x,y)
Out[32]:
array([[ 28., 64.],
[ 67., 181.]])
一個二維陣列跟一個大小合適的一維陣列的矩陣點積運算之後將會得到一個一維陣列:
In [33]: np.dot(x, np.ones(3))
Out[33]: array([ 6., 15.])
@符(類似Python 3.5)也可以用作中綴運算子,進行矩陣乘法:
In [34]: x @ np.ones(3)
Out[34]: array([ 6., 15.])
numpy.linalg中有一組標準的矩陣分解運算以及諸如求逆和行列式之類的東西。它們跟MATLAB和R等語言所使用的是相同的行業標準線性代數庫,如BLAS、LAPACK、Intel MKL(Math Kernel Library,可能有,取決於你的NumPy版本)等:
In [36]: from numpy.linalg import inv, qr
In [37]: x = np.random.rand(5, 5)
In [38]: mat = x.T.dot(x)
In [39]: inv(mat)
Out[39]:
array([[ 24.25773408, 3.09904598, -10.49682692, -23.85575513,
7.05938574],
[ 3.09904598, 2.63122057, -2.15738514, -2.84277908,
-0.4158934 ],
[-10.49682692, -2.15738514, 6.91840105, 9.17259942,
-4.10985733],
[-23.85575513, -2.84277908, 9.17259942, 25.96915031,
-7.82584247],
[ 7.05938574, -0.4158934 , -4.10985733, -7.82584247,
6.09047827]])
In [40]: mat.dot(inv(mat))
Out[40]:
array([[ 1.00000000e+00, 1.08490052e-15, 8.63057225e-17,
1.50049454e-16, -8.01708725e-16],
[ 4.76736930e-15, 1.00000000e+00, -7.01519917e-17,
-5.44503590e-15, 3.71342593e-16],
[ 5.85014959e-15, 1.23737453e-15, 1.00000000e+00,
-3.13033967e-15, 9.67269400e-16],
[-3.18522853e-15, 6.30828711e-16, 3.05358331e-15,
1.00000000e+00, -4.94376604e-16],
[ 2.91234983e-15, 1.16323349e-15, -6.18529051e-16,
-2.43961692e-15, 1.00000000e+00]])
In [41]: q, r = qr(mat)
In [42]: r
Out[42]:
array([[-3.1693019 , -2.59359093, -3.57833789, -2.67421668, -2.40708289],
[ 0. , -0.88414288, -0.52716641, -0.06335676, -0.54377552],
[ 0. , 0. , -0.28652111, 0.08813383, -0.10099342],
[ 0. , 0. , 0. , -0.08465008, -0.23351305],
[ 0. , 0. , 0. , 0. , 0.07779415]])
表示式X.T.dot(X)計算X和它的轉置X.T的點積。
下表列出了一些最常用的線性代數函式:
4.6 偽隨機數生成
numpy.random模組對Python內建的random進行了補充,增加了一些用於高效生成多種概率分佈的樣本值的函式。例如,你可以用normal來得到一個標準正態分佈的4×4樣本陣列:
In [43]: samples = np.random.normal(size =(4, 4))
In [44]: samples
Out[44]:
array([[ 0.26031823, 0.98180027, -0.12338395, -1.66632337],
[ 0.71714765, 0.3356353 , 0.40195087, 0.24449305],
[ 1.55741521, 0.73821615, -0.635539 , 1.01435757],
[ 0.97447286, 0.32439164, 0.12732344, -0.36458816]])
而Python內建的random模組則只能一次生成一個樣本值。從下面的測試結果中可以看出,如果需要產生大量樣本值,numpy.random快了不止一個數量級:
In [45]: from random import normalvariate
In [46]: N = 1000000
In [47]: %timeit smples = [normalvariate(0, 1) for _ in range(N)]
617 ms ± 1.57 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
In [49]: %timeit np.random.normal(size = N)
25 ms ± 209 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
我們說這些都是偽隨機數,是因為它們都是通過演算法基於隨機數生成器種子,在確定性的條件下生成的。你可以用NumPy的np.random.seed更改隨機數生成種子:
In [50]: np.random.seed(1234)
numpy.random的資料生成函式使用了全域性的隨機種子。要避免全域性狀態,你可以使用numpy.random.RandomState,建立一個與其它隔離的隨機數生成器:
In [51]: rng = np.random.RandomState(1234)
In [52]: rng.randn(10)
Out[52]:
array([ 0.47143516, -1.19097569, 1.43270697, -0.3126519 , -0.72058873,
0.88716294, 0.85958841, -0.6365235 , 0.01569637, -2.24268495])
下表列出了numpy.random中的部分函式。在下一節中,我將給出一些利用這些函式一次性生成大量樣本值的範例。
4.7 隨機漫步
4.7.1 單個隨機漫步
我們通過模擬隨機漫步來說明如何運用陣列運算。先來看一個簡單的隨機漫步的例子:從0開始,步長1和-1出現的概率相等。
下面是一個通過內建的random模組以純Python的方式實現1000步的隨機漫步:
In [56]:
import random
import matplotlib.pyplot as plt
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)
plt .plot(walk[:100])
不難看出,這其實就是隨機漫步中各步的累計和,可以用一個數組運算來實現。因此,我用np.random模組一次性隨機產生1000個“擲硬幣”結果(即兩個數中任選一個),將其分別設定為1或-1,然後計算累計和:
In [57]: nsteps = 1000
In [58]: draws = np.random.randint(0, 2, size = nsteps)
In [59]: steps = np.where(draws > 0, 1, -1)
In [60]: walk = steps.cumsum()
有了這些資料之後,我們就可以沿著漫步路徑做一些統計工作了,比如求取最大值和最小值:
In [61]: walk.min()
Out[61]: -9
In [62]: walk.max()
Out[62]: 60
現在來看一個複雜點的統計任務——首次穿越時間,即隨機漫步過程中第一次到達某個特定值的時間。假設我們想要知道本次隨機漫步需要多久才能距離初始0點至少10步遠(任一方向均可)。np.abs(walk)>=10可以得到一個布林型陣列,它表示的是距離是否達到或超過10,而我們想要知道的是第一個10或-10的索引。可以用argmax來解決這個問題,它返回的是該布林型陣列第一個最大值的索引(True就是最大值):
In [63]: (np.abs(walk) >= 10).argmax()
Out[63]: 297
注意,這裡使用argmax並不是很高效,因為它無論如何都會對陣列進行完全掃描。在本例中,只要發現了一個True,那我們就知道它是個最大值了。
4.7.2 一次模擬多個隨機漫步
如果你希望模擬多個隨機漫步過程(比如5000個),只需對上面的程式碼做一點點修改即可生成所有的隨機漫步過程。只要給numpy.random的函式傳入一個二元元組就可以產生一個二維陣列,然後我們就可以一次性計算5000個隨機漫步過程(一行一個)的累計和了:
In [64]: nwalks = 5000
In [65]: nsteps = 1000
In [66]: draws = np.random.randint(0, 2, size = (nwalks, nsteps)) # 0 或者 1
In [68]: steps = np.where(draws > 0, 1, -1)
In [69]: walks = steps.cumsum(1)
In [70]: walks
Out[70]:
array([[ 1, 2, 3, ..., 46, 47, 46],
[ 1, 0, 1, ..., 40, 41, 42],
[ 1, 2, 3, ..., -26, -27, -28],
...,
[ 1, 0, 1, ..., 64, 65, 66],
[ 1, 2, 1, ..., 2, 1, 0],
[ -1, -2, -3, ..., 32, 33, 34]], dtype=int32)
現在,我們來計算所有隨機漫步過程的最大值和最小值:
In [71]: walks.max()
Out[71]: 122
In [72]: walks.min()
Out[72]: -128
得到這些資料之後,我們來計算30或-30的最小穿越時間。這裡稍微複雜些,因為不是5000個過程都到達了30。我們可以用any方法來對此進行檢查:
In [73]: hits30 = (np.abs(walks) >= 30).any(1)
In [74]: hits30
Out[74]: array