1. 程式人生 > >第四篇 NumPy基礎:陣列和⽮量計算

第四篇 NumPy基礎:陣列和⽮量計算

NumPy(Numerical Python的簡稱)是Python數值計算最重要的基礎包。⼤多數提供科學計算的包都是⽤NumPy的陣列作為構建基礎。
NumPy的部分功能如下:
  ndarray,⼀個具有⽮量算術運算和複雜⼴播能⼒的快速且節省空間的多維陣列。
  ⽤於對整組資料進⾏快速運算的標準數學函式(⽆需編寫迴圈)。
  ⽤於讀寫磁碟資料的⼯具以及⽤於操作記憶體對映⽂件的⼯具。
  線性代數、隨機數⽣成以及傅⾥葉變換功能。
  ⽤於整合由C、C++、Fortran等語⾔編寫的程式碼的A C API。
由於NumPy提供了⼀個簡單易⽤的C API,因此很容易將資料傳遞給由低階語⾔編寫的外部庫,外部庫也能以NumPy陣列的形式將資料返回給Python。這個功能使Python成為⼀種包裝C/C++/Fortran歷史程式碼庫的選擇,並使被包裝庫擁有⼀個動態的、易⽤的接⼝。

⼤部分資料分析應⽤,最關注的功能主要集中在:
  ⽤於資料整理和清理、⼦集構造和過濾、轉換等快速的⽮量化陣列運算。
  常⽤的陣列演算法,如排序、唯⼀化、集合運算等。
  ⾼效的描述統計和資料聚合/摘要運算。
  ⽤於異構資料集的合併/連線運算的資料對⻬和關係型資料運算。
  將條件邏輯表述為陣列表示式(⽽不是帶有if-elif-else分⽀的迴圈)。
  資料的分組運算(聚合、轉換、函式應⽤等)。
Python的⾯向陣列計算起源1995年,JimHugunin建立了Numeric庫。在2005年,Travis Oliphant從Numeric和Numarray項⽬整合出了NumPy項⽬,進⽽所有社群都集合到了這個框架下。

NumPy之於數值計算特別重要的原因之⼀,是因為它可以⾼效處理⼤陣列的資料。這是因為:NumPy是在⼀個連續的記憶體塊中儲存資料,獨⽴於其他Python內建物件。
  NumPy的C語⾔編寫的演算法庫可以操作記憶體,⽽不必進⾏型別檢查或其它前期⼯作。
  ⽐起Python的內建序列,NumPy陣列使⽤的記憶體更少。
NumPy可以在整個陣列上執⾏複雜的計算,⽽不需要Python的for迴圈。要搞明⽩具體的效能差距,考察⼀個包含⼀百萬整數的陣列,和⼀個等價的Python列表:
import numpy as np
my_arr = np.arange(1000000)
my_list = list(range(1000000))


%time for _ in range(10): my_arr2 = my_arr * 2         # 輸出:Wall time: 30 ms
%time for _ in range(10): my_list2 = [x * 2 for x in my_list]   # 輸出:Wall time: 1.24 s
基於NumPy的演算法要⽐純Python快10到100倍(甚⾄更快),並且使⽤的記憶體更少。

 

一、NumPy的ndarray:⼀種多維陣列物件
  NumPy最重要的⼀個特點就是其N維陣列物件(即ndarray),該物件是⼀個快速⽽靈活的⼤資料集容器。你可以利⽤這種陣列對整塊資料執⾏⼀些數學運算,其語法跟標量元素之間的運算⼀樣。要明⽩Python是如何利⽤與標量值類似的語法進⾏批次計算,我先引⼊NumPy,然後⽣成⼀個包含隨機資料的⼩陣列:
import numpy as np
data = np.random.randn(2, 3)  # Generate some random data
data               # 輸出如下:
array([[-0.77451155, 0.44650751, -0.04153134],
   [ 2.94286194, -1.39901653, -0.07859503]])
data * 10           # 輸出如下:,該操作不會改變原陣列data,所有的元素都乘以10
array([[ -7.7451155 , 4.46507507, -0.41531343],
   [ 29.4286194 , -13.99016527, -0.78595028]])
data + data          # 輸出如下:,該操作不會改變原陣列data,每個元素都與⾃身相加。
array([[-1.5490231 , 0.89301501, -0.08306269],
   [ 5.88572388, -2.79803305, -0.15719006]])
ndarray是⼀個通⽤的同構資料多維容器,也就是說,其中的所有元素必須是相同型別的。
每個陣列都有⼀個shape(⼀個表示各維度⼤⼩的元組)和⼀個dtype(⼀個⽤於說明陣列資料型別的物件)
data.shape     # 輸出:(2,3),表示2行3列
data.dtype     # 輸出:dtype('float64')
下面看到的“陣列”、“NumPy陣列”、"ndarray"時,都指的ndarray物件。

1、建立ndarray

建立陣列最簡單的辦法就是使⽤array函式。它接受⼀切序列型的物件(包括其他陣列),然後產⽣⼀個新的含有傳⼊資料的NumPy陣列。以⼀個列表的轉換為例:
data1 = [6, 7.5, 8, 0, 1]
arr1 = np.array(data1)     # python列表轉換為Numpy陣列
arr1             # 輸出:array([ 6. , 7.5, 8. , 0. , 1. ])

巢狀序列(⽐如由⼀組等⻓列表組成的列表)將會被轉換為⼀個多維陣列:
data2 = [[1, 2, 3, 4], [5, 6, 7, 8]]
arr2 = np.array(data2)
In [24]: arr2         # 輸出如下:
array([[1, 2, 3, 4],
   [5, 6, 7, 8]])
因為data2是列表的列表,NumPy陣列arr2的兩個維度的shape是從data2引⼊的。可以⽤屬性ndimshape驗證:
arr2.ndim          # 輸出:2
arr2.shape          # 輸出:(2, 4)
除⾮特別說明,np.array會嘗試為新建的這個陣列推斷出⼀個較為合適的資料型別。資料型別儲存在⼀個特殊的dtype物件中。⽐如說,在上⾯的兩個例⼦中,我們有:
arr1.dtype          # 輸出:dtype('float64')
arr2.dtype            # 輸出:dtype('int32')
除np.array之外,還有⼀些函式也可以新建陣列。⽐如,zerosones分別可以建立指定⻓度或形狀的全0或全1陣列。empty可以建立⼀個沒有任何具體值的陣列。要⽤這些⽅法建立多維陣列,只需傳⼊⼀個表示形狀的元組即可:
np.zeros(10)       # 輸出:array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
np.zeros((3, 6))        # 輸出如下:
array([[0., 0., 0., 0., 0., 0.],
   [0., 0., 0., 0., 0., 0.],
   [0., 0., 0., 0., 0., 0.]])
np.empty((2, 3, 2))     # 輸出如下:
array([[[0., 0.],
   [0., 0.],
   [0., 0.]],
     [[0., 0.],
   [0., 0.],
   [0., 0.]]])
np.empty會返回的都是⼀些未初始化的垃圾值。認為返回全0是不安全的。
arange是Python內建函式range的陣列版:
np.arange(15)     # 輸出:array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14])
Numpy陣列建立常用函式如下表4-1所示:

 

2、ndarray的資料型別
dtype(資料型別)是⼀個特殊的物件,它含有ndarray將⼀塊記憶體解釋為特定資料型別所需的資訊:
arr1 = np.array([1, 2, 3], dtype=np.float64)
arr2 = np.array([1, 2, 3], dtype=np.int32)
arr1.dtype                # 輸出:dtype('float64')
arr2.dtype                # 輸出:dtype('int32')
dtype是NumPy靈活互動其它系統的源泉之⼀。多數情況下,它們直接對映到相應的機器表示,這使得“讀寫磁碟上的⼆進位制資料流”以及“整合低階語⾔程式碼(如C、Fortran)”等⼯作變得更加簡單。

數值型dtype的命名⽅式相同:⼀個型別名(如float或int),後⾯跟⼀個⽤於表示各元素位⻓的數字。
標準的雙精度浮點值(即Python中的float物件)需要佔⽤8位元組(即64位)。因此,該型別在NumPy中就記作float64。
NumPy的資料型別如下表4-2所示:

 

 

可通過ndarray的astype⽅法明確地將⼀個數組從⼀個dtype轉換成另⼀個dtype:
arr = np.array([1, 2, 3, 4, 5])
arr.dtype          # dtype('int64')
float_arr = arr.astype(np.float64)
float_arr.dtype       # dtype('float64')

如果將浮點數轉換成整數,則⼩數部分將會被擷取刪除:
arr = np.array([3.7, -1.2, -2.6, 0.5, 12.9, 10.1])
arr.astype(np.int32)      # array([ 3, -1, -2, 0, 12, 10], dtype=int32)

如果某字串陣列表示的全是數字,也可以⽤astype將其轉換為數值形式:
numeric_strings = np.array(['1.25', '-9.6', '42'], dtype=np.string_)
numeric_strings.astype(float)   # array([ 1.25, -9.6 , 42. ])
注意:使⽤numpy.string_型別時,⼀定要⼩⼼,因為NumPy的字串資料是⼤⼩固定的,發⽣擷取時,不會發出警告。轉換過程因為某種原因⽽失敗了(⽐如某個不能被轉換為float64的字串),就會引發⼀個ValueError。

陣列的dtype還有另⼀個屬性:將本陣列的dtype設定為另一個數組的dtype
int_array = np.arange(10)
calibers = np.array([.22, .270, .357, .380, .44, .50], dtype=np.float64)
int_array.astype(calibers.dtype)           # array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9.])
int_array.astype(calibers.dtype).dtype == 'float64'   # 輸出:True

還可以⽤簡潔的型別程式碼來表示dtype:

empty_uint32 = np.empty(8, dtype='u4')
empty_uint32     # 輸出如下:
array([ 0, 1075314688, 0, 1075707904, 0,
   1075838976, 0, 1072693248], dtype=uint32)

調⽤astype總會建立⼀個新的陣列(⼀個數據的備份),即使新的dtype與舊的dtype相同。

3、NumPy陣列的運算
陣列不⽤編寫迴圈即可對資料執⾏批量運算。稱為⽮量化(vectorization)。⼤⼩相等的陣列之間的任何算術運算都會將運算應⽤到元素級:
arr = np.array([[1., 2., 3.], [4., 5., 6.]])   # arr = np.arange(1,7).reshape(2,3),用該方法建立同等shape的整數陣列
arr * arr     # 輸出如下:
array([[ 1., 4., 9.],
   [ 16., 25., 36.]])
arr - arr     # 輸出如下:
array([[ 0., 0., 0.],
   [ 0., 0., 0.]])
陣列與標量的算術運算會將標量值傳播到各個元素:
1 / arr     # 輸出如下:
array([[ 1. , 0.5 , 0.3333],
   [ 0.25 , 0.2 , 0.1667]])
arr ** 0.5   # 輸出如下:
array([[ 1. , 1.4142, 1.7321],
   [ 2. , 2.2361, 2.4495]])
⼤⼩相同的陣列之間的⽐較會⽣成布林值陣列:
arr2 = np.array([[0., 4., 1.], [7., 2., 12.]])
arr2 > arr     # 輸出如下:元素級比較
array([[False, True, False],
   [ True, False, True]], dtype=bool)
不同⼤⼩的陣列之間的運算叫做⼴播(broadcasting)

4、基本的索引和切⽚

⼀維陣列很簡單。從表⾯上看,它們跟Python列表的功能差不多。
arr = np.arange(10)
arr     # 輸出:array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
arr[5:8] = 12   # 切片賦值,自動傳播到整個選區
arr       # 輸出:array([ 0, 1, 2, 3, 4, 12, 12, 12, 8, 9])
當你將⼀個標量值賦值給⼀個切⽚時(如arr[5:8]=12),該值會⾃動傳播到整個選區。跟列表最重要的區別在於,陣列切⽚是原始陣列的檢視。這意味著資料不會被複制,檢視上的任何修改都會直接反映到源陣列上。作為例⼦,先建立⼀個arr的切⽚:

arr_slice = arr[5:8]
當修改arr_slice中的值,變動也會體現在原始陣列arr中:
arr_slice[1] = 12345      # 原始陣列中的值也會改變
arr             # 輸出:array([ 0, 1, 2, 3, 4, 12, 12345, 12, 8, 9])
切⽚[ : ]會給陣列中的所有值賦值:
arr_slice[:] = 64
arr       # array([ 0, 1, 2, 3, 4, 64, 64, 64, 8, 9])
注意:如果要得到的是ndarray切⽚的⼀份副本⽽⾮檢視,就需要明確地進⾏複製操作,例如 arr[5:8].copy()

對於⾼維度陣列,能做的事情更多。在⼀個⼆維陣列中,各索引位置上的元素不再是標量⽽是⼀維陣列:
arr2d = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
arr2d[2]   # array([7, 8, 9])
因此,可以對各個元素進⾏遞迴訪問,但這樣需要做的事情有點多。你可以傳⼊⼀個以逗號隔開的索引列表來選取單個元素。也就是說,下⾯兩種⽅式是等價的:
arr2d[0][2]    # 輸出:3
arr2d[0, 2]    # 輸出:3
arr2d的索引方式中,軸0代表行,軸2代表列

在多維陣列中,如果省略了後⾯的索引,則返回物件會是⼀個維度低⼀點的ndarray(它含有⾼⼀級維度上的所有資料)。因此,在2×2×3陣列arr3d中:
arr3d = np.array([[[1, 2, 3], [4, 5, 6]], [[7, 8, 9], [10, 11, 12]]])
arr3d[0]是⼀個2×3陣列:
arr3d[0]    # 輸出如下:
array([[1, 2, 3],
   [4, 5, 6]])

標量值和陣列都可以被賦值給 arr3d[0]:
old_values = arr3d[0].copy()
arr3d[0] = 42
arr3d     # 輸出如下:
array([[[42, 42, 42],
    [42, 42, 42]],
   [[ 7, 8, 9],
    [10, 11, 12]]])
arr3d[0] = old_values    # 賦回原值
arr3d          # 輸出如下:
array([[[ 1, 2, 3],
    [ 4, 5, 6]],
   [[ 7, 8, 9],
   [10, 11, 12]]])

相似的,arr3d[1,0]可以訪問索引以(1,0)開頭的那些值(以⼀維陣列的形式返回):
arr3d[1, 0]      # 輸出:array([7, 8, 9])
雖然是⽤兩步進⾏索引的,表示式是相同的:
x = arr3d[1]     # 先賦值給x
x          # 輸出如下:
array([[ 7, 8, 9],
   [10, 11, 12]])
x[0]         # 輸出:array([7, 8, 9])
在上⾯所有這些選取陣列⼦集的例⼦中,返回的陣列都是檢視。

 

5、切⽚索引

ndarray的切⽚語法跟Python列表這樣的⼀維物件差不多:對於之前的⼆維陣列arr2d,其切⽚⽅式稍顯不同
arr2d = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
arr2d[:2]     # 輸出如下:
array([[1, 2, 3],
   [4, 5, 6]])
切⽚是沿著⼀個軸向選取元素的。表示式arr2d[:2]可以被認為是“選取arr2d的前兩⾏”。
可以⼀次傳⼊多個切⽚,就像傳⼊多個索引那樣:
arr2d[:2, 1:]   # 輸出如下:
array([[2, 3],
   [5, 6]])
像這樣進⾏切⽚時,只能得到相同維數的陣列檢視。通過將整數索引和切⽚混合,可以得到低維度的切⽚。
選取第⼆⾏的前兩列:
arr2d[1, :2]     # 輸出:array([4, 5])
相似的,還可以選擇第三列的前兩⾏:
arr2d[:2, 2]     # 輸出:array([3, 6])

只有冒號”表示選取整個軸,因此你可以像下⾯這樣只對⾼維軸進⾏切⽚:
arr2d[:, :1]      # 輸出如下:
array([[1],
   [4],
   [7]])

對切⽚表示式的賦值操作也會被擴散到整個選區:
arr2d[:2, 1:] = 0
arr2d        # 輸出如下:
array([[1, 0, 0],
   [4, 0, 0],
   [7, 8, 9]])

 

6、布林型索引,看一個例子:
假設我們有⼀個⽤於儲存資料的陣列以及⼀個儲存姓名的陣列(含有重複項)。先使⽤numpy.random中的randn函式⽣成⼀些正態分佈的隨機資料
names = np.array(['Bob', 'Joe', 'Will', 'Bob', 'Will', 'Joe', 'Joe'])
data = np.random.randn(7, 4)
names     # 輸出:array(['Bob', 'Joe', 'Will', 'Bob', 'Will', 'Joe', 'Joe'], dtype='<U4')
data      # 輸出如下:
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"的⽐較運算將會產⽣⼀個布林型陣列:
names == 'Bob'     # 輸出:array([ True, False, False, True, False, False, False])
這個布林型陣列可⽤於陣列索引:
data[names == 'Bob']   # 輸出如下:(如果要取相應的列,可這樣做:data[:,names == 'Bob'],前提是長度要一致)
array([[ 0.0929, 0.2817, 0.769 , 1.2464],
   [ 1.669 , -0.4386, -0.5397, 0.477 ]])

布林型陣列的⻓度必須跟被索引的軸⻓度⼀致。此外,還可以將布林型陣列跟切⽚、整數(或整數序列)混合使⽤。
注意:如果布林型陣列的⻓度不對,布林型選擇就會出錯,因此⼀定要⼩⼼。

下⾯的例⼦,選取了names == 'Bob'的⾏,並索引了列:布林型陣列與切片混合使用
data[names == 'Bob', 2:]  # 輸出如下:
array([[ 0.769 , 1.2464],
   [-0.5397, 0.477 ]])
data[names == 'Bob', 3]   # 輸出:array([ 1.2464, 0.477 ])

要選擇除"Bob"以外的其他值,既可以使⽤不等於符號(!=),也可以通過(~)對條件進⾏否定
names != 'Bob'       # 輸出:array([False, True, True, False, True, True, True])
data[~(names == 'Bob')]   # 輸出如下:
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]])
(~)操作符⽤來反轉條件很好⽤:
cond = names == 'Bob'
data[~cond]     # 輸出與前面的一樣

選取這三個名字中的兩個需要組合應⽤多個布林條件,使⽤&(與)、|(或)之類的布林算術運算子即可:
mask = (names == 'Bob') | (names == 'Will')
mask       # 輸出:array([ True, False, True, True, True, False, False])
data[mask]     # 輸出如下:
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,我們只需:
data[data < 0] = 0
data       # 輸出如下:
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. ]])

通過⼀維布林陣列設定整⾏或列的值也很簡單:
data[names != 'Joe'] = 7
data         # 輸出如下:
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. ]])

 

7、花式索引(複雜索引,Fancy indexing)

花式索引(Fancy indexing)是⼀個NumPy術語,它指的是利⽤整數陣列進⾏索引。假設有⼀個8×4陣列:
arr = np.empty((8, 4))
for i in range(8):
  arr[i] = i
arr     # 輸出如下:
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即可:
arr[[4, 3, 0, 6]]
array([[4., 4., 4., 4.],
  [3., 3., 3., 3.],
  [0., 0., 0., 0.],
  [6., 6., 6., 6.]])
這段程式碼確實達到我們的要求了!使⽤負數索引將會從末尾開始選取⾏
arr[[-3, -5, -7]]
array([[5., 5., 5., 5.],
  [3., 3., 3., 3.],
  [1., 1., 1., 1.]])
⼀次傳⼊多個索引陣列會有⼀點特別。它返回的是⼀個⼀維陣列,其中的元素對應各個索引元組
arr = np.arange(32).reshape((8, 4))
arr     # 輸出如下:
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]])
arr[[1, 5, 7, 2], [0, 3, 1, 2]]     #輸出:array([ 4, 23, 29, 10])
最終選出的是元素(1,0)、(5,3)、(7,1)和(2,2)。⽆論陣列是多少維的,花式索引總是⼀維的。
花式索引的⾏為可能會跟預期不⼀樣,選取矩陣的⾏列⼦集應該是矩形區域的形式才對。下⾯是得到該結果的⼀個辦法:
arr[[1, 5, 7, 2]][:, [0, 3, 1, 2]]   # 輸出如下:
array([[ 4, 7, 5, 6],      # 該行選取的元素是:(1, 0),(1, 3),(1, 1),(1, 2)
  [20, 23, 21, 22],     # 該行選取的元素是:(5, 0),(5, 3),(5, 1),(5, 2),以此類推
  [28, 31, 29, 30],
  [ 8, 11, 9, 10]])
花式索引跟切⽚不⼀樣,它總是將資料複製到新陣列中。

 

8、陣列轉置和軸對換
轉置是重塑的⼀種特殊形式,它返回的是源資料的檢視(不會進⾏任何複製操作)。
陣列不僅有transpose⽅法,還有⼀個特殊的T屬性
arr = np.arange(15).reshape((3,5))
arr         # 輸出如下:
array([[ 0, 1, 2, 3, 4],
  [ 5, 6, 7, 8, 9],
  [10, 11, 12, 13, 14]])
arr.T       # 輸出如下:使用T屬性轉置
array([[ 0, 5, 10],
  [ 1, 6, 11],
  [ 2, 7, 12],
  [ 3, 8, 13],
  [ 4, 9, 14]])
在進⾏矩陣計算時,經常需要⽤到該操作,⽐如利⽤np.dot計算矩陣內積
arr = np.random.randn(6, 3)
arr       # 輸出如下:
array([[-2.13468628, 1.06582408, -0.72489158],
  [-1.17032574, -0.82376495, 0.07941442],
  [-0.29019922, -1.98512296, 0.12906019],
  [ 0.14798076, 0.9989999 , -0.83814468],
  [-1.61693841, -0.23485318, -0.38723287],
  [ 0.09921124, -1.21635853, -2.44492694]])
np.dot(arr.T, arr)   # 輸出如下:
array([[ 8.65699441, -0.32814613, 1.67656039],
  [-0.32814613, 8.28796773, 1.1333181 ],
  [ 1.67656039, 1.1333181 , 7.37853451]])

對於⾼維陣列,transpose需要得到⼀個由軸編號組成的元組才能對這些軸進⾏轉置
arr = np.arange(16).reshape((2, 2, 4))
arr      # 輸出如下:
array([[[ 0, 1, 2, 3],
  [ 4, 5, 6, 7]],
  [[ 8, 9, 10, 11],
  [12, 13, 14, 15]]])
arr.transpose((1, 0, 2))   # 輸出如下:
array([[[ 0, 1, 2, 3],
  [ 8, 9, 10, 11]],
  [[ 4, 5, 6, 7],
  [12, 13, 14, 15]]])
這⾥,第⼀個軸被換成了第⼆個,第⼆個軸被換成了第⼀個,最後⼀個軸不變。

簡單的轉置可以使⽤.T,它其實就是進⾏軸對換⽽已。ndarray還有⼀個swapaxes⽅法,它需要接受⼀對軸編號:
arr.swapaxes(1, 2)   # 輸出如下:
array([[[ 0, 4],
   [ 1, 5],
   [ 2, 6],
   [ 3, 7]],

   [[ 8, 12],
   [ 9, 13],
   [10, 14],
   [11, 15]]])
swapaxes也是返回源資料的檢視(不會進⾏任何複製操作)。

 

二、通⽤函式:快速的元素級陣列函式

 

通⽤函式(即ufunc)是⼀種對ndarray中的資料執⾏元素級運算的函式。你可以將其看做簡單函式(接受⼀個或多個標量值,併產⽣⼀個或多個標量值)的⽮量化包裝器。許多ufunc都是簡單的元素級變體,如sqrt和exp
arr = np.arange(10)
arr       # 輸出:array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
np.sqrt(arr)   # 輸出:(相當於arr ** 0.5)
array([0. , 1. , 1.41421356, 1.73205081, 2. ,
    2.23606798, 2.44948974, 2.64575131, 2.82842712, 3. ])
這些都是⼀元(unary)ufunc。另外⼀些(如add或maximum)接受2個數組(因此也叫⼆元(binary)ufunc),並返回⼀個結果陣列:
x = np.random.randn(8)
y = np.random.randn(8)
x     # 輸出如下:
array([ 1.20051688, -1.65710862, -0.53963283, -1.21187257, 0.63507697,
  -0.37209111, -0.93943269, -1.00426802])
y     # 輸出如下:
array([ 1.05467508, 1.4219286 , -0.92969225, -1.04814943, -0.82260389,
  0.91427391, 0.2504812 , 0.97563071])
np.maximum(x, y)   # 輸出如下:實際比較的是(x[0], y[0]), (x[1], y[1])...
array([ 1.20051688, 1.4219286 , -0.53963283, -1.04814943, 0.63507697,
  0.91427391, 0.2504812 , 0.97563071])
numpy.maximum計算了x和y中元素級別最⼤的元素。

modf就是Python內建函式divmod的⽮量化版本,它會返回浮點數陣列的⼩數和整數部分
arr = np.random.randn(7) * 5
arr     # 原始資料如下:
array([-3.52389072, -3.4832977 , 2.56810837, 1.2339455 , 6.10703696,
  7.61798681, 3.05483746])
remainder, whole_part = np.modf(arr)
remainder      # 輸出小數部分
array([-0.52389072, -0.4832977 , 0.56810837, 0.2339455 , 0.10703696,
  0.61798681, 0.05483746])
whole_part     # 輸出整數部分
array([-3., -3., 2., 1., 6., 7., 3.])

Ufuncs接受out選項引數,可以讓它們在陣列的原地進⾏操作:
arr       # 輸出如下:
array([ 3.43072847, 10.82286699, -3.63291583, 4.07716536, -1.60187118, -1.24558177, -5.24618548])
np.sqrt(arr)   # 輸出如下:
array([1.85222258, 3.28981261, nan, 2.01919919, nan, nan, nan])
arr       # 輸出如下:
array([1.85222258, 3.28981261, nan, 2.01919919, nan, nan, nan])

一元ufunc(一元通用函式)如下表4-3所示:

二元ufunc(二元通用函式)如下表4-4所示:

 

三、利⽤陣列進⾏資料處理


⽤陣列表示式代替迴圈的做法,通常被稱為⽮量化。⽮量化陣列運算要⽐等價的純Python⽅式快上⼀兩個數量級(甚⾄更多),尤其是各種數值計算。在後⾯將介紹⼴播,這是⼀種針對⽮量化計算的強⼤⼿段。看一個例子:假設我們想要在⼀組值(⽹格型)上計算函式sqrt(x^2+y^2)。

np.meshgrid函式接受兩個⼀維陣列,併產⽣兩個⼆維矩陣(對應於兩個陣列中所有的(x,y)對):
points = np.arange(-5, 5, 0.01)    # 1000 equally spaced points
xs, ys = np.meshgrid(points, points)
ys               # 輸出如下:
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]])
現在,把這兩個陣列當做兩個浮點數那樣編寫表示式即可:
z = np.sqrt(xs ** 2 + ys ** 2)
z           # 輸出如下:
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
import matplotlib.pyplot as plt
plt.imshow(z, cmap=plt.cm.gray); plt.colorbar()   # 輸出:<matplotlib.colorbar.Colorbar at 0x1f488261e80>
plt.title("Image plot of $\sqrt{x^2 + y^2}$ for a grid of values")    # 輸出:Text(0.5,1,'Image plot of $\\sqrt{x^2 + y^2}$ for a grid of values')

 

1、將條件邏輯表述為陣列運算
numpy.where函式是三元表示式x if condition else y的⽮量化版本。假設我們有⼀個布林陣列和兩個值陣列:
xarr = np.array([1.1, 1.2, 1.3, 1.4, 1.5])
yarr = np.array([2.1, 2.2, 2.3, 2.4, 2.5])
cond = np.array([True, False, True, True, False])
根據cond中的值選取xarr和yarr的值:當cond中的值為True時,選取xarr的值,否則從yarr中選取。列表推導式的寫法應該如下所示:
result = [(x if c else y) for x, y, c in zip(xarr, yarr, cond)]
result     # 輸出:[1.1, 2.2, 1.3, 1.4, 2.5]
第⼀,它對⼤陣列的處理速度不是很快(因為所有⼯作都是由純Python完成的)。
第⼆,⽆法⽤於多維陣列。若使⽤np.where,則可以將該功能寫得⾮常簡潔:
result1 = np.where(cond, xarr, yarr)
result1     # 輸出:array([1.1, 2.2, 1.3, 1.4, 2.5])
np.where的第⼆個和第三個引數不必是陣列,它們都可以是標量值。在資料分析⼯作中,where通常⽤於根據另⼀個數組⽽產⽣⼀個新的陣列。

假設有⼀個由隨機資料組成的矩陣,你希望將所有正值替換為2,將所有負值替換為-2。若利⽤np.where,則會⾮常簡單:
arr = np.random.randn(4,4)
arr       # 輸出如下:
array([[ 0.80426609, 1.0974934 , 1.52641017, 0.26308459],
  [ 0.39494499, -0.09829609, 2.93032116, 0.77040912],
  [-1.35818382, 0.78914328, 1.38073343, 0.98920804],
  [ 0.93070469, 0.32393877, -0.76748938, -0.94859212]])
arr > 0     # 輸出如下:
array([[ True, True, True, True],
  [ True, False, True, True],
  [False, True, True, True],
  [ True, True, False, False]])
np.where(arr > 0, 2, -2)   # 輸出如下:
array([[ 2, 2, 2, 2],
  [ 2, -2, 2, 2],
  [-2, 2, 2, 2],
  [ 2, 2, -2, -2]])
使⽤np.where,可以將標量和陣列結合起來。例如,我可⽤常數2替換arr中所有正的值:
np.where(arr > 0, 2, arr)     # 大於0的值設定為2
array([[ 2. , 2. , 2. , 2. ],
  [ 2. , -0.09829609, 2. , 2. ],
  [-1.35818382, 2. , 2. , 2. ],
  [ 2. , 2. , -0.76748938, -0.94859212]])
傳遞給where的陣列⼤⼩可以不相等,甚⾄可以是標量值。

 

2、數學和統計⽅法

可以通過陣列上的⼀組數學函式對整個陣列或某個軸向的資料進⾏統計計算。sum、mean以及標準差std等聚合計算(aggregation,通常叫做約簡(reduction))既可以當做陣列的例項⽅法調⽤,也可以當做頂級NumPy函式使⽤。下面例子,先生成一些正態分佈隨機資料,然後做了聚類統計:
arr = np.random.randn(5,4)
arr        # 輸出如下:
array([[-0.1133632 , 0.58949215, 0.275524 , -0.3855296 ],
  [ 2.55070219, 0.63003176, -0.34253814, -0.01309244],
  [ 1.04023817, -1.9999175 , -0.40581791, -0.36883092],
  [-0.71064938, -1.40922625, -1.18840513, -0.92986452],
  [-0.05223665, 0.5629013 , 0.69139223, -1.45967786]])
arr.mean()      # 輸出:-0.15194338563310547
np.mean(arr)      # 輸出:-0.15194338563310547
arr.sum()      # 輸出:-3.0388677126621095

mean和sum這類的函式可以接受⼀個axis選項引數,⽤於計算該軸向上的統計值,最終結果是⼀個少⼀維的陣列:
arr.mean(axis=1)   # 輸出:array([ 0.09153084, 0.70627584, -0.43358204, -1.05953632, -0.06440525])
arr.sum(axis=0)    # 輸出:array([ 2.71469112, -1.62671854, -0.96984495, -3.15699534])
這⾥,arr.mean(1)是“計算⾏的平均值”,arr.sum(0)是“計算每列的和”。

其他如cumsum和cumprod之類的⽅法則不聚合,⽽是產⽣⼀個由中間結果組成的陣列:
arr = np.array([0, 1, 2, 3, 4, 5, 6, 7])
arr.cumsum()     # 輸出:array([ 0, 1, 3, 6, 10, 15, 21, 28], dtype=int32),累加

在多維陣列中,累加函式(如cumsum)返回的是同樣⼤⼩的陣列,但是會根據每個低維的切⽚沿著標記軸計算部分聚類:
arr = np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])
arr.cumsum(axis=0)      # 縱向累計求和,輸出如下:
array([[ 0, 1, 2],
  [ 3, 5, 7],
  [ 9, 12, 15]], dtype=int32)
arr.cumprod(axis=1)      # 橫向累計求積,輸出如下:
array([[ 0, 0, 0],
  [ 3, 12, 60],
  [ 6, 42, 336]], dtype=int32)

全部的基本陣列統計方法如下表4-5所示:

 

3、⽤於布林型陣列的⽅法
在上⾯這些⽅法中,布林值會被強制轉換為1(True)和0(False)。因此,sum經常被⽤來對布林型陣列中的True值計數:
arr = np.random.randn(100)
(arr > 0).sum()     # 輸出:一個隨機數,大於0的個數
另外還有兩個⽅法any和all,它們對布林型陣列⾮常有⽤。any⽤於測試陣列中是否存在⼀個或多個True,⽽all則檢查陣列中所有值是否都是True:
bools = np.array([False, False, True, False])
bools.any()       # 輸出:True
bools.all()        # 輸出:False
這兩個⽅法也能⽤於⾮布林型陣列,所有⾮0元素將會被當做True。

 

4、排序

NumPy陣列也可以通過sort⽅法就地排序:要改變原始陣列
arr = np.random.randn(6)
arr          # 輸出如下:
array([ 1.37750254, 0.36255115, -0.87588226, -0.80512566, -1.33751372, 0.18521148])
arr.sort()        # 就地排序
arr            # 輸出如下:
array([-1.33751372, -0.87588226, -0.80512566, 0.18521148, 0.36255115, 1.37750254])

多維陣列可以在任何⼀個軸向上進⾏排序,只需將軸編號傳給sort即可:
arr = np.random.randn(5, 3)
arr        # 輸出如下:
array([[-0.70497864, 0.75979962, -0.04453997],
  [ 1.44800768, -1.67338973, 1.1078619 ],
  [ 0.40954131, 1.03518218, -0.49360008],
  [ 0.30501957, 1.29319448, 1.02856608],
  [-1.16897745, -1.51113507, 0.9717796 ]])
arr.sort(1)      # 對每一行的數進行排序
arr
array([[-0.70497864, -0.04453997, 0.75979962],
  [-1.67338973, 1.1078619 , 1.44800768],
  [-0.49360008, 0.40954131, 1.03518218],
  [ 0.30501957, 1.02856608, 1.29319448],
  [-1.51113507, -1.16897745, 0.9717796 ]])
頂級⽅法np.sort返回的是陣列的已排序副本,⽽就地排序則會修改陣列本身。計算陣列分位數最簡單的辦法是對其進⾏排序,然後選取特定位置的值:
large_arr = np.random.randn(1000)     # 生成1000個隨機數
large_arr.sort()              # 從小到大排序
large_arr[int(0.05 * len(large_arr))]      # 5% 分位數(quantile)
# 輸出:-1.5311513550102103

 

5、唯⼀化以及其它的集合邏輯
NumPy提供了⼀些針對⼀維ndarray的基本集合運算。最常⽤的可能要數np.unique了,
它⽤於找出陣列中的唯⼀值並返回已排序的結果:
names = np.array(['Bob', 'Joe', 'Will', 'Bob', 'Will', 'Joe', 'Joe'])
# 去掉重複元素,並按字母排序,
np.unique(names) # 輸出:array(['Bob', 'Joe', 'Will'], dtype='<U4')
ints = np.array([3, 3, 3, 2, 2, 1, 1, 4, 4]) # 對數字排序
np.unique(ints) # 輸出:array([1, 2, 3, 4])
拿np.unique跟等價的純Python程式碼來對⽐⼀下:
sorted(set(names)) # 輸出:['Bob', 'Joe', 'Will']

另⼀個函式np.in1d⽤於測試⼀個數組中的值在另⼀個數組中的成員資格,返回⼀個布林型陣列:
values = np.array([6, 0, 0, 3, 2, 5, 6])
np.in1d(values, [2, 3, 6]) # 輸出:array([ True, False, False, True, True, False, True])

NumPy陣列的集合運算(表4-6):


四、⽤於陣列的⽂件輸⼊輸出


NumPy能夠讀寫磁碟上的⽂本資料或⼆進位制資料。先討論NumPy的內建⼆進位制格式,因為更多的⽤戶會使⽤pandas或其它⼯具載入⽂本或表格資料。

np.save和np.load是讀寫磁碟陣列資料的兩個主要函式。預設情況下,陣列是以未壓縮的原始⼆進位制格式儲存在副檔名為.npy的⽂件中的:例如:
arr = np.arange(10)
np.save('some_array', arr)
如果⽂件路徑末尾沒有副檔名.npy,則該副檔名會被⾃動加上。然後就可以通過np.load讀取磁碟上的陣列:
np.load('some_array.npy')   # 輸出:array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

通過np.savez可以將多個數組儲存到⼀個未壓縮⽂件中,將陣列以關鍵字引數的形式傳⼊即可:
np.savez('array_archive.npz', a=arr, b=arr)
載入.npz⽂件時,你會得到⼀個類似字典的物件,該物件會對各個陣列進⾏延遲載入:
arch = np.load('array_archive.npz')
arch['b']      # 輸出:array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

如果資料壓縮的很好,就可以使⽤numpy.savez_compressed
np.savez_compressed('arrays_compressed.npz', a=arr, b=arr)

 

五、線性代數


線性代數(如矩陣乘法、矩陣分解、⾏列式以及其他⽅陣數學等)是任何陣列庫的重要組成部分。MATLAB通過 對兩個⼆維陣列相乘得到的是⼀個元素級的積,⽽不是⼀個矩陣點積。NumPy提供了⼀個⽤於矩陣乘法的dot函式(既是⼀個數組⽅法也是numpy名稱空間中的⼀個函式):
x = np.array([[1., 2., 3.], [4., 5., 6.]])
y = np.array([[6., 23.], [-1, 7], [8, 9]])
x.dot(y)        # 輸出如下:得到的結果是行數等於x的行,列數等於y的列
array([[ 28., 64.],
  [ 67., 181.]])
x.dot(y)等價於np.dot(x, y):
np.dot(x, y)     # 輸出同上

⼀個⼆維陣列跟⼀個⼤⼩合適的⼀維陣列的矩陣點積運算之後將會得到⼀個⼀維陣列:
np.dot(x, np.ones(3))      # 輸出:array([ 6., 15.])
x陣列的每個元素乘以第二個陣列的每個元素後相加的結果

@符(類似Python 3.5)也可以⽤作中綴運算子,進⾏矩陣乘法
x @ np.ones(3)        # 輸出:array([ 6., 15.])

numpy.linalg中有⼀組標準的矩陣分解運算以及諸如求逆和⾏列式之類的東⻄。它們跟MATLAB和R等語⾔所使⽤的是相同的⾏業標準線性代數庫,如BLAS、LAPACK、Intel MKL(MathKernel Library,可能有,取決於你的NumPy版本)等:
from numpy.linalg import inv, qr
X = np.random.randn(5, 5)
mat = X.T.dot(X)
inv(mat)          # 輸出如下:
array([[ 0.28133188, -0.0977005 , -0.47041506, 0.01876554, 0.32981926],
  [-0.0977005 , 0.31194322, 0.15049957, -0.27441188, -0.58712427],
  [-0.47041506, 0.15049957, 1.1281112 , -0.07145714, -0.59076183],
  [ 0.01876554, -0.27441188, -0.07145714, 0.56844553, 0.74076293],
  [ 0.32981926, -0.58712427, -0.59076183, 0.74076293, 1.68251546]])
mat.dot(inv(mat))     # 輸出如下:
array([[ 1.00000000e+00, -7.03456949e-17, -2.26629202e-16, -3.58528577e-17, 2.46988638e-16],
  [-4.69384500e-16, 1.00000000e+00, -1.58143616e-17, -3.57503814e-16, -5.44789041e-16],
  [-2.31542198e-16, 1.52236389e-16, 1.00000000e+00,-1.57566693e-16, -4.74249967e-16],
  [ 2.37545108e-16, 1.70674192e-16, 1.02588110e-16, 1.00000000e+00, -5.48276527e-16],
  [-2.69371029e-16, 1.04852032e-16, -5.72195720e-17, -7.39959808e-17, 1.00000000e+00]])
q, r = qr(mat)
r         # 輸出如下:
array([[-16.80324568, -0.94736212, -5.62018513, -5.74117222, 3.61353809],
   [ 0., -10.27117765, -0.92342637, 0.583159 , -4.36651116],
   [ 0.  , 0. , -1.83415284, 4.48040783, -2.86906256],
   [ 0. , 0. , 0. , -3.01189889, 1.35175951],
   [ 0. , 0. , 0. , 0. , 0.48899415]])
表示式X.T.dot(X)計算X和它的轉置X.T的點積。

常用的numpy.linalg函式(線性代數函式)表4-7:

 

六、偽隨機數⽣成


numpy.random模組對Python內建的random進⾏了補充,增加了⼀些⽤於⾼效⽣成多種概率分佈的樣本值的函式。例如,可⽤normal來得到⼀個標準正態分佈的4×4樣本陣列:
samples = np.random.normal(size=(4, 4))
samples      # 輸出如下:
array([[-1.22523734, 1.0020195 , 0.62450895, -0.30984473],
  [ 0.17302941, -1.1195615 , -1.4692486 , 0.47059667],
  [ 1.12656905, 0.66808849, -0.69090706, 0.13304373],
  [ 0.46783217, -1.33587931, 1.19216589, 0.13784556]])
⽽Python內建的random模組則只能⼀次⽣成⼀個樣本值。從下⾯的測試結果中可以看出,如果需要產⽣⼤量樣本值,numpy.random快了不⽌⼀個數量級:
from random import normalvariate
N = 1000000
%timeit samples = [normalvariate(0, 1) for _ in range(N)]   # 輸出如下:
886 ms ± 5.22 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit np.random.normal(size=N)              # 輸出如下:
40.4 ms ± 229 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

說它是偽隨機數,是因為它們都是通過演算法基於隨機數⽣成器種⼦,在確定性的條件下⽣成的
你可以⽤NumPy的np.random.seed更改隨機數⽣成種⼦:
np.random.seed(1234)

numpy.random的資料⽣成函式使⽤了全域性的隨機種⼦。要避免全域性狀態,可以使⽤numpy.random.RandomState,建立⼀個與其它隔離的隨機數⽣成器:
rng = np.random.RandomState(1234)
rng.randn(10)             # 輸出如下:
array([ 0.47143516, -1.19097569, 1.43270697, -0.3126519 , -0.72058873,
  0.88716294, 0.85958841, -0.6365235 , 0.01569637, -2.24268495])

下面列出了numpy.random中的部分函式(表4-8):

 

七、示例:隨機漫步


通過模擬隨機漫步來說明如何運⽤陣列運算。⼀個簡單的隨機漫步的例⼦:從0開始,步⻓1和-1出現的概率相等。
下⾯是⼀個通過內建的random模組以純Python的⽅式實現1000步的隨機漫步:
import random
position = 0
walk = [ position ]
steps = 1000
for i in range(steps):
  step = 1 if random.randint(0, 1) else -1      # random.randint(0, 1)隨機選取0或1
  position += step
  walk.append(position)
根據前100個隨機漫步值⽣成的折線圖:
%matplotlib
import matplotlib.pyplot as plt
plt.plot(walk[:100])    # 輸出如下圖照片所示

從輸出照片可以看出,這其實就是隨機漫步中各步的累計和,可以⽤⼀個數組運算來實現。⽤np.random模組⼀次性隨機產⽣1000個“擲硬幣”結果(即兩個數中任選⼀個),將其分別設定為1或 -1,然後計算累計和:
nsteps = 1000
draws = np.random.randint(0, 2, size=nsteps)
steps = np.where(draws > 0, 1, -1)
walk = steps.cumsum()    # cumsum()方法是累加

有了這些資料就可以沿著漫步路徑做⼀些統計⼯作,⽐如求取最⼤值和最⼩值:
walk.min()      # 輸出:-9,每次結果不一樣
walk.max()     # 輸出:60

現在來看⼀個複雜點的統計任務——⾸次穿越時間,即隨機漫步過程中第⼀次到達某個特定值的時間。假設想要知道本次隨機漫步需要多久才能距離初始0點⾄少10步遠(任⼀⽅向均可)。np.abs(walk)>=10可以得到⼀個布林型陣列,它表示的是距離是否達到或超過10,⽽我們想要知道的是第⼀個10或 -10的索引。可以⽤argmax來解決這個問題,它返回的是該布林型陣列第⼀個最⼤值的索引(True就是最⼤值):
(np.abs(walk) >= 10).argmax()     # 輸出:20
使⽤argmax並不是很⾼效,它會對陣列進⾏完全掃描。在本例中,只要發現⼀個True,就知道它是個最⼤值。

 

1、⼀次模擬多個隨機漫步
模擬多個隨機漫步過程(⽐如5000個),只需對上⾯的程式碼做⼀點點修改即可⽣成所有的隨機漫步過程。只要給numpy.random的函式傳⼊⼀個⼆元元組就可以產⽣⼀個⼆維陣列,然後我們就可以⼀次性計算5000個隨機漫步過程(⼀⾏⼀個)的累計和了:
nwalks = 5000
nsteps = 1000
draws = np.random.randint(0, 2, size=(nwalks, nsteps))    # 在[0,1]之間選取,size引數指定要生成的矩陣行數和列數
steps = np.where(draws > 0, 1, -1)
walks = steps.cumsum(1)       # 輸出如下:
array([[ 1, 0, 1, ..., 42, 43, 44],
  [ -1, -2, -1, ..., 40, 41, 42],
  [ -1, -2, -1, ..., -42, -43, -44],
  ...,
  [ -1, 0, 1, ..., 66, 65, 64],
  [ -1, -2, -1, ..., 0, -1, 0],
  [ -1, -2, -1, ..., 34, 35, 34]], dtype=int32)
現在,來計算所有隨機漫步過程的最⼤值和最⼩值:
walks.max()       # 輸出:118
walks.min()       # 輸出:-130
得到這些資料之後,來計算30或 -30的最⼩穿越時間。這⾥稍微複雜些,因為不是5000個過程都到達了30。可以⽤any⽅法來對此進⾏檢查:
hits30 = (np.abs(walks) >= 30).any(1)  # 返回每一行的“滿足條件的第一個True,均不滿足條件就返回False”,組成的一維陣列
hits30         # array([ True, True, True, ..., True, False, True])
hits30.sum()      # 輸出:3381,達到絕對值30次以上的個數
然後利⽤這個布林型陣列選出那些穿越了30(絕對值)的隨機漫步(⾏),並調⽤argmax在軸1上獲取穿越時間:
crossing_times = (np.abs(walks[hits30]) >= 30).argmax(1)
crossing_times.mean()   # 輸出:502.03756285122745

嘗試⽤其他分佈⽅式得到漫步資料。只需使⽤不同的隨機數⽣成函式即可,如normal⽤於⽣成指定均值和標準差的正態分佈資料:
steps = np.random.normal(loc=0, scale=0.25, size=(nwalks, nsteps))