python數據分析 Numpy基礎 數組和矢量計算


NumPy(Numerical Python的簡稱)是Python數值計算最重要的基礎包。大多數提供科學計算的包都是用NumPy的數組作為構建基礎。

NumPy的部分功能如下:

  • ndarray,一個具有矢量算術運算和復雜廣播能力的快速且節省空間的多維數組。
  • 用於對整組數據進行快速運算的標准數學函數(無需編寫循環)。
  • 用於讀寫磁盤數據的工具以及用於操作內存映射文件的工具。
  • 線性代數、隨機數生成以及傅里葉變換功能。
  • 用於集成由C、C++、Fortran等語言編寫的代碼的A C API。

由於NumPy提供了一個簡單易用的C API,因此很容易將數據傳遞給由低級語言編寫的外部庫,外部庫也能以NumPy數組的形式將數據返回給Python。這個功能使Python成為一種包裝C/C++/Fortran歷史代碼庫的選擇,並使被包裝庫擁有一個動態的、易用的接口。

NumPy本身並沒有提供多么高級的數據分析功能,理解NumPy數組以及面向數組的計算將有助於你更加高效地使用諸如pandas之類的工具。因為NumPy是一個很大的題目,我會在附錄A中介紹更多NumPy高級功能,比如廣播。

對於大部分數據分析應用而言,我最關注的功能主要集中在:

  • 用於數據整理和清理、子集構造和過濾、轉換等快速的矢量化數組運算。
  • 常用的數組算法,如排序、唯一化、集合運算等。
  • 高效的描述統計和數據聚合/摘要運算。
  • 用於異構數據集的合並/連接運算的數據對齊和關系型數據運算。
  • 將條件邏輯表述為數組表達式(而不是帶有if-elif-else分支的循環)。
  • 數據的分組運算(聚合、轉換、函數應用等)。。

雖然NumPy提供了通用的數值數據處理的計算基礎,但大多數讀者可能還是想將pandas作為統計和分析工作的基礎,尤其是處理表格數據時。pandas還提供了一些NumPy所沒有的領域特定的功能,如時間序列處理等。

筆記:Python的面向數組計算可以追溯到1995年,Jim Hugunin創建了Numeric庫。接下來的10年,許多科學編程社區紛紛開始使用Python的數組編程,但是進入21世紀,庫的生態系統變得碎片化了。2005年,Travis Oliphant從Numeric和Numarray項目整合出了NumPy項目,進而所有社區都集合到了這個框架下。

NumPy之於數值計算特別重要的原因之一,是因為它可以高效處理大數組的數據。這是因為:

  • NumPy是在一個連續的內存塊中存儲數據,獨立於其他Python內置對象。NumPy的C語言編寫的算法庫可以操作內存,而不必進行類型檢查或其它前期工作。比起Python的內置序列,NumPy數組使用的內存更少。
  • NumPy可以在整個數組上執行復雜的計算,而不需要Python的for循環。

要搞明白具體的性能差距,考察一個包含一百萬整數的數組,和一個等價的Python列表:

In [7]: import numpy as np In [8]: my_arr = np.arange(1000000) In [9]: my_list = list(range(1000000))

各個序列分別乘以2:

In [10]: %time for _ in range(10): my_arr2 = my_arr * 2 CPU times: user 20 ms, sys: 50 ms, total: 70 ms Wall time: 72.4 ms In [11]: %time for _ in range(10): my_list2 = [x * 2 for x in my_list] CPU times: user 760 ms, sys: 290 ms, total: 1.05 s Wall time: 1.05 s

基於NumPy的算法要比純Python快10到100倍(甚至更快),並且使用的內存更少。

 

4.1 NumPy的ndarray:一種多維數組對象

NumPy最重要的一個特點就是其N維數組對象(即ndarray),該對象是一個快速而靈活的大數據集容器。你可以利用這種數組對整塊數據執行一些數學運算,其語法跟標量元素之間的運算一樣。

要明白Python是如何利用與標量值類似的語法進行批次計算,我先引入NumPy,然后生成一個包含隨機數據的小數組:

In [12]: import numpy as np # Generate some random data In [13]: data = np.random.randn(2, 3) In [14]: data Out[14]: array([[-0.2047, 0.4789, -0.5194], [-0.5557, 1.9658, 1.3934]])

然后進行數學運算:

In [15]: data * 10 Out[15]: array([[ -2.0471, 4.7894, -5.1944], [ -5.5573, 19.6578, 13.9341]]) In [16]: data + data Out[16]: array([[-0.4094, 0.9579, -1.0389], [-1.1115, 3.9316, 2.7868]])

第一個例子中,所有的元素都乘以10。第二個例子中,每個元素都與自身相加。

筆記:在本章及全書中,我會使用標准的NumPy慣用法import numpy as np。你當然也可以在代碼中使用from numpy import *,但不建議這么做。numpy的命名空間很大,包含許多函數,其中一些的名字與Python的內置函數重名(比如min和max)。

ndarray是一個通用的同構數據多維容器,也就是說,其中的所有元素必須是相同類型的。每個數組都有一個shape(一個表示各維度大小的元組)和一個dtype(一個用於說明數組數據類型的對象):

In [17]: data.shape Out[17]: (2, 3) In [18]: data.dtype Out[18]: dtype('float64')

本章將會介紹NumPy數組的基本用法,這對於本書后面各章的理解基本夠用。雖然大多數數據分析工作不需要深入理解NumPy,但是精通面向數組的編程和思維方式是成為Python科學計算牛人的一大關鍵步驟。

筆記:當你在本書中看到“數組”、“NumPy數組”、"ndarray"時,基本上都指的是同一樣東西,即ndarray對象。

 

創建ndarray

創建數組最簡單的辦法就是使用array函數。它接受一切序列型的對象(包括其他數組),然后產生一個新的含有傳入數據的NumPy數組。以一個列表的轉換為例:

In [19]: data1 = [6, 7.5, 8, 0, 1] In [20]: arr1 = np.array(data1) In [21]: arr1 Out[21]: array([ 6. , 7.5, 8. , 0. , 1. ])

嵌套序列(比如由一組等長列表組成的列表)將會被轉換為一個多維數組:

In [22]: data2 = [[1, 2, 3, 4], [5, 6, 7, 8]] In [23]: arr2 = np.array(data2) In [24]: arr2 Out[24]: array([[1, 2, 3, 4], [5, 6, 7, 8]])

因為data2是列表的列表,NumPy數組arr2的兩個維度的shape是從data2引入的。可以用屬性ndim和shape驗證:

In [25]: arr2.ndim Out[25]: 2 In [26]: arr2.shape Out[26]: (2, 4)

除非特別說明(稍后將會詳細介紹),np.array會嘗試為新建的這個數組推斷出一個較為合適的數據類型。數據類型保存在一個特殊的dtype對象中。比如說,在上面的兩個例子中,我們有:

In [27]: arr1.dtype Out[27]: dtype('float64') In [28]: arr2.dtype Out[28]: dtype('int64')

除np.array之外,還有一些函數也可以新建數組。比如,zeros和ones分別可以創建指定長度或形狀的全0或全1數組。empty可以創建一個沒有任何具體值的數組。要用這些方法創建多維數組,只需傳入一個表示形狀的元組即可:

In [29]: np.zeros(10) Out[29]: array([ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]) In [30]: np.zeros((3, 6)) Out[30]: array([[ 0., 0., 0., 0., 0., 0.], [ 0., 0., 0., 0., 0., 0.], [ 0., 0., 0., 0., 0., 0.]]) In [31]: np.empty((2, 3, 2)) Out[31]: array([[[ 0., 0.], [ 0., 0.], [ 0., 0.]], [[ 0., 0.], [ 0., 0.], [ 0., 0.]]])

注意:認為np.empty會返回全0數組的想法是不安全的。很多情況下(如前所示),它返回的都是一些未初始化的垃圾值。

arange是Python內置函數range的數組版:

In [32]: np.arange(15) Out[32]: array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14])

表4-1列出了一些數組創建函數。由於NumPy關注的是數值計算,因此,如果沒有特別指定,數據類型基本都是float64(浮點數)。

表4-1 數組創建函數

 

ndarray的數據類型

dtype(數據類型)是一個特殊的對象,它含有ndarray將一塊內存解釋為特定數據類型所需的信息:

In [33]: arr1 = np.array([1, 2, 3], dtype=np.float64) In [34]: arr2 = np.array([1, 2, 3], dtype=np.int32) In [35]: arr1.dtype Out[35]: dtype('float64') In [36]: arr2.dtype Out[36]: dtype('int32')

dtype是NumPy靈活交互其它系統的源泉之一。多數情況下,它們直接映射到相應的機器表示,這使得“讀寫磁盤上的二進制數據流”以及“集成低級語言代碼(如C、Fortran)”等工作變得更加簡單。數值型dtype的命名方式相同:一個類型名(如float或int),后面跟一個用於表示各元素位長的數字。標准的雙精度浮點值(即Python中的float對象)需要占用8字節(即64位)。因此,該類型在NumPy中就記作float64。表4-2列出了NumPy所支持的全部數據類型。

筆記:記不住這些NumPy的dtype也沒關系,新手更是如此。通常只需要知道你所處理的數據的大致類型是浮點數、復數、整數、布爾值、字符串,還是普通的Python對象即可。當你需要控制數據在內存和磁盤中的存儲方式時(尤其是對大數據集),那就得了解如何控制存儲類型。

你可以通過ndarray的astype方法明確地將一個數組從一個dtype轉換成另一個dtype:

In [37]: arr = np.array([1, 2, 3, 4, 5]) In [38]: arr.dtype Out[38]: dtype('int64') In [39]: float_arr = arr.astype(np.float64) In [40]: float_arr.dtype Out[40]: dtype('float64')

在本例中,整數被轉換成了浮點數。如果將浮點數轉換成整數,則小數部分將會被截取刪除:

In [41]: arr = np.array([3.7, -1.2, -2.6, 0.5, 12.9, 10.1]) In [42]: arr Out[42]: array([ 3.7, -1.2, -2.6, 0.5, 12.9, 10.1]) In [43]: arr.astype(np.int32) Out[43]: array([ 3, -1, -2, 0, 12, 10], dtype=int32)

如果某字符串數組表示的全是數字,也可以用astype將其轉換為數值形式:

In [44]: numeric_strings = np.array(['1.25', '-9.6', '42'], dtype=np.string_) In [45]: numeric_strings.astype(float) Out[45]: array([ 1.25, -9.6 , 42. ])

注意:使用numpy.string_類型時,一定要小心,因為NumPy的字符串數據是大小固定的,發生截取時,不會發出警告。pandas提供了更多非數值數據的便利的處理方法。

如果轉換過程因為某種原因而失敗了(比如某個不能被轉換為float64的字符串),就會引發一個ValueError。這里,我比較懶,寫的是float而不是np.float64;NumPy很聰明,它會將Python類型映射到等價的dtype上。

數組的dtype還有另一個屬性:

In [46]: int_array = np.arange(10) In [47]: calibers = np.array([.22, .270, .357, .380, .44, .50], dtype=np.float64) In [48]: int_array.astype(calibers.dtype) Out[48]: array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9.])

你還可以用簡潔的類型代碼來表示dtype:

In [49]: empty_uint32 = np.empty(8, dtype='u4') In [50]: empty_uint32 Out[50]: array([ 0, 1075314688, 0, 1075707904, 0, 1075838976, 0, 1072693248], dtype=uint32)

筆記:調用astype總會創建一個新的數組(一個數據的備份),即使新的dtype與舊的dtype相同。

 

NumPy數組的運算

數組很重要,因為它使你不用編寫循環即可對數據執行批量運算。NumPy用戶稱其為矢量化(vectorization)。大小相等的數組之間的任何算術運算都會將運算應用到元素級:

In [51]: arr = np.array([[1., 2., 3.], [4., 5., 6.]]) In [52]: arr Out[52]: array([[ 1., 2., 3.], [ 4., 5., 6.]]) In [53]: arr * arr Out[53]: array([[ 1., 4., 9.], [ 16., 25., 36.]]) In [54]: arr - arr Out[54]: array([[ 0., 0., 0.], [ 0., 0., 0.]])

數組與標量的算術運算會將標量值傳播到各個元素:

In [55]: 1 / arr Out[55]: array([[ 1. , 0.5 , 0.3333], [ 0.25 , 0.2 , 0.1667]]) In [56]: arr ** 0.5 Out[56]: array([[ 1. , 1.4142, 1.7321], [ 2. , 2.2361, 2.4495]])

大小相同的數組之間的比較會生成布爾值數組:

In [57]: arr2 = np.array([[0., 4., 1.], [7., 2., 12.]]) In [58]: arr2 Out[58]: array([[ 0., 4., 1.], [ 7., 2., 12.]]) In [59]: arr2 > arr Out[59]: array([[False, True, False], [ True, False, True]], dtype=bool)

不同大小的數組之間的運算叫做廣播(broadcasting),將在附錄A中對其進行詳細討論。本書的內容不需要對廣播機制有多深的理解。

 

基本的索引和切片

NumPy數組的索引是一個內容豐富的主題,因為選取數據子集或單個元素的方式有很多。一維數組很簡單。從表面上看,它們跟Python列表的功能差不多:

In [60]: arr = np.arange(10) In [61]: arr Out[61]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]) In [62]: arr[5] Out[62]: 5 In [63]: arr[5:8] Out[63]: array([5, 6, 7]) In [64]: arr[5:8] = 12 In [65]: arr Out[65]: array([ 0, 1, 2, 3, 4, 12, 12, 12, 8, 9])

如上所示,當你將一個標量值賦值給一個切片時(如arr[5:8]=12),該值會自動傳播(也就說后面將會講到的“廣播”)到整個選區。跟列表最重要的區別在於,數組切片是原始數組的視圖。這意味着數據不會被復制,視圖上的任何修改都會直接反映到源數組上。

作為例子,先創建一個arr的切片:

In [66]: arr_slice = arr[5:8] In [67]: arr_slice Out[67]: array([12, 12, 12])

現在,當我修稿arr_slice中的值,變動也會體現在原始數組arr中:

In [68]: arr_slice[1] = 12345 In [69]: arr Out[69]: array([ 0, 1, 2, 3, 4, 12, 12345, 12, 8, 9])

切片[ : ]會給數組中的所有值賦值:

In [70]: arr_slice[:] = 64 In [71]: arr Out[71]: array([ 0, 1, 2, 3, 4, 64, 64, 64, 8, 9])

如果你剛開始接觸NumPy,可能會對此感到驚訝(尤其是當你曾經用過其他熱衷於復制數組數據的編程語言)。由於NumPy的設計目的是處理大數據,所以你可以想象一下,假如NumPy堅持要將數據復制來復制去的話會產生何等的性能和內存問題。

注意:如果你想要得到的是ndarray切片的一份副本而非視圖,就需要明確地進行復制操作,例如arr[5:8].copy()

對於高維度數組,能做的事情更多。在一個二維數組中,各索引位置上的元素不再是標量而是一維數組:

In [72]: arr2d = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) In [73]: arr2d[2] Out[73]: array([7, 8, 9])

因此,可以對各個元素進行遞歸訪問,但這樣需要做的事情有點多。你可以傳入一個以逗號隔開的索引列表來選取單個元素。也就是說,下面兩種方式是等價的:

In [74]: arr2d[0][2] Out[74]: 3 In [75]: arr2d[0, 2] Out[75]: 3

圖4-1說明了二維數組的索引方式。軸0作為行,軸1作為列。

圖4-1 NumPy數組中的元素索引

在多維數組中,如果省略了后面的索引,則返回對象會是一個維度低一點的ndarray(它含有高一級維度上的所有數據)。因此,在2×2×3數組arr3d中:

In [76]: arr3d = np.array([[[1, 2, 3], [4, 5, 6]], [[7, 8, 9], [10, 11, 12]]]) In [77]: arr3d Out[77]: array([[[ 1, 2, 3], [ 4, 5, 6]], [[ 7, 8, 9], [10, 11, 12]]])

arr3d[0]是一個2×3數組:

In [78]: arr3d[0] Out[78]: array([[1, 2, 3], [4, 5, 6]])

標量值和數組都可以被賦值給arr3d[0]:

In [79]: old_values = arr3d[0].copy() In [80]: arr3d[0] = 42 In [81]: arr3d Out[81]: array([[[42, 42, 42], [42, 42, 42]], [[ 7, 8, 9], [10, 11, 12]]]) In [82]: arr3d[0] = old_values In [83]: arr3d Out[83]: array([[[ 1, 2, 3], [ 4, 5, 6]], [[ 7, 8, 9], [10, 11, 12]]])

相似的,arr3d[1,0]可以訪問索引以(1,0)開頭的那些值(以一維數組的形式返回):

In [84]: arr3d[1, 0] Out[84]: array([7, 8, 9])

雖然是用兩步進行索引的,表達式是相同的:

In [85]: x = arr3d[1] In [86]: x Out[86]: array([[ 7, 8, 9], [10, 11, 12]]) In [87]: x[0] Out[87]: array([7, 8, 9])

注意,在上面所有這些選取數組子集的例子中,返回的數組都是視圖。

 

切片索引

ndarray的切片語法跟Python列表這樣的一維對象差不多:

In [88]: arr Out[88]: array([ 0, 1, 2, 3, 4, 64, 64, 64, 8, 9]) In [89]: arr[1:6] Out[89]: array([ 1, 2, 3, 4, 64])

對於之前的二維數組arr2d,其切片方式稍顯不同:

In [90]: arr2d Out[90]: array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) In [91]: arr2d[:2] Out[91]: array([[1, 2, 3], [4, 5, 6]])

可以看出,它是沿着第0軸(即第一個軸)切片的。也就是說,切片是沿着一個軸向選取元素的。表達式arr2d[:2]可以被認為是“選取arr2d的前兩行”。

你可以一次傳入多個切片,就像傳入多個索引那樣:

In [92]: arr2d[:2, 1:] Out[92]: array([[2, 3], [5, 6]])

像這樣進行切片時,只能得到相同維數的數組視圖。通過將整數索引和切片混合,可以得到低維度的切片。

例如,我可以選取第二行的前兩列:

In [93]: arr2d[1, :2] Out[93]: array([4, 5])

相似的,還可以選擇第三列的前兩行:

In [94]: arr2d[:2, 2] Out[94]: array([3, 6])

圖4-2對此進行了說明。注意,“只有冒號”表示選取整個軸,因此你可以像下面這樣只對高維軸進行切片:

In [95]: arr2d[:, :1] Out[95]: array([[1], [4], [7]])

圖4-2 二維數組切片

自然,對切片表達式的賦值操作也會被擴散到整個選區:

In [96]: arr2d[:2, 1:] = 0 In [97]: arr2d Out[97]: array([[1, 0, 0], [4, 0, 0], [7, 8, 9]])

 

布爾型索引

來看這樣一個例子,假設我們有一個用於存儲數據的數組以及一個存儲姓名的數組(含有重復項)。在這里,我將使用numpy.random中的randn函數生成一些正態分布的隨機數據:

In [98]: names = np.array(['Bob', 'Joe', 'Will', 'Bob', 'Will', 'Joe', 'Joe']) In [99]: data = np.random.randn(7, 4) In [100]: names Out[100]: array(['Bob', 'Joe', 'Will', 'Bob', 'Will', 'Joe', 'Joe'], dtype='<U4') In [101]: data Out[101]: array([[ 0.0929, 0.2817, 0.769 , 1.2464], [ 1.0072, -1.2962, 0.275 , 0.2289], [ 1.3529, 0.8864, -2.0016, -0.3718], [ 1.669 , -0.4386, -0.5397, 0.477 ], [ 3.2489, -1.0212, -0.5771, 0.1241], [ 0.3026, 0.5238, 0.0009, 1.3438], [-0.7135, -0.8312, -2.3702, -1.8608]])

假設每個名字都對應data數組中的一行,而我們想要選出對應於名字"Bob"的所有行。跟算術運算一樣,數組的比較運算(如==)也是矢量化的。因此,對names和字符串"Bob"的比較運算將會產生一個布爾型數組:

In [102]: names == 'Bob' Out[102]: array([ True, False, False, True, False, False, False], dtype=bool)

這個布爾型數組可用於數組索引:

In [103]: data[names == 'Bob'] Out[103]: array([[ 0.0929, 0.2817, 0.769 , 1.2464], [ 1.669 , -0.4386, -0.5397, 0.477 ]])

布爾型數組的長度必須跟被索引的軸長度一致。此外,還可以將布爾型數組跟切片、整數(或整數序列,稍后將對此進行詳細講解)混合使用:

In [103]: data[names == 'Bob'] Out[103]: array([[ 0.0929, 0.2817, 0.769 , 1.2464], [ 1.669 , -0.4386, -0.5397, 0.477 ]])

注意:如果布爾型數組的長度不對,布爾型選擇就會出錯,因此一定要小心。

下面的例子,我選取了names == 'Bob'的行,並索引了列:

In [104]: data[names == 'Bob', 2:] Out[104]: array([[ 0.769 , 1.2464], [-0.5397, 0.477 ]]) In [105]: data[names == 'Bob', 3] Out[105]: array([ 1.2464, 0.477 ])

要選擇除"Bob"以外的其他值,既可以使用不等於符號(!=),也可以通過~對條件進行否定:

In [106]: names != 'Bob' Out[106]: array([False, True, True, False, True, True, True], dtype=bool) In [107]: data[~(names == 'Bob')] Out[107]: array([[ 1.0072, -1.2962, 0.275 , 0.2289], [ 1.3529, 0.8864, -2.0016, -0.3718], [ 3.2489, -1.0212, -0.5771, 0.1241], [ 0.3026, 0.5238, 0.0009, 1.3438], [-0.7135, -0.8312, -2.3702, -1.8608]])

~操作符用來反轉條件很好用:

In [108]: cond = names == 'Bob' In [109]: data[~cond] Out[109]: array([[ 1.0072, -1.2962, 0.275 , 0.2289], [ 1.3529, 0.8864, -2.0016, -0.3718], [ 3.2489, -1.0212, -0.5771, 0.1241], [ 0.3026, 0.5238, 0.0009, 1.3438], [-0.7135, -0.8312, -2.3702, -1.8608]])

選取這三個名字中的兩個需要組合應用多個布爾條件,使用&(和)、|(或)之類的布爾算術運算符即可:

In [110]: mask = (names == 'Bob') | (names == 'Will') In [111]: mask Out[111]: array([ True, False, True, True, True, False, False], dtype=bool) In [112]: data[mask] Out[112]: array([[ 0.0929, 0.2817, 0.769 , 1.2464], [ 1.3529, 0.8864, -2.0016, -0.3718], [ 1.669 , -0.4386, -0.5397, 0.477 ], [ 3.2489, -1.0212, -0.5771, 0.1241]])

通過布爾型索引選取數組中的數據,將總是創建數據的副本,即使返回一模一樣的數組也是如此。

注意:Python關鍵字and和or在布爾型數組中無效。要使用&與|。

通過布爾型數組設置值是一種經常用到的手段。為了將data中的所有負值都設置為0,我們只需:

In [113]: data[data < 0] = 0 In [114]: data Out[114]: array([[ 0.0929, 0.2817, 0.769 , 1.2464], [ 1.0072, 0. , 0.275 , 0.2289], [ 1.3529, 0.8864, 0. , 0. ], [ 1.669 , 0. , 0. , 0.477 ], [ 3.2489, 0. , 0. , 0.1241], [ 0.3026, 0.5238, 0.0009, 1.3438], [ 0. , 0. , 0. , 0. ]])

通過一維布爾數組設置整行或列的值也很簡單:

In [115]: data[names != 'Joe'] = 7 In [116]: data Out[116]: array([[ 7. , 7. , 7. , 7. ], [ 1.0072, 0. , 0.275 , 0.2289], [ 7. , 7. , 7. , 7. ], [ 7. , 7. , 7. , 7. ], [ 7. , 7. , 7. , 7. ], [ 0.3026, 0.5238, 0.0009, 1.3438], [ 0. , 0. , 0. , 0. ]])

后面會看到,這類二維數據的操作也可以用pandas方便的來做。

 

花式索引

花式索引(Fancy indexing)是一個NumPy術語,它指的是利用整數數組進行索引。假設我們有一個8×4數組:

In [117]: arr = np.empty((8, 4)) In [118]: for i in range(8): .....: arr[i] = i In [119]: arr Out[119]: 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 [120]: arr[[4, 3, 0, 6]] Out[120]: array([[ 4., 4., 4., 4.], [ 3., 3., 3., 3.], [ 0., 0., 0., 0.], [ 6., 6., 6., 6.]])

這段代碼確實達到我們的要求了!使用負數索引將會從末尾開始選取行:

In [121]: arr[[-3, -5, -7]] Out[121]: array([[ 5., 5., 5., 5.], [ 3., 3., 3., 3.], [ 1., 1., 1., 1.]])

一次傳入多個索引數組會有一點特別。它返回的是一個一維數組,其中的元素對應各個索引元組:

In [122]: arr = np.arange(32).reshape((8, 4)) In [123]: arr Out[123]: 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 [124]: arr[[1, 5, 7, 2], [0, 3, 1, 2]] Out[124]: array([ 4, 23, 29, 10])

附錄A中會詳細介紹reshape方法。

最終選出的是元素(1,0)、(5,3)、(7,1)和(2,2)。無論數組是多少維的,花式索引總是一維的。

這個花式索引的行為可能會跟某些用戶的預期不一樣(包括我在內),選取矩陣的行列子集應該是矩形區域的形式才對。下面是得到該結果的一個辦法:

In [125]: arr[[1, 5, 7, 2]][:, [0, 3, 1, 2]] Out[125]: array([[ 4, 7, 5, 6], [20, 23, 21, 22], [28, 31, 29, 30], [ 8, 11, 9, 10]])

記住,花式索引跟切片不一樣,它總是將數據復制到新數組中。

 

數組轉置和軸對換

轉置是重塑的一種特殊形式,它返回的是源數據的視圖(不會進行任何復制操作)。數組不僅有transpose方法,還有一個特殊的T屬性:

In [126]: arr = np.arange(15).reshape((3, 5)) In [127]: arr Out[127]: array([[ 0, 1, 2, 3, 4], [ 5, 6, 7, 8, 9], [10, 11, 12, 13, 14]]) In [128]: arr.T Out[128]: array([[ 0, 5, 10], [ 1, 6, 11], [ 2, 7, 12], [ 3, 8, 13], [ 4, 9, 14]])

在進行矩陣計算時,經常需要用到該操作,比如利用np.dot計算矩陣內積:

In [129]: arr = np.random.randn(6, 3) In [130]: arr Out[130]: array([[-0.8608, 0.5601, -1.2659], [ 0.1198, -1.0635, 0.3329], [-2.3594, -0.1995, -1.542 ], [-0.9707, -1.307 , 0.2863], [ 0.378 , -0.7539, 0.3313], [ 1.3497, 0.0699, 0.2467]]) In [131]: np.dot(arr.T, arr) Out[131]: array([[ 9.2291, 0.9394, 4.948 ], [ 0.9394, 3.7662, -1.3622], [ 4.948 , -1.3622, 4.3437]])

對於高維數組,transpose需要得到一個由軸編號組成的元組才能對這些軸進行轉置(比較費腦子):

In [132]: arr = np.arange(16).reshape((2, 2, 4)) In [133]: arr Out[133]: array([[[ 0, 1, 2, 3], [ 4, 5, 6, 7]], [[ 8, 9, 10, 11], [12, 13, 14, 15]]]) In [134]: arr.transpose((1, 0, 2)) Out[134]: array([[[ 0, 1, 2, 3], [ 8, 9, 10, 11]], [[ 4, 5, 6, 7], [12, 13, 14, 15]]])

這里,第一個軸被換成了第二個,第二個軸被換成了第一個,最后一個軸不變。

簡單的轉置可以使用.T,它其實就是進行軸對換而已。ndarray還有一個swapaxes方法,它需要接受一對軸編號:

In [135]: arr Out[135]: array([[[ 0, 1, 2, 3], [ 4, 5, 6, 7]], [[ 8, 9, 10, 11], [12, 13, 14, 15]]]) In [136]: arr.swapaxes(1, 2) Out[136]: 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 [137]: arr = np.arange(10) In [138]: arr Out[138]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]) In [139]: np.sqrt(arr) Out[139]: array([ 0. , 1. , 1.4142, 1.7321, 2. , 2.2361, 2.4495, 2.6458, 2.8284, 3. ]) In [140]: np.exp(arr) Out[140]: array([ 1. , 2.7183, 7.3891, 20.0855, 54.5982, 148.4132, 403.4288, 1096.6332, 2980.958 , 8103.0839])

這些都是一元(unary)ufunc。另外一些(如add或maximum)接受2個數組(因此也叫二元(binary)ufunc),並返回一個結果數組:

In [141]: x = np.random.randn(8) In [142]: y = np.random.randn(8) In [143]: x Out[143]: array([-0.0119, 1.0048, 1.3272, -0.9193, -1.5491, 0.0222, 0.7584, -0.6605]) In [144]: y Out[144]: array([ 0.8626, -0.01 , 0.05 , 0.6702, 0.853 , -0.9559, -0.0235, -2.3042]) In [145]: np.maximum(x, y) Out[145]: array([ 0.8626, 1.0048, 1.3272, 0.6702, 0.853 , 0.0222, 0.7584, -0.6605])

這里,numpy.maximum計算了x和y中元素級別最大的元素。

雖然並不常見,但有些ufunc的確可以返回多個數組。modf就是一個例子,它是Python內置函數divmod的矢量化版本,它會返回浮點數數組的小數和整數部分:

In [146]: arr = np.random.randn(7) * 5 In [147]: arr Out[147]: array([-3.2623, -6.0915, -6.663 , 5.3731, 3.6182, 3.45 , 5.0077]) In [148]: remainder, whole_part = np.modf(arr) In [149]: remainder Out[149]: array([-0.2623, -0.0915, -0.663 , 0.3731, 0.6182, 0.45 , 0.0077]) In [150]: whole_part Out[150]: array([-3., -6., -6., 5., 3., 3., 5.])

Ufuncs可以接受一個out可選參數,這樣就能在數組原地進行操作:

In [151]: arr Out[151]: array([-3.2623, -6.0915, -6.663 , 5.3731, 3.6182, 3.45 , 5.0077]) In [152]: np.sqrt(arr) Out[152]: array([ nan, nan, nan, 2.318 , 1.9022, 1.8574, 2.2378]) In [153]: np.sqrt(arr, arr) Out[153]: array([ nan, nan, nan, 2.318 , 1.9022, 1.8574, 2.2378]) In [154]: arr Out[154]: array([ nan, nan, nan, 2.318 , 1.9022, 1.8574, 2.2378])

表4-3和表4-4分別列出了一些一元和二元ufunc。

 

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 values Out[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))


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM