數(shù)學(xué)建模三劍客之Numpy

不管是不是巴薩的球迷，只要你喜歡足球，就一定聽說過梅西(Messi)、蘇亞雷斯(Suarez)和內(nèi)馬爾(Neymar)這個(gè)MSN組合。在眾多的數(shù)學(xué)建模輔助工具中，也有一個(gè)犀利無比的MSN組合，他們就是python麾下大名鼎鼎的 Matplotlib + Scipy + Numpy三劍客。

本文是我整理的MSN學(xué)習(xí)筆記，有些理解可能比較膚淺，甚至是錯(cuò)誤的。如果因此誤導(dǎo)了某位看官，在工作中造成重大失誤或損失，我頂多只能賠償一頓飯——還得是我們樓下的十元盒飯。特此聲明。

文中代碼均從我的這臺(tái)時(shí)不時(shí)出點(diǎn)問題、鬧個(gè)情緒的Yoga 3 pro上復(fù)制而來，這意味著所有的代碼均可在下面的運(yùn)行環(huán)境中順利運(yùn)行：

pyhton 2.7.8

numpy 1.11.1

scipy 0.16.1

matplotlib 1.5.1

三劍客之Numpy

numpy是一個(gè)開源的python科學(xué)計(jì)算庫，包含了很多實(shí)用的數(shù)學(xué)函數(shù)，涵蓋線性代數(shù)、傅里葉變換和隨機(jī)數(shù)生成等功能。最初的numpy其實(shí)是scipy的一部分，后來才從scipy中分離出來。

numpy不是python的標(biāo)準(zhǔn)庫，需要單獨(dú)安裝。假定你的運(yùn)行環(huán)境已經(jīng)安裝了python包管理工具pip，numpy的安裝就非常簡單：

pip install numpy

數(shù)組對(duì)象

ndarray是多維數(shù)組對(duì)象，也是numpy最核心的對(duì)象。在numpy中，數(shù)組的維度(dimensions)叫做軸(axes)，軸的個(gè)數(shù)叫做秩(rank)。通常，一個(gè)numpy數(shù)組的所有元素都是同一種類型的數(shù)據(jù)，而這些數(shù)據(jù)的存儲(chǔ)和數(shù)組的形式無關(guān)。

下面的例子，創(chuàng)建了一個(gè)三維的數(shù)組（在導(dǎo)入numpy時(shí)，一般都簡寫成np）。

import numpy asnp

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

數(shù)據(jù)類型

numpy支持的數(shù)據(jù)類型主要有布爾型(bool)、整型(integrate)、浮點(diǎn)型(float)和復(fù)數(shù)型(complex)，每一種數(shù)據(jù)類型根據(jù)占用內(nèi)存的字節(jié)數(shù)又分為多個(gè)不同的子類型。常見的數(shù)據(jù)類型見下表。

創(chuàng)建數(shù)組

通常，我們用np.array()創(chuàng)建數(shù)組。如果僅僅是創(chuàng)建一維數(shù)組，也可以使用np.arange()或者np.linspace()的方法。np.zeros()、np.ones()、np.eye()則可以構(gòu)造特殊的數(shù)據(jù)。np.random.randint()和np.random.random()則可以構(gòu)造隨機(jī)數(shù)數(shù)組。

構(gòu)造復(fù)雜數(shù)組

很多時(shí)候，我們需要從簡單的數(shù)據(jù)結(jié)構(gòu)，構(gòu)造出復(fù)雜的數(shù)組。例如，用一維的數(shù)據(jù)生成二維格點(diǎn)。

重復(fù)數(shù)組: tile

一維數(shù)組網(wǎng)格化: meshgrid

指定范圍和分割方式的網(wǎng)格化: mgrid

上面的例子中用到了虛數(shù)。構(gòu)造虛數(shù)的方法如下：

>>> complex(2,5)

(2+5j)

數(shù)組的屬性

numpy的數(shù)組對(duì)象除了一些常規(guī)的屬性外，也有幾個(gè)類似轉(zhuǎn)置、扁平迭代器等看起來更像是方法的屬性。扁平迭代器也許是遍歷多維數(shù)組的一個(gè)簡明方法，下面的代碼給出了一個(gè)例子。

改變數(shù)組維度

numpy數(shù)組的存儲(chǔ)順序和數(shù)組的維度是不相干的，因此改變數(shù)組的維度是非常便捷的操作，除resize()外，這一類操作不會(huì)改變所操作的數(shù)組本身的存儲(chǔ)順序。

索引和切片

對(duì)于一維數(shù)組的索引和切片，numpy和python的list一樣，甚至更靈活。

假設(shè)有一棟2層樓，每層樓內(nèi)的房間都是3排4列，那我們可以用一個(gè)三維數(shù)組來保存每個(gè)房間的居住人數(shù)（當(dāng)然，也可以是房間面積等其他數(shù)值信息）。

數(shù)組合并

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

數(shù)組拆分

拆分是合并的逆過程，概念是一樣的，但稍微有一點(diǎn)不同：

數(shù)組運(yùn)算

數(shù)組和常數(shù)的四則運(yùn)算，是數(shù)組的每一個(gè)元素分別和常數(shù)運(yùn)算；數(shù)組和數(shù)組的四則運(yùn)算則是兩個(gè)數(shù)組對(duì)應(yīng)元素的運(yùn)算（兩個(gè)數(shù)組有相同的shape，否則拋出異常）。

特別提示：如果想對(duì)數(shù)組內(nèi)符合特定條件的元素做特殊處理，下面的代碼也許有用。

數(shù)組方法和常用函數(shù)

數(shù)組對(duì)象本身提供了計(jì)算算數(shù)平均值、求最大最小值等內(nèi)置方法，numpy也提供了很多實(shí)用的函數(shù)。為了縮減篇幅，下面的代碼僅以一維數(shù)組為例，展示了這些方法和函數(shù)用法。事實(shí)上，大多數(shù)情況下這些方法和函數(shù)對(duì)于多維數(shù)組同樣有效，只有少數(shù)例外，比如compress函數(shù)。

矩陣對(duì)象

matrix是矩陣對(duì)象，繼承自ndarray類型，因此含有ndarray的所有數(shù)據(jù)屬性和方法。不過，當(dāng)你把矩陣對(duì)象當(dāng)數(shù)組操作時(shí)，需要注意以下幾點(diǎn)：

matrix對(duì)象總是二維的，即使是展平（ravel函數(shù)）操作或是成員選擇，返回值也是二維的

matrix對(duì)象和ndarray對(duì)象混合的運(yùn)算總是返回matrix對(duì)象

創(chuàng)建矩陣

matrix對(duì)象可以使用一個(gè)Matlab風(fēng)格的字符串來創(chuàng)建（以空格分隔列，以分號(hào)分隔行的字符串），也可以用數(shù)組來創(chuàng)建。

矩陣的特有屬性

矩陣有幾個(gè)特有的屬性使得計(jì)算更加容易，這些屬性有：

矩陣乘法

對(duì)ndarray對(duì)象而言，星號(hào)是按元素相乘，dot()函數(shù)則當(dāng)作矩陣相乘。對(duì)于matrix對(duì)象來說，星號(hào)和dot()函數(shù)都是矩陣相乘。特別的，對(duì)于一維數(shù)組，dot()函數(shù)實(shí)現(xiàn)的是向量點(diǎn)乘（結(jié)果是標(biāo)量），但星號(hào)實(shí)現(xiàn)的卻不是差乘。

線性代數(shù)模塊

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

計(jì)算逆矩陣

盡管matrix對(duì)象本身有逆矩陣的屬性，但用numpy.linalg模塊求解矩陣的逆，也是非常簡單的。

計(jì)算行列式

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

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

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

2.0

計(jì)算特征值和特征向量

截至目前，我的工作和特征值、特征向量還有沒任何關(guān)聯(lián)。記錄這一節(jié)，純粹是為了我女兒，她正在讀數(shù)學(xué)專業(yè)。

求解線性方程組

有線性方程組如下：

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

求解過程如下：

三劍客之Matplotlib

matplotlib 是python最著名的繪圖庫，它提供了一整套和Matlab相似的命令A(yù)PI，十分適合交互式地進(jìn)行制圖。而且也可以方便地將它作為繪圖控件，嵌入GUI應(yīng)用程序中。matplotlib 可以繪制多種形式的圖形包括普通的線圖，直方圖，餅圖，散點(diǎn)圖以及誤差線圖等；可以比較方便的定制圖形的各種屬性比如圖線的類型，顏色，粗細(xì)，字體的大小等；它能夠很好地支持一部分 TeX 排版命令，可以比較美觀地顯示圖形中的數(shù)學(xué)公式。

pyplot介紹

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

中文顯示問題的解決方案

有很多方法可以解決此問題，但下面的方法恐怕是最簡單的解決方案了（我只在windows平臺(tái)上測(cè)試過，其他平臺(tái)請(qǐng)看官自測(cè)）。

繪制最簡單的圖形

>>> import numpy asnp

>>> import matplotlib.pyplot asplt

>>> x = np.arange(0,2*np.pi,0.01)

>>> y = np.sin(x)

>>> plt.plot(x,y)

>>> plt.show()

設(shè)置標(biāo)題、坐標(biāo)軸名稱、坐標(biāo)軸范圍

如果你在python的shell中運(yùn)行下面的代碼，而shell的默認(rèn)編碼又不是utf-8的話，中文可能仍然會(huì)顯示為亂碼。你可以嘗試著把 u’正弦曲線’ 寫成 ‘正弦曲線’.decode(‘gbk’)或者‘正弦曲線’.decode(‘utf-8’)

設(shè)置點(diǎn)和線的樣式、寬度、顏色

plt.plot函數(shù)的調(diào)用形式如下：

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”。顏色的選項(xiàng)為：

藍(lán)色： ‘b’ (blue)

綠色： ‘g’ (green)

紅色： ‘r’ (red)

墨綠： ‘c’ (cyan)

洋紅： ‘m’ (magenta)

黃色： ‘y’ (yellow)

黑色： ‘k’ (black)

白色： ‘w’ (white)

灰度表示： e.g. 0.75 ([0,1]內(nèi)任意浮點(diǎn)數(shù))

RGB表示法： e.g. ‘#2F4F4F’ 或 (0.18, 0.31, 0.31)

linestyle指定線型，可簡寫為“l(fā)s”。線型的選項(xiàng)為：

實(shí)線： ‘-’ (solid line)

虛線： ‘–’ (dashed line)

虛點(diǎn)線： ‘-.’ (dash-dot line)

點(diǎn)線： ‘:’ (dotted line)

無： ”或’ ‘或’None’

linewidth指定線寬，可簡寫為“l(fā)w”。

marker描述數(shù)據(jù)點(diǎn)的形狀

點(diǎn)線： ‘.’

點(diǎn)線： ‘o’

加號(hào)： ‘+

叉號(hào)： ‘x’

上三角： ‘^’

上三角： ‘v’

markerfacecolor指定數(shù)據(jù)點(diǎn)標(biāo)記的表面顏色，可 簡寫為“ mfc”。

markersize指定數(shù)據(jù)點(diǎn)標(biāo)記的大小，可 簡寫為“ ms”。

文本標(biāo)注和圖例

我們分別使用不同的線型、顏色來繪制以10、e、2為基的一組冪函數(shù)曲線，演示文本標(biāo)注和圖例的使用。

在繪制圖例時(shí)，loc用于指定圖例的位置，可用的選項(xiàng)有：

best

upper right

upper left

lower left

lower right

繪制多軸圖

在介紹如何將多幅子圖繪制在同一畫板的同時(shí)，順便演示如何繪制直線和矩形。我們可以使用subplot函數(shù)快速繪制有多個(gè)軸的圖表。subplot函數(shù)的調(diào)用形式如下：

subplot(numRows, numCols, plotNum)

subplot將整個(gè)繪圖區(qū)域等分為numRows行 * numCols列個(gè)子區(qū)域，然后按照從左到右，從上到下的順序?qū)γ總€(gè)子區(qū)域進(jìn)行編號(hào)，左上的子區(qū)域的編號(hào)為1。如果numRows，numCols和plotNum這三個(gè)數(shù)都小于10的話，可以把它們縮寫為一個(gè)整數(shù)，例如subplot(323)和subplot(3,2,3)是相同的。subplot在plotNum指定的區(qū)域中創(chuàng)建一個(gè)軸對(duì)象。如果新創(chuàng)建的軸和之前創(chuàng)建的軸重疊的話，之前的軸將被刪除。

三劍客之Scipy

前面已經(jīng)說過，最初的numpy其實(shí)是scipy的一部分，后來才從scipy中分離出來。scipy函數(shù)庫在numpy庫的基礎(chǔ)上增加了眾多的數(shù)學(xué)、科學(xué)以及工程計(jì)算中常用的庫函數(shù)。例如線性代數(shù)、常微分方程數(shù)值求解、信號(hào)處理、圖像處理、稀疏矩陣等等。由于其涉及的領(lǐng)域眾多，我之于scipy，就像盲人摸大象，只能是摸到哪兒算哪兒。

插值

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

一維插值和樣條插值

下面的例子清楚地展示了線性插值和樣條插值之后的數(shù)據(jù)形態(tài)。

將原始數(shù)據(jù)以及線性插值和樣條插值之后的數(shù)據(jù)繪制在一起，效果會(huì)比較明顯：

代碼如下：

特別說明：樣條插值附帶了很多默認(rèn)參數(shù)，下面是簡單的說明。詳情請(qǐng)自行搜索。

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)

# 參數(shù)s用來確定平滑點(diǎn)數(shù)，通常是m-SQRT(2m),m是曲線點(diǎn)數(shù)。如果在插值中不需要平滑應(yīng)該設(shè)定s=0。splrep()輸出的是一個(gè)3元素的元胞數(shù)組（t,c,k）,其中t是曲線點(diǎn)，c是計(jì)算出來的系數(shù)，k是樣條階數(shù)，通常是3階，但可以通過k改變。

scipy.interpolate.splev(x,tck,der=0)

# 其中的der是進(jìn)行樣條計(jì)算是需要實(shí)際計(jì)算到的階數(shù)，必須滿足條件der<=k

擬合

在工作中，我們常常需要在圖中描繪某些實(shí)際數(shù)據(jù)觀察的同時(shí)，使用一個(gè)曲線來擬合這些實(shí)際數(shù)據(jù)。所謂擬合，就是找出符合數(shù)據(jù)變化趨勢(shì)的曲線方程，或者直接繪制出擬合曲線。

使用numpy.polyfit擬合

下面這段代碼，基于Numpy模塊，可以直接繪制出擬合曲線，但我無法得到曲線方程（盡管輸出了一堆曲線參數(shù)）。這是一個(gè)值得繼續(xù)深入研究的問題。

3個(gè)擬合結(jié)果顯示在下圖中。

使用scipy.optimize.optimize.curve_fit擬合

scipy提供的擬合，貌似需要先確定帶參數(shù)的曲線方程，然后由scipy求解方程，返回曲線參數(shù)。我們還是以上面的一組數(shù)據(jù)為例使用scipy擬合曲線。

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

求解非線性方程（組）

在數(shù)學(xué)建模中，需要對(duì)一些稀奇古怪的方程（組）求解，Matlab自然是首選，但Matlab不是免費(fèi)的，scipy則為我們提供了免費(fèi)的午餐！scipy.optimize庫中的fsolve函數(shù)可以用來對(duì)非線性方程（組）進(jìn)行求解。它的基本調(diào)用形式如下：

fsolve(func, x0)

func(x)是計(jì)算方程組誤差的函數(shù)，它的參數(shù)x是一個(gè)矢量，表示方程組的各個(gè)未知數(shù)的一組可能解，func返回將x代入方程組之后得到的誤差；x0為未知數(shù)矢量的初始值。

我們先來求解一個(gè)簡單的方程：sin(x)?cos(x)=0.2

>>> from scipy.optimize import fsolve

>>> import numpy asnp

>>> deff(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]

>>> printf(result)

[2.7977428707082197e-09]

哈哈，易如反掌！再來一個(gè)方程組：

4x2?2sin(yz)=0

5y+3=0

yz?1.5=0

圖像處理

在scipy.misc模塊中，有一個(gè)函數(shù)可以載入Lena圖像——這副圖像是被用作圖像處理的經(jīng)典示范圖像。我只是簡單展示一下在該圖像上的幾個(gè)操作。

載入Lena圖像，并顯示灰度圖像

應(yīng)用中值濾波掃描信號(hào)的每一個(gè)數(shù)據(jù)點(diǎn)，并替換為相鄰數(shù)據(jù)點(diǎn)的中值

旋轉(zhuǎn)圖像

應(yīng)用Prewitt濾波器（基于圖像強(qiáng)度的梯度計(jì)算）

后記

這篇博文自2016年9月初動(dòng)筆，斷斷續(xù)續(xù)寫了5個(gè)多月。延宕這么久，除了自身懶惰的原因外，主要是因?yàn)镸SN這個(gè)主題涉及的內(nèi)容太過繁雜，又極其晦澀，無論怎么努力，總怕掛一漏萬、貽笑大方。