• {{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=min⁡S(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。根据多元函数极值的必要条件，

    ∂I∂ak=2∑i=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,j≠k,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=∣∣f∣∣22−∑k=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