科學(xué)計算包SciPY及應(yīng)用,第11講,Scipy簡介,解決python科學(xué)計算而編寫的一組程序包 快速實現(xiàn)相關(guān)的數(shù)據(jù)處理 如以前的課程中的積分,Scipy提供的數(shù)據(jù)I/O,相比numpy，scipy提供了更傻瓜式的操作方式 二進(jìn)制存儲 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) datamat1,Scipy的IO,datamat1 array( 1., 1., 1., 1., 1., 1.) datamat2 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)計假設(shè)與檢驗 stats包,stats提供了產(chǎn)生連續(xù)性分布的函數(shù)， 均勻分布(uniform)、 正態(tài)分布（norm）、 貝塔分布（beta）； 產(chǎn)生離散分布的函數(shù)， 伯努利分布（bernoulli）、 幾何分布（geom） 泊松分布 poisson 使用時，調(diào)用分布的rvs方法即可,統(tǒng)計假設(shè)與檢驗 stats包,import scipy.stats as stats x=stats.uniform.rvs(size=20) #產(chǎn)生20個在0，1均勻分布的隨機(jī)數(shù) y=stats.beta.rvs(size=20,a=3,b=4) 產(chǎn)生20個服從參數(shù)a=3，b=4的貝塔分布隨機(jī)數(shù) z=stats.norm.rvs(size=20,loc=0,scale=1) 產(chǎn)生了20個服從0，1正態(tài)分布的隨機(jī)數(shù) x=stats.poisson.rvs(0.6,loc=0,size=20) 產(chǎn)生poisson分布,假設(shè)檢驗,假設(shè)給定的樣本服從某種分布，如何驗證？ 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ù)據(jù)，如何驗證它是否是正態(tài)分布的？,假設(shè)檢驗,假設(shè)給定的樣本服從某種分布，如何驗證？ statVal, pVal = stats.kstest(z,norm,(mean,dev) print(pVal=,pVal) 計算得到： pVal= 0.667359687999 結(jié)論：我們接受假設(shè)，既z數(shù)據(jù)是服從正態(tài)分布的,信號特征,低頻的原始信號，加高頻的噪聲 如何去掉噪聲？,快速傅里葉變換 FFT,應(yīng)用范圍非常廣，在物理學(xué)、化學(xué)、電子通訊、信號處理、概率論、統(tǒng)計學(xué)、密碼學(xué)、聲學(xué)、光學(xué)、海洋學(xué)、結(jié)構(gòu)動力學(xué)等 原理是將時域中的測量信號轉(zhuǎn)換到頻域，然后再將得到的頻域信號合成為時域中的信號 時域信號轉(zhuǎn)換為頻域信號時，將信號分解成幅值譜，顯示與頻率對應(yīng)的幅值大小 一個信號是由多種頻率的信號合成的 將這些信號分離，就能去掉信號中的噪聲,快速傅里葉變換 FFT,Scipy類庫中提供了快速傅里葉變換包fftpack,FFT應(yīng)用案例產(chǎn)生帶噪聲的信號,import numpy as np import matplotlib.pyplot as plt from scipy import fftpack as fft timeStep = 0.02 # 定義采樣點間隔 timeVec = np.arange(0, 20, timeStep) # 生成采樣點 # 模擬帶高頻噪聲的信號sig sig = np.sin( np.pi / 5 * timeVec) + 0.1 * np.random.randn(timeVec.size) # 表達(dá)式中后面一項為噪聲 plt.plot(timeVec, sig) plt.show(),信號特征,低頻的原始信號，加高頻的噪聲 如何去掉噪聲？,FFT信號變換 sig已知,n=sig.size sig_fft = fft.fft(sig) # 正變換后的結(jié)果保存在 sig_fft中 sampleFreq = fft.fftfreq(n, d=timeStep) # 獲得每個采樣點頻率幅值 pidxs = np.where(sampleFreq 0) # 結(jié)果是對稱的，僅需要使用譜的正值部分來找出信號頻率 freqs = sampleFreqpidxs # 取得所有正值部分 power = np.abs(sig_fft)pidxs # 找到對應(yīng)的頻率振幅值 freq = freqspower.argmax() # 假設(shè)信噪比足夠強(qiáng)，則最大幅值對應(yīng)的頻率，就是信號的頻率 sig_fftnp.abs(sampleFreq) freq = 0 # 舍棄所有非信號頻率 main_sig = fft.ifft(sig_fft) # 用傅立葉反變換，重構(gòu)去除噪聲的信號 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)起始點設(shè)置為0時，它找到了0附近的最小極值點，該解也是全局最優(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)起始點設(shè)置為5時，它找到了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)，采用步長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 # 不同的初值，會怎么樣？,非線性方程組求解,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ù)可以用下面微分方程描述 方程定義了三維空間中各個坐標(biāo)點上的速度矢量。從某個坐標(biāo)開始沿著速度矢量進(jìn)行積分，就可以計算出無質(zhì)量點在此空間中的運(yùn)動軌跡 ,為三個常數(shù)，x,y,z為點的坐標(biāo),常微分方程組求解,Scipy中提供了函數(shù)odeint，負(fù)責(zé)微分方程組的求解 是一個參數(shù)非常復(fù)雜的函數(shù)，但通常我們關(guān)心的只是該函數(shù)的前3項，因此，函數(shù)的調(diào)用格式可以縮寫為： odeint(func, y0, t) func是有關(guān)微分方程組的函數(shù)， y0是一個元組，記錄每個變量的初值， t則是一個時間序列。 使用時請注意，oedint函數(shù)，要求微分方程必須化為標(biāo)準(zhǔn)形式，即dy/dt=f(y,t)。,常微分方程組求解 lorenz函數(shù),def lorenz(w, t): # 給出位置矢量w，和三個參數(shù)r,p, b計算出 r=10.0 p=28.0 b=3.0 # w是 x,y,z的值 x, y, z = w # 直接與lorenz的計算公式對應(yīng) return np.array(r*(y-x), x*(p-z)-y, x*y-b*z) # 三個微分方程，每個作為一項，寫進(jìn)一個列表中,常微分方程組求解 loren函數(shù),import numpy as np t = np.arange(0, 30, 0.01) # 創(chuàng)建時間點 # 調(diào)用odeint對lorenz進(jìn)行求解, 用兩個不同的初始值 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ù)上的一系列采樣點，并在這些采樣點上增加一些噪聲，然后利用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ù)點 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-fit,popt= 1.99022068 0.34002638 pcov= 6.14619911e-03 -1.53673628e-11 -1.53673628e-11 2.19002498e-01 popt列表包含每個參數(shù)的擬合值，此例求得的a為1.99022068，b為0.34002638。pcov列表的對角元素是每個參數(shù)的方差。通過方差，可以評判擬合的質(zhì)量，方差越小，擬合越可靠,插值,根據(jù)現(xiàn)有的試驗點值，去預(yù)測中間的點值 采用線性、兩次樣條、三次樣條插值,插值-案例,首先在0,10區(qū)間里等間距產(chǎn)生了20個采樣點，并計算其sin值，模擬試驗采集得到的20個點 采用線性、兩次樣條、三次樣條插值函數(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個點 y=np.sin(x) # x,y 現(xiàn)在是假設(shè)的采樣點,插值-案例,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)生插值點 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ù)據(jù)，進(jìn)行樣條插值并自動過濾部分噪聲，使用UnivariateSpline函數(shù)，并啟動其s參數(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(),多項式擬合處理,import numpy as np import matplotlib.pyplot as plt from scipy import signal # 引入信號處理包 from pylab import * mpl.rcParamsfont.sans-serif = SimHei X=np.mafromtxt(r“E:teach教改項目教材墨翠樣品拉曼光譜墨翠墨綠四季豆.txt“) X=X.data x=X:,0 # 文件的第一列為拉曼測量的波數(shù) y=X:,1 # 第二列為拉曼響應(yīng)信號 plt.ylabel(u拉曼響應(yīng)) plt.xlabel(u波數(shù)) plt.plot(x,y,r,label=orig) # 畫帶有基線的信號 plt.plot(x, signal.detrend(y), b,label=detrend) # 畫去掉基線后的信號 plt.legend() plt.show(),多項式擬合處理,模式聚類,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)生第一個正態(tài)分布類，基礎(chǔ)抬高10 class2=np.random.randn(40,2)-10 # 產(chǎn)生第二個正態(tài)分布類，基礎(chǔ)降低10 class3=np.random.randn(50,2) # 產(chǎn)生第三個正態(tài)分布類 data=np.vstack(class1,class2,class3) #將數(shù)據(jù)疊合到一起，形成一個矩陣,模式聚類,centroid, var=vq.kmeans(data,3) # 用k均值聚類法聚類，指定按3個類別聚類，獲取類中心和方差 key,distance=vq.vq(data,centroid) # 根據(jù)聚類中心，將不同的樣本分類 vqclass1=datakey=0 # 取出類別0 vqclass2=datakey=1 vqclass3=datakey=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(),模式聚類,