利用leastsq()函数对数据进行最小二乘算法拟合。

先来看个简单的线性的例子：

#假设要拟合的数据点(xdata,ydata)如下

import numpy as np

from scipy import optimize

import matplotlib.pyplot as plt

xdata=np.array([8.19,2.72,6.39,8.71,4.7,2.66,3.78])

ydata=np.array([7.01,2.78,6.47,6.71,4.1,4.23,4.05])

#先作图看下大概的(xdata,ydata)的形状

plt.plot(xdata, ydata, 'b.', label='data')

#plt.show()

#可得形状如下

#通过散点(xdata,ydata)的形状，可使用线性函数y=k*x+b来拟合数据点。首先定义残差函数y-f(x)

def residual(p):

k,b=p

return ydata-(k*xdata+b)

result=optimize.leastsq(residual,[0,0])

#leastsq函数的调用形式：scipy.optimize.leastsq(func, x0, args=(), Dfun=None, full_output=0, col_deriv=0, ftol=1.49012e-08, xtol=1.49012e-08, gtol=0.0, maxfev=0, epsfcn=None, factor=100, diag=None)

#leastsq函数中的func和x0为必须的参数，其中func为残差函数y-f(x)(该例为residual函数，即ydata-(k*xdata+b))；x0为待定参数估计的初始值，可以验证，leastsq函数对初始估计值不敏感，但是最好还是基于经验取一个大概的值，不然结果会失真，这一点在后面会继续讲到。

#由于在该例的residual函数中，p为两个待定参数k和b的值，所以该例中leastsq函数的x0位置要传入一个二维的列表，而且不能写成r=optimize.leastsq(residual,0,0)。如果写成r=optimize.leastsq(residual,0,0)，那么后面的0会传入leastsq函数的args这个参数中。这里需要注意位置的对应关系。

#所以无论residual函数中有多少个待定参数(该例为k和b两个)，都写成一个类似于该例一样的数组p，然后在调用leastsq函数时，用一个待定参数初始值的估计值组成的数组(该例为[0,0])传入leastsq函数中。

print (result[0])

#最后将数据点和拟合曲线放在一张图中

ydatafit=result[0][0]*xdata+result[0][1]

#k=result[0][0]; b=result[0][1]

plt.plot(xdata,ydatafit, 'r-', label='fit: a=%5.3f, b=%5.3f' % tuple(result[0]))

plt.show()

#再来看个数据稍微复杂的例子：

#假设要拟合的数据点(xdata,ydata)如下，这里通过函数加入噪声来形成要拟合的数据点

import numpy as np

from scipy import optimize

import matplotlib.pyplot as plt

def func(x):

return 2.5 * np.exp(-1.3 * x) + 0.5

xdata = np.linspace(0, 4, 50)

y=func(xdata)

ydata=y+0.2*np.random.normal(size=xdata.size)

#至此，得到了要拟合的数据点(xdata,ydata)，如下图所示

plt.plot(xdata, ydata, 'b.', label='data')

#plt.show()

'''

#这里我们首先以先知者的角色来以形式为y=a* np.exp(-b * x) + c的函数来拟合数据：

def residual(p):

a,b,c=p

return ydata-(a* np.exp(-b * xdata) + c)

#注意这里不要写成returnydata-a* np.exp(-b * xdata) + c

result=optimize.leastsq(residual,[1,1,1])

print (result[0])

#最后将数据点和拟合曲线放在一张图中

a,b,c=result[0]

ydatafit=a* np.exp(-b * xdata) + c

plt.plot(xdata,ydatafit, 'r-', label='fit: a=%5.3f, b=%5.3f, c=%5.3f' % tuple(result[0]))

plt.show()

'''

#以上，我们是以一个先知者的角色来选取函数y=a* np.exp(-b * x) + c来拟合，因为数据点是通过这个函数再加入一些噪声形成的。正常地，我们应该通过观测数据点的分布图，然后想到可以利用类似于一个指数函数的形式y=a* x**b + c来拟合。过程如下

def residual(p):

a,b,c=p

return ydata-(a* xdata**b + c)

result=optimize.leastsq(residual,[1,1,1])

print (result[0])

#最后将数据点和拟合曲线放在一张图中

a,b,c=result[0]

ydatafit=a*xdata**b+c

plt.plot(xdata,ydatafit, 'r-', label='fit: a=%5.3f, b=%5.3f, c=%5.3f' % tuple(result[0]))

plt.show()

#这里，如果我们改变待定参数的初始估计值，如改为result=optimize.leastsq(residual,[1,100,1])，则形状如下

#显然不是我们想要的结果，所以说虽然leastsq()函数对待定参数的初始估计值不敏感，但是估计值也应该在合理范围中。如果出现类似于上图这种不合理的结果，则再次调整初始估计值则可。

总结一下leastsq()函数的使用：1.首先对于给定数据点的散点形状先进行一个分析，大致判断一个函数的形式。2.写成这个函数的残差函数y-f(x)。3.传入两个必须的参数，即残差函数和待定参数的初始估计值进入leastsq()函数中进行拟合，注意待定参数的初始值要取得较为合理。

#讲到将残差函数传入某函数中进行分析，optimize.fsolve()函数也是需要传入残差函数来计算的，该函数主要是来计算非线性方程组

#optimize.fsolve()函数的调用形式为：scipy.optimize.fsolve(func, x0, args=(), fprime=None, full_output=0, col_deriv=0, xtol=1.49012e-08, maxfev=0, band=None, epsfcn=None, factor=100, diag=None)，其中前两个参数func和x0是必须的参数，func为残差函数，x0为待求参数的初始估计值。

#例子如下，要求的非线性方程组如下图

import numpy as np

from scipy import optimize

#定义残差函数

def residual(p):

x0,x1,x2 = p

return[5*np.cos(x1)+3*x0,4*x0**2-2*np.sin(x1*x2),x1*np.sin(x2)-1.5]

result = optimize.fsolve(residual,[1,1,1])

print (result)

#result=[0.06717167 1.61111025 1.94435388]

#此外，optimize.fslove()函数提供了fprime参数(见该函数的调用形式的第三个参数)来传递未知数的雅各比矩阵从而加速计算，传递的雅各比矩阵每i行是第i个方程分别为x0、x1、x2...的导数。如求解上例的代码可以改写为：

import numpy as np

from scipy import optimize

def residual(p):

x0,x1,x2 = p

return[5*np.cos(x1)+3*x0,4*x0**2-2*np.sin(x1*x2),x1*np.sin(x2)-1.5]

def func(p):

x0,x1,x2 = p

return[[3, -5*np.sin(x1), 0], [8*x0, -2*x2*np.cos(x1*x2), -2*x1*np.cos(x1*x2)], [0, np.sin(x2), x1*np.cos(x2)]]

#如[3, -5*np.sin(x1), 0]的三个数则为第一个方程5*np.cos(x1)+3*x0分别对x0、x1、x2的导数，其余类似

result = optimize.fsolve(residual,[1,1,1], fprime=func)

print (result)

#result=[0.06717167 1.61111025 1.94435388]