5.Python科学计算与数据处理

1SciPy—数值计算库目录常数和特殊函数优化—optimize最小二乘拟合函数最小值非线性方程组求解插值—interpolateB样条曲线插值外推和Spline拟合二维插值2目录数值积分—integrate球的体积解常微分方程组统计—stats连续和离散概率分布二项、泊松、伽玛分布34SciPy在NumPy的基础上增加了众多的数学、科学以及工程计算中常用的模块，例如线性代数、常微分方程数值求解、信号处理、图像处理、稀疏矩阵，等等。在本节将通过实例介绍SciPy中常用的一些模块。在实例程序中会使用matplotlib和Mayavi绘制二维和三维图表，在后续章节中这些绘图库进行详细介绍。常数和特殊函数SciPy的constants模块包含了众多的物理常数:在字典physical_constants中，以物理常量名为键，对应的值是一个含有三个元素的元组，分别为常数值、单位及误差，例如下面的程序可以查看电子质量：5fromscipyimportconstantsasCC.c#真空中的光速299792458.0C.h#普朗克常数6.62606930800080626-34C.physical_constants[electronmass](9.1093825999999998e-31,’kg',1.5999999999999999e-37)常数和特殊函数除了物理常数之外，constants模块中还包括许多单位信息，例如:6C.mile#1英里等于多少米1609.3439999999998C.inch#1英寸等于多少米0.025399999999999999C.gram#1克等于多少千克0.001C.pound#1磅等于多少千克0.45359236999999997常数和特殊函数SciPy的special模块是一个非常完整的函数库，其中包含了基本数学函数、特殊数学函数以及NumPy中出现的所有函数。由于函数数量众多，本节仅对其进行简要介绍。至于具体包含的函数列表，请参考SciPy的帮助文档。伽玛函数是概率统计学中经常出现的一个特殊函数，它的计算公式如下：710()ztztedt常数和特殊函数显然，通过此公式计算函数的值是比较麻烦的，可以用special模块中的gamma()进行计算：8importscipy.specialasSS.gamma(4)6.0S.gamma(0.5)1.772458509055159S.gamma(1+1j)#gamma函数支持复数(0.49801566811835629-0.15494982836181106j)S.gamma(1000)inf常数和特殊函数函数是阶乘函数在实数和复数范围上的扩展，它的增长速度非常快，因为1000的阶乘已经超过了双精度浮点数的表示范围，因此结果是无穷大。为了计算更大的范围，可以使用S.gammaln()：S.gammaln(x)计算ln(||)的值，它使用特殊的算法，直接计算函数的对数值，因此可以表示更大的范围。9S.gammaln(1000)5905.2204232091817()x常数和特殊函数special模块中的某些函数并不是数学意义上的特殊函数，例如log1p(x)计算log(1+x)的值。这是由于浮点数的精度有限，无法很精确地表示十分接近1的实数。例如无法用浮点数表示”1+1e-20”的值，因此“log(1+1e-20)”的值为0,而当使用log1p()时，则可以很精确地计算：实际上当x非常小时，log1p(x)约等于x。101+1e-201.0log(1+1e-20)0.0S.log1p(1e-20)9.9999999999999995e-21优化—optimize11最小二乘拟合假设有一组实验数据(xi，yi),事先知道它们之间应该满足某函数关系yi=f(xi)，通过这些已知信息，需要确定函数f的一些参数。例如，如果函数f是线性函数f(x)=kx+b,那么参数k和b就是需要确定的值。优化—optimize如果用p表示函数中需要确定的参数，那么目标就是找到一组p,使得下面的函数S的值最小：这种算法被称为最小二乘拟合(Least-squarefitting)。在optimize模块中，可以使用leastsq()对数据进行最小二乘拟合计算。leastsq()的用法很简单，只需要将计算误差的函数和待确定参数的初始值传递给它即可。下面是用leastsq()对线性函数进行拟合的程序。1221()[(,)]MiiiSpyfxp优化—optimize用最小二乘法拟合直线，并显示误差曲面(scipy_least_square_line.py)13importnumpyasnpfromscipy.optimizeimportleastsqX=np.array([8.19,2.72,6.39,8.71,4.7,2.66,3.78])Y=np.array([7.01,2.78,6.47,6.71,4.1,4.23,4.05])defresiduals(p):计算以p为参数的直线和原始数据之间的误差k,b=preturnY-(k*X+b)#leastsq使得residuals()的输出数组的平方和最小，参数的初始值为[1,0]r=leastsq(residuals,[1,0])k,b=r[0]printk=,k,b=,b优化—optimizeleastsq()函数传入误差计算函数和初始值[1,0]，该初始值将作为误差计算函数的第一个参数传入；计算的结果r是一个包含两个元素的元组，第一个元素是一个数组，表示拟合后的参数k、b；第二个元素如果等于1、2、3、4中的其中一个整数，则拟合成功，否则将会返回mesg.程序的输出为：residuals()的参数p是拟合直线的参数，函数返回的是原始数据和拟合直线之间的误差。14k=0.61349S346194b=1.79489255555优化—optimize再看一个对正弦波数据进行拟合的例子:使用最小二乘法对带噪声的正玄波数据进行拟合(scipy_least_square_sin.py).15使用leastsq()对带噪声的正弦波数据进行拟合。拟合所得到的参数虽然和实际的参数有可能完全不同，但是由于正弦函数具有周期性，实际上拟合的结果和实际的函数是一致的。importnumpyasnpfromscipy.optimizeimportleastsqdeffunc(x,p):数据拟合所用的函数:A*sin(2*pi*k*x+theta)A,k,theta=preturnA*np.sin(2*np.pi*k*x+theta)优化—optimize16defresiduals(p,y,x):实验数据x,y和拟合函数之间的差，p为拟合需要找到的系数returny-func(x,p)x=np.linspace(-2*np.pi,0,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]#第一次猜测的函数拟合参数优化—optimize17#调用leastsq进行数据拟合,residuals为计算误差的函数#p0为拟合参数的初始值,#args为需要拟合的实验数据plsq=leastsq(residuals,p0,args=(y1,x))#除了初始值之外，还调用了args参数，用于指定residuals中使用到的其他参数（直线拟合时直接使用了X,Y的全局变量）,同样也返回一个元组，第一个元素为拟合后的参数数组；printu真实参数:,[A,k,theta]printu“拟合参数”,plsq[0]#实验数据拟合后的参数importpylabasplpl.plot(x,y0,label=u真实数据)pl.plot(x,y1,label=u带噪声的实验数据)pl.plot(x,func(x,plsq[0]),label=u拟合数据)pl.legend()pl.show()优化—optimize程序中，要拟合的func()是一个正弦函数，它的参数p是一个数组，包含决定正弦波的三个参数：A、k、theta,分别对应正弦函数的振幅、频率和相角。一组包含噪声的数据:(x,y1),其中数组y1在标准正弦波数据y0之上添加了随机噪声。用leastsq()对带噪声的实验数据(x,y1)进行数据拟合，它可以找到数组x和真实数据y0之间的正弦关系，即确定A、k、theta等参数。这里将(y1,x)传递给args参数。Leastsq()会将这两个额外的参数传递给residuals()。因此residuals()有三个参数，p是正弦函数的参数，y和x是表示实验数据的数组。18优化—optimize下面是程序的输出：看到拟合参数虽然和真实参数完全不同，但是由于正弦函数具有周期性，实际上拟合参数得到的函数和真实参数对应的函数是一致的。19真实参数:[10,0.34000000000000002,0.52359877559829882]拟合参数[-9.841527750.33829767-2.68899335]优化—optimize函数最小值optimize模块还提供了许多求函数最小值的算法：fmin、fmin_powell、fmin_cg、fmin_bfgs等。下面用一个实例观察这些“fmin*()”是如何找到函数的最小值的。在本例中，要计算最小值的函数f(x,y)为：f(x,y)=(1-x)2+100(y-x2)2为了提高运算速度和精度，有些“fmin()”带有一个fprime参数，它是计算目标函数f对各个自变量的偏导数的函数。20优化—optimizef(x,y)对变量x和y的偏导函数为：这个函数叫做Rosenbrock(罗森布罗克)函数，它经常用来测试最小化算法的收敛速度。它有一个十分平坦的山谷区域，收敛到此山谷区域比较容易，但是在山谷区域搜索到最小点则比较困难。根据函数的计算公式不难看出此函数的最小值是0,在(1,1)处。下面的程序计算f(x,y)的最小值，并且绘制出f所表示的曲面和寻找最小值时的搜索路径。212222400(),200200ffxxyxyxxy优化—optimize观察fmin*函数计算最小值时的路径(scipy_fmin_demo.py)22使用fmin()计算函数最小值，并用matplotlib绘制搜索最小值的路径。“”importscipy.optimizeasoptimportnumpyasnpimportsyspoints=[]deff(p):x,y=pz=(1-x)**2+100*(y-x**2)**2points.append((x,y,z))returnz优化—optimize23deffprime(p):x,y=pdx=-2+2*x-400*x*(y-x**2)dy=200*y-200*x**2returnnp.array([dx,dy])init_point=(-2,-2)try:method=sys.argv[1]except:method=fmin_bfgs”优化—optimize24fmin_func=opt.__dict__[method]ifmethodin[fmin,fmin_powell]:result=fmin_func(f,init_point)#参数为目标函数和初始值elifmethodin[fmin_cg,fmin_bfgs,fmin_l_bfgs_b,fmin_tnc]:result=fmin_func(f,init_point,fprime)#参数为目标函数、初始值和导函数elifmethodin[fmin_cobyla]:r