4.1簡介
4.2統計假設與檢驗 stats包
4.3信號特徵
4.4尋優
4.5求解
4.6曲線擬合 curve-fit
4.7插值
4.8模式聚類
Scipy是一個高級的科學計算庫，它和Numpy聯繫很密切，Scipy一般都是操控Numpy數組來進行科學計算，Scipy讓Python成為了半個MATLAB。Scipy包含的功能有最優化、線性代數、積分、插值、擬合、特殊函數、快速傅立葉變換、信號處理和圖像處理、常微分方程求解和其他科學與工程中常用的計算，而這些功能都是我們在之後進行數據分析需要的。
Scipy依賴於Numpy庫，因此安裝Scipy時應先安裝Numpy庫，Scipy安裝與其他庫一樣，可通過pip install Scipy安裝，也可以自行下載原始碼，然後用pip install 路徑+文件名全稱（包括.後綴文件名）進行安裝，源碼下載連結：https://pypi.python.org/pypi/scipy/1.0.0，選擇對應版本下載即可。
scipy.io模塊提供了多種功能來解決不同格式的文件輸入和輸出，包括Matlab，Wave，Arff，Matrix Market等等，最常見的是Matlab格式的。
（1）導入相關依賴包
from scipy import io as spioimport numpy as np（2）創建一個數組並保存到文件中
a = np.ones((3,3))spio.savemat('file.mat',{'a':a})data = spio.loadmat('file.mat', struct_as_record=True)data['a']載入txt文件：
numpy.loadtxt()/numpy.savetxt()
智能導入文本/csv文件：
numpy.genfromtxt()/numpy.recfromcsv()
高速，有效率但numpy特有的二進位格式：
numpy.save()/numpy.load()
stats提供了產生連續性分布的函數：均勻分布(uniform)、正態分布（norm）、貝塔分布（beta）；
產生離散分布的函數：伯努利分布（bernoulli）、幾何分布（geom）泊松分布 poisson。
使用時，調用分布的rvs方法即可。
例如：
import scipy.stats as statsx=stats.uniform.rvs(size=20)
y=stats.beta.rvs(size=20,a=3,b=4)
z=stats.norm.rvs(size=20,loc=0,scale=1)
x=stats.poisson.rvs(0.6,loc=0,size=20)import numpy as npimport scipy.stats as statsnormDist=stats.norm(loc=2.5,scale=0.5)z=normDist.rvs(size=400)mean=np.mean(z)med=np.median(z)dev=np.std(z)print('mean=',mean,' med=',med,' dev=',dev)設z是實驗獲得的數據，如何驗證它是否是正態分布的？
statVal, pVal = stats.kstest(z,'norm',(mean,dev))print('pVal=',pVal)計算得到
mean= 2.4788577873926516
med= 2.473711775463392
dev= 0.5343868193748538
pVal= 0.8908805240503633
結論：我們接受假設，即z數據是服從正態分布的
低頻的原始信號，加高頻的噪聲，如何去掉噪聲？
快速傅立葉變換 FFT應用範圍非常廣，在物理學、化學、電子通訊、信號處理、概率論、統計學、密碼學、聲學、光學、海洋學、結構動力學等。原理是將時域中的測量信號轉換到頻域，然後再將得到的頻域信號合成為時域中的信號；時域信號轉換為頻域信號時，將信號分解成幅值譜，顯示與頻率對應的幅值大小。一個信號是由多種頻率的信號合成的，將這些信號分離，就能去掉信號中的噪聲。
Scipy類庫中提供了快速傅立葉變換包fftpack。
接下來我們通過一個FFT應用案例—產生帶噪聲的信號來進一步學習這個包吧~
import numpy as npimport matplotlib.pyplot as pltfrom scipy import fftpack as ffttimeStep = 0.02 timeVec = np.arange(0, 20, timeStep) sig = np.sin( np.pi / 5 * timeVec) + 0.1 * np.random.randn(timeVec.size) plt.plot(timeVec, sig)plt.show()n=sig.sizesig_fft = fft.fft(sig) sampleFreq = fft.fftfreq(n, d=timeStep) pidxs = np.where(sampleFreq > 0) freqs = sampleFreq[pidxs] power = np.abs(sig_fft)[pidxs] freq = freqs[power.argmax()] sig_fft[np.abs(sampleFreq) > freq] = 0 main_sig = fft.ifft(sig_fft) plt.plot(timeVec, main_sig, linewidth=3)f(x)=x2-4*x+8 在x=2的位置，函數有最小值4。
這個題可以通過scipy.optimize包提供的求極值的函數fmin來計算：
from scipy.optimize import fminimport numpy as npdef f(x): return x**2-4*x+8print (fmin(f, 0))Optimization terminated successfully.
Current function value: 4.000000
Iterations: 27
Function evaluations: 54
這個函數不僅可以進行單維尋優，也可以進行多維尋優哦~例如，求解z=x2+y2+8的最小值
from scipy.optimize import fmindef myfunc(p): x,y=p return x**2+y**2+8p=(1,1)print (fmin(myfunc, p ))Optimization terminated successfully.
Current function value: 8.000000
Iterations: 38
Function evaluations: 69
[-2.10235293e-05 2.54845649e-05]
from scipy import optimizeimport numpy as npdef f(x): return x**2 + 20 * np.sin(x)ans = optimize.fsolve(f, -4)print(ans)print(f(ans))結果為
[-2.75294663]
[1.687539e-14]
求解非線性方程組
scipy. optimize的fsolve函數也可以方便用於求解非線性方程組。
求解原則：將方程組的右邊全部規劃為0，將方程組定義為如下格式的python函數：
def f(x): x1,x2,…, xn=x return [f1(x1, x2,…, xn), f2(x1,x2,…, xn),….]from scipy.optimize import fsolvefrom math import *def f(x): x0, x1,x2 = x return [2*x1+3, 4*x0*x0 + sin(x1*x2), x1*x2/2 - 3]ans = fsolve(f, [1.0,1.0,1.0])print (ans)print (f(ans))結果為
[-0.26429884 -1.5 -4. ]
[0.0, 1.1482453876610066e-10, 6.4002136923591024e-12]
給定的y=ax+b函數上的一系列採樣點，並在這些採樣點上增加一些噪聲，然後利用scipy optimize包中提供的curve_fit方法，求解係數a和b
from scipy import optimizeimport matplotlib.pyplot as pltimport numpy as npdef f(x,a,b): return a*x + bx = np.linspace(-10, 10, 30)y = f(x,2,1)ynew= y+ 3*np.random.normal(size=x.size) popt, pcov = optimize.curve_fit(f,x,ynew)print(popt)print (pcov)plt.plot(x,y,color='r',label='orig')plt.plot(x,popt[0]*x+popt[1],color='b',label='fitting')plt.legend(loc='upper left')plt.scatter(x,ynew)plt.show()結果為
popt=[2.06246704 1.65798483]
pcov=[[6.31189436e-03 1.45358283e-10]
[1.45358283e-10 2.24906575e-01]]
popt列表包含每個參數的擬合值，此例求得的a為2.06246704，b為1.65798483。pcov列表的對角元素是每個參數的方差。通過方差，可以評判擬合的質量，方差越小，擬合越可靠。
根據現有的試驗點值，去預測中間的點值
採用線性、兩次樣條、三次樣條插值
我們同樣通過一個案例來學習
首先在[0,10π]區間裡等間距產生了20個採樣點，並計算其sin值，模擬試驗採集得到的20個點，採用線性、兩次樣條、三次樣條插值函數插值擬合原函數sin(x)
import numpy as npfrom scipy.interpolate import interp1dimport matplotlib.pyplot as plt
x=np.linspace(0,10*np.pi,20) y=np.sin(x)
fl=interp1d(x,y,kind='linear') fq=interp1d(x,y,kind='quadratic') fc=interp1d(x,y,kind='cubic') xint=np.linspace(x.min(),x.max(),1000) yintl=fl(xint) yintq=fq(xint) yintc=fc(xint) plt.plot(xint,yintl,color='r', linestyle='--',label='linear')plt.plot(xint,yintq,color='b' ,label='quadr')plt.plot(xint,yintc,color='b' ,linestyle='-.',label='cubic')plt.legend()plt.show()Scipy聚類分析中主要提供了矢量量化分析和系統聚類法
import numpy as npfrom scipy.cluster import vqimport matplotlib.pyplot as pltclass1=np.random.randn(30,2)+10 class2=np.random.randn(40,2)-10 class3=np.random.randn(50,2) data=np.vstack([class1,class2,class3]) centroid, var=vq.kmeans(data,3) key,distance=vq.vq(data,centroid) vqclass1=data[key==0] vqclass2=data[key==1]vqclass3=data[key==2]
plt.scatter(vqclass1[:,0],vqclass1[:,1],marker='o', color="r",label='c1') # 為類0 製圖plt.scatter(vqclass2[:,0],vqclass2[:,1],marker='1', color="g" ,label='c2')plt.scatter(vqclass3[:,0],vqclass3[:,1],marker='2', color="b",label='c3')plt.legend(loc='upper left')plt.show()Scipy的簡單介紹及應用到這裡就結束啦，能看到最後的大家都棒棒噠！總的來說，Scipy庫的強大不是這一篇推文可以概括完的，但還是希望這篇推文能給正在學習python的你帶來幫助！本期內容結束，感想閱讀和關注~
本期作者：張怡
本期編輯校對：李嘉楠
長按關注 「 數據皮皮俠 」
Python的大門為你敞開