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

科學(xué)計(jì)算包SciPY及應(yīng)用第11講Scipy簡(jiǎn)介解決python科學(xué)計(jì)算而編寫(xiě)的一組程序包快速實(shí)現(xiàn)相關(guān)的數(shù)據(jù)處理如以前的課程中的積分Scipy提供的數(shù)據(jù)I/O相比numpy，scipy提供了更傻瓜式的操作方式二進(jìn)制存儲(chǔ)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)計(jì)假設(shè)與檢驗(yàn)stats包stats提供了產(chǎn)生連續(xù)性分布的函數(shù)，均勻分布(uniform)、正態(tài)分布〔norm〕、貝塔分布〔beta〕；產(chǎn)生離散分布的函數(shù)，伯努利分布〔bernoulli〕、幾何分布〔geom〕泊松分布poisson使用時(shí)，調(diào)用分布的rvs方法即可統(tǒng)計(jì)假設(shè)與檢驗(yàn)stats包importscipy.statsasstatsx=stats.uniform.rvs(size=20)#產(chǎn)生20個(gè)在[0，1]均勻分布的隨機(jī)數(shù)y=stats.beta.rvs(size=20,a=3,b=4)產(chǎn)生20個(gè)服從參數(shù)a=3，b=4的貝塔分布隨機(jī)數(shù)z=stats.norm.rvs(size=20,loc=0,scale=1)產(chǎn)生了20個(gè)服從[0，1]正態(tài)分布的隨機(jī)數(shù)x=stats.poisson.rvs(0.6,loc=0,size=20)產(chǎn)生poisson分布假設(shè)檢驗(yàn)假設(shè)給定的樣本服從某種分布，如何驗(yàn)證？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)設(shè)z是實(shí)驗(yàn)獲得的數(shù)據(jù)，如何驗(yàn)證它是否是正態(tài)分布的？假設(shè)檢驗(yàn)假設(shè)給定的樣本服從某種分布，如何驗(yàn)證？statVal,pVal=stats.kstest(z,'norm',(mean,dev))print('pVal=',pVal)計(jì)算得到：pVal=0.667359687999結(jié)論：我們接受假設(shè)，既z數(shù)據(jù)是服從正態(tài)分布的信號(hào)特征低頻的原始信號(hào)，加高頻的噪聲如何去掉噪聲？快速傅里葉變換FFT應(yīng)用范圍非常廣，在物理學(xué)、化學(xué)、電子通訊、信號(hào)處理、概率論、統(tǒng)計(jì)學(xué)、密碼學(xué)、聲學(xué)、光學(xué)、海洋學(xué)、結(jié)構(gòu)動(dòng)力學(xué)等原理是將時(shí)域中的測(cè)量信號(hào)轉(zhuǎn)換到頻域，然后再將得到的頻域信號(hào)合成為時(shí)域中的信號(hào)時(shí)域信號(hào)轉(zhuǎn)換為頻域信號(hào)時(shí)，將信號(hào)分解成幅值譜，顯示與頻率對(duì)應(yīng)的幅值大小一個(gè)信號(hào)是由多種頻率的信號(hào)合成的將這些信號(hào)別離，就能去掉信號(hào)中的噪聲快速傅里葉變換FFTScipy類庫(kù)中提供了快速傅里葉變換包fftpackFFT應(yīng)用案例—產(chǎn)生帶噪聲的信號(hào)importnumpyasnpimportmatplotlib.pyplotaspltfromscipyimportfftpackasffttimeStep=0.02#定義采樣點(diǎn)間隔timeVec=np.arange(0,20,timeStep)#生成采樣點(diǎn)#模擬帶高頻噪聲的信號(hào)sigsig=np.sin(np.pi/5*timeVec)+0.1*np.random.randn(timeVec.size)#表達(dá)式中后面一項(xiàng)為噪聲plt.plot(timeVec,sig)plt.show()信號(hào)特征低頻的原始信號(hào)，加高頻的噪聲如何去掉噪聲？FFT信號(hào)變換sign=sig.sizesig_fft=fft.fft(sig)#正變換后的結(jié)果保存在sig_fft中sampleFreq=fft.fftfreq(n,d=timeStep)#獲得每個(gè)采樣點(diǎn)頻率幅值pidxs=np.where(sampleFreq>0)#結(jié)果是對(duì)稱的，僅需要使用譜的正值局部來(lái)找出信號(hào)頻率freqs=sampleFreq[pidxs]#取得所有正值局部power=np.abs(sig_fft)[pidxs]#找到對(duì)應(yīng)的頻率振幅值freq=freqs[power.argmax()]#假設(shè)信噪比足夠強(qiáng)，那么最大幅值對(duì)應(yīng)的頻率，就是信號(hào)的頻率sig_fft[np.abs(sampleFreq)>freq]=0#舍棄所有非信號(hào)頻率main_sig=fft.ifft(sig_fft)#用傅立葉反變換，重構(gòu)去除噪聲的信號(hào)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開(kāi)始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當(dāng)起始點(diǎn)設(shè)置為0時(shí)，它找到了0附近的最小極值點(diǎn)，該解也是全局最優(yōu)解全局尋優(yōu)---5開(kāi)始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]當(dāng)起始點(diǎn)設(shè)置為5時(shí)，它找到了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高維網(wǎng)格尋優(yōu)deff(x,y):z=(np.sin(x)+0.05*x**2)+np.sin(y)+0.05*y**2returnz高維網(wǎng)格尋優(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ū)間內(nèi)，采用步長(zhǎng)0.1進(jìn)行網(wǎng)格搜索求最優(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]#不同的初值，會(huì)怎么樣？非線性方程組求解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ù)可以用下面微分方程描述方程定義了三維空間中各個(gè)坐標(biāo)點(diǎn)上的速度矢量。從某個(gè)坐標(biāo)開(kāi)始沿著速度矢量進(jìn)行積分，就可以計(jì)算出無(wú)質(zhì)量點(diǎn)在此空間中的運(yùn)動(dòng)軌跡σ,ρ,β為三個(gè)常數(shù)，x,y,z為點(diǎn)的坐標(biāo)常微分方程組求解Scipy中提供了函數(shù)odeint，負(fù)責(zé)微分方程組的求解是一個(gè)參數(shù)非常復(fù)雜的函數(shù)，但通常我們關(guān)心的只是該函數(shù)的前3項(xiàng)，因此，函數(shù)的調(diào)用格式可以縮寫(xiě)為：odeint(func,y0,t)func是有關(guān)微分方程組的函數(shù)，y0是一個(gè)元組，記錄每個(gè)變量的初值，t那么是一個(gè)時(shí)間序列。使用時(shí)請(qǐng)注意，oedint函數(shù)，要求微分方程必須化為標(biāo)準(zhǔn)形式，即dy/dt=f(y,t)。常微分方程組求解lorenz函數(shù)deflorenz(w,t):#給出位置矢量w，和三個(gè)參數(shù)r,p,b計(jì)算出

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

x,y,z=w#直接與lorenz的計(jì)算公式對(duì)應(yīng)

returnnp.array([r*(y-x),x*(p-z)-y,x*y-b*z])#三個(gè)微分方程，每個(gè)作為一項(xiàng)，寫(xiě)進(jìn)一個(gè)列表中常微分方程組求解loren函數(shù)importnumpyasnpt=np.arange(0,30,0.01)#創(chuàng)立時(shí)間點(diǎn)#調(diào)用odeint對(duì)lorenz進(jìn)行求解,用兩個(gè)不同的初始值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ù)上的一系列采樣點(diǎn)，并在這些采樣點(diǎn)上增加一些噪聲，然后利用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)#產(chǎn)生帶噪聲的數(shù)據(jù)點(diǎn)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列表包含每個(gè)參數(shù)的擬合值，此例求得的a為1.99022068，b為0.34002638。pcov列表的對(duì)角元素是每個(gè)參數(shù)的方差。通過(guò)方差，可以評(píng)判擬合的質(zhì)量，方差越小，擬合越可靠插值根據(jù)現(xiàn)有的試驗(yàn)點(diǎn)值，去預(yù)測(cè)中間的點(diǎn)值采用線性、兩次樣條、三次樣條插值插值---案例首先在[0,10π]區(qū)間里等間距產(chǎn)生了20個(gè)采樣點(diǎn)，并計(jì)算其sin值，模擬試驗(yàn)采集得到的20個(gè)點(diǎn)采用線性、兩次樣條、三次樣條插值函數(shù)插值擬合原函數(shù)sin(x)插值---案例importnumpyasnpfromerpolateimportinterp1dimportmatplotlib.pyplotaspltx=np.linspace(0,10*np.pi,20)#產(chǎn)生20個(gè)點(diǎn)y=np.sin(x)#x,y現(xiàn)在是假設(shè)的采樣點(diǎ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)#產(chǎn)生插值點(diǎn)yintl=fl(xint)#由線性插值得到的函數(shù)值yintq=fq(xint)#由二次樣條插值函數(shù)計(jì)算得到的函數(shù)值yintc=fc(xint)#由三次樣條插值函數(shù)計(jì)算得到的函數(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()插值---案例插值---模擬帶噪聲的問(wèn)題Scipy還可以對(duì)含有噪聲的數(shù)據(jù)，進(jìn)行樣條插值并自動(dòng)過(guò)濾局部噪聲，使用UnivariateSpline函數(shù)，并啟動(dòng)其s參數(shù)即可實(shí)現(xiàn)該功能fromerpolateimportUnivariateSpline插值---模擬帶噪聲的問(wèn)題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)分布的隨機(jī)噪聲插值---模擬帶噪聲的問(wèn)題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()多項(xiàng)式擬合處理importnumpyasnpimportmatplotlib.pyplotaspltfromscipyimportsignal#引入信號(hào)處理包frompylabimport*mpl.rcParams['font.sans-serif']=['SimHei']X=np.mafromtxt(r"E:\teach\教改工程教材\墨翠樣品拉曼光譜\墨翠墨綠四季豆.txt")X=X.datax=X[:,0]#文件的第一列為拉曼測(cè)量的波數(shù)y=X[:,1]#第二列為拉曼響應(yīng)信號(hào)plt.ylabel(u