3.2.2 最小二乘拟合

假设有一组实验数据(xi ,yi )，我们事先知道它们之间应该满足某函数关系：yi =f(xi )。通过这些已知信息，需要确定函数f( )的一些参数。例如，如果函数f( )是线性函数f(x)=kx+b，那么参数k和b就是需要确定的值。

如果用p表示函数中需要确定的参数，则目标是找到一组p使得函数S的值最小：

这种算法被称为最小二乘拟合(least-square fitting)。在optimize模块中，可以使用leastsq( )对数据进行最小二乘拟合计算。leastsq( )的用法很简单，只需要将计算误差的函数和待确定参数的初始值传递给它即可。下面是用leastsq( )对线性函数进行拟合的程序：

    import numpy as np
    from scipy import optimize
    
    X = 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])
    
    def residuals(p): ❶
        "计算以p为参数的直线和原始数据之间的误差"
        k, b = p
        return Y - (k*
X + b)
    
    # leastsq使得residuals( )的输出数组的平方和最小，参数的初始值为[1,0]
    r = optimize.leastsq(residuals, [1, 0]) ❷
    k, b = r[0]
    print "k =",k, "b =",b
    k = 0.613495349193 b = 1.79409254326

图3-2（左）直观地显示了原始数据、拟合直线以及它们之间的误差。❶residuals( )的参数p是拟合直线的参数，函数返回的是原始数据和拟合直线之间的误差。图中用数据点到拟合直线在Y轴上的距离表示误差。❷leastsq( )使得这些误差的平方和最小，即图中所有正方形的面积之和最小。

图3-2 最小化正方形面积之和（左），误差曲面（右）

由前面的函数S的公式可知，对于直线拟合来说，误差的平方和是直线参数k和b的二次多项式函数，因此可以用如图3-2（右）所示的曲面直观地显示误差平方和与两个参数之间的关系。图中用红色圆点表示曲面的最小点，它的X-Y轴的坐标就是leastsq( )的拟合结果。

接下来，让我们再看一个对正弦波数据进行拟合的例子：

    def func(x, p): ❶
        """
    数据拟合所用的函数: A*
sin(2*
pi*
k*
x + theta)
        """
        A, k, theta = p
        return A*
np.sin(2*
np.pi*
k*
x+theta)   
    
    def residuals(p, y, x): ❷
        """
    实验数据x, y和拟合函数之间的差，p为拟合需要找到的系数
        """
        return y - func(x, p)
    
    x = np.linspace(0, 2*
np.pi, 100)
    A, k, theta = 10, 0.34, np.pi/6 # 真实数据的函数参数
    y0 = func(x, [A, k, theta]) # 真实数据
    # 加入噪声之后的实验数据
    np.random.seed(0)
    y1 = y0 + 2 *
 np.random.randn(len(x)) ❸
    
    p0 = [7, 0.40, 0] # 第一次猜测的函数拟合参数
    
    # 调用leastsq进行数据拟合
    # residuals为计算误差的函数
    # p0为拟合参数的初始值
    # args为需要拟合的实验数据
    plsq = optimize.leastsq(residuals, p0, args=(y1, x)) ❹
    
    print u"真实参数:", [A, k, theta] 
    print u"拟合参数", plsq[0] # 实验数据拟合后的参数
    
    pl.plot(x, y1, "o", label=u"带噪声的实验数据")
    pl.plot(x, y0, label=u"真实数据")
    pl.plot(x, func(x, plsq[0]), label=u"拟合数据")
    pl.legend(loc="best")
    真实参数: [10, 0.34, 0.5235987755982988]
    拟合参数 [ 10.25218748   0.3423992    0.50817424]

图3-3显示了带噪声的正弦波拟合。

图3-3 带噪声的正弦波拟合

程序中，❶要拟合的目标函数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是表示实验数据的数组。

对于这种一维曲线拟合，optimize库还提供了一个curve_fit( )函数，下面使用此函数对正弦波数据进行拟合。它的目标函数与leastsq( )稍有不同，各个待优化参数直接作为函数的参数传入。

    def func2(x, A, k, theta):
        return A*
np.sin(2*
np.pi*
k*
x+theta)   
    
    popt, _ = optimize.curve_fit(func2, x, y1, p0=p0)
    print popt
    [ 10.25218748   0.3423992    0.50817424]

如果频率的初值和真实值的差别较大，拟合结果中的频率参数可能无法收敛于实际的频率。在下面的例子中，由于频率初值的选择不当，导致curve_fit( )未能收敛到真实的参数。这时可以通过其他方法先估算一个频率的近似值，或者使用全局优化算法。在后面的例子中，我们会使用全局优化算法重新对正弦波数据进行拟合。

    popt, _ = optimize.curve_fit(func2, x, y1, p0=[10, 1, 0])
    print u"真实参数:", [A, k, theta]
    print u"拟合参数", popt
    真实参数: [10, 0.34, 0.5235987755982988]
    拟合参数 [ 0.71093473  1.02074599 -0.1277666 ]
上一章目录下一章