第11講--scipy-數(shù)據(jù)處理應用

1、科學計算包SciPY及應用第11講Scipy簡介解決python科學計算而編寫的一組程序包快速實現(xiàn)相關的數(shù)據(jù)處理如以前的課程中的積分Scipy提供的數(shù)據(jù)I/O相比numpy，scipy提供了更傻瓜式的操作方式二進制存儲from scipy import io as fioimport numpy as npx=np.ones(3,2)y=np.ones(5,5)fio.savemat(rd:111.mat,mat1:x,mat2:y)data=fio.loadmat(rd:111.mat,struct_as_record=True)datamat1Scipy的IOdatamat1array(

2、1., 1., 1., 1., 1., 1.) datamat2array( 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.)統(tǒng)計假設與統(tǒng)計假設與檢驗檢驗 stats包包stats提供了產(chǎn)生連續(xù)性分布的函數(shù)提供了產(chǎn)生連續(xù)性分布的函數(shù)，均勻分布均勻分布(uniform)、正態(tài)分布正態(tài)分布（norm）、貝塔貝塔分布（分布（beta）；產(chǎn)生產(chǎn)生離散分布的函數(shù)離散分布的函數(shù)，伯努利分布（伯努利分布（bernoulli）、幾何分布幾何分布（geom）泊松分

3、布泊松分布 poisson使用時，調用分布的使用時，調用分布的rvs方法即可統(tǒng)計假設與統(tǒng)計假設與檢驗檢驗 stats包包import scipy.stats as statsx=stats.uniform.rvs(size=20)#產(chǎn)生產(chǎn)生20個在個在0，1均勻分布的隨機數(shù)均勻分布的隨機數(shù)y=stats.beta.rvs(size=20,a=3,b=4)產(chǎn)生產(chǎn)生20個服從參數(shù)個服從參數(shù)a=3，b=4的貝塔分布的貝塔分布隨機數(shù)隨機數(shù)z=stats.norm.rvs(size=20,loc=0,scale=1)產(chǎn)生了產(chǎn)生了20個服從個服從0，1正態(tài)分布的正態(tài)分布的隨機數(shù)隨機數(shù)x=stats.poi

4、sson.rvs(0.6,loc=0,size=20)產(chǎn)生產(chǎn)生poisson分布分布假設檢驗假設檢驗假設給定的樣本服從某種假設給定的樣本服從某種分布分布，如何驗證？，如何驗證？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是實驗獲得的數(shù)據(jù)，如何驗證它是否是正態(tài)分布

5、的？是實驗獲得的數(shù)據(jù)，如何驗證它是否是正態(tài)分布的？假設檢驗假設檢驗假設給定的樣本服從某種假設給定的樣本服從某種分布分布，如何驗證？，如何驗證？statVal, pVal = stats.kstest(z,norm,(mean,dev)print(pVal=,pVal)計算得到：計算得到：pVal= 0.667359687999結論：結論：我們我們接受假設，既接受假設，既z數(shù)據(jù)是服從正態(tài)分布的數(shù)據(jù)是服從正態(tài)分布的信號特征低頻的原始信號，加高頻的噪聲低頻的原始信號，加高頻的噪聲如何去掉噪聲？如何去掉噪聲？快速傅里葉變換快速傅里葉變換 FFT應用范圍非常廣，在物理學、化學、電子通訊、信號處理、概率論

6、、統(tǒng)計學、密碼學、聲學、光學、海洋學、結構動力學等原理是將時域中的測量信號轉換到頻域，然后再將得到的頻域信號合成為時域中的信號時域信號轉換為頻域信號時，將信號分解成幅值譜，顯示與頻率對應的幅值大小一個信號是由多種頻率的信號合成的將這些信號分離，就能去掉信號中的噪聲快速傅里葉變換快速傅里葉變換 FFTScipy類庫中提供了快速傅里葉變換包fftpackFFT應用案例產(chǎn)生帶噪聲的信號import numpy as npimport matplotlib.pyplot as pltfrom scipy import fftpack as ffttimeStep = 0.02 # 定義采樣點間隔tim

7、eVec = np.arange(0, 20, timeStep) # 生成采樣點# 模擬帶高頻噪聲的信號sigsig = np.sin( np.pi / 5 * timeVec) + 0.1 * np.random.randn(timeVec.size) # 表達式中后面一項為噪聲 plt.plot(timeVec, sig)plt.show()信號特征低頻的原始信號，加高頻的噪聲低頻的原始信號，加高頻的噪聲如何去掉噪聲？如何去掉噪聲？FFT信號變換 sig已知n=sig.sizesig_fft = fft.fft(sig) # 正正變換后的結果保存在變換后的結果保存在 sig_fft中中s

8、ampleFreq = fft.fftfreq(n, d=timeStep) # 獲得每個采樣點頻率幅值獲得每個采樣點頻率幅值pidxs = np.where(sampleFreq 0) # 結果結果是對稱的是對稱的，僅，僅需要使用譜的正需要使用譜的正值部分來找出信號頻率值部分來找出信號頻率freqs = sampleFreqpidxs # 取得所有正值部分取得所有正值部分power = np.abs(sig_fft)pidxs # 找到對應的頻率振幅值找到對應的頻率振幅值freq = freqspower.argmax() # 假設信噪比足夠強，則最大幅值對應的假設信噪比足夠強，則最大幅值對

9、應的頻率，就是信號的頻率頻率，就是信號的頻率sig_fftnp.abs(sampleFreq) freq = 0 # 舍棄所有非信號頻率舍棄所有非信號頻率main_sig = fft.ifft(sig_fft) # 用傅立葉反變換，重構去除噪聲的信號用傅立葉反變換，重構去除噪聲的信號plt.plot(timeVec, main_sig, linewidth=3)尋優(yōu) f(x)=x2-4*x+8 在在x=2的位置，函數(shù)有最小值的位置，函數(shù)有最小值4尋優(yōu) scipy.optimize包提供了求極值的函數(shù)包提供了求極值的函數(shù)fminfrom scipy.optimize import fminimp

10、ort 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多維尋優(yōu)z=x2+y2+8from scipy.optimize import fmindef myfunc(p): # 注意定義 x,y=p return x*2+y*2+8p=(1,1)print (fmin(myfunc, p )多維尋優(yōu)Optimizat

11、ion terminated successfully. Current function value: 8.000000 Iterations: 38 Function evaluations: 69 -2.10235293e-05 2.54845649e-05全局尋優(yōu)全局尋優(yōu)y=x2 + 20 sin(x)全局尋優(yōu)-0開始from scipy import optimizeimport matplotlib.pyplot as pltimport numpy as npdef f(x): return x*2 + 20 * np.sin(x)ans=optimize.fmin_bfgs(f

12、, 0)print(ans)全局最優(yōu)求解Optimization terminated successfully. Current function value: -17.757257 Iterations: 5 Function evaluations: 24 Gradient evaluations: 8當起始點設置為當起始點設置為0時，它找到了時，它找到了0附近的最小極值點，該解附近的最小極值點，該解也是全局最優(yōu)解也是全局最優(yōu)解全局尋優(yōu)-5開始from scipy import optimizeimport matplotlib.pyplot as pltimport numpy as

13、npdef f(x): return x*2 + 20 * np.sin(x)ans=optimize.fmin_bfgs(f, 5)print(ans)全局最優(yōu)求解Optimization terminated successfully. Current function value: 0.158258 Iterations: 5 Function evaluations: 24 Gradient evaluations: 8 4.27109534當當起始點設置起始點設置為為5時時，它找到，它找到了了5附近的附近的局部最優(yōu)局部最優(yōu)全局最優(yōu)求解代替方案optimize.fminboundfro

14、m scipy import optimizeimport numpy as npdef f(x): return x*2 + 20 * np.sin(x)ans = optimize.fminbound(f, -100, 100)print(ans)print(f(ans)-1.42755181859-17.7572565315高維網(wǎng)格尋優(yōu)def f(x,y): z=(np.sin(x)+0.05*x*2) + np.sin(y)+0.05*y*2 return z高維網(wǎng)格尋優(yōu)import numpy as npdef f(p): x,y=p ans=(np.sin(x)+0.05*x*2)

15、 + np.sin(y)+0.05*y*2 return ansimport scipy.optimize as optrranges = (slice(-10, 10, 0.1), slice(-10, 10, 0.1)res=opt.brute(f,rranges)print(res)print(f(res)x和和y都是在都是在-10,10區(qū)間內(nèi)，采用步長區(qū)間內(nèi)，采用步長0.1進行網(wǎng)格搜索求進行網(wǎng)格搜索求最優(yōu)解最優(yōu)解-1.42755002 -1.42749423-1.77572565134求解一元高次方程的根求解一元高次方程的根from scipy import optimizeimpor

16、t 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.68753900e-14# 不同的初值，會怎么樣？不同的初值，會怎么樣？非線性方程組求解scipy. optimize的的fsolve函數(shù)也可以方便用于求解非函數(shù)也可以方便用于求解非線性線性方程組方程組求解原則：求解原則：將方程組的右邊全部規(guī)劃為將方程組的右邊全部規(guī)劃為0將方程組將方程組定義為如下格式的定義為如下格式的python函數(shù)函數(shù)def f(x): x1,x2,

17、 xn=x return f1(x1, x2, xn), f2(x1,x2, xn),.非線性方程組求解 例子2*x1+3=04*x02 + sin(x1*x2)=0 x1*x2/2 3=0非線性方程組求解 例子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 - 3ans = fsolve(f, 1.0,1.0,1.0)print (ans)print (f(ans)-0.26429884 -1.5 -4

18、.0.0, 1.1482453876610066e-10, 6.4002136923591024e-12常微分方程組求解洛侖茲函數(shù)可以用下面微分方程洛侖茲函數(shù)可以用下面微分方程描述描述方程定義了三維空間中各個坐標點上的速度矢量。從某個坐標開始沿方程定義了三維空間中各個坐標點上的速度矢量。從某個坐標開始沿著速度矢量進行積分，就可以計算出無質量點在此空間中的運動著速度矢量進行積分，就可以計算出無質量點在此空間中的運動軌跡軌跡,為三個常數(shù)，為三個常數(shù)，x,y,z為點的坐標為點的坐標常微分方程組求解Scipy中提供了函數(shù)中提供了函數(shù)odeint，負責微分方程組的，負責微分方程組的求求解解是是一個參數(shù)非

19、常復雜的函數(shù)，但通常我們關心的一個參數(shù)非常復雜的函數(shù)，但通常我們關心的只是該函數(shù)的前只是該函數(shù)的前3項，因此，函數(shù)的調用格式可項，因此，函數(shù)的調用格式可以縮寫為：以縮寫為：odeint(func, y0, t)func是有關微分方程組的函數(shù)是有關微分方程組的函數(shù)，y0是一個元組，記錄每個變量的初值是一個元組，記錄每個變量的初值，t則是一個時間序列則是一個時間序列。使用使用時請注意，時請注意，oedint函數(shù)，要求微分方程必須函數(shù)，要求微分方程必須化為標準形式，即化為標準形式，即dy/dt=f(y,t)。常微分方程組求解 lorenz函數(shù)def lorenz(w, t): # 給出位置矢量給出位

20、置矢量w，和三個參數(shù)，和三個參數(shù)r,p, b計算出計算出 r=10.0 p=28.0 b=3.0 # w是是 x,y,z的的值值 x, y, z = w # 直接與直接與lorenz的計算公式對應的計算公式對應 return np.array(r*(y-x), x*(p-z)-y, x*y-b*z) # 三個微分方程，每個作為一項，寫進一個列表中三個微分方程，每個作為一項，寫進一個列表中常微分方程組求解 loren函數(shù)import numpy as np t = np.arange(0, 30, 0.01) # 創(chuàng)建時間點創(chuàng)建時間點 # 調用調用odeint對對lorenz進行求解進行求解,

21、用兩個不同的初始值用兩個不同的初始值 track1 = odeint(lorenz, (0.0, 1.00, 0.0), t) track2 = odeint(lorenz, (0.0, 1.01, 0.0), t) # 繪圖繪圖from mpl_toolkits.mplot3d import Axes3Dimport matplotlib.pyplot as plt fig = plt.figure()ax = Axes3D(fig)ax.plot(track1:,0, track1:,1, track1:,2)ax.plot(track2:,0, track2:,1, track2:,2)

22、plt.show()print(track1)曲線擬合 curve-fit給定的給定的y=ax+b函數(shù)上的一系列采樣點，并在這些函數(shù)上的一系列采樣點，并在這些采樣點上增加一些噪聲，然后利用采樣點上增加一些噪聲，然后利用scipy optimize包中提供的包中提供的curve_fit方法，求解系數(shù)方法，求解系數(shù)a和和b曲線擬合 curve-fitfrom scipy import optimizeimport matplotlib.pyplot as pltimport numpy as npdef f(x,a,b): return a*x + b曲線擬合 curve-fitx = np.li

23、nspace(-10, 10, 30)y = f(x,2,1)ynew= y+ 3*np.random.normal(size=x.size) # 產(chǎn)生帶噪聲產(chǎn)生帶噪聲的數(shù)據(jù)點的數(shù)據(jù)點popt, pcov = optimize.curve_fit(f,x,ynew)print(popt)print (pcov)plt.plot(x,y,color=r,label=orig)plt.plot(x,popt0*x+popt1,color=b,label=fitting)plt.legend(loc=upper left)plt.scatter(x,ynew)plt.show()曲線擬合 curve

24、-fitpopt= 1.99022068 0.34002638pcov= 6.14619911e-03 -1.53673628e-11 -1.53673628e-11 2.19002498e-01popt列表包含每個參數(shù)的擬合值，此例求得的列表包含每個參數(shù)的擬合值，此例求得的a為為1.99022068，b為為0.34002638。pcov列表的對角元素是每個列表的對角元素是每個參數(shù)的方差。通過方差，可以評判擬合的質量，方差越小參數(shù)的方差。通過方差，可以評判擬合的質量，方差越小，擬合越可靠，擬合越可靠插值根據(jù)現(xiàn)有的試驗點值，去預測中間的點根據(jù)現(xiàn)有的試驗點值，去預測中間的點值值采用線性、兩次樣條、

25、三次樣條插值采用線性、兩次樣條、三次樣條插值插值-案例首先在首先在0,10區(qū)間里等間距產(chǎn)生了區(qū)間里等間距產(chǎn)生了20個采樣點，個采樣點，并計算其并計算其sin值，模擬試驗采集得到的值，模擬試驗采集得到的20個個點點采用線性、兩次樣條、三次樣條插值函數(shù)插值擬采用線性、兩次樣條、三次樣條插值函數(shù)插值擬合原函數(shù)合原函數(shù)sin(x)插值-案例import numpy as npfrom erpolate import interp1dimport matplotlib.pyplot as pltx=np.linspace(0,10*np.pi,20) #產(chǎn)生產(chǎn)生20個點個點y=np.s

26、in(x)# x,y 現(xiàn)在是假設的采樣點現(xiàn)在是假設的采樣點插值-案例fl=interp1d(x,y,kind=linear) # 線性插值函數(shù)線性插值函數(shù)fq=interp1d(x,y,kind=quadratic) # 二次樣條插值二次樣條插值fc=interp1d(x,y,kind=cubic) # 三次樣條插值三次樣條插值xint=np.linspace(x.min(),x.max(),1000) # 產(chǎn)生插值產(chǎn)生插值點點yintl=fl(xint) # 由線性插值得到的函數(shù)值由線性插值得到的函數(shù)值yintq=fq(xint) # 由二次樣條插值函數(shù)計算得到的函數(shù)值由二次樣條插值函數(shù)計算

27、得到的函數(shù)值yintc=fc(xint) # 由三次樣條插值函數(shù)計算得到的函數(shù)由三次樣條插值函數(shù)計算得到的函數(shù)值值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還可以對含有噪聲還可以對含有噪聲的的數(shù)據(jù)，數(shù)據(jù)，進行進行樣條插值并自動過濾樣條插值并自動過濾部分噪聲，

28、使用部分噪聲，使用UnivariateSpline函數(shù)，并啟動其函數(shù)，并啟動其s參數(shù)即可參數(shù)即可實現(xiàn)該實現(xiàn)該功能功能from erpolate import UnivariateSpline插值-模擬帶噪聲的問題import numpy as npfrom erpolate import UnivariateSplineimport matplotlib.pyplot as pltsample=50 x=np.linspace(1,20*np.pi,sample)y=np.sin(x) + np.log(x) + np.random.randn(sample

29、)/10#在函數(shù)取值上增加了正態(tài)分布的隨機噪聲在函數(shù)取值上增加了正態(tài)分布的隨機噪聲插值-模擬帶噪聲的問題f=UnivariateSpline(x,y,s=1) # s=1 啟用啟用s參數(shù)，生成行條函參數(shù)，生成行條函數(shù)數(shù)xint=np.linspace(x.min(),x.max(),1000)yint=f(xint)plt.plot(xint,yint,color=r, linestyle=-,label=interpolation)plt.plot(x,y,color=b ,label=orig)plt.legend(loc=upper left)plt.show()多項式擬合處理impor

30、t numpy as npimport matplotlib.pyplot as pltfrom scipy import signal # 引入信號處理包引入信號處理包from pylab import *mpl.rcParamsfont.sans-serif = SimHeiX=np.mafromtxt(rE:teach教改項目教材教改項目教材墨翠樣品拉曼光譜墨翠樣品拉曼光譜墨翠墨綠四季豆墨翠墨綠四季豆.txt)X=X.datax=X:,0 # 文件的第一列為拉曼測量的波數(shù)文件的第一列為拉曼測量的波數(shù)y=X:,1 # 第二列為拉曼響應信號第二列為拉曼響應信號plt.ylabel(u拉曼響應拉曼響應) plt.xlabel(u