• {{article.mytitle}}
  • 1. 曲线拟合的最小二乘法及python实现

    1.1. 返回首页

    首页地址https://kz16.top/na
    本文来源https://kz16.top/na/leastSquare

    1.2. 最小二乘法理论

    1.2.1. 一般最小二乘逼近

    最小二乘法一般提法: 对于给定的一组数据(xi,yi),i=0,1,,m(x_i,y_i),i=0,1,\cdots,m,要求在函数空间Φ=span{ϕ0,ϕ1,,ϕn}\Phi = \text{span}\{ \phi_0,\phi_1,\cdots,\phi_n \}找到一个函数y=S(x)y=S^*(x),使误差加权平方和

    δ22=i=0mω(xi)[S(xi)yi]2=minS(x)ϕi=0mω(xi)[S(xi)yi]2.|| \bm \delta ||_2^2 = \sum_{i=0}^m \omega(x_i)[S^*(x_i)-y_i]^2 = \min_{S(x) \in \phi} \sum_{i=0}^m \omega(x_i) [S(x_i)-y_i]^2 .

    这里ω(x)>0\omega(x)>0是权函数,这里的

    S(x)=a0ϕ0+a1ϕ1++anϕn(x),(n<m).(1)\tag{1} S(x) = a_0\phi_0 + a_1 \phi_1 + \cdots + a_n\phi_n(x) ,\quad(n<m).

    其中a0,a1,,ana_0,a_1,\cdots, a_n为待求系数。

    由(1)式可知,最小二乘法可以进一步转化为求有限维多元函数

    I(a0,a1,,an)=i=0mω(xi)[j=0najϕj(xi)f(xi)]2(2)\tag{2} I(a_0,a_1,\cdots,a_n) = \sum_{i=0}^m \omega(x_i) \left [\sum_{j=0}^n a_j \phi_j (x_i) -f(x_i) \right]^{\color{red}2}

    的极小点(a0,a1,,an)(a_0^*, a_1^*,\cdots, a_n^*)问题。这里f(xi)=yif(x_i) = y_i。根据多元函数极值的必要条件,

    Iak=2i=0mω(xi)[j=0najϕj(xi)f(xi)]ϕk=0,k=0,1,,n.\frac{\partial I}{\partial a_k} = 2 \sum_{i=0}^m \omega(x_i) \left [\sum_{j=0}^n a_j \phi_j (x_i) -f(x_i) \right] \phi_k=0, \quad k=0,1,\cdots,n.

    上式其实就是一个关于a=(a0,a1,,an)T\bm a = (a_0,a_1, \cdots, a_n)^T的线性方程组,矩阵形式为

    Ga=d,(3)\tag{3} \mathbb G \bm a = \bm d ,

    其中

    G=[gkj]k=0,1,,nj=0,1,,n,gkj=i=0mω(xi)ϕk(xi)ϕj(xi)\mathbb G=[g_{kj}]_{k=0,1,\cdots,n}^{j=0,1,\cdots,n}, \qquad g_{kj} = \sum_{i=0}^{m}\omega(x_i)\phi_k(x_i)\phi_j(x_i)

    以及

    d=[dk]k=0,1,,n,dk=i=0mω(xi)ϕk(xi)f(xi)\bm d = [d_k]_{k=0,1,\cdots,n}, \quad d_k = \sum_{i=0}^{m}\omega(x_i)\phi_k(x_i)f(x_i)

    因为基函数ϕ0,ϕ1,,ϕn\phi_0,\phi_1,\cdots,\phi_n线性无关,因此矩阵G\mathbb G是一个满秩矩阵,从而线性方程组(3)存在唯一解。 根据gkj=gjkg_{kj}=g_{jk},因此矩阵G\mathbb G是对称矩阵。

    1.2.2. 用正交函数作最小二乘法

    我们当然最希望线性方程组(3)中的矩阵G\mathbb G是一个对角矩阵,也就是需要满足

    gkj=i=0mω(xi)ϕk(xi)ϕj(xi)={0,jk,Ak>0,j=k.(4)\tag{4} g_{kj} = \sum_{i=0}^{m}\omega(x_i)\phi_k(x_i)\phi_j(x_i) = \begin{cases} 0, & j\neq k,\\ A_k>0, & j=k. \end{cases}
    我们称满足上面(4)式基函数族{ϕj}j=0,1,,n\{\phi_j\}_{j=0,1,\cdots,n}为关于点集{xi}i=0,1,,m\{x_i\}_{i=0,1,\cdots, m}带权{w(xi)}i=0,1,,m\{w(x_i)\}_{i=0,1,\cdots, m}的正交基函数族。

    如果采用这种带权正交基函数,则线性方程组(3)的解为

    ak=dkAk=i=0mω(xi)ϕk(xi)f(xi)i=0mω(xi)ϕk2(xi)a_k^* = \dfrac{d_k}{A_k} = \dfrac{\sum_{i=0}^{m}\omega(x_i)\phi_k(x_i)f(x_i)}{ \sum_{i=0}^{m}\omega(x_i)\phi^2_k(x_i)}

    且平方误差为

    δ22=f22k=0kAk(ak)2.|| \bm \delta ||_2^2 = || f ||_2^2 - \sum_{k=0}^k A_k(a_k^*)^2.

    1.3. 最小二乘法数值例子

    1.3.1. 问题

    已知一组实验数据如下表,求它的拟合曲线

    xix_i 1 2 3 4 5
    f(xi)f(x_i) 4 4.5 6 8 8.5
    频数 2 1 3 1 1

    1.3.2. 分析

    很显然,这里的频数可以当成权,这组实验数据的点分布在一条直线附近,因此可以选择线性拟合,基函数取

    ϕ0=1,ϕ1=x\phi_0 = 1, \quad \phi_1= x

    1.3.3. python代码

    import numpy as np
    
    def phi(x):
        return np.array([1.0+0*x, x])
    
    x = np.array([1, 2, 3, 4, 5])
    y = np.array([4, 4.5, 6, 8, 8.5])
    w = np.array([2, 1, 3, 1, 1])
    
    phixsw = phi(x)*np.sqrt(w)
    leftG = phixsw.dot(phixsw.T)
    rightd = phi(x).dot(y*w)
    
    ans = np.linalg.solve(leftG, rightd)
    print("Left matrix G = ")
    print(leftG)
    print("Right vector G = ")
    print(rightd)
    print("Solver ans = ")
    print(ans)
    

    结果

    Left matrix G =
    [[ 8. 22.]
     [22. 74.]]
    Right vector G =
    [ 47.  145.5]
    Solver ans =
    [2.56481481 1.2037037 ]
    

    因此拟合曲线为
    f(x)=2.56481481+1.2037037xf(x) = 2.56481481 + 1.2037037x

    注1: 用python或matlab运算为了提高效率,可以把问题改成矩阵向量的运算,从而可以直接使用内置函数

    1.4. 返回首页

    首页地址https://kz16.top/na
    本文来源https://kz16.top/na/leastSquare


    仅供参考,部分漏洞在所难免
    欢迎转载,转载请指明来源,
    请勿用于商业
    作者:Zhou Bingzhen
    2019.11.21