《scipy-數(shù)據(jù)處理應(yīng)用.ppt》由會(huì)員分享，可在線閱讀，更多相關(guān)《scipy-數(shù)據(jù)處理應(yīng)用.ppt（52頁(yè)珍藏版）》請(qǐng)?jiān)谘b配圖網(wǎng)上搜索。

科學(xué)計(jì)算包SciPY及應(yīng)用,第11講,Scipy簡(jiǎn)介,解決python科學(xué)計(jì)算而編寫的一組程序包 快速實(shí)現(xiàn)相關(guān)的數(shù)據(jù)處理 如以前的課程中的積分,Scipy提供的數(shù)據(jù)I/O,相比numpy，scipy提供了更傻瓜式的操作方式 二進(jìn)制存儲(chǔ) from scipy import io as fio import numpy as np x=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) data[mat1],Scipy的IO,data[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包,import scipy.stats as stats x=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)證？ import numpy as np import scipy.stats as stats normDist=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)中的噪聲,快速傅里葉變換 FFT,Scipy類庫(kù)中提供了快速傅里葉變換包fftpack,FFT應(yīng)用案例—產(chǎn)生帶噪聲的信號(hào),import numpy as np import matplotlib.pyplot as plt from scipy import fftpack as fft timeStep = 0.02 # 定義采樣點(diǎn)間隔 timeVec = np.arange(0, 20, timeStep) # 生成采樣點(diǎn) # 模擬帶高頻噪聲的信號(hào)sig sig = 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)變換 sig已知,n=sig.size sig_fft = fft.fft(sig) # 正變換后的結(jié)果保存在 sig_fft中 sampleFreq = fft.fftfreq(n, d=timeStep) # 獲得每個(gè)采樣點(diǎn)頻率幅值 pidxs = np.where(sampleFreq 0) # 結(jié)果是對(duì)稱的，僅需要使用譜的正值部分來找出信號(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ù)fmin from scipy.optimize import fmin import numpy as np def f(x): return x**2-4*x+8 print (fmin(f, 0)),Optimization terminated successfully. Current function value: 4.000000 Iterations: 27 Function evaluations: 54,多維尋優(yōu),z=x2+y2+8 from scipy.optimize import fmin def myfunc(p): # 注意定義 x,y=p return x**2+y**2+8 p=(1,1) print (fmin(myfunc, p )),多維尋優(yōu),Optimization terminated successfully. Current function value: 8.000000 Iterations: 38 Function evaluations: 69 [ -2.10235293e-05 2.54845649e-05],全局尋優(yōu),y=x2 + 20 sin(x),全局尋優(yōu)---0開始,from scipy import optimize import matplotlib.pyplot as plt import numpy as np def f(x): return x**2 + 20 * np.sin(x) ans=optimize.fmin_bfgs(f, 0) print(ans),全局最優(yōu)求解,Optimization terminated successfully. Current function value: -17.757257 Iterations: 5 Function evaluations: 24 Gradient evaluations: 8 當(dāng)起始點(diǎn)設(shè)置為0時(shí)，它找到了0附近的最小極值點(diǎn)，該解也是全局最優(yōu)解,全局尋優(yōu)---5開始,from scipy import optimize import matplotlib.pyplot as plt import numpy as np def 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] 當(dāng)起始點(diǎn)設(shè)置為5時(shí)，它找到了5附近的局部最優(yōu),全局最優(yōu)求解—代替方案optimize.fminbound,from scipy import optimize import numpy as np def 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 np def f(p): x,y=p ans=(np.sin(x)+0.05*x**2) + np.sin(y)+0.05*y**2 return ans import scipy.optimize as opt rranges = (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,求解一元高次方程的根,from scipy import optimize import numpy as np def f(x): return x**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ù) def f(x): x1,x2,…, xn=x return [f1(x1, x2,…, xn), f2(x1,x2,…, xn),….],非線性方程組求解 例子,2*x1+3=0 4*x02 + sin(x1*x2)=0 x1*x2/2 – 3=0,非線性方程組求解 例子,from scipy.optimize import fsolve from 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],常微分方程組求解,洛侖茲函數(shù)可以用下面微分方程描述 方程定義了三維空間中各個(gè)坐標(biāo)點(diǎn)上的速度矢量。從某個(gè)坐標(biāo)開始沿著速度矢量進(jìn)行積分，就可以計(jì)算出無質(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)用格式可以縮寫為： 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ù),def lorenz(w, t): # 給出位置矢量w，和三個(gè)參數(shù)r,p, b計(jì)算出 r=10.0 p=28.0 b=3.0 # w是 x,y,z的值 x, y, z = w # 直接與lorenz的計(jì)算公式對(duì)應(yīng) return np.array([r*(y-x), x*(p-z)-y, x*y-b*z]) # 三個(gè)微分方程，每個(gè)作為一項(xiàng)，寫進(jìn)一個(gè)列表中,常微分方程組求解 loren函數(shù),import numpy as np t = 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) # 繪圖 from mpl_toolkits.mplot3d import Axes3D import 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]) plt.show() print(track1),曲線擬合 curve-fit,給定的y=ax+b函數(shù)上的一系列采樣點(diǎn)，并在這些采樣點(diǎn)上增加一些噪聲，然后利用scipy optimize包中提供的curve_fit方法，求解系數(shù)a和b,曲線擬合 curve-fit,from scipy import optimize import matplotlib.pyplot as plt import numpy as np def f(x,a,b): return a*x + b,曲線擬合 curve-fit,x = 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=upper left) plt.scatter(x,ynew) plt.show(),曲線擬合 curve-fit,popt= [ 1.99022068 0.34002638] pcov= [[ 6.14619911e-03 -1.53673628e-11] [ -1.53673628e-11 2.19002498e-01]] popt列表包含每個(gè)參數(shù)的擬合值，此例求得的a為1.99022068，b為0.34002638。pcov列表的對(duì)角元素是每個(gè)參數(shù)的方差。通過方差，可以評(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),插值---案例,import numpy as np from scipy.interpolate import interp1d import matplotlib.pyplot as plt x=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(),插值---案例,插值---模擬帶噪聲的問題,Scipy還可以對(duì)含有噪聲的數(shù)據(jù)，進(jìn)行樣條插值并自動(dòng)過濾部分噪聲，使用UnivariateSpline函數(shù)，并啟動(dòng)其s參數(shù)即可實(shí)現(xiàn)該功能 from scipy.interpolate import UnivariateSpline,插值---模擬帶噪聲的問題,import numpy as np from scipy.interpolate import UnivariateSpline import matplotlib.pyplot as plt sample=50 x=np.linspace(1,20*np.pi,sample) y=np.sin(x) + np.log(x) + np.random.randn(sample)/10 #在函數(shù)取值上增加了正態(tài)分布的隨機(jī)噪聲,插值---模擬帶噪聲的問題,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=upper left) plt.show(),多項(xiàng)式擬合處理,import numpy as np import matplotlib.pyplot as plt from scipy import signal # 引入信號(hào)處理包 from pylab import * mpl.rcParams[font.sans-serif] = [SimHei] X=np.mafromtxt(r“E:\teach\教改項(xiàng)目教材\墨翠樣品拉曼光譜\墨翠墨綠四季豆.txt“) X=X.data x=X[:,0] # 文件的第一列為拉曼測(cè)量的波數(shù) y=X[:,1] # 第二列為拉曼響應(yīng)信號(hào) plt.ylabel(u拉曼響應(yīng)) plt.xlabel(u波數(shù)) plt.plot(x,y,r,label=orig) # 畫帶有基線的信號(hào) plt.plot(x, signal.detrend(y), b,label=detrend) # 畫去掉基線后的信號(hào) plt.legend() plt.show(),多項(xiàng)式擬合處理,模式聚類,Scipy的聚類分析中主要提供了矢量量化分析和系統(tǒng)聚類法,模式聚類,import numpy as np from scipy.cluster import vq import matplotlib.pyplot as plt class1=np.random.randn(30,2)+10 # 產(chǎn)生第一個(gè)正態(tài)分布類，基礎(chǔ)抬高10 class2=np.random.randn(40,2)-10 # 產(chǎn)生第二個(gè)正態(tài)分布類，基礎(chǔ)降低10 class3=np.random.randn(50,2) # 產(chǎn)生第三個(gè)正態(tài)分布類 data=np.vstack([class1,class2,class3]) #將數(shù)據(jù)疊合到一起，形成一個(gè)矩陣,模式聚類,centroid, var=vq.kmeans(data,3) # 用k均值聚類法聚類，指定按3個(gè)類別聚類，獲取類中心和方差 key,distance=vq.vq(data,centroid) # 根據(jù)聚類中心，將不同的樣本分類 vqclass1=data[key==0] # 取出類別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(),模式聚類,