SciPy函数库在NumPy库的基础上增加了众多的数学、科学以及工程计算中常用的库函数。例如线性代数、常微分方程数值求解、信号处理、图像处理、稀疏矩阵等。
最小二乘拟合
假设有一组实验数据(x[i],y[i])，我们知道他们之间的函数关系：y=f(x),通过这些已知信息，需要确定函数中的一些参数项。例如：如果f是一个线形函数f(x)=k*x+b,那么参数k和b就是我们需要确定的值。如果将这些参数用p表示的话，那么我们就要找到一组p值使得如下公式的S函数最小：
其中y为原函数，f为拟合函数。拟合函数可通过给定值找出表达式，不要求多精确，之要求能反映数据变化的趋势。
这种算法被称为最小二乘拟合。
scipy中的子函数库optimize已经提供了实现最小二乘拟合算法的函数leastsq。
#-*- coding: utf-8 -*-import numpy as npfrom scipy.optimize import leastsqimport matplotlib.pylab as pldef func(x,p):"""数据拟合所用的函数：A*sin(2*pi*k*x+theta)"""A,k,theta=preturn A*np.sin(2*np.pi*k*x+theta)def residuals(p,y,x):"""实验数据x,y和拟合函数之间的差，p为拟合需要找到的系数"""return y-func(x,p)#在0～-2*np.pi均匀返回100个数字x=np.linspace(0,-2*np.pi,100)A,k,theta=10,0.34,np.pi/6 #真实数据的函数参数y0=func(x,[A,k,theta])y1=y0+2*np.random.randn(len(x)) #加入噪声之后的实验数据p0=[7,0.2,0] #第一次猜测的函数拟合参数#调用leastsq进行数据拟合#residuals为计算误差的函数#p0为拟合参数的初始值#args为需要拟合的实验数据plsq=leastsq(residuals,p0,args=(y1,x))print("Real_parameter：",[A,k,theta])print("Fitting_parameters" ,plsq[0])pl.plot(x,y0,label="RealData")pl.plot(x,y1,label="NoiseData")pl.plot(x,func(x,plsq[0]),label="FittingData")pl.legend()pl.show()
这里我们要拟合的函数是一个正弦波函数，它有参个参数A，k，theta分别对应振幅，频率，相角。假设实验数据是一组包含噪声的数据x,y1，其中y1是在真实数据的基础上加入噪声得到的。
通过leastsq函数对带噪声的实验数据x,y1进行数据拟合，可以找到x和真实数据y0之间的正弦关系的参数：A,k,theta。
Real_parameter： [10, 0.34, 0.5235987755982988]Fitting_parameters [10.72946651 0.34093584 0.58900944]
函数最小值
optimize库提供了几个求函数最小值的算法：fmin，fmin_powell，fmin_cg,fim_bfgs。下面的程序通过求解卷积的逆运算岩石fmin的功能。
对于一个离散的线性时不变系统h，如果他的输入是x，那么其输出也可以用x和h卷积表示，即y=x*h。
现在的问题是如果已知系统的输入x和输出y，如何计算系统的传递函数h;或者如果已知系统的传递函数h和系统的输出y，如何计算系统的输入x。这种运算被称为反卷积运算。
#coding:utf-8import scipy.optimize as optimport numpy as npdef test_fmin_convolve(fminfunc,x,h,y,yn,x0):"""x(*)h=y卷积yn为在y的基础上添加一些干扰噪声的结果x0为求解x的初始值"""def convolve_func(h):"""计算yn-x(*)h的powerfmin将通过计算使得此power最小"""return np.sum((yn-np.convolve(x,h))**2)#调用fmin函数，以x0为初始值h0=fminfunc(convolve_func,x0)print(fminfunc.__name__)print("----------------------")#输出x(*)h0和y之间的相对误差print("error of y:",float(np.sum(np.convolve(x,h0)-y)**2)/np.sum(y**2))#输出h0和h之间的相对误差print("error of h:",float(np.sum((h0-h)**2)/np.sum(h**2)))print('/n')def test_n(m,n,nscale):"""随机产生x,h,y,yn,x0等数列，调用个和宗fmin函数求解bm为x的长度，n为h的长度，nscale为干扰强度"""x=np.random.rand(m)h=np.random.rand(n)y=np.convolve(x,h)yn=y+np.random.rand(len(y))*nscalex0=np.random.rand(n)test_fmin_convolve(opt.fmin,x,h,y,yn,x0)test_fmin_convolve(opt.fmin_powell,x,h,y,yn,x0)test_fmin_convolve(opt.fmin_cg,x,h,y,yn,x0)test_fmin_convolve(opt.fmin_bfgs,x,h,y,yn,x0)if __name__=="__main__":test_n(200,20,0.1)
fmin----------------------error of y: 0.03524512840813301error of h: 0.10872448809534953/nOptimization terminated successfully.Current function value: 0.214475Iterations: 32Function evaluations: 5999fmin_powell----------------------error of y: 0.04085253762222568error of h: 0.00025751019243286225/nOptimization terminated successfully.Current function value: 0.214441Iterations: 19Function evaluations: 858Gradient evaluations: 39fmin_cg----------------------error of y: 0.04084548316927832error of h: 0.00025921974956644503/nOptimization terminated successfully.Current function value: 0.214441Iterations: 32Function evaluations: 1012Gradient evaluations: 46fmin_bfgs----------------------error of y: 0.0408454823459306error of h: 0.00025922110838010043/n
非线性方程组求解
optimize库中的fsolve函数可以用来对非线性方程组进行进行求解。
func(x)是计算方程组误差的函数，它的参数x是一个矢量，表示方程组的各个未知数的一组可能解，func返回将x带入方程组之后得到的误差；x0为未知数矢量的初始值。
实例：
求解
5*x1+3=04x0x0-2sin(x1x2)=0x1*x2-1.5=0
#coding:utf-8from scipy.optimize import fsolvefrom math import sin,cosdef f(x):x0=float(x[0])x1=float(x[1])x2=float(x[2])return [5*x1+3,4*x0*x0-2*sin(x1*x2),x1*x2-1.5]result=fsolve(f,[1,1,1]) #输入一组可能解print(result)print(f(result)) #返回误差
[-0.70622057 -0.6 -2.5 ][0.0, -9.126033262418787e-14, 5.32907051851e-15]
在对方程组进行求解时，fsolve会自动计算方程组的雅可比矩阵，如果方程组中的未知数很多，而与每个方程有关的未知数较少时，传递一个计算雅可比矩阵的函数将能大幅度提高运算速度。
什么是雅可比矩阵？
雅可比矩阵是一阶偏导数以一定方式排列的矩阵，它给出可微分方程与给定点的最优线性逼近，因此类似于多元函数的导数。
使用雅可比矩阵的fsolve，计算雅可比矩阵的函数j通过fprime参数传递给fsolve，函数j和函数f一样，有一个未知数的解矢量参数x，函数j计算非线性方程组在矢量x点上的雅可比矩阵。
#coding:utf-8from scipy.optimize import fsolvefrom math import sin,cosdef f(x):x0=float(x[0])x1=float(x[1])x2=float(x[2])return [5*x1+3,4*x0*x0-2*sin(x1*x2),x1*x2-1.5]def j(x):x0 = float(x[0])x1 = float(x[1])x2 = float(x[2])#求每个方程的偏导数，生成雅可比矩阵return [[0,5,0],[8*x0, -2*x2*cos(x1*x2), -2*x1*cos(x1*x2)],[0,x2,x1]]#传入生成的雅可比矩阵result=fsolve(f,[1,1,1],fprime=j) #输入一组可能解print(result)print(f(result)) #返回误差
运行可以感觉到效率明显提高。
B-Spline样条曲线
interpolate库提供了许多对数据进行插值运算的函数。下面使用直线和B-Spline对正弦波上的点进行插值。
#coding:utf-8import numpy as npimport matplotlib.pylab as plimport matplotlib.pyplot as pltfrom scipy import interpolateplt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签x=np.linspace(0,2*np.pi+np.pi/4,10)y=np.sin(x)x_new=np.linspace(0,2*np.pi+np.pi/4,100)#得到一个新的线性插值函数f_linear=interpolate.interp1d(x,y)#计算出B-Spline曲线的参数tck=interpolate.splrep(x,y)#将参数传递给splev函数计算出各个取样点的插值结果y_bspline=interpolate.splev(x_new,tck)pl.plot(x,y,"o",label="orginallyData")pl.plot(x_new,f_linear(x_new),label="LineInsert")pl.plot(x_new,y_bspline,label="B-splineInsert")pl.legend()pl.show()
数值积分
数值积分是对定积分的数值求解，例如可以利用数值积分计算某个形状的面积。
实例：
计算半径为1的半圆的面积。
单位半圆曲线可以用下面函数表示。
#coding:utf-8import numpy as npdef half_circle(x):return (1-x**2)**0.5#使用经典的分小矩形计算面积总和，计算出单位面积。N=1000x=np.linspace(-1,1,N)y=half_circle(x)#使用分割法计算dx=2.0/NS=dx*np.sum(y[:-1]+y[1:])#进行数值积分S1=np.trapz(y,x) *2#此函数计算的是以x,y为顶点坐标的折线与X所夹的面积。同样的分割点数，结果更精确#使用quad函数积分from scipy import integratepi_half,err=integrate.quad(half_circle,-1,1)S2=pi_half*2print("Split_run: %s" %S)print("Trapz: %s" %S1)print("Quad: %s" %S2)
Split_run: 3.138345831727239Trapz: 3.141487319046285Quad: 3.1415926535897967
可以看出quad函数是非常精确的。
多重定积分的求值可以通过多次调用quad函数实现，为了方便，integrate库提供了dblquad函数进行二重积分，tplquad函数进行三重积分。
实例：
计算单位半球体积。
符合x**2+y**2+z**2=1
from scipy import integrateimport numpy as npdef half_circle(x):return (1-x**2)**0.5def half_sphere(x,y):return (1-x**2-y**2)**0.5# X-Y轴平面与此球体交线为一个单位圆，因此积分区间为此单位圆，可以考虑# 为X轴坐标从-1到1进行积分，而Y轴从-half_circle(x)到half_circle(x)#进行积分V=integrate.dblquad(half_sphere,-1,1,lambda x:-half_circle(x),lambda x:half_circle(x))print(V)#通过球体体积公式计算print(np.pi*4/3/2)
(2.0943951023931984, 1.0002354500215915e-09)2.0943951023931953
dblquad函数的调用方式：
dblquad(func2d,a,b,gfun,hfun)
对于func2d(x,y)函数进行二重积分，其中a,b为变量x的积分期间。
解常微分方程组
scipy.integrate库提供了数值积分和常微分方程组求解算法odeint。
实例：
计算洛伦兹吸引子的轨迹，洛伦兹吸引子有三个微分方程定义：
dx/dt=σ(y-x)dy/dt=x(ρ-z)-ydz/dt=xy-βz
这三个方程定义了三维空间中各个坐标点上的速度矢量。从某个坐标开始沿着速度矢量进行积分，就可以计算出无质量点在此空间中的运动轨迹。其中 σ, ρ, β 为三个常数,不同的参数可以计算出不同的
运动轨迹: x(t), y(t), z(t)。 当参数为某些值时,轨迹出现馄饨现象:即微小的初值差别也会显著地影响运动轨迹。下面是洛仑兹吸引子的轨迹计算和绘制程序:
#coding:utf-8from scipy.integrate import odeintimport numpy as npdef lorenz(w,t,p,r,b):# 给出位置矢量w和三个参数p,r,b计算出#dx/dtt,dy/dt/dz/dt的值x,y,z=w#直接用lorenz的计算公式对应return np.array([p*(y-x),x*(r-z)-y,x*y-b*z])t=np.arange(0,30,0.01) #创建时间点#调用ode对lorenz进行求解，用两个不同的初始值track1=odeint(lorenz,(0.0,1.00,0.0),t,args=(10.0,28.0,3.0))track2=odeint(lorenz,(0.0,1.01,0.0),t,args=(10.0,28.0,3.0))#绘图from mpl_toolkits.mplot3d import Axes3Dimport matplotlib.pyplot as pltfig=plt.figure()ax=Axes3D(fig)ax.plot(track1[:,0],track1[:,1],track1[:,2])ax.plot(track2[:,0],track1[:,1],track1[:,2])plt.show()
在程序中先定义一个lozren函数，它的任务是计算出某个位置的各个方向的微分值，这个计算直接根据洛伦兹吸引子公式得出。但后调用odeint，对微分方程求解：
参数祥解：
lorenz：计算某个位移上的各个方向的速度；(0.0,1.0,0.0)，位移初始值。计算常微分方程所需的各个变量的初始值。t:表示时间的数组，odeint对于此数组中的每个时间点进行求解，的出所有时间点的位置。args：这些参数直接传递给lorenz函数，因此他们都是常量。
滤波器设计
scipy.signal库提供了西多信号处理反面的函数。如何利用signal库设计滤波器，查看滤波器的频率响应，以及如何使用滤波器对信号进行滤波。
设计一个带通滤波器：
import scipy.signal as signalb,a=signal.iirdesign([0.2,0.5],[0.1,0.6],2,40)#参数解析：[0.2,0.5]通带范围# [0.1,0.6]阻带小于0.1*f0大于0.6*f0# 通带的最大增益衰减为2dB，阻带的最小增益衰减为40dB，# 其中f0为1/2的信号取样频率# iirdesign返回的两个数组b和a,他们分别是IIR滤波器的分子和分母部分的系数.# 其中a[0]恒等于1
通过调用freqz计算所得到的滤波器的频率响应:
w,h=signal.freqz(b,a)
frez返回两个数组w和h,其中w是圆频率数组,通过w/pi*f0可以计算出其对应的实际频率.h是w中的对应频率点的响应.其幅值为滤波器的增益,相角为滤波器的相位特性.
计算h的增益特性,并转换为dB度量,由于h中存在幅值几乎为0的值,因此先用clip函数对其进行裁剪,在调用对数函数,避免计算出错。
#绘制出滤波器的增益特性图pl.plot(w/np.pi*4000,power)pl.show()
测量未知系统的频率特性
将频率扫描到系统中，观察系统的输出，从而计算其频率特性。
# 为了调用chirp函数产生频率扫描波形的数据,首先产生一个等差数组代表取样时间,#产生2秒钟取样频率为8khz的取样时间数组.t=np.arange(0,2,1/8000.0)#调用chirp得到2秒钟的频率扫描波形数据sweep=signal.chirp(t,f0=0,t1=2,f1=4000.0)#频率扫描波的开始频率f0为0HZ，结束频率为f1为4kHz，到达4kHz的时间为2秒#使用数组t作为取样时间点。#调用lfilter函数计算sweep波形经过带通滤波器值后的效果out=signal.lfilter(b,a,sweep)
lfilter内部通过如下算是计算IIR滤波器的输出：
数组x代表输入信号，y代表输出信号：
y[n]=b[0]x[n]+b[1]x[n-1]+...+b[P]x[n-P]-a[1]y[n-1]-a[2]y[n-2]-...-a[Q]y[n-Q]
获取输出波形的包络，所以先将输出波形数据转换为能量值：
# 为了调用chirp函数产生频率扫描波形的数据,首先产生一个等差数组代表取样时间,#产生2秒钟取样频率为8khz的取样时间数组.t=np.arange(0,2,1/8000.0)#调用chirp得到2秒钟的频率扫描波形数据sweep=signal.chirp(t,f0=0,t1=2,f1=4000.0)#频率扫描波的开始频率f0为0HZ，结束频率为f1为4kHz，到达4kHz的时间为2秒#使用数组t作为取样时间点。#调用lfilter函数计算sweep波形经过带通滤波器值后的效果out=signal.lfilter(b,a,sweep)out=20*np.log10(np.abs(out))#找到所有能量大于前后两个取样点的下标index=np.where(np.logical_and(out[1:-1]>out[:-2],out[1:-1]>out[2:]))[0]+1#将时间转换为对用频率，绘制所有局部最大点的能量值：pl.plot(t[index]/2.0*4000,out[index])pl.show()
用Weave嵌入C语言
Python作为动态语言其功能虽然强大,但是在数值计算方面有一个最大的缺点:速度不够快。在Python级别的循环和计算的速度只有C语言程序的百分之一。因此才有了NumPy, SciPy这样的函数库,将高度优化的C、Fortran的函数库进行包装,以供Python程序调用。如果这些高度优化的函数库无法实现我们的算法,必须从头开始写循环、计算的话,那么用Python来做显然是不合适的。因此SciPy提供了快速调用C++语言程序的方法-- Weave。
如果觉得《《用Python进行科学计算》——SciPy数值计算库》对你有帮助，请点赞、收藏，并留下你的观点哦！
