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

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

[1.,1.],

[1.,1.]])

data['mat2']array([[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)計假設與檢驗stats包stats提供了產生連續(xù)性分布的函數(shù)，均勻分布(uniform)、正態(tài)分布〔norm〕、貝塔分布〔beta〕；產生離散分布的函數(shù)，伯努利分布〔bernoulli〕、幾何分布〔geom〕泊松分布poisson使用時，調用分布的rvs方法即可統(tǒng)計假設與檢驗stats包importscipy.statsasstatsx=stats.uniform.rvs(size=20)#產生20個在[0，1]均勻分布的隨機數(shù)y=stats.beta.rvs(size=20,a=3,b=4)產生20個服從參數(shù)a=3，b=4的貝塔分布隨機數(shù)z=stats.norm.rvs(size=20,loc=0,scale=1)產生了20個服從[0，1]正態(tài)分布的隨機數(shù)x=stats.poisson.rvs(0.6,loc=0,size=20)產生poisson分布假設檢驗假設給定的樣本服從某種分布，如何驗證？importnumpyasnpimportscipy.statsasstatsnormDist=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ù)據，如何驗證它是否是正態(tài)分布的？假設檢驗假設給定的樣本服從某種分布，如何驗證？statVal,pVal=stats.kstest(z,'norm',(mean,dev))print('pVal=',pVal)計算得到：pVal=0.667359687999結論：我們接受假設，既z數(shù)據是服從正態(tài)分布的信號特征低頻的原始信號，加高頻的噪聲如何去掉噪聲？快速傅里葉變換FFT應用范圍非常廣，在物理學、化學、電子通訊、信號處理、概率論、統(tǒng)計學、密碼學、聲學、光學、海洋學、結構動力學等原理是將時域中的測量信號轉換到頻域，然后再將得到的頻域信號合成為時域中的信號時域信號轉換為頻域信號時，將信號分解成幅值譜，顯示與頻率對應的幅值大小一個信號是由多種頻率的信號合成的將這些信號別離，就能去掉信號中的噪聲快速傅里葉變換FFTScipy類庫中提供了快速傅里葉變換包fftpackFFT應用案例—產生帶噪聲的信號importnumpyasnpimportmatplotlib.pyplotaspltfromscipyimportfftpackasffttimeStep=0.02#定義采樣點間隔timeVec=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信號變換sign=sig.sizesig_fft=fft.fft(sig)#正變換后的結果保存在sig_fft中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)尋優(yōu)

f(x)=x2-4*x+8在x=2的位置，函數(shù)有最小值4尋優(yōu)

scipy.optimize包提供了求極值的函數(shù)fminfromscipy.optimizeimportfminimportnumpyasnpdeff(x):returnx**2-4*x+8print(fmin(f,0))

Optimizationterminatedsuccessfully.Currentfunctionvalue:4.000000Iterations:27Functionevaluations:54多維尋優(yōu)z=x2+y2+8fromscipy.optimizeimportfmindefmyfunc(p):#注意定義x,y=preturnx**2+y**2+8p=(1,1)print(fmin(myfunc,p))多維尋優(yōu)Optimizationterminatedsuccessfully.Currentfunctionvalue:8.000000Iterations:38Functionevaluations:69[-2.10235293e-052.54845649e-05]全局尋優(yōu)y=x2+20sin(x)全局尋優(yōu)---0開始fromscipyimportoptimizeimportmatplotlib.pyplotaspltimportnumpyasnpdeff(x):returnx**2+20*np.sin(x)ans=optimize.fmin_bfgs(f,0)print(ans)全局最優(yōu)求解Optimizationterminatedsuccessfully.Currentfunctionvalue:-17.757257Iterations:5Functionevaluations:24Gradientevaluations:8當起始點設置為0時，它找到了0附近的最小極值點，該解也是全局最優(yōu)解全局尋優(yōu)---5開始fromscipyimportoptimizeimportmatplotlib.pyplotaspltimportnumpyasnpdeff(x):returnx**2+20*np.sin(x)ans=optimize.fmin_bfgs(f,5)print(ans)全局最優(yōu)求解Optimizationterminatedsuccessfully.Currentfunctionvalue:0.158258Iterations:5Functionevaluations:24Gradientevaluations:8[4.27109534]當起始點設置為5時，它找到了5附近的局部最優(yōu)全局最優(yōu)求解—代替方案optimize.fminboundfromscipyimportoptimizeimportnumpyasnpdeff(x):returnx**2+20*np.sin(x)ans=optimize.fminbound(f,-100,100)print(ans)print(f(ans))-1.42755181859-17.7572565315高維網格尋優(yōu)deff(x,y):z=(np.sin(x)+0.05*x**2)+np.sin(y)+0.05*y**2returnz高維網格尋優(yōu)importnumpyasnpdeff(p):x,y=pans=(np.sin(x)+0.05*x**2)+np.sin(y)+0.05*y**2returnansimportscipy.optimizeasoptrranges=(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ū)間內，采用步長0.1進行網格搜索求最優(yōu)解[-1.42755002-1.42749423]-1.77572565134求解一元高次方程的根fromscipyimportoptimizeimportnumpyasnpdeff(x):returnx**2+20*np.sin(x)ans=optimize.fsolve(f,-4)print(ans)print(f(ans))[-2.75294663][1.68753900e-14]#不同的初值，會怎么樣？非線性方程組求解scipy.optimize的fsolve函數(shù)也可以方便用于求解非線性方程組求解原那么：將方程組的右邊全部規(guī)劃為0將方程組定義為如下格式的python函數(shù)deff(x):x1,x2,…,xn=xreturn[f1(x1,x2,…,xn),f2(x1,x2,…,xn),….]非線性方程組求解

例子2*x1+3=04*x02+sin(x1*x2)=0x1*x2/2–3=0非線性方程組求解

例子fromscipy.optimizeimportfsolvefrommathimport*deff(x):x0,x1,x2=xreturn[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.][常微分方程組求解洛侖茲函數(shù)可以用下面微分方程描述方程定義了三維空間中各個坐標點上的速度矢量。從某個坐標開始沿著速度矢量進行積分，就可以計算出無質量點在此空間中的運動軌跡σ,ρ,β為三個常數(shù)，x,y,z為點的坐標常微分方程組求解Scipy中提供了函數(shù)odeint，負責微分方程組的求解是一個參數(shù)非常復雜的函數(shù)，但通常我們關心的只是該函數(shù)的前3項，因此，函數(shù)的調用格式可以縮寫為：odeint(func,y0,t)func是有關微分方程組的函數(shù)，y0是一個元組，記錄每個變量的初值，t那么是一個時間序列。使用時請注意，oedint函數(shù)，要求微分方程必須化為標準形式，即dy/dt=f(y,t)。常微分方程組求解lorenz函數(shù)deflorenz(w,t):#給出位置矢量w，和三個參數(shù)r,p,b計算出

r=10.0p=28.0b=3.0#w是x,y,z的值

x,y,z=w#直接與lorenz的計算公式對應

returnnp.array([r*(y-x),x*(p-z)-y,x*y-b*z])#三個微分方程，每個作為一項，寫進一個列表中常微分方程組求解loren函數(shù)importnumpyasnpt=np.arange(0,30,0.01)#創(chuàng)立時間點#調用odeint對lorenz進行求解,用兩個不同的初始值track1=odeint(lorenz,(0.0,1.00,0.0),t)track2=odeint(lorenz,(0.0,1.01,0.0),t)#繪圖frommpl_toolkits.mplot3dimportAxes3Dimportmatplotlib.pyplotaspltfig=plt.figure()ax=Axes3D(fig)ax.plot(track1[:,0],track1[:,1],track1[:,2])ax.plot(track2[:,0],track2[:,1],track2[:,2])plt.show()print(track1)曲線擬合curve-fit給定的y=ax+b函數(shù)上的一系列采樣點，并在這些采樣點上增加一些噪聲，然后利用scipyoptimize包中提供的curve_fit方法，求解系數(shù)a和b曲線擬合curve-fitfromscipyimportoptimizeimportmatplotlib.pyplotaspltimportnumpyasnpdeff(x,a,b):returna*x+b曲線擬合curve-fitx=np.linspace(-10,10,30)y=f(x,2,1)ynew=y+3*np.random.normal(size=x.size)#產生帶噪聲的數(shù)據點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='upperleft')plt.scatter(x,ynew)plt.show()曲線擬合curve-fitpopt=[1.990220680.34002638]pcov=[[6.14619911e-03-1.53673628e-11][-1.53673628e-112.19002498e-01]]popt列表包含每個參數(shù)的擬合值，此例求得的a為1.99022068，b為0.34002638。pcov列表的對角元素是每個參數(shù)的方差。通過方差，可以評判擬合的質量，方差越小，擬合越可靠插值根據現(xiàn)有的試驗點值，去預測中間的點值采用線性、兩次樣條、三次樣條插值插值---案例首先在[0,10π]區(qū)間里等間距產生了20個采樣點，并計算其sin值，模擬試驗采集得到的20個點采用線性、兩次樣條、三次樣條插值函數(shù)插值擬合原函數(shù)sin(x)插值---案例importnumpyasnpfromerpolateimportinterp1dimportmatplotlib.pyplotaspltx=np.linspace(0,10*np.pi,20)#產生20個點y=np.sin(x)#x,y現(xiàn)在是假設的采樣點插值---案例fl=interp1d(x,y,kind='linear')#線性插值函數(shù)fq=interp1d(x,y,kind='quadratic')#二次樣條插值fc=interp1d(x,y,kind='cubic')#三次樣條插值xint=np.linspace(x.min(),x.max(),1000)#產生插值點yintl=fl(xint)#由線性插值得到的函數(shù)值yintq=fq(xint)#由二次樣條插值函數(shù)計算得到的函數(shù)值yintc=fc(xint)#由三次樣條插值函數(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ù)據，進行樣條插值并自動過濾局部噪聲，使用UnivariateSpline函數(shù)，并啟動其s參數(shù)即可實現(xiàn)該功能fromerpolateimportUnivariateSpline插值---模擬帶噪聲的問題importnumpyasnpfromerpolateimportUnivariateSplineimportmatplotlib.pyplotaspltsample=50x=np.linspace(1,20*np.pi,sample)y=np.sin(x)+np.log(x)+np.random.randn(sample)/10#在函數(shù)取值上增加了正態(tài)分布的隨機噪聲插值---模擬帶噪聲的問題f=UnivariateSpline(x,y,s=1)#s=1啟用s參數(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='upperleft')plt.show()多項式擬合處理importnumpyasnpimportmatplotlib.pyplotaspltfromscipyimportsignal#引入信號處理包frompylabimport*mpl.rcParams['font.sans-serif']=['SimHei']X=np.mafromtxt(r"E:\teach\教改工程教材\墨翠樣品拉曼光譜\墨翠墨綠四季豆.txt")X=X.datax=X[:,0]#文件的第一列為拉曼測量的波數(shù)y=X[:,1]#第二列為拉曼響應信號plt.ylabel(u