數學建模三劍客 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),每一種數據類型根據佔用內存的位元組數又分為多個不同的子類型。常見的數據類型見下表。
創建數組
通常,我們用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
])
>>>
result
[
0.92729522
]
>>>
f
(
result
)
[
2.7977428707082197e
-
09
]
哈哈,易如反掌!再來一個方程組:
4x2?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/
③ 最後請附上您的個人簡介哈~
看完本文有收穫?請轉
發分享給更多人
關注「P
ython開發者」,提升Python技能
![](https://pic.pimg.tw/zzuyanan/1488615166-1259157397.png)
![](https://pic.pimg.tw/zzuyanan/1482887990-2595557020.jpg)
※主辦方很無奈:Linus 搞錯內核維護者峰會的時間地點
※API Star:一個 Python 3 的 API 框架
TAG:Python開發者 |