(點選上方公號，快速關註我們)

作者：許向武

前言

不管是不是巴薩的球迷，只要你喜歡足球，就一定聽說過梅西(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)，每一種資料型別根據佔用記憶體的位元組數又分為多個不同的子型別。常見的資料型別見下表。

建立陣列

通常，我們用np.array()建立陣列。如果僅僅是建立一維陣列，也可以使用np.arange()或者np.linspace()的方法。np.zeros()、np.ones()、np.eye()則可以構造特殊的資料。np.random.randint()和np.random.random()則可以構造隨機數陣列。

構造複雜陣列

很多時候，我們需要從簡單的資料結構，構造出複雜的陣列。例如，用一維的資料生成二維格點。

重覆陣列: tile

一維陣列網格化: meshgrid

指定範圍和分割方式的網格化: mgrid

上面的例子中用到了虛數。構造虛數的方法如下：

>>> complex(2,5)

(2+5j)

陣列的屬性

numpy的陣列物件除了一些常規的屬性外，也有幾個類似轉置、扁平迭代器等看起來更像是方法的屬性。扁平迭代器也許是遍歷多維陣列的一個簡明方法，下麵的程式碼給出了一個例子。

改變陣列維度

numpy陣列的儲存順序和陣列的維度是不相干的，因此改變陣列的維度是非常便捷的操作，除resize()外，這一類操作不會改變所操作的陣列本身的儲存順序。

索引和切片

對於一維陣列的索引和切片，numpy和python的list一樣，甚至更靈活。

假設有一棟2層樓，每層樓內的房間都是3排4列，那我們可以用一個三維陣列來儲存每個房間的居住人數（當然，也可以是房間面積等其他數值資訊）。

陣列合併

陣列合併除了下麵介紹的水平合併、垂直合併、深度合併外，還有行合併、列合併，以及concatenate()等方式。假如你比我還懶，那就只瞭解前三種方法吧，足夠用了。

陣列拆分

拆分是合併的逆過程，概念是一樣的，但稍微有一點不同：

陣列運算

陣列和常數的四則運算，是陣列的每一個元素分別和常數運算；陣列和陣列的四則運算則是兩個陣列對應元素的運算（兩個陣列有相同的shape，否則丟擲異常）。

特別提示：如果想對陣列內符合特定條件的元素做特殊處理，下麵的程式碼也許有用。

陣列方法和常用函式

陣列物件本身提供了計算算數平均值、求最大最小值等內建方法，numpy也提供了很多實用的函式。為了縮減篇幅，下麵的程式碼僅以一維陣列為例，展示了這些方法和函式用法。事實上，大多數情況下這些方法和函式對於多維陣列同樣有效，只有少數例外，比如compress函式。

矩陣物件

matrix是矩陣物件，繼承自ndarray型別，因此含有ndarray的所有資料屬性和方法。不過，當你把矩陣物件當陣列操作時，需要註意以下幾點：

matrix物件總是二維的，即使是展平（ravel函式）操作或是成員選擇，傳回值也是二維的
matrix物件和ndarray物件混合的運算總是傳回matrix物件

建立矩陣

matrix物件可以使用一個Matlab風格的字串來建立（以空格分隔列，以分號分隔行的字串），也可以用陣列來建立。

矩陣的特有屬性

矩陣有幾個特有的屬性使得計算更加容易，這些屬性有：

矩陣乘法

對ndarray物件而言，星號是按元素相乘，dot()函式則當作矩陣相乘。對於matrix物件來說，星號和dot()函式都是矩陣相乘。特別的，對於一維陣列，dot()函式實現的是向量點乘（結果是標量），但星號實現的卻不是差乘。

線性代數模組

numpy.linalg 是numpy的線性代數模組，可以用來解決逆矩陣、特徵值、線性方程組以及行列式等問題。

計算逆矩陣

儘管matrix物件本身有逆矩陣的屬性，但用numpy.linalg模組求解矩陣的逆，也是非常簡單的。

計算行列式

如何計算行列式，我早已經不記得了，但手工計算行列式的痛苦，我依然記憶猶新。現在好了，你在手機上都可以用numpy輕鬆搞定（前提是你的手機上安裝了python + numpy）。

m = np.mat(‘0 1 2; 1 0 3; 4 -3 8’)

np.linalg.det(m) # 什麼？這就成了？

2.0

計算特徵值和特徵向量

截至目前，我的工作和特徵值、特徵向量還有沒任何關聯。記錄這一節，純粹是為了我女兒，她正在讀數學專業。

求解線性方程組

有線性方程組如下：

x – 2y + z = 0
2y -8z = 8
-4x + 5y + 9z = -9

求解過程如下：

三劍客之Matplotlib

matplotlib 是python最著名的繪相簿，它提供了一整套和Matlab相似的命令API，十分適合互動式地進行製圖。而且也可以方便地將它作為繪圖控制元件，嵌入GUI應用程式中。matplotlib 可以繪製多種形式的圖形包括普通的線圖，直方圖，餅圖，散點圖以及誤差線圖等；可以比較方便的定製圖形的各種屬性比如圖線的型別，顏色，粗細，字型的大小等；它能夠很好地支援一部分 TeX 排版命令，可以比較美觀地顯示圖形中的數學公式。

pyplot介紹

Matplotlib 包含了幾十個不同的模組，如 matlab、mathtext、finance、dates 等，而 pyplot 則是我們最常用的繪圖模組，這也是本文介紹的重點。

中文顯示問題的解決方案

有很多方法可以解決此問題，但下麵的方法恐怕是最簡單的解決方案了（我只在windows平臺上測試過，其他平臺請看官自測）。

繪製最簡單的圖形

>>> 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’)

設定點和線的樣式、寬度、顏色

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為基的一組冪函式曲線，演示文字標註和圖例的使用。

在繪製圖例時，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指定的區域中建立一個軸物件。如果新建立的軸和之前建立的軸重疊的話，之前的軸將被刪除。

三劍客之Scipy

前面已經說過，最初的numpy其實是scipy的一部分，後來才從scipy中分離出來。scipy函式庫在numpy庫的基礎上增加了眾多的數學、科學以及工程計算中常用的庫函式。例如線性代數、常微分方程數值求解、訊號處理、影象處理、稀疏矩陣等等。由於其涉及的領域眾多，我之於scipy，就像盲人摸大象，只能是摸到哪兒算哪兒。

插值

一維插值和二維插值，是我最常用的scipy的功能之一，也是最容易上手的。

一維插值和樣條插值

下麵的例子清楚地展示了線性插值和樣條插值之後的資料形態。

將原始資料以及線性插值和樣條插值之後的資料繪製在一起，效果會比較明顯：

程式碼如下：

特別說明：樣條插值附帶了很多預設引數，下麵是簡單的說明。詳情請自行搜尋。

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

擬合

在工作中，我們常常需要在圖中描繪某些實際資料觀察的同時，使用一個曲線來擬合這些實際資料。所謂擬合，就是找出符合資料變化趨勢的曲線方程，或者直接繪製出擬合曲線。

使用numpy.polyfit擬合

下麵這段程式碼，基於Numpy模組，可以直接繪製出擬合曲線，但我無法得到曲線方程（儘管輸出了一堆曲線引數）。這是一個值得繼續深入研究的問題。

3個擬合結果顯示在下圖中。

使用scipy.optimize.optimize.curve_fit擬合

scipy提供的擬合，貌似需要先確定帶引數的曲線方程，然後由scipy求解方程，傳回曲線引數。我們還是以上面的一組資料為例使用scipy擬合曲線。

可以看出，曲線近似正弦函式。構建函式y=a*sin(x*pi/6+b)+c，使用scipy的optimize.curve_fit函式求出a、b、c的值：

求解非線性方程（組）

在數學建模中，需要對一些稀奇古怪的方程（組）求解，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×2−2sin(yz)=0

5y+3=0

yz−1.5=0

影象處理

在scipy.misc模組中，有一個函式可以載入Lena影象——這副影象是被用作影象處理的經典示範影象。我只是簡單展示一下在該影象上的幾個操作。

載入Lena影象，並顯示灰度影象
應用中值濾波掃描訊號的每一個資料點，並替換為相鄰資料點的中值
旋轉影象
應用Prewitt濾波器（基於影象強度的梯度計算）

後記

這篇博文自2016年9月初動筆，斷斷續續寫了5個多月。延宕這麼久，除了自身懶惰的原因外，主要是因為MSN這個主題涉及的內容太過繁雜，又極其晦澀，無論怎麼努力，總怕掛一漏萬、貽笑大方。

現在好了，終於寫完了。倘若哪位看官發現了謬誤，請自行修改，順便通知我一聲；若因此文受益而想約飯、約酒，請發郵件至：

xufive@gmail.com

【關於作者】

許向武：山東遠思資訊科技有限公司CEO，網名牧碼人，齊國土著，太公之後。少小離家，獨闖江湖，後歸隱於華不註山。素以敲擊鍵盤為業，偶爾遊戲於各網路對局室，擅長送財送分，深為眾棋友所喜聞樂見。

【關於投稿】

如果大家有原創好文投稿，請直接給公號傳送留言。

① 留言格式：
【投稿】+《文章標題》+ 文章連結

② 示例：
【投稿】《不要自稱是程式員，我十多年的 IT 職場總結》：http://blog.jobbole.com/94148/

③ 最後請附上您的個人簡介哈~

看完本文有收穫？請轉發分享給更多人

關註「資料分析與開發」，提升資料技能

數學建模三劍客 MSN

前言