Python可以這樣玩(10)：分析與科學計算
學習完基本的GUI程式與繪圖功能之後，我們更進一步的來學習數學運算，並且將數學運算的結果用圖形來表示出來。很多人聽到數學就害怕，在這裡真的不需要，相反的，Python 反而可以幫你解決數學問題，絕對可以相輔相成。
是樂於資料分析與科學計算視覺化的模組相當多，我們就從一定要會的幾個模組開始。
擴展庫 numpy
如果你還沒有安裝 numpy 擴展庫，請透過 pip3 insatll numpy 安裝，這是一個基本的擴展庫，很多其他的擴展庫在安裝之前都需要先安裝它，所以應該不會有什麼問題。
在數學裡面有一樣東西叫做矩陣，電腦上稱為陣列，Python 的基本庫所提供的列表雖然很像陣列，但是並沒有計算功能。Numpy 提供了陣列，可以讓我們完成數學上的矩陣運算。
先看一下矩陣在數學上的定義，把行與列弄清楚，電腦也是一樣的定義：
矩陣的基本名詞：
(a)
元(element)：矩陣中列出來的每個數稱為矩陣的元。
(b)
列(row)：同一水平線各元合稱此矩陣的一列。
(c)
行(column)：同一鉛直線各元合稱此矩陣的一行。
(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.],
[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)
32
>>> np.dot(a,b)
32
>>>
我們再看一下二維與一維的運算：
>>> 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]
2
>>>
陣列元素還支援多個元素同時存取：
>>> 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)
45
>>> 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
留言
張貼留言