Python可以這樣玩(10)：分析與科學計算

- 4月 09, 2018

學習完基本的GUI程式與繪圖功能之後，我們更進一步的來學習數學運算，並且將數學運算的結果用圖形來表示出來。很多人聽到數學就害怕，在這裡真的不需要，相反的，Python 反而可以幫你解決數學問題，絕對可以相輔相成。

是樂於資料分析與科學計算視覺化的模組相當多，我們就從一定要會的幾個模組開始。

擴展庫 numpy

如果你還沒有安裝 numpy 擴展庫，請透過 pip3 insatll numpy 安裝，這是一個基本的擴展庫，很多其他的擴展庫在安裝之前都需要先安裝它，所以應該不會有什麼問題。

在數學裡面有一樣東西叫做矩陣，電腦上稱為陣列，Python 的基本庫所提供的列表雖然很像陣列，但是並沒有計算功能。Numpy 提供了陣列，可以讓我們完成數學上的矩陣運算。

先看一下矩陣在數學上的定義，把行與列弄清楚，電腦也是一樣的定義：

矩陣的基本名詞：

(a) 元(element)：矩陣中列出來的每個數稱為矩陣的元。

(b) 列(row)：同一水平線各元合稱此矩陣的一列。

(d) 位於第 i 列，第 j 行的元稱為(i,j)元。

(e) 當一個矩陣 M 有 n 列 m 行時，我們稱 M 為 n×m 階的矩陣。

(f) 當一個矩陣M有n列n行時，我們稱M為n階方陣。

產生陣列(Array)

電腦語言對矩陣的行、列定義都一樣，只有一點，就是數學的 (i,j) 從1開始，Python的 [i][j] 從0開始，關於這點到後面的運算便可以體會。下面直接看範例：

>>> import numpy as np

>>> np.array((1,2,3,4,5))

array([1, 2, 3, 4, 5])

>>> np.array([1,2,3,4,5]) # 元組或列表都可以轉成 array

array([1, 2, 3, 4, 5])

>>> np.array(range(1,6))

array([1, 2, 3, 4, 5])

>>> np.array([[1,2,3],[4,5,6]])

array([[1, 2, 3],

[4, 5, 6]])

>>> np.linspace(0,10,11) # 0 到 10 之間產生十個等差數列

array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.])

>>> np.linspace(0,8,11)

array([0. , 0.8, 1.6, 2.4, 3.2, 4. , 4.8, 5.6, 6.4, 7.2, 8. ])

>>> np.logspace(0,100,10) # 對數數列

array([1.00000000e+000, 1.29154967e+011, 1.66810054e+022, 2.15443469e+033,

2.78255940e+044, 3.59381366e+055, 4.64158883e+066, 5.99484250e+077,

7.74263683e+088, 1.00000000e+100])

>>> np.zeros((3,4)) # 有三列，四行，同數學定義

array([[0., 0., 0., 0.],

[0., 0., 0., 0.],

[0., 0., 0., 0.]])

>>> np.ones((4,3))

array([[1., 1., 1.],

[1., 1., 1.],

[1., 1., 1.]])

>>> np.identity(3) # 單位距離

array([[1., 0., 0.],

[0., 1., 0.],

[0., 0., 1.]])

>>> np.empty((2,5)) # 空陣列，只給空間不初始化，元素值不確定

array([[1.00000000e+000, 1.29154967e+011, 1.66810054e+022,

2.15443469e+033, 2.78255940e+044],

[3.59381366e+055, 4.64158883e+066, 5.99484250e+077,

7.74263683e+088, 1.00000000e+100]])

>>> np.empty((2,3))

array([[0., 0., 0.],

[0., 0., 0.]])

>>>

陣列與數值的演算

先看下面的數學定義：

>>> import numpy as np

>>> A = np.array([[1,2,3],[4,5,6]])

>>> A

array([[1, 2, 3],

[4, 5, 6]])

>>> A * 2

array([[ 2, 4, 6],

[ 8, 10, 12]])

>>> A * (-1)

array([[-1, -2, -3],

[-4, -5, -6]])

>>> A / 2

array([[0.5, 1. , 1.5],

[2. , 2.5, 3. ]])

>>> A // 2

array([[0, 1, 1],

[2, 2, 3]], dtype=int32)

>>> A + 2

array([[3, 4, 5],

[6, 7, 8]])

>>> A ** 2

array([[ 1, 4, 9],

[16, 25, 36]], dtype=int32)

>>> A % 3

array([[1, 2, 0],

[1, 2, 0]], dtype=int32)

>>>

關於數值與矩陣的運算，數學上其實只有定義 rA，並沒有定義 r + A，所以前面出現的加法減法，單純屬於 Python 運算子功能的範疇，只能用電腦程式語言的觀點來看，不能用數學來思考。

陣列與陣列的運算

數學上定義的矩陣加減法，與Python 相同，沒有問題：

>>> A = np.array([[9,4,0],[-5,6,7]])

>>> B = np.array([[-1,2,3],[4,5,6]])

>>> A + B # 與上面答案相同

array([[ 8, 6, 3],

[-1, 11, 13]])

>>> A - B

array([[10, 2, -3],

[-9, 1, 1]])

>>> C = A * B # 數學上乘法的定義並非如此

>>> C

array([[ -9, 8, 0],

[-20, 30, 42]])

>>> C / A # 分母為0 報錯

Warning (from warnings module):

File "__main__", line 1

RuntimeWarning: invalid value encountered in true_divide

array([[-1., 2., nan],

[ 4., 5., 6.]])

>>> C / B

array([[ 9., 4., 0.],

[-5., 6., 7.]])

>>> B * B

array([[ 1, 4, 9],

[16, 25, 36]])

>>> B / B

array([[1., 1., 1.],

[1., 1., 1.]])

>>> B - B

array([[0, 0, 0],

[0, 0, 0]])

Python 的 * 運算子作用在兩個 array 上面的時候，就是很單純的把兩陣列每個元素相乘，這在數學上並不叫做矩陣相乘，而是叫做 “阿達馬乘積”，定義如下：

>>> x = np.array([[1,3,2],[1,0,0,],[1,2,2]])

>>> x = np.array([[1,3,2],[1,0,0],[1,2,2]])

>>> x

array([[1, 3, 2],

[1, 0, 0],

[1, 2, 2]])

>>> y = np.array([[0,0,2],[7,5,0],[2,1,1]])

>>> y

array([[0, 0, 2],

[7, 5, 0],

[2, 1, 1]])

>>> x * y

array([[0, 0, 4],

[7, 0, 0],

[2, 2, 2]])

>>>

但是真正的矩陣相乘的定義，稱為一般矩陣乘積，定義如下：

直接由定義計算：

這種矩陣乘積亦可由稍微不同的觀點來思考：把向量和各係數相乘後相加起來。設A和B是兩個給定如下的矩陣：

舉個例子來說明：

>>> a = np.array([[1,0,2],[-1,3,1]])

>>> a

array([[ 1, 0, 2],

[-1, 3, 1]])

>>> b = np.array([[3,1],[2,1],[1,0]])

>>> b

array([[3, 1],

[2, 1],

[1, 0]])

>>> a * b

Traceback (most recent call last):

File "<pyshell#45>", line 1, in <module>

a * b

ValueError: operands could not be broadcast together with shapes (2,3) (3,2)

>>>

當然報錯啦，因為它們的維度不同，前面說過維度不同的兩個矩陣式不可以做加法運算的，這裡 array 的乘法比照矩陣加法運算原則。那麼數學上的矩陣相乘在 Python 要如何表呢？先別著急，我們後面對詳細說明。由於 array 與數學的矩陣有很多共通處，所以我利用數學觀念來談 array，順便讓讀者溫習一下數學，對後面的大數據分析會有幫助。

二維陣列轉置

上面的例子，因為維度不同，我們不能運算，一個是 3x2，一個是 2x3，如果們將其中一個 array 透過轉置，就可以運算了。所以，轉置的定義很簡單，就是將 3x2 變成 2x3。數學上的定義是這樣：

不過如果是一維陣列，轉置出來的結果是一樣的，其實這點與數學的矩陣是不同的，這點要記住：

>>> a = np.array([1,2,3,4,5])

>>> a

array([1, 2, 3, 4, 5])

>>> a.T

array([1, 2, 3, 4, 5])

>>>

向量積

又稱為向量內積，還記得數學上矩陣的乘法運算嗎？寫成 A x B = C，乘積C中的每一個元素，就是A的列元素與B的行元素的向量積，也就是把元素相乘之後再加總。

>>> a = np.array([1,2,3])

>>> b = np.array([4,5,6])

>>> a.dot(b)

>>> np.dot(a,b)

>>>

我們再看一下二維與一維的運算：

>>> c = np.array([[1,2,3],[4,5,6],[7,8,9]])

>>> c

array([[1, 2, 3],

[4, 5, 6],

[7, 8, 9]])

>>> a

array([1, 2, 3])

>>> c.dot(a)

array([14, 32, 50])

>>> ct = c.T

>>> ct

array([[1, 4, 7],

[2, 5, 8],

[3, 6, 9]])

>>> ct.dot(a)

array([30, 36, 42])

有趣的事情來了，兩個一維陣列的內積就是兩個向量的內積會等於一個數字，但是兩個多維陣列的內積則會是一個多維陣列，這點跟前面談到的數學上矩陣乘積有點類似，所以我們再看前面的例子：

>>> a = np.array([[1,0,2],[-1,3,1]])

>>> b = np.array([[3,1],[2,1],[1,0]])

>>> a

array([[ 1, 0, 2],

[-1, 3, 1]])

>>> b

array([[3, 1],

[2, 1],

[1, 0]])

>>> a.dot(b)

array([[5, 1],

[4, 2]])

答案就對了，所以說，我們可以下個結論，Python 的陣列A * 陣列B，其結果是數學矩陣的阿達馬乘積；Python 的陣列A.dot(陣列B)、或是 dot(陣列A, 陣列B)才是數學矩陣的乘積。

特別延伸說明，A * B = B * A，乘法有交換率，但是向量內積的交換率不存在：

不過只有一個例外，就是當 A 或 B 其中一個為單位方陣的時候才成立。請參考數學範例：

什麼是單位方陣呢？就是一開始我們提到的 np.identity(3)，數學定義如下：

看完這樣的說明，相信對 identity() 的用處就更清楚了：

>>> a

array([[ 1, 0, 2],

[-1, 3, 1]])

>>> b

array([[3, 1],

[2, 1],

[1, 0]])

>>> a.dot(b)

array([[5, 1],

[4, 2]])

>>> b.dot(a)

array([[2, 3, 7],

[1, 3, 5],

[1, 0, 2]])

>>> c=b.dot(a)

>>> c

array([[2, 3, 7],

[1, 3, 5],

[1, 0, 2]])

>>> d = np.identity(3)

>>> d

array([[1., 0., 0.],

[0., 1., 0.],

[0., 0., 1.]])

>>> c.dot(d)

array([[2., 3., 7.],

[1., 3., 5.],

[1., 0., 2.]])

>>> d.dot(c)

array([[2., 3., 7.],

[1., 3., 5.],

[1., 0., 2.]])

>>>

存取陣列元素

數學的部分談完了，接下來就是 Python 語法部分，如何透過 index 取出陣列的元素呢？

>>> c = np.array([[1,2,3],[4,5,6],[7,8,9]])

>>> c

array([[1, 2, 3],

[4, 5, 6],

[7, 8, 9]])

>>> c[0]

array([1, 2, 3])

>>> c[0][1]

>>>

陣列元素還支援多個元素同時存取：

>>> x = np.linspace(0,100,11)

>>> x

array([ 0., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.])

>>> index = np.random.randint(0, len(x), 5)

>>> index

array([ 9, 10, 4, 9, 7])

>>> x[index]

array([ 90., 100., 40., 90., 70.])

>>>

對陣列進行函數運算

np.arange() 的作用就跟之前談過去 range() 類似：

>>> x = np.arange(0,100,10)

>>> x

array([ 0, 10, 20, 30, 40, 50, 60, 70, 80, 90])

>>> np.sin(x)

array([ 0. , -0.54402111, 0.91294525, -0.98803162, 0.74511316,

-0.26237485, -0.30481062, 0.77389068, -0.99388865, 0.89399666])

>>> x = np.arange(0,100,9)

>>> x

array([ 0, 9, 18, 27, 36, 45, 54, 63, 72, 81, 90, 99])

對陣列不同維度上的元素進行運算

直接看例子：

>>> x = np.arange(0,10).reshape(2,5)

>>> x

array([[0, 1, 2, 3, 4],

[5, 6, 7, 8, 9]])

>>> np.sum(x)

>>> np.sum(x, axis=0)

array([ 5, 7, 9, 11, 13])

>>> np.sum(x, axis=1)

array([10, 35])

>>>

改變陣列大小

>>> a = np.arange(1,11,1)

>>> a

array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])

>>> a.reshape(2,5)

array([[ 1, 2, 3, 4, 5],

[ 6, 7, 8, 9, 10]])

>>> a

array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])

>>> a.shape = 2,5

>>> a

array([[ 1, 2, 3, 4, 5],

[ 6, 7, 8, 9, 10]])

剛才一維陣列無法轉置問題，可以這樣解決，就是用 shape 先轉置一次，以後就可以做一維的轉置了，不同處在紅色處：

>>> a = np.arange(1,11,1)

>>> a

array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])

>>> a.shape = 10,1

>>> a

array([[ 1],

[ 2],

[ 3],

[ 4],

[ 5],

[ 6],

[ 7],

[ 8],

[ 9],

[10]])

>>> a.T

array([[ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]])

>>> b = a.T

>>> b

array([[ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]])

>>> b.T

array([[ 1],

[ 2],

[ 3],

[ 4],

[ 5],

[ 6],

[ 7],

[ 8],

[ 9],

[10]])

>>>

這種雙括號的一列陣列，其實就是 Python 的矩陣了，我們來比較一下下面三種的不同：

>>> x = np.matrix(range(1,11,1))

>>> x

matrix([[ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]])

>>> y = np.arange(1,11,1)

>>> y

array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])

>>> z = list(range(1,11,1))

>>> z

[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

>>> x.T

matrix([[ 1],

[ 2],

[ 3],

[ 4],

[ 5],

[ 6],

[ 7],

[ 8],

[ 9],

[10]])

>>> y.T

array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])

>>> z.T

Traceback (most recent call last):

File "<pyshell#141>", line 1, in <module>

z.T

AttributeError: 'list' object has no attribute 'T'

>>> y

array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])

>>> y.shape = 10,1

>>> y

array([[ 1],

[ 2],

[ 3],

[ 4],

[ 5],

[ 6],

[ 7],

[ 8],

[ 9],

[10]])

>>> y = y.T

>>> y

array([[ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]])

>>>

科學計算擴展庫 scipy

透過 pip install scipy 就可以安裝此擴展庫，此擴展庫在 numpy 的基礎上面增加大量應用於科學計算與工程計算的模組。SciPy包含的模組有最佳化、線性代數、積分、插值、特殊函式、快速傅立葉變換、訊號處理和圖像處理、常微分方程求解和其他科學與工程中常用的計算。

由於這個擴展庫太過於複雜，因為科學的領域太大，我就只針對最基本的淺談，其他深入的部分，就搭配像是 matplotlib 擴展庫使用的時候，再行討論。

數學、物理常用常數與單位模組 constants

例如，下面的方法可以用來存取預先定義的常數：

>>> from scipy import constants as C

>>> C.pi

3.141592653589793

>>> C.golden

1.618033988749895

>>> C.c

299792458.0

>>> C.mile

1609.3439999999998

>>> C.inch

0.0254

>>> C.oz

0.028349523124999998

>>> C.lb

0.45359236999999997

特殊函數模組

一些特別的數學函數，包含計算立方根，排列組合等。

>>> from scipy import special as S

>>> S.cbrt(8)

2.0

>>> S.exp10(3)

1000.0

>>> S.sindg(90)

1.0

>>> S.round(3.5)

4.0

>>> S.round(3.4)

3.0

>>> S.comb(5,3)

10.0

>>> S.perm(5,3)

60.0

搜尋此網誌

AndrewPythonQuasar