最小二乘法拟合
在函数拟合中，如果用p表示函数中需要确定的参数，那么目标就是找到一组p，使得下面函数S的值最小：
最小二乘
这种算法称为最小二乘法拟合。Python的Scipy数值计算库中的optimize模块提供了leastsq()函数，可以对数据进行最小二乘拟合计算。
此处利用该函数对一段弧线使用圆方程进行了拟合，并通过Matplotlib模块进行了作图，程序内容如下：
#使用最小二乘法拟合圆方程
import numpy as np
from scipy.optimize import leastsq
from pylab import *
x,y=np.loadtxt('num.dat',delimiter=' ',dtype=float,usecols=(0,1),unpack=True)
#for i in x:
# print i
#for j in y:
# print j
def residuals(p):
a,b,r=p
return r**2-(y-b)**2-(x-a)**2
result=leastsq(residuals,[1,1,1])
a,b,r=result[0]
print "a=",a,"b=",b,"r=",r
plot(x,y,color="red",label="origin",linewidth=2)
yfit=sqrt(r**2-(x-a)**2)+b
plot(x,yfit,"b--",label="fit",linewidth=2)
legend(loc="upper right",frameon=False)
annotate('(x-112.892403261)$^2$+(y+58.8238235027)$^2$=110.123575696$^2$',xy=(40,15))
annotate('R=110.123575696',xy=(100,10))
show()
拟合过程说明
导入模块
Python的使用中需要导入相应的模块，此处首先用import语句
import numpy as np
from scipy.optimize import leastsq
from pylab import *
分别导入了numpy, leastsq与pylab模块，其中numpy模块常用用与数组类型的建立，读入等过程。leastsq则为最小二乘法拟合函数。pylab是绘图模块。
读入数据
接下来我们需要读入需要进行拟合的数据，这里使用了numpy.loadtxt()函数：
x,y=np.loadtxt('num.dat',delimiter=' ',dtype=float,usecols=(0,1),unpack=True)
其参数有：
fname：读入文件名，也可省略参数直接写'filename'。
dtype：读入后数据的存放类型，默认为float。
comments：注释符号，读入时忽略以该符号开头的行，默认为'#'。
delimiter：分隔符，存在多列数据时以该符号进行分割，默认为空格。
converters：将列号映射到所需值，可以为遗失的数据提供默认值，默认为不开启。
skiprows：读入时跳过前多少行，默认为0。
usecols：读取第多少列，从0开始计数，可跳跃读取，如usecols=(0,3,8)。1.11.0版更新，取单独一列时可使用usecols=3，并非必须用元组的形式。
unpack：布尔值，如果为真，则对返回的数组进行置换，以便可以使用x,y,z = loadtxt(…)来解压缩数组。当与结构化数据类型一起使用时，将为每个字段返回数组。默认为False。
ndmin：返回的数组将至少具有ndmin维度。否则，一维轴将被压缩。合法值:0(默认值)、1或2。1.6.0版本中的新功能。
encoding：用于解码。1.14.0版本中的新功能。（不是很理解英文原文的意思）
max_rows：skiprows后读入的最大行数，默认读取全部。1.16.0版本中的新功能。
进行拟合
进行拟合时，首先我们需要定义一个目标函数。对于圆的方程，我们需要圆心坐标(a,b)以及半径r三个参数，方便起见用p来存储：
def residuals(p):
a,b,r=p
return r**2-(y-b)**2-(x-a)**2
紧接着就可以进行拟合了，leastsq()函数需要至少提供拟合的函数名与参数的初始值：
result=leastsq(residuals,[1,1,1])
a,b,r=result[0]
print "a=",a,"b=",b,"r=",r
返回的结果为一数组，分别为拟合得到的参数与其误差值等，这里只取拟合参数值。
leastsq()的参数具体有：
func：函数名，该函数应至少具有一个参数并返回多个浮点数，不能返回NaNs
x0：需拟合参数的初始值
Dfun：用行间导数计算func的雅可比行列的函数或方法。如果为空，则估计雅可比行列式。
full_output：布尔值，非零以返回所有可选输出。
col_deriv：布尔值，非零以指定雅可比函数计算列中的导数（更快，因为没有转置操作）。
ftol：平方和中所需的相对误差。
xtol：近似解中所需的相对误差。
gtol：函数向量和雅可比行列之间所需的正交性。
maxfev：对函数的最大调用次数。如果提供了Dfun，则默认maxfev为100*(N+1)，其中N是x0中元素的数量，否则默认maxfev为200*(N+1)。
epsfn：用于确定雅可比行列式的前向误差合适步长的变量（对于Dfun=None）。通常，实际步长将为sqrt(epsfcn)*x如果epsfcn小于机器精度，则假定相对误差为机器精度的数量级。
factor：确定初始步长的参数(factor * || diag * x||)，应在(0.1,100)的区间内。
diag：N个正条目，用作变量的比例因子。
输出选项有：
x：解（或拟合不成功时的最后一次迭代的结果）。
cov_x：使用fjac和ipvt可选输出来构建解决方案周围雅可比的估计值。如果遇到奇异矩阵则表示无（表示在某个方向上的曲率非常平坦）。该矩阵必须乘以残差方差，得到参数估计的协方差.详见curve_fit。
infodict：密钥s的可选输出字典：
- nfev：函数调用的次数。
- fvec：在输出处评估函数。
- fjac：列存储最终近似雅可比矩阵的QR分解的R矩阵的排列。 与
ipvt一起，估计的协方差的可能近似。- ipvt：长度为N的整数数组，它定义了一个置换矩阵p，使得
fjac*p=q*r，其中r是上三角形，对角元素的幅度不增大。 p的列j是单位矩阵的列ipvt(j)。- qtf：向量（转置
(q)*fvec）。
mesg：一个字符串消息，提供有关失败原因的信息。
ies：整数标志。如果它等于1,2,3或4，则发现解决方案。否则，找不到解决方案。在任何一种情况下，可选的输出变量mesg都会提供更多信息。
对结果作图
最后我们可以将原数据与拟合结果一同做成线状图，可采用pylab.plot()函数：
plot(x,y,color="red",label="origin",linewidth=2)
plot(x,yfit,"b--",label="fit",linewidth=2)
legend(loc="upper right",frameon=False)
annotate('R=110.123575696',xy=(100,10))
show()
pylab.plot()函数需提供两列数组作为输入，其他参数可调控线条颜色，形状，粗细以及对应名称等性质。视需求而定，此处不做详解。
pylab.legend()函数可以调控图像标签的位置，有无边框等性质。
pylab.annotate()函数设置注释，需至少提供注释内容与放置位置坐标的参数。
pylab.show()函数用于显示图像。
最终结果如下图所示：
拟合结果
网友评论