第11讲--scipy-数据处理应用

科学计算包SciPY及应用第11讲Scipy简介解决python科学计算而编写的一组程序包快速实现相关的数据处理如以前的课程中的积分Scipy提供的数据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.]])统计假设与检验stats包stats提供了产生连续性分布的函数，均匀分布(uniform)、正态分布（norm）、贝塔分布（beta）；产生离散分布的函数，伯努利分布（bernoulli）、几何分布（geom）泊松分布poisson使用时，调用分布的rvs方法即可统计假设与检验stats包importscipy.statsasstatsx=stats.uniform.rvs(size=20)#产生20个在[0，1]均匀分布的随机数y=stats.beta.rvs(size=20,a=3,b=4)产生20个服从参数a=3，b=4的贝塔分布随机数z=stats.norm.rvs(size=20,loc=0,scale=1)产生了20个服从[0，1]正态分布的随机数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是实验获得的数据，如何验证它是否是正态分布的？假设检验假设给定的样本服从某种分布，如何验证？statVal,pVal=stats.kstest(z,'norm',(mean,dev))print('pVal=',pVal)计算得到：pVal=0.667359687999结论：我们接受假设，既z数据是服从正态分布的信号特征低频的原始信号，加高频的噪声如何去掉噪声？快速傅里叶变换FFT应用范围非常广，在物理学、化学、电子通讯、信号处理、概率论、统计学、密码学、声学、光学、海洋学、结构动力学等原理是将时域中的测量信号转换到频域，然后再将得到的频域信号合成为时域中的信号时域信号转换为频域信号时，将信号分解成幅值谱，显示与频率对应的幅值大小一个信号是由多种频率的信号合成的将这些信号分离，就能去掉信号中的噪声快速傅里叶变换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信号变换sig已知n=sig.sizesig_fft=fft.fft(sig)#正变换后的结果保存在sig_fft中sampleFreq=fft.fftfreq(n,d=timeStep)#获得每个采样点频率幅值pidxs=np.where(sampleFreq0)#结果是对称的，仅需要使用谱的正值部分来找出信号频率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)寻优f(x)=x2-4*x+8在x=2的位置，函数有最小值4寻优scipy.optimize包提供了求极值的函数fminfromscipy.optimizeimportfminimportnumpyasnpdeff(x):returnx**2-4*x+8print(fmin(f,0))Optimizationterminatedsuccessfully.Currentfunctionvalue:4.000000Iterations:27Functionevaluations:54多维寻优z=x2+y2+8fromscipy.optimizeimportfmindefmyfunc(p):#注意定义x,y=preturnx**2+y**2+8p=(1,1)print(fmin(myfunc,p))多维寻优Optimizationterminatedsuccessfully.Currentfunctionvalue:8.000000Iterations:38Functionevaluations:69[-2.10235293e-052.54845649e-05]全局寻优y=x2+20sin(x)全局寻优---0开始fromscipyimportoptimizeimportmatplotlib.pyplotaspltimportnumpyasnpdeff(x):returnx**2+20*np.sin(x)ans=optimize.fmin_bfgs(f,0)print(ans)全局最优求解Optimizationterminatedsuccessfully.Currentfunctionvalue:-17.757257Iterations:5Functionevaluations:24Gradientevaluations:8当起始点设置为0时，它找到了0附近的最小极值点，该解也是全局最优解全局寻优---5开始fromscipyimportoptimizeimportmatplotlib.pyplotaspltimportnumpyasnpdeff(x):returnx**2+20*np.sin(x)ans=optimize.fmin_bfgs(f,5)print(ans)全局最优求解Optimizationterminatedsuccessfully.Currentfunctionvalue:0.158258Iterations:5Functionevaluations:24Gradientevaluations:8[4.27109534]当起始点设置为5时，它找到了5附近的局部最优全局最优求解—代替方案optimize.fminboundfromscipyimportoptimizeimportnumpyasnpdeff(x):returnx**2+20*np.sin(x)ans=optimize.fminbound(f,-100,100)print(ans)print(f(ans))-1.42755181859-17.7572565315高维网格寻优deff(x,y):z=(np.sin(x)+0.05*x**2)+np.sin(y)+0.05*y**2returnz高维网格寻优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区间内，采用步长0.1进行网格搜索求最优解[-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函数也可以方便用于求解非线性方程组求解原则：将方程组的右边全部规划为0将方程组定义为如下格式的python函数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.][0.0,1.1482453876610066e-10,6.4002136923591024e-12]常微分方程组求解洛仑兹函数可以用下面微分方程描述方程定义了三维空间中各个坐标点上的速度矢量。从某个坐标开始沿着速度矢量进行积分，就可以计算出无质量点在此空间中的运动轨迹σ,ρ,β为三个常数，x,y,z为点的坐标常微分方程组求解Scipy中提供了函数odeint，负责微分方程组的求解是一个参数非常复杂的函数，但通常我们关心的只是该函数的前3项，因此，函数的调用格式可以缩写为：odeint(func,y0,t)func是有关微分方程组的函数，y0是一个元组，记录每个变量的初值，t则是一个时间序列。使用时请注意，oedint函数，要求微分方程必须化为标准形式，即dy/dt=f(y,t)。常微分方程组求解lorenz函数deflorenz(w,t):#给出位置矢量w，和三个参数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函数importnumpyasnpt=np.arange(0,30,0.01)#创建时间点#调用odeint对lorenz进行求解,用两个不同的初始值track1=odeint(lorenz,(0.0,1.00,0.0),t)track2=odeint(loren