python 最小二乘法

scipy中最小二乘法函数leastsq的用法

好久没有写Blog了，最近都没有啥好写的。 今天我研究了一下scipy里面的那个最小二乘法的函数的用法，一开始，没弄懂那个函数是怎么调用了，只知道敲进示例程序能用，自己写的程序却报错，后来搜索了一下，看了看别人的代码，搞明白了一点。

基本原理

首先介绍一下最小二乘法的一个原理吧。最小二乘法是一种数学优化方法，通过最小化误差的平方和寻找最佳的匹配函数。一般常用于曲线的拟合，我最早接触到最小二乘法也是在高中的数学课上用来拟合一次函数曲线。

最小二乘法的基本公式如下: (来自百度百科)

∂ϕ∂a0=−2∑(Yi−a0−aiXi)=0

∂ϕ∂ai=−2Xi∑(Yi−a0−aiXi)=0

简单的一次函数拟合

假设一组数据的分布呈线性，我们就可以用一个函数$y=k\times x +b$ 去拟合它，那么那两个系数 k和b怎么求呢？

直接给出公式解: (参考维基百科)

 和  x0=y¯−x1t¯

其中 t¯=1n∑ni=1ti 为t值的算术平均值，也可以解得如下形式:

x1=∑ni=1(ti−t¯)(yi−y¯)∑ni=1(ti−t¯)2

示例程序

我用python根据以上算法写了一个简单的一次函数最小二乘法拟合的程序:

 1 #!/usr/bin/env python
 2 # *-* encoding: utf-8 *-*
 3 #
 4 # =============================================
 5 #      Author   : Andy Scout
 6 #    Homepage   : http://andyhuzhill.github.com
 7 #    E-mail     : andyhuzhill@gmail.com
 8 #
 9 #  Description  :
10 #  Revision     :
11 #
12 # =============================================
13 
14 from __future__ import division
15 
16 def leastsq(x,y):
17     """
18     x,y分别是要拟合的数据的自变量列表和因变量列表
19     """
20     meanx = sum(x) / len(x)   #求x的平均值
21     meany = sum(y) / len(y)   #求y的平均值
22 
23     xsum = 0.0
24     ysum = 0.0
25 
26     for i in range(len(x)):
27         xsum += (x[i] - meanx)*(y[i]-meany)
28         ysum += (x[i] - meanx)**2
29 
30     k = xsum/ysum
31     b = meany - k*meanx
32 
33     return k,b   #返回拟合的两个参数值

使用scipy提供的最小二乘法函数

以上的方法理解起来也比较容易，不过如果需要拟合的函数不是一次函数，就比较麻烦了。python的科学计算包scipy的里面提供了一个函数，可以求出任意的想要拟合的函数的参数。那就是scipy.optimize包里面的leastsq函数。 函数原型是:

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=0.0, factor=100, diag=None, warning=True)

一般我们只要指定前三个参数就可以了。

func 是我们自己定义的一个计算误差的函数，

x0 是计算的初始参数值

args 是指定func的其他参数

下面用一个示例程序来解释他的用法,同样也是使用上面的求一次函数的拟合参数

 1 #!/usr/bin/env python
 2 # *-* encoding: utf-8 *-*
 3 #
 4 # =============================================
 5 #      Author   : Andy Scout
 6 #    Homepage   : http://andyhuzhill.github.com
 7 #    E-mail     : andyhuzhill@gmail.com
 8 #
 9 #  Description  :
10 #  Revision     :
11 #
12 # =============================================
13 
14 import numpy as np
15 from scipy.optimize import leastsq
16 
17 from data import x,y  #从外部导入想要拟合的数据, x和y都是list类型
18 
19 def fun(p, x):
20     """
21     定义想要拟合的函数
22     """
23     k, b = p  #从参数p获得拟合的参数
24     return k*x+b
25 
26 def err(p, x, y):
27     """
28     定义误差函数
29     """
30     return fun(p,x) -y
31 
32 #定义起始的参数 即从 y = 1*x+1 开始，其实这个值可以随便设，只不过会影响到找到最优解的时间
33 p0 = [1,1]   
34 
35 #将list类型转换为 numpy.ndarray 类型，最初我直接使用
36 #list 类型,结果 leastsq函数报错，后来在别的blog上看到了，原来要将类型转
37 #换为numpy的类型
38 
39 x1 = np.array(x)  
40 y1 = np.array(y)
41 
42 xishu = leastsq(err, p0, args=(x1,y1))
43 
44 print xishu[0]
45 
46 # xishu[0]，即为获得的参数