數學建模三劍客 MSN
前言
不管是不是巴薩的球迷，只要你喜歡足球，就一定聽說過梅西(Messi)、蘇亞雷斯(Suarez)和內馬爾(Neymar)這個MSN組合。在眾多的數學建模輔助工具中，也有一個犀利無比的MSN組合，他們就是python麾下大名鼎鼎的 Matplotlib + Scipy + Numpy三劍客。
本文是我整理的MSN學習筆記，有些理解可能比較膚淺，甚至是錯誤的。如果因此誤導了某位看官，在工作中造成重大失誤或損失，我頂多只能賠償一頓飯——還得是我們樓下的十元盒飯。特此宣告。
文中程式碼均從我的這臺時不時出點問題、鬧個情緒的Yoga 3 pro上覆制而來，這意味著所有的程式碼均可在下面的執行環境中順利執行：
- pyhton 2.7.8
- numpy 1.11.1
- scipy 0.16.1
- matplotlib 1.5.1
三劍客之 Numpy
numpy是一個開源的python科學計算庫，包含了很多實用的數學函式，涵蓋線性代數、傅立葉變換和隨機數生成等功能。最初的numpy其實是scipy的一部分，後來才從scipy中分離出來。
numpy不是python的標準庫，需要單獨安裝。假定你的執行環境已經安裝了python包管理工具pip，numpy的安裝就非常簡單：
pip install numpy
陣列物件
ndarray是多維陣列物件，也是numpy最核心的物件。在numpy中，陣列的維度(dimensions)叫做軸(axes)，軸的個數叫做秩(rank)。通常，一個numpy陣列的所有元素都是同一種類型的資料，而這些資料的儲存和陣列的形式無關。
下面的例子，建立了一個三維的陣列（在匯入numpy時，一般都簡寫成np）。
import numpy as np a = np.array([[1,2,3],[4,5,6],[7,8,9]])
資料型別
numpy支援的資料型別主要有布林型(bool)、整型(integrate)、浮點型(float)和複數型(complex)，每一種資料型別根據佔用記憶體的位元組數又分為多個不同的子型別。常見的資料型別見下表。
| 型別 | 描述 |
|---|---|
| bool | 用1位儲存的布林型別（值為TRUE或FALSE） |
| inti | 由所在平臺決定其精度的整數（一般為int32或int64） |
| int8 | 1位元組整數 |
| int16 | 2位元組整數 |
| int32 | 4位元組整數 |
| int64 | 8位元組整數 |
| uint8 | 1位元組無符號整數 |
| uint16 | 2位元組無符號整數 |
| uint32 | 4位元組無符號整數 |
| uint64 | 8位元組無符號整數 |
| float16 | 半精度浮點數（16位），1位符號，5位指數，10位尾數 |
| float32 | 單精度浮點數（32位），1位符號，8位指數，23位尾數 |
| float64/float | 雙精度浮點數（64位），1位符號，11位指數，52位尾數 |
| complex64 | 複數，分別用32位表示實部和虛部 |
| complex128/complex | 複數，分別用64位表示實部和虛部 |
建立陣列
通常，我們用np.array()建立陣列。如果僅僅是建立一維陣列，也可以使用np.arange()或者np.linspace()的方法。np.zeros()、np.ones()、np.eye()則可以構造特殊的資料。np.random.randint()和np.random.random()則可以構造隨機數陣列。
>>> np.array([[1,2,3],[4,5,6]]) # 預設元素型別為int32 array([[1, 2, 3], [4, 5, 6]]) >>> np.array([[1,2,3],[4,5,6]], dtype=np.int8) # 指定元素型別為int8 array([[1, 2, 3], [4, 5, 6]], dtype=int8) >>> np.arange(5) # 預設元素型別為int32 array([0, 1, 2, 3, 4]) >>> np.arange(3,8, dtype=np.int8) # 指定元素型別為int8 array([3, 4, 5, 6, 7], dtype=int8) >>> np.arange(12).reshape(3,4) # 改變shape array([[ 0, 1, 2, 3], [ 4, 5, 6, 7], [ 8, 9, 10, 11]]) >>> np.linspace(1,2,5) # 從1到2生成5個浮點數 array([ 1. , 1.25, 1.5 , 1.75, 2. ]) >>> np.zeros((2,3)) # 全0陣列 array([[ 0., 0., 0.], [ 0., 0., 0.]]) >>> np.ones((2,3)) # 全1陣列 array([[ 1., 1., 1.], [ 1., 1., 1.]]) >>> np.eye(3) # 主對角線元素為1其他元素為0 array([[ 1., 0., 0.], [ 0., 1., 0.], [ 0., 0., 1.]]) >>> np.random.random((2,3)) # 生成[0,1)之間的隨機浮點數 array([[ 0.84731148, 0.8222318 , 0.85799278], [ 0.59371558, 0.92330741, 0.04518351]]) >>> np.random.randint(0,10,(3,2)) # 生成[0,10)之間的隨機整數 array([[2, 4], [8, 3], [8, 5]])
構造複雜陣列
很多時候，我們需要從簡單的資料結構，構造出複雜的陣列。例如，用一維的資料生成二維格點。
重複陣列: tile
>>> a = np.arange(5) >>> a array([0, 1, 2, 3, 4]) >>> np.tile(a, 2) array([0, 1, 2, 3, 4, 0, 1, 2, 3, 4]) >>> np.tile(a, (3,2)) array([[0, 1, 2, 3, 4, 0, 1, 2, 3, 4], [0, 1, 2, 3, 4, 0, 1, 2, 3, 4], [0, 1, 2, 3, 4, 0, 1, 2, 3, 4]])
一維陣列網格化: meshgrid
>>> a = np.arange(5) >>> b = np.arange(5,10) >>> np.meshgrid(a,b) [array([[0, 1, 2, 3, 4], [0, 1, 2, 3, 4], [0, 1, 2, 3, 4], [0, 1, 2, 3, 4], [0, 1, 2, 3, 4]]), array([[5, 5, 5, 5, 5], [6, 6, 6, 6, 6], [7, 7, 7, 7, 7], [8, 8, 8, 8, 8], [9, 9, 9, 9, 9]])] >>>
指定範圍和分割方式的網格化: mgrid
>>> np.mgrid[0:1:2j, 1:2:3j] array([[[ 0. , 0. , 0. ], [ 1. , 1. , 1. ]], [[ 1. , 1.5, 2. ], [ 1. , 1.5, 2. ]]]) >>> np.mgrid[0:1:0.3, 1:2:0.4] array([[[ 0. , 0. , 0. ], [ 0.3, 0.3, 0.3], [ 0.6, 0.6, 0.6], [ 0.9, 0.9, 0.9]], [[ 1. , 1.4, 1.8], [ 1. , 1.4, 1.8], [ 1. , 1.4, 1.8], [ 1. , 1.4, 1.8]]])
上面的例子中用到了虛數。構造虛數的方法如下：
>>> complex(2,5) (2+5j)
陣列的屬性
numpy的陣列物件除了一些常規的屬性外，也有幾個類似轉置、扁平迭代器等看起來更像是方法的屬性。扁平迭代器也許是遍歷多維陣列的一個簡明方法，下面的程式碼給出了一個例子。
>>> a = np.array([[1,2,3],[4,5,6]])
>>> a.dtype # 陣列元素的資料型別
dtype('int32')
>>> a.dtype.itemsize # 陣列元素佔據的記憶體位元組數
4
>>> a.itemsize # 陣列元素佔據的記憶體位元組數
4
>>> a.shape # 陣列的維度
(2, 3)
>>> a.size # 陣列元素個數
6
>>> a.T # 陣列行變列，類似於transpose()
array([[1, 4],
[2, 5],
[3, 6]])
>>> a.flat # 返回一個扁平迭代器，用於遍歷多維陣列
<numpy.flatiter object at 0x037188F0>
>>> for item in a.flat:
print item
1
2
...
改變陣列維度
numpy陣列的儲存順序和陣列的維度是不相干的，因此改變陣列的維度是非常便捷的操作，除resize()外，這一類操作不會改變所操作的陣列本身的儲存順序。
>>> a = np.array([[1,2,3],[4,5,6]]) >>> a.shape # 檢視陣列維度 (2, 3) >>> a.reshape(3,2) # 返回3行2列的陣列 array([[1, 2], [3, 4], [5, 6]]) >>> a.ravel() # 返回一維陣列 array([1, 2, 3, 4, 5, 6]) >>> a.transpose() # 行變列（類似於矩陣轉置） array([[1, 4], [2, 5], [3, 6]]) >>> a.resize((3,2)) # 類似於reshape，但會改變所操作的陣列 >>> a array([[1, 2], [3, 4], [5, 6]])
索引和切片
對於一維陣列的索引和切片，numpy和python的list一樣，甚至更靈活。
a = np.arange(9) >>> a[-1] # 最後一個元素 8 >>> a[2:5] # 返回第2到第5個元素 array([2, 3, 4]) >>> a[:7:3] # 返回第0到第7個元素，步長為3 array([0, 3, 6]) >>> a[::-1] # 返回逆序的陣列 array([8, 7, 6, 5, 4, 3, 2, 1, 0])
假設有一棟2層樓，每層樓內的房間都是3排4列，那我們可以用一個三維陣列來儲存每個房間的居住人數（當然，也可以是房間面積等其他數值資訊）。
>>> a = np.arange(24).reshape(2,3,4) # 2層3排4列 >>> a 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]]]) >>> a[1][2][3] # 雖然可以這樣 23 >>> a[1,2,3] # 但這才是規範的用法 23 >>> a[:,0,0] # 所有樓層的第1排第1列 array([ 0, 12]) >>> a[0,:,:] # 1樓的所有房間，等價與a[0]或a[0,...] array([[ 0, 1, 2, 3], [ 4, 5, 6, 7], [ 8, 9, 10, 11]]) >>> a[:,:,1:3] # 所有樓層所有排的第2到4列 array([[[ 1, 2], [ 5, 6], [ 9, 10]], [[13, 14], [17, 18], [21, 22]]]) >>> a[1,:,-1] # 2層每一排的最後一個房間 array([15, 19, 23])
數組合並
數組合併除了下面介紹的水平合併、垂直合併、深度合併外，還有行合併、列合併，以及concatenate()等方式。假如你比我還懶，那就只瞭解前三種方法吧，足夠用了。
>>> a = np.arange(9).reshape(3,3) >>> b = np.arange(9,18).reshape(3,3) >>> a array([[0, 1, 2], [3, 4, 5], [6, 7, 8]]) >>> b array([[ 9, 10, 11], [12, 13, 14], [15, 16, 17]]) >>> np.hstack((a,b)) # 水平合併 array([[ 0, 1, 2, 9, 10, 11], [ 3, 4, 5, 12, 13, 14], [ 6, 7, 8, 15, 16, 17]]) >>> np.vstack((a,b)) # 垂直合併 array([[ 0, 1, 2], [ 3, 4, 5], [ 6, 7, 8], [ 9, 10, 11], [12, 13, 14], [15, 16, 17]]) >>> np.dstack((a,b)) # 深度合併 array([[[ 0, 9], [ 1, 10], [ 2, 11]], [[ 3, 12], [ 4, 13], [ 5, 14]], [[ 6, 15], [ 7, 16], [ 8, 17]]])
陣列拆分
拆分是合併的逆過程，概念是一樣的，但稍微有一點不同：
>>> a = np.arange(9).reshape(3,3) >>> np.hsplit(a, 3) # 水平拆分，返回list [array([[0], [3], [6]]), array([[1], [4], [7]]), array([[2], [5], [8]])] >>> np.vsplit(a, 3) # 垂直拆分，返回list [array([[0, 1, 2]]), array([[3, 4, 5]]), array([[6, 7, 8]])] >>> a = np.arange(27).reshape(3,3,3) >>> np.dsplit(a, 3) # 深度拆分，返回list [array([[[ 0], [ 3], [ 6]], [[ 9], [12], [15]], [[18], [21], [24]]]), array([[[ 1], [ 4], [ 7]], [[10], [13], [16]], [[19], [22], [25]]]), array([[[ 2], [ 5], [ 8]], [[11], [14], [17]], [[20], [23], [26]]])]
陣列運算
陣列和常數的四則運算，是陣列的每一個元素分別和常數運算；陣列和陣列的四則運算則是兩個陣列對應元素的運算（兩個陣列有相同的shape，否則丟擲異常）。
>>> a = np.arange(4, dtype=np.float32).reshape(2,2) >>> b = np.arange(4, 8, dtype=np.float32).reshape(2,2) >>> a+2 # 陣列和常數可以進行四則運算 array([[ 2., 3.], [ 4., 5.]], dtype=float32) >>> a/b # 陣列和陣列可以進行四則運算 array([[ 0. , 0.2 ], [ 0.33333334, 0.42857143]], dtype=float32) >>> a == b # 最神奇的是，陣列可以判斷對應元素是否相等 array([[False, False], [False, False]], dtype=bool) >>> (a == b).all() # 判斷陣列是否相等 False
特別提示：如果想對陣列內符合特定條件的元素做特殊處理，下面的程式碼也許有用。
>>> a = np.arange(6).reshape((2,3)) >>> a array([[0, 1, 2], [3, 4, 5]]) >>> (a>2)&(a<=4) array([[False, False, False], [ True, True, False]], dtype=bool) >>> a[(a>2)&(a<=4)] array([3, 4]) >>> a[(a>2)&((a<=4))] += 10 >>> a array([[ 0, 1, 2], [13, 14, 5]])
陣列方法和常用函式
陣列物件本身提供了計算算數平均值、求最大最小值等內建方法，numpy也提供了很多實用的函式。為了縮減篇幅，下面的程式碼僅以一維陣列為例，展示了這些方法和函式用法。事實上，大多數情況下這些方法和函式對於多維陣列同樣有效，只有少數例外，比如compress函式。
>>> a = np.array([3,2,4]) >>> a.sum() # 所有元素的和 9 >>> a.prod() # 所有元素的乘積 24 >>> a.mean() # 所有元素的算數平均值 3.0 >>> a.max() # 所有元素的最大值 4 >>> a.min() # 所有元素的最小值 2 >>> a.clip(3,4) # 小於3的元素替換為3，大於4的元素替換為4 array([3, 3, 4]) >>> a.compress(a>2) # 返回大於2的元素組成的陣列 array([3, 4]) >>> a.tolist() # 返回python的list [3, 2, 4] >>> a.var() # 計算方差（元素與均值之差的平方的均值） 0.66666666666666663 >>> a.std() # 計算標準差（方差的算術平方根） 0.81649658092772603 >>> a.ptp() # 返回陣列的最大值和最小值之差 2 >>> a.argmin() # 返回最小值在扁平陣列中的索引 1 >>> a.argmax() # 返回最大值在扁平陣列中的索引 2 >>> np.where(a == 2) # 返回所有值為2的元素的索引 (array([1]),) >>> np.diff(a) # 返回相鄰元素的差 array([-1, 2]) >>> np.log(a) # 返回對數陣列 array([ 1.09861229, 0.69314718, 1.38629436]) >>> np.exp(a) # 返回指數陣列 array([ 20.08553692, 7.3890561 , 54.59815003]) >>> np.sqrt(a) # 返回開方陣列 array([ 1.73205081, 1.41421356, 2. ]) >>> np.msort(a) # 陣列排序 array([2, 3, 4]) >>> a = np.array([1,4,7]) >>> b = np.array([8,5,2]) >>> np.maximum(a, b) # 返回多個數組中對應位置元素的最大值陣列 array([8, 5, 7]) >>> np.minimum(a, b) # 返回多個數組中對應位置元素的最小值陣列 array([1, 4, 2]) >>> np.true_divide(a, b) # 對整數實現真正的數學除法運算 array([ 0.125, 0.8 , 3.5 ])
矩陣物件
matrix是矩陣物件，繼承自ndarray型別，因此含有ndarray的所有資料屬性和方法。不過，當你把矩陣物件當陣列操作時，需要注意以下幾點：
- matrix物件總是二維的，即使是展平（ravel函式）操作或是成員選擇，返回值也是二維的
- matrix物件和ndarray物件混合的運算總是返回matrix物件
建立矩陣
matrix物件可以使用一個Matlab風格的字串來建立（以空格分隔列，以分號分隔行的字串），也可以用陣列來建立。
>>> np.mat('1 4 7; 2 5 8; 3 6 9')
matrix([[1, 4, 7],
[2, 5, 8],
[3, 6, 9]])
>>> np.mat(np.arange(1,10).reshape(3,3))
matrix([[1, 2, 3],
[4, 5, 6],
[7, 8, 9]])
矩陣的特有屬性
矩陣有幾個特有的屬性使得計算更加容易，這些屬性有：
>>> m = np.mat(np.arange(1,10).reshape(3,3)) >>> m matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) >>> m.T # 返回自身的轉置 matrix([[1, 4, 7], [2, 5, 8], [3, 6, 9]]) >>> m.H # 返回自身的共軛轉置 matrix([[1, 4, 7], [2, 5, 8], [3, 6, 9]]) >>> m.I # 返回自身的逆矩陣 matrix([[ -4.50359963e+15, 9.00719925e+15, -4.50359963e+15], [ 9.00719925e+15, -1.80143985e+16, 9.00719925e+15], [ -4.50359963e+15, 9.00719925e+15, -4.50359963e+15]]) >>> m.A # 返回自身資料的二維陣列的一個檢視 array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
矩陣乘法
對ndarray物件而言，星號是按元素相乘，dot()函式則當作矩陣相乘。對於matrix物件來說，星號和dot()函式都是矩陣相乘。特別的，對於一維陣列，dot()函式實現的是向量點乘（結果是標量），但星號實現的卻不是差乘。
>>> a = np.array([1,2,3]) >>> b = np.array([4,5,6]) >>> a*b # 一維陣列，元素相乘 array([ 4, 10, 18]) >>> np.dot(a,b) # 一維陣列，元素相乘再求和 32 >>> a = np.array([[1,2],[3,4]]) >>> b = np.array([[5,6],[7,8]]) >>> a*b # 多維陣列，元素相乘 array([[ 5, 12], [21, 32]]) >>> np.dot(a,b) # 多維陣列，實現的是矩陣相乘 array([[19, 22], [43, 50]]) >>> m = np.mat(a) >>> n = np.mat(b) >>> np.dot(m,n) # 矩陣相乘 matrix([[19, 22], [43, 50]]) >>> m*n # 矩陣相乘 matrix([[19, 22], [43, 50]])
線性代數模組
numpy.linalg 是numpy的線性代數模組，可以用來解決逆矩陣、特徵值、線性方程組以及行列式等問題。
計算逆矩陣
儘管matrix物件本身有逆矩陣的屬性，但用numpy.linalg模組求解矩陣的逆，也是非常簡單的。
m = np.mat('0 1 2; 1 0 3; 4 -3 8')
mi = np.linalg.inv(m) # mi即為m的逆矩陣。何以證明？
m * mi # 矩陣與其逆矩陣相乘，結果為單位矩陣
matrix([[ 1., 0., 0.],
[ 0., 1., 0.],
[ 0., 0., 1.]])
計算行列式
如何計算行列式，我早已經不記得了，但手工計算行列式的痛苦，我依然記憶猶新。現在好了，你在手機上都可以用numpy輕鬆搞定（前提是你的手機上安裝了python + numpy）。
m = np.mat('0 1 2; 1 0 3; 4 -3 8')
np.linalg.det(m) # 什麼？這就成了？
2.0
計算特徵值和特徵向量
截至目前，我的工作和特徵值、特徵向量還有沒任何關聯。記錄這一節，純粹是為了我女兒，她正在讀數學專業。
m = np.mat('0 1 2; 1 0 3; 4 -3 8')
>>> np.linalg.eigvals(m) # 計算特徵值
array([ 7.96850246, -0.48548592, 0.51698346])
>>> np.linalg.eig(m) # 返回特徵值及其對應特徵向量的元組
(array([ 7.96850246, -0.48548592, 0.51698346]), matrix([[ 0.26955165, 0.90772191, -0.74373492],
[ 0.36874217, 0.24316331, -0.65468206],
[ 0.88959042, -0.34192476, 0.13509171]]))
求解線性方程組
有線性方程組如下：
x – 2y + z = 02y -8z = 8-4x + 5y + 9z = -9
求解過程如下：
>>> A = np.mat('1 -2 1; 0 2 -8; -4 5 9')
>>> b = np.array([0, 8, -9])
>>> np.linalg.solve(A, b)
array([ 29., 16., 3.]) # x = 29, y = 16, z = 3
三劍客之 Matplotlib
matplotlib 是python最著名的繪相簿，它提供了一整套和Matlab相似的命令API，十分適合互動式地進行製圖。而且也可以方便地將它作為繪圖控制元件，嵌入GUI應用程式中。matplotlib 可以繪製多種形式的圖形包括普通的線圖，直方圖，餅圖，散點圖以及誤差線圖等；可以比較方便的定製圖形的各種屬性比如圖線的型別，顏色，粗細，字型的大小等；它能夠很好地支援一部分 TeX 排版命令，可以比較美觀地顯示圖形中的數學公式。
pylot介紹
Matplotlib 包含了幾十個不同的模組， 如 matlab、mathtext、finance、dates 等，而 pylot 則是我們最常用的繪圖模組，這也是本文介紹的重點。
中文顯示問題的解決方案
有很多方法可以解決此問題，但下面的方法恐怕是最簡單的解決方案了（我只在windows平臺上測試過，其他平臺請看官自測）。如果想了解更多，也可以參考我N年前的一片博文： ofollow,noindex">matplotlib顯示中文的解決方案
>>> from pylab import mpl >>> mpl.rcParams['font.sans-serif'] = ['FangSong'] # 指定預設字型 >>> mpl.rcParams['axes.unicode_minus'] = False # 解決儲存影象時'-'顯示為方塊的問題
繪製最簡單的圖形
>>> import numpy as np >>> import matplotlib.pyplot as plt >>> x = np.arange(0, 2*np.pi, 0.01) >>> y = np.sin(x) >>> plt.plot(x, y) >>> plt.show()
設定標題、座標軸名稱、座標軸範圍
如果你在python的shell中執行下面的程式碼，而shell的預設編碼又不是utf-8的話，中文可能仍然會顯示為亂碼。你可以嘗試著把 u’正弦曲線’ 寫成 ‘正弦曲線’.decode(‘gbk’) 或者 ‘正弦曲線’.decode(‘utf-8’)
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from pylab import mpl
>>> mpl.rcParams['font.sans-serif'] = ['FangSong']
>>> mpl.rcParams['axes.unicode_minus'] = False
>>> x = np.arange(0, 2*np.pi, 0.01)
>>> y = np.sin(x)
>>> plt.plot(x, y)
>>> plt.title(u'正弦曲線', fontdict={'size':20}) # 設定標題
>>> plt.xlabel(u'弧度', fontdict={'size':16}) # 顯示橫軸名稱
>>> plt.ylabel(u'正弦值', fontdict={'size':16}) # 顯示縱軸名稱
>>> plt.axis([-0.1*np.pi, 2.1*np.pi, -1.1, 1.1]) # 設定座標軸範圍
>>> plt.show()
設定點和線的樣式、寬度、顏色
plt.plot函式的呼叫形式如下：
plot(x, y, color='green', linestyle='dashed', linewidth=1, marker='o', markerfacecolor='blue', markersize=6) plot(x, y, c='g', ls='--', lw=1, marker='o', mfc='blue', ms=6)
- color指定線的顏色，可簡寫為“c”。顏色的選項為：
- 藍色： ‘b’ (blue)
- 綠色： ‘g’ (green)
- 紅色： ‘r’ (red)
- 墨綠： ‘c’ (cyan)
- 洋紅： ‘m’ (magenta)
- 黃色： ‘y’ (yellow)
- 黑色： ‘k’ (black)
- 白色： ‘w’ (white)
- 灰度表示： e.g. 0.75 ([0,1]內任意浮點數)
- RGB表示法： e.g. ‘#2F4F4F’ 或 (0.18, 0.31, 0.31)
- linestyle指定線型，可簡寫為“ls”。線型的選項為：
- 實線： ‘-’ (solid line)
- 虛線： ‘–’ (dashed line)
- 虛點線： ‘-.’ (dash-dot line)
- 點線： ‘:’ (dotted line)
- 無： ”或’ ‘或’None’
- linewidth指定線寬，可簡寫為“lw”。
- marker描述資料點的形狀
- 點線： ‘.’
- 點線： ‘o’
- 加號： ‘+
- 叉號： ‘x’
- 上三角： ‘^’
- 上三角： ‘v’
- markerfacecolor指定資料點標記的表面顏色，可 簡寫為“ mfc”。
- markersize指定資料點標記的大小，可 簡寫為“ ms”。
文字標註和圖例
我們分別使用不同的線型、顏色來繪製以10、e、2為基的一組冪函式曲線，演示文字標註和圖例的使用。
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from pylab import mpl
>>> mpl.rcParams['font.sans-serif'] = ['FangSong']
>>> mpl.rcParams['axes.unicode_minus'] = False
>>> x = np.linspace(-4, 4, 200)
>>> f1 = np.power(10, x)
>>> f2 = np.power(np.e, x)
>>> f3 = np.power(2, x)
>>> plt.plot(x, f1, 'r', ls='-', linewidth=2, label='$10^x$')
>>> plt.plot(x, f2, 'b', ls='--', linewidth=2, label='$e^x$')
>>> plt.plot(x, f3, 'g', ls=':', linewidth=2, label='$2^x$')
>>> plt.axis([-4, 4, -0.5, 8])
>>> plt.text(1, 7.5, r'$10^x$', fontsize=16)
>>> plt.text(2.2, 7.5, r'$e^x$', fontsize=16)
>>> plt.text(3.2, 7.5, r'$2^x$', fontsize=16)
>>> plt.title('冪函式曲線', fontsize=16)
>>> plt.legend(loc='upper left')
>>> plt.show()
在繪製圖例時，loc用於指定圖例的位置，可用的選項有：
- best
- upper right
- upper left
- lower left
- lower right
繪製多軸圖
在介紹如何將多幅子圖繪製在同一畫板的同時，順便演示如何繪製直線和矩形。我們可以使用subplot函式快速繪製有多個軸的圖表。subplot函式的呼叫形式如下：
subplot(numRows, numCols, plotNum)
subplot將整個繪圖區域等分為numRows行 * numCols列個子區域，然後按照從左到右，從上到下的順序對每個子區域進行編號，左上的子區域的編號為1。如果numRows，numCols和plotNum這三個數都小於10的話，可以把它們縮寫為一個整數，例如subplot(323)和subplot(3,2,3)是相同的。subplot在plotNum指定的區域中建立一個軸物件。如果新建立的軸和之前建立的軸重疊的話，之前的軸將被刪除。
>>> import matplotlib.pyplot as plt >>> plt.subplot(221) # 兩行兩列的第1個位置 >>> plt.axis([-1, 2, -1, 2]) >>> plt.axhline(y=0.5, color='b') >>> plt.axhline(y=0.5, xmin=0.25, xmax=0.75, color='r') >>> plt.subplot(222) # 兩行兩列的第2個位置 >>> plt.axis([-1, 2, -1, 2]) >>> plt.axvline(x=0, ymin=0, linewidth=4, color='r') >>> plt.axvline(x=1.0, ymin=-0.5, ymax=0.5, linewidth=4, color='g') >>> plt.subplot(212) # 兩行一列的第2個位置 >>> plt.axis([-1, 2, -1, 2]) >>> plt.axvspan(1.25, 1.55, facecolor='g', alpha=0.5) >>> plt.axhspan(0.25, 0.75, facecolor='0.5', alpha=0.5) >>> plt.show()
常用繪圖型別
直方圖
用numpy隨機生成一個符合正態分佈的資料集，統計分段區域內資料的個數。
>>> import numpy as np >>> import matplotlib.pyplot as plt >>> data = np.random.normal(5.0, 3.0, 1000) >>> plt.hist(data) >>> bins = np.arange(-5, 16, 1) >>> plt.hist(data, bins) # 使用自定義的分段區域 >>> plt.show()
散點圖
使用plot()繪圖時，如果指定樣式引數為僅繪製資料點（linestyle=’None’），那麼所繪製的就是一幅雜湊圖。這種方法所繪製的點無法單獨指定資料點的顏色和大小，而使用scatter()繪製雜湊圖就可以指定每個點的顏色和大小。
plt.scatter函式的呼叫形式如下：
scatter(x, y, s=None, c=None, marker=None, cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, verts=None, edgecolors=None, hold=None, data=None, **kwargs)
scatter()的前兩個引數是陣列，分別指定每個點的X軸和Y軸的座標。s引數指定點的大 小，值和點的面積成正比，它可以是一個數，指定所有點的大小，也可以是陣列，分別對每個點指定大小。c引數指定每個點的顏色，可以是數值或陣列。這裡使用一維陣列為每個點指定了一個數值。通過顏色對映表，每個數值都會與一個顏色相對應。預設的顏色對映表中藍色與最小值對應，紅色與最大值對應。當c引數是形狀為(N,3)或(N,4)的二維陣列時，則直接表示每個點的RGB顏色。marker引數設定點的形狀，可以是個表示形狀的字串，也可以是表示多邊形的兩個元素的元組，第一個元素表示多邊形的邊數，第二個元素表示多邊形的樣式，取值範圍為0、1、2、3。0表示多邊形，1表示星形，2表示放射形，3表示忽略邊數而顯示為圓形。alpha引數設定點的透明度。facecolors引數為“none”時，表示雜湊點沒有填充色。
>>> import numpy as np >>> import matplotlib.pyplot as plt >>> x = np.random.rand(50) >>> y = np.random.rand(50) >>> area = np.pi * (15 * np.random.rand(50))**2 >>> color = 2 * np.pi * np.random.rand(50) >>> plt.scatter(x, y, s=area, c=color, alpha=0.5, cmap=plt.cm.hsv) >>> plt.show()
梯形圖、柱狀圖、填充圖
>>> import numpy as np >>> import matplotlib.pyplot as plt >>> n = np.array([0,1,2,3,4,5]) >>> x = np.linspace(-0.75, 1., 100) >>> plt.subplot(131) >>> plt.step(n, n**2, lw=2) >>> plt.subplot(132) >>> plt.bar(n, n**2, align="center", width=0.5, alpha=0.5) >>> plt.subplot(133) >>> plt.fill_between(x, x**2, x**3, color="green", alpha=0.5) >>> plt.show()
對數座標
plot()所繪製圖表的X-Y軸座標都是算術座標。繪製對數座標圖的函式有三個：semilogx()、semilogy()和loglog(),它們分別繪製X軸為對數座標、Y軸為對數座標以及兩個軸都為對數座標時的圖表。
>>> import numpy as np >>> import matplotlib.pyplot as plt >>> x = np.linspace(0, 3, 100) >>> y = np.power(2, np.power(2,x)) >>> plt.subplot(121) >>> plt.semilogy(x, y , '-r') >>> plt.subplot(122) >>> plt.plot(x,y, '--g') >>> plt.show()
極座標繪圖
極座標系是和笛卡爾(X-Y)座標系完全不同的座標系，極座標系中的點由一個夾角和一段相對中心點的距離來表示。polar(theta, r, **kwargs)可以直接創建極座標子圖並在其中繪製曲線。也可以使用程式中呼叫subplot()建立子圖時通過設 polar引數為True,建立一個極座標子圖，然後呼叫plot()在極座標子圖中繪圖。
>>> import numpy as np >>> import matplotlib.pyplot as plt >>> theta = np.arange(0, 2*np.pi, 0.02) >>> plt.polar(theta, 1.4*np.cos(5*theta), "--", linewidth=2) >>> plt.polar(theta, 1.8*np.cos(4*theta), linewidth=2) >>> plt.rgrids(np.arange(0.5, 2, 0.5), angle=45) >>> plt.thetagrids([0, 45])thetagridlabel objects>) >>> plt.show() >>>
2D繪圖
等值線圖
所謂等值線，是指由函式值相等的各點連成的平滑曲線。等值線可以直觀地表示二元函式值的變化趨勢，例如等值線密集的地方表示函式值在此處的變化較大。matplotlib中可以使用contour()和contourf()描繪等值線，它們的區別是：contourf()所得到的是帶填充效果的等值線。
>>> import numpy as np >>> import matplotlib.pyplot as plt >>> y, x = np.ogrid[-2:2:200j, -3:3:300j] >>> z = x * np.exp( - x**2 - y**2) >>> extent = [np.min(x), np.max(x), np.min(y), np.max(y)] >>> plt.subplot(121) >>> cs = plt.contour(z, 10, extent=extent) >>> plt.clabel(cs) <a><a list of 8 text.Text objects> >>> plt.subplot(122) >>> plt.contourf(x.reshape(-1), y.reshape(-1), z, 20) >>> plt.show()</a>
為了更淸楚地區分X軸和Y軸，這裡讓它們的取值範圍和等分次數均不相同.這樣得 到的陣列z的形狀為(200, 300),它的第0軸對應Y軸、第1軸對應X軸。
呼叫contour()繪製陣列z的等值線圖，第二個引數為10,表示將整個函式的取值範圍等分為10個區間,即顯示的等值線圖中將有9條等值線。可以使用extent引數指定等值線圖的X軸和Y軸的資料範圍。
contour()所返回的是一個QuadContourSet物件， 將它傳遞給clabel()，為其中的等值線標上對應的值。
呼叫contourf(),繪製將取值範圍等分為20份、帶填充效果的等值線圖。這裡演示了另外一種設定X、Y軸取值範圍的方法，它的前兩個引數分別是計算陣列z時所使用的X軸和Y軸上的取樣點,這兩個陣列必須是一維的。
二維資料的平面色彩顯示
>>> import numpy as np >>> import matplotlib.pyplot as plt >>> data=np.clip(np.random.randn(5,5),-1,1) >>> plt.subplot(221) >>> plt.imshow(data) >>> plt.subplot(222) >>> plt.imshow(data,cmap=plt.cm.cool) >>> plt.subplot(223) >>> plt.imshow(data,cmap=plt.cm.hot) >>> plt.colorbar() >>> plt.subplot(224) >>> im = plt.imshow(data,cmap=plt.cm.winter) >>> plt.colorbar(im, cmap=plt.cm.winter, ticks=[-1,0,1]) >>> plt.show()
3D繪圖
雖然matplotlib主要專注於繪圖，並且主要是二維的圖形，但是它也有一些不同的擴充套件，能讓我們在地理圖上繪圖，讓我們把Excel和3D圖表結合起來。在matplotlib的世界裡，這些擴充套件叫做工具包（toolkits）。工具包是一些關注在某個話題（如3D繪圖）的特定函式的集合。
比較流行的工具包有Basemap、GTK 工具、Excel工具、Natgrid、AxesGrid和mplot3d。
mpl_toolkits.mplot3工具包提供了一些基本的3D繪圖功能，其支援的圖表型別包括散點圖（scatter）、曲面圖（surf）、線圖（line）和網格圖（mesh）。雖然mplot3d不是一個最好的3D圖形繪製庫，但是它是伴隨著matplotlib產生的，因此我們對其介面已經很熟悉了。
下面是一個使用plot_surface繪製3d曲面圖的例子。
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> import mpl_toolkits.mplot3d
>>> x, y = np.mgrid[-2:2:50j,-2:2:50j]
>>> z = x*np.exp(-x**2-y**2)
>>> ax = plt.subplot(111,projection='3d')
>>> ax.plot_surface(x,y,z,rstride=2,cstride=1,cmap=plt.cm.coolwarm,alpha=0.8)
>>> ax.set_xlabel('x')
>>> ax.set_ylabel('y')
>>> ax.set_zlabel('z')
>>> plt.show()
三劍客之 Scipy
前面已經說過，最初的numpy其實是scipy的一部分，後來才從scipy中分離出來。scipy函式庫在numpy庫的基礎上增加了眾多的數學、科學以及工程計算中常用的庫函式。例如線性代數、常微分方程數值求解、訊號處理、影象處理、稀疏矩陣等等。由於其涉及的領域眾多，我之於scipy，就像盲人摸大象，只能是摸到哪兒算哪兒。
插值
一維插值和二維插值，是我最常用的scipy的功能之一，也是最容易上手的。
一維插值和樣條插值
下面的例子清楚地展示了線性插值和樣條插值之後的資料形態。
>>> import numpy as np >>> from scipy import interpolate >>> x = np.arange(0,10) >>> y = np.exp(-x/3.0) >>> x array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]) >>> y array([ 1. , 0.71653131, 0.51341712, 0.36787944, 0.26359714, 0.1888756 , 0.13533528, 0.09697197, 0.06948345, 0.04978707]) >>> f_line = interpolate.interp1d(x, y) >>> x_new = np.arange(0,9,0.2) >>> f_line(x_new) array([ 1. , 0.94330626, 0.88661252, 0.82991879, 0.77322505, 0.71653131, 0.67590847, 0.63528563, 0.5946628 , 0.55403996, 0.51341712, 0.48430958, 0.45520205, 0.42609451, 0.39698698, 0.36787944, 0.34702298, 0.32616652, 0.30531006, 0.2844536 , 0.26359714, 0.24865283, 0.23370852, 0.21876422, 0.20381991, 0.1888756 , 0.17816754, 0.16745947, 0.15675141, 0.14604335, 0.13533528, 0.12766262, 0.11998996, 0.11231729, 0.10464463, 0.09697197, 0.09147426, 0.08597656, 0.08047886, 0.07498115, 0.06948345, 0.06554417, 0.0616049 , 0.05766562, 0.05372634]) >>> bs = interpolate.splrep(x, y) >>> interpolate.splev(x_new, bs) array([ 1. , 0.93571489, 0.8754193 , 0.8189194 , 0.76602135, 0.71653131, 0.67025545, 0.62699994, 0.58657094, 0.54877461, 0.51341712, 0.48031625, 0.44933621, 0.42035284, 0.39324198, 0.36787944, 0.34414628, 0.32194438, 0.30118082, 0.28176271, 0.26359714, 0.24659576, 0.23068846, 0.2158097 , 0.20189393, 0.1888756 , 0.17669225, 0.16529366, 0.15463272, 0.1446623 , 0.13533528, 0.1266067 , 0.11844022, 0.11080172, 0.10365702, 0.09697197, 0.09071432, 0.0848594 , 0.07938446, 0.07426673, 0.06948345, 0.06501185, 0.06082918, 0.05691267, 0.05323955])
將原始資料以及線性插值和樣條插值之後的資料繪製在一起，效果會比較明顯：
程式碼如下：
import numpy as np from scipy import interpolate import matplotlib.pyplot as plt from pylab import mpl mpl.rcParams['font.sans-serif'] = ['FangSong'] # 指定預設字型 mpl.rcParams['axes.unicode_minus'] = False # 解決儲存影象時'-'顯示為方塊的問題 x = np.arange(0,10) y = np.exp(-x/3.0) x_new = np.arange(0,9,0.2) f_linear = interpolate.interp1d(x, y) bs = interpolate.splrep(x, y) plt.plot(x, y, "o", label=u"原始資料") plt.plot(x_new, f_linear(x_new), label=u"線性插值") plt.plot(x_new, interpolate.splev(x_new, bs), label=u"B-spline插值") plt.legend() plt.show()
特別說明：樣條插值附帶了很多預設引數，下面是簡單的說明。詳情請自行搜尋。
scipy.interpolate.splrep(x, y, w=None, xb=None, xe=None, k=3, task=0, s=None, t=None, full_output=0, per=0, quiet=1) # 引數s用來確定平滑點數，通常是m-SQRT(2m),m是曲線點數。如果在插值中不需要平滑應該設定s=0。splrep()輸出的是一個3元素的元胞陣列（t,c,k）,其中t是曲線點，c是計算出來的係數，k是樣條階數，通常是3階，但可以通過k改變。 scipy.interpolate.splev(x, tck, der=0) # 其中的der是進行樣條計算是需要實際計算到的階數，必須滿足條件der<=k
二維插值
在一個房間的地板上按九宮格的位置放置9個溫度感測器，測得溫度如下：
>>> import numpy as np >>> temp = np.random.randint(20,30,(3,3)) >>> temp array([[21, 29, 21], [24, 27, 20], [25, 28, 22]])
在不增加感測器的前提下，我們採用二維插值的方法，可以使得資料變化較為平滑：
>>> import numpy as np >>> temp = np.random.randint(20,30,(3,3)) >>> temp array([[21, 29, 21], [24, 27, 20], [25, 28, 22]]) >>> x = np.arange(3) >>> y = np.arange(3) >>> ip = interpolate.interp2d(x,y,temp) >>> x_new = np.linspace(0,2,9) >>> x_new = np.linspace(0,2,5) >>> y_new = np.linspace(0,2,5) >>> temp_new = ip(x_new,y_new) >>> temp_new array([[ 21. , 25. , 29. , 25. , 21. ], [ 22.5 , 25.25, 28. , 24.25, 20.5 ], [ 24. , 25.5 , 27. , 23.5 , 20. ], [ 24.5 , 26. , 27.5 , 24.25, 21. ], [ 25. , 26.5 , 28. , 25. , 22. ]])
下圖是根據原始資料和插值資料繪製的該房間溫度平面圖。
擬合
在工作中，我們常常需要在圖中描繪某些實際資料觀察的同時，使用一個曲線來擬合這些實際資料。所謂擬合，就是找出符合資料變化趨勢的曲線方程，或者直接繪製出擬合曲線。
使用numpy.polyfit擬合
下面這段程式碼，基於Numpy模組，可以直接繪製出擬合曲線，但我無法得到曲線方程（儘管輸出了一堆曲線引數）。這是一個值得繼續深入研究的問題。
#! python2 # coding: utf-8 import numpy import pylab def plot_polynomail_fit(y, *deg): ''''' 這個函式一次只擬合一組資料。但是可以對這一組資料同時擬合多條曲線並顯示 y是一個list，存放的是需要擬合的資料 *deg是一個元組，長度不定,裡面存放擬合的次數，可以對一組資料擬合出多條直線進行比較 ''' x = xrange(len(y)) COLOR = ['c','m','y','k','r','p','o','g','b'] temp = [] numOfLineToFit = len(deg) # 需要擬合的次數列表 for index,item in enumerate (deg): param = numpy.polyfit(x,y,item) # 曲線的引數 equation = numpy.poly1d(param) # 曲線方程 temp.extend(param[:]) # 提取曲線引數 print param pylab.subplot(numOfLineToFit,1,index+1) pylab.plot(x, equation(x),'%s--'%COLOR[index]) pylab.plot(x,y,'b--') pylab.show() return temp if __name__ == '__main__': y = [17, 19, 21, 28, 33, 38, 37, 37, 31, 23, 19, 18] plot_polynomail_fit(y, 2, 3, 4)
3個擬合結果顯示在下圖中。
使用scipy.optimize.optimize.curve_fit擬合
scipy提供的擬合，貌似需要先確定帶引數的曲線方程，然後由scipy求解方程，返回曲線引數。我們還是以上面的一組資料為例使用scipy擬合曲線。
>>> import numpy as np >>> import matplotlib.pyplot as plt >>> from scipy import optimize >>> x = np.arange(1,13,1) >>> y = np.array([17,19,21,28,33,38,37,37,31,23,19,18 ]) >>> plt.plot(x, y) [<matplotlib.lines.Line2D object at 0x04799D10>] >>> plt.show()
可以看出，曲線近似正弦函式。構建函式y=a*sin(x*pi/6+b)+c，使用scipy的optimize.curve_fit函式求出a、b、c的值：
>>> def fmax(x,a,b,c): return a*np.sin(x*np.pi/6+b)+c >>> fita, fitb = optimize.curve_fit(fmax, x, y, [1,1,1]) >>> print fita [ 10.93254951 -1.9496096 26.75 ] >>> xn = np.arange(1,13,0.1) >>> plt.plot(x, y) [<matplotlib.lines.Line2D object at 0x04B160B0>] >>> plt.plot(xn, fmax(xn, fita[0],fita[1],fita[2])) [<matplotlib.lines.Line2D object at 0x04B16510>] >>> plt.show()
求解非線性方程（組）
在數學建模中，需要對一些稀奇古怪的方程（組）求解，Matlab自然是首選，但Matlab不是免費的，scipy則為我們提供了免費的午餐！scipy.optimize庫中的fsolve函式可以用來對非線性方程（組）進行求解。它的基本呼叫形式如下：
fsolve(func, x0)
func(x)是計算方程組誤差的函式，它的引數x是一個向量，表示方程組的各個未知數的一組可能解，func返回將x代入方程組之後得到的誤差；x0為未知數向量的初始值。
我們先來求解一個簡單的方程： sin ( x ) − cos ( x ) = 0.2
>>> from scipy.optimize import fsolve >>> import numpy as np >>> def f(A): x = float(A[0]) return [np.sin(x) - np.cos(x) - 0.2] >>> result = fsolve(f, [1]) array([ 0.92729522]) >>> print result [0.92729522] >>> print f(result) [2.7977428707082197e-09]
哈哈，易如反掌！再來一個方程組：
4 x 2 − 2 sin ( y z ) = 0
5 y + 3 = 0
y z − 1.5 = 0
>>> from scipy.optimize import fsolve >>> import numpy as np >>> def f(A): x = float(A[0]) y = float(A[1]) z = float(A[2]) return [4*x*x - 2*np.sin(y*z), 5*y + 3, y*z - 1.5] >>> result = fsolve(f, [1, 1, 1]) >>> print result [-0.70622057 -0.6 -2.5 ] >>> print f(result) [-9.1260332624187868e-14, 0.0, 5.329070518200751e-15]
數值積分
數值積分是對定積分的數值求解，例如可以利用數值積分計算某個形狀的面積。我們知道，半徑為1的圓的方程可寫成：
x 2 + y 2 = 1
下面讓我們來考慮一下如何計算半徑為1的半圓的面積，根據圓的面積公式，其面積應該等於PI/2。單位半圓曲線可以用下面的函式表示：
y = 1 − x 2 −−−−−√
我們先定義一個計算根據x計算y的函式：
>>> def half_circle(x): return (1-x**2)**0.5
經典微分法
下面的程式使用經典的分小矩形計算面積總和的方式，計算出單位半圓的面積：
>>> N = 10000 >>> x = np.linspace(-1, 1, N) >>> dx = 2.0/N >>> y = half_circle(x) >>> dx * np.sum(y[:-1] + y[1:]) # 面積的兩倍 3.1412751679988937
使用定積分求解函式
如果我們呼叫scipy.integrate庫中的quad函式的話，將會得到非常精確的結果：
>>> from scipy import integrate >>> pi_half, err = integrate.quad(half_circle, -1, 1) >>> pi_half*2 3.1415926535897984
影象處理
在scipy.misc模組中，有一個函式可以載入Lena影象——這副影象是被用作影象處理的經典示範影象。我只是簡單展示一下在該影象上的幾個操作。
- 載入Lena影象，並顯示灰度影象
- 應用中值濾波掃描訊號的每一個數據點，並替換為相鄰資料點的中值
- 旋轉影象
- 應用Prewitt濾波器（基於影象強度的梯度計算）
>>> from scipy import misc
>>> from scipy import ndimage
>>> img = misc.lena().astype(np.float32)
>>> plt.subplot(221)
>>> plt.title('Original Image')
>>> plt.imshow(img, cmap=plt.cm.gray)
>>> plt.subplot(222)
>>> plt.title('Median Filter')
>>> filtered = ndimage.median_filter(img, size=(42,42))
>>> plt.imshow(filtered, cmap=plt.cm.gray)
>>> plt.subplot(223)
>>> plt.title('Rotated')
>>> rotated = ndimage.rotate(img, 90)
>>> plt.imshow(rotated, cmap=plt.cm.gray)
>>> plt.subplot(224)
>>> plt.title('Prewitt Filter')
>>> filtered = ndimage.prewitt(img)
>>> plt.imshow(filtered, cmap=plt.cm.gray)
>>> plt.show()
後記
這篇博文自2016年9月初動筆，斷斷續續寫了5個多月。延宕這麼久，除了自身懶惰的原因外，主要是因為MSN這個主題涉及的內容太過繁雜，又極其晦澀，無論怎麼努力，總怕掛一漏萬、貽笑大方。
現在好了，終於寫完了。倘若哪位看官發現了謬誤，請自行修改，順便通知我一聲；若因此文受益而想約飯、約酒，請發郵件至：[email protected]