文章目录
- 1 什么是最小二乘
- 2 最小二乘原理
- 3 最小二乘应用示例
- 4 法方程到底是什么
1 什么是最小二乘
在科学实验的统计方法研究中,往往要从一组实验数据 ( x i , y i ) ( i = 0 , 1 , 2 , … , m ) (x_i,y_i)(i=0,1,2,…,m) (xi,yi)(i=0,1,2,…,m) 中寻找自变量 x x x 与因变量 y y y 之间的函数关系 y = F ( x ) y=F(x) y=F(x). 由于观测数据往往不准确,因此不要求 y = F ( x ) y=F(x) y=F(x) 经过所有点 ( x i , y i ) (x_i,y_i) (xi,yi),而只要求在给定点 x i x_i xi 上误差 δ i = F ( x i ) − y i ( i = 0 , 1 , 2 , … , m ) δ_i=F(x_i )-y_i (i=0,1,2,…,m) δi=F(xi)−yi(i=0,1,2,…,m)按某种标准最小。若记 δ = ( δ 0 , δ 1 , … , δ m ) T δ=(δ_0,δ_1,…,δ_m)^T δ=(δ0,δ1,…,δm)T,就是要求向量 δ δ δ 的范数最小,通常采用计算较为简单的欧式范数 ‖ δ ‖ 2 ‖δ‖_2 ‖δ‖2 作为误差衡量的标准。
关于最小二乘的一般提法是:对给定的一组数据 ( x i , y i ) ( i = 0 , 1 , 2 , … , m ) (x_i,y_i)(i=0,1,2,…,m) (xi,yi)(i=0,1,2,…,m),要求在函数类 φ { φ 0 , φ 1 , … , φ n } φ\{φ_0,φ_1,…,φ_n\} φ{φ0,φ1,…,φn}中找到一个函数 y = S ∗ ( x ) y=S^* (x) y=S∗(x),使误差平方和最小,即
∥ δ ∥ 2 2 = ∑ i = 0 m δ i 2 = ∑ i = 0 m [ S ∗ ( x i ) − y i ] 2 = m i n S ( x ) ∈ φ ∑ i = 0 m [ S ( x i ) − y i ] 2 − − − − ( F o r m u l a . 1 ) \Vert δ \Vert_2^2=∑_{i=0}^m{δ_i^2}=∑_{i=0}^m[S^* (x_i )-y_i ]^2=min_{S(x)∈φ}∑_{i=0}^m[S(x_i )-y_i ]^2----(Formula.1) ∥δ∥22=i=0∑mδi2=i=0∑m[S∗(xi)−yi]2=minS(x)∈φi=0∑m[S(xi)−yi]2−−−−(Formula.1)
其中,
S ( x ) = a 0 φ 0 ( x ) + a 1 φ 1 ( x ) + ⋯ + a 1 φ 1 ( x ) , ( n < m ) − − − − ( F o r m u l a . 2 ) S(x)=a_0 φ_0 (x)+a_1 φ_1 (x)+⋯+a_1 φ_1 (x) ,(n<m)----(Formula.2) S(x)=a0φ0(x)+a1φ1(x)+⋯+a1φ1(x),(n<m)−−−−(Formula.2)
这就是一般的最小二乘逼近,用几何语言说,就称为曲线拟合的最小二乘法。
2 最小二乘原理
用最小二乘法求拟合曲线时,首先要确定 S ( x ) S(x) S(x) 的形式。这不单纯是数学问题,还与所研究问题的运动规律及所得观测数据 ( x i , y i ) (x_i,y_i) (xi,yi) 有关;通常要从问题的运动规律及给定数据描图,确定 S ( x ) S(x) S(x) 的形式,并通过实际计算选出较好的结果。 S ( x ) S(x) S(x)的一般表达式为 F o r m u l a . 2 Formula.2 Formula.2 式表示的线性形式。为了使问题的提法更有一般性,通常把最小二乘法中 ∥ δ ∥ 2 2 \Vert δ \Vert_2^2 ∥δ∥22都考虑为加权平方和
∥ δ ∥ 2 2 = ∑ i = 0 m ω ( x i ) [ S ( x i ) − f ( x i ) ] 2 − − − − ( F o r m u l a . 3 ) \Vert δ \Vert_2^2=\sum_{i=0}^m\omega(x_i)\left[S(x_i)-f(x_i)\right]^2----(Formula.3) ∥δ∥22=i=0∑mω(xi)[S(xi)−f(xi)]2−−−−(Formula.3)
这里, ω ( x ) ≥ 0 ω(x)≥0 ω(x)≥0是 [ a , b ] [a ,b] [a,b]上的权函数,它表示不同点 ( x i , f ( x i ) ) (x_i,f(x_i)) (xi,f(xi))处的数据比重不同,通常表示在点 ( x i , f ( x i ) ) (x_i,f(x_i)) (xi,f(xi)) 处重复观测的次数。用最小二乘法求拟合曲线的问题,就是在形如 ( F o r m u l a . 2 ) (Formula.2) (Formula.2) 式的 S ( x ) S(x) S(x) 中求一函数 y = S ∗ ( x ) y=S^* (x) y=S∗(x),使 F o r m u l a . 3 Formula.3 Formula.3 式取得最小。可以将这个问题转化为求多元函数
I ( a 0 , a 1 , … , a n ) = ∑ i = 0 m ω ( x i ) [ ∑ j = 0 n a j φ j ( x i ) − f ( x i ) ] 2 − − − − − ( F o r m u l a . 4 ) I(a_0,a_1,…,a_n )=∑_{i=0}^mω(x_i)\left[∑_{j=0}^na_j φ_j (x_i )-f(x_i)\right]^2-----(Formula.4) I(a0,a1,…,an)=i=0∑mω(xi)[j=0∑najφj(xi)−f(xi)]2−−−−−(Formula.4)
的极小值点 ( a 0 ∗ , a 1 ∗ , … , a n ∗ ) (a_0^*,a_1^*,…,a_n^*) (a0∗,a1∗,…,an∗) 的问题。
由求多元函数极值的必要条件,有
∂ I ∂ a k = 2 ∑ i = 0 m ω ( x i ) [ ∑ j = 0 n a j φ j ( x i ) − f ( x i ) ] φ k ( x i ) = 0 , ( k = 0 , 1 , … , n ) − − − − ( F o r m u l a . 5 ) \cfrac{∂I}{∂a_k}=2∑_{i=0}^mω(x_i )\left[∑_{j=0}^na_jφ_j(x_i)-f(x_i)\right] φ_k (x_i )=0,(k=0,1,…,n)----(Formula.5) ∂ak∂I=2i=0∑mω(xi)[j=0∑najφj(xi)−f(xi)]φk(xi)=0,(k=0,1,…,n)−−−−(Formula.5)
若记
( φ j , φ k ) = ∑ i = 0 m ω ( x i ) φ j ( x i ) φ k ( x i ) − − − − ( F o r m u l a . 6 ) (φ_j,φ_k )= ∑_{i=0}^mω(x_i ) φ_j (x_i ) φ_k (x_i )----(Formula.6) (φj,φk)=i=0∑mω(xi)φj(xi)φk(xi)−−−−(Formula.6)
( f , φ k ) = ∑ i = 0 m ω ( x i ) f ( x i ) φ k ( x i ) ≡ d k , ( k = 0 , 1 , … , n ) − − − − ( F o r m u l a . 7 ) (f,φ_k )=∑_{i=0}^mω(x_i ) f(x_i ) φ_k (x_i )≡d_k,(k=0,1,…,n)----(Formula.7) (f,φk)=i=0∑mω(xi)f(xi)φk(xi)≡dk,(k=0,1,…,n)−−−−(Formula.7)
则 F o r m u l a . 4 Formula.4 Formula.4 可改写为
∑ j = 0 n ( φ k , φ j ) a j = d k , ( k = 0 , 1 , … , n ) − − − − ( F o r m u l a . 8 ) ∑_{j=0}^n(φ_k,φ_j)a_j=d_k,(k=0,1,…,n)----(Formula.8) j=0∑n(φk,φj)aj=dk,(k=0,1,…,n)−−−−(Formula.8)
F o r m u l a . 8 Formula.8 Formula.8 式称为法方程,矩阵形式为
G a = d − − − − ( F o r m u l a . 9 ) Ga=d----(Formula.9) Ga=d−−−−(Formula.9)
其中, a = ( a 0 , a 1 , . . . , a n ) T a=(a_0,a_1,...,a_n)^T a=(a0,a1,...,an)T, d = ( d 0 , d 1 , . . . , d n ) T d=(d_0,d_1,...,d_n)^T d=(d0,d1,...,dn)T
G = [ ( φ 0 , φ 0 ) ( φ 0 , φ 1 ) ⋯ ( φ 0 , φ n ) ( φ 1 , φ 0 ) ( φ 1 , φ 1 ) ⋯ ( φ 1 , φ n ) ⋮ ⋮ ⋱ ⋮ ( φ n , φ 0 ) ( φ n , φ 1 ) ⋯ ( φ n , φ n ) ] G=\begin{bmatrix} {(φ_0,φ_0)}&{(φ_0,φ_1)}&{\cdots}&{(φ_0,φ_n)}\\ {(φ_1,φ_0)}&{(φ_1,φ_1)}&{\cdots}&{(φ_1,φ_n)}\\ {\vdots}&{\vdots}&{\ddots}&{\vdots}\\ {(φ_n,φ_0)}&{(φ_n,φ_1)}&{\cdots}&{(φ_n,φ_n)}\\ \end{bmatrix} G=⎣⎢⎢⎢⎡(φ0,φ0)(φ1,φ0)⋮(φn,φ0)(φ0,φ1)(φ1,φ1)⋮(φn,φ1)⋯⋯⋱⋯(φ0,φn)(φ1,φn)⋮(φn,φn)⎦⎥⎥⎥⎤
由于 φ 0 , φ 1 , . . . , φ n φ_0,φ_1,...,φ_n φ0,φ1,...,φn 线性无关,因此 ∣ G ∣ ≠ 0 |G|≠0 ∣G∣=0,方程组 F o r m u l a . 8 Formula.8 Formula.8 存在唯一解
a k = a k ∗ , ( k = 0 , 1 , . . . , n ) − − − − ( F o r m u l a . 10 ) a_k=a_k^*,(k=0,1,...,n)----(Formula.10) ak=ak∗,(k=0,1,...,n)−−−−(Formula.10)
从而得到函数 f ( x ) f(x) f(x) 的最小二乘解为
S ∗ ( x ) = a 0 ∗ φ 0 ( x ) + a 1 ∗ φ 1 ( x ) + . . . + a n ∗ φ n ( x ) − − − − ( F o r m u l a . 11 ) S^*(x)=a_0^*φ_0(x)+a_1^*φ_1(x)+...+a_n^*φ_n(x)----(Formula.11) S∗(x)=a0∗φ0(x)+a1∗φ1(x)+...+an∗φn(x)−−−−(Formula.11)
可以证明 S ∗ ( x ) S^*(x) S∗(x) 对于任何形如 F o r m u l a . 2 Formula.2 Formula.2 式的 S ( x ) S(x) S(x) ,都有
∑ i = 0 m ω ( x i ) [ S ∗ ( x i ) − f ( x i ) ] 2 ≤ ∑ i = 0 m ω ( x i ) [ S ( x i ) − f ( x i ) ] 2 − − − − ( F o r m u l a . 12 ) ∑_{i=0}^mω(x_i ) [S^* (x_i )-f(x_i )]^2≤∑_{i=0}^mω(x_i ) [S(x_i)-f(x_i )]^2 ----(Formula.12) i=0∑mω(xi)[S∗(xi)−f(xi)]2≤i=0∑mω(xi)[S(xi)−f(xi)]2−−−−(Formula.12)
也就是说,只要对一组数据的法方程进行求解,就可以得到唯一一组多项式的系数解。如何求解方程组,将会在后续的博客中给出。
3 最小二乘应用示例
下面通过一个例题进一步理解曲线的最小二乘
4 法方程到底是什么
相信不少人对于法方程 G a = d Ga=d Ga=d 中 G G G 的元素到底是什么存在疑问,那么 G G G 中的 φ φ φ 到底是什么呢?下面通过分析一阶、二阶、三阶多项式拟合的法方程,帮助大家理解这个问题。
回顾一下 G G G
G = [ ( φ 0 , φ 0 ) ( φ 0 , φ 1 ) ⋯ ( φ 0 , φ n ) ( φ 1 , φ 0 ) ( φ 1 , φ 1 ) ⋯ ( φ 1 , φ n ) ⋮ ⋮ ⋱ ⋮ ( φ n , φ 0 ) ( φ n , φ 1 ) ⋯ ( φ n , φ n ) ] G=\begin{bmatrix} {(φ_0,φ_0)}&{(φ_0,φ_1)}&{\cdots}&{(φ_0,φ_n)}\\ {(φ_1,φ_0)}&{(φ_1,φ_1)}&{\cdots}&{(φ_1,φ_n)}\\ {\vdots}&{\vdots}&{\ddots}&{\vdots}\\ {(φ_n,φ_0)}&{(φ_n,φ_1)}&{\cdots}&{(φ_n,φ_n)}\\ \end{bmatrix} G=⎣⎢⎢⎢⎡(φ0,φ0)(φ1,φ0)⋮(φn,φ0)(φ0,φ1)(φ1,φ1)⋮(φn,φ1)⋯⋯⋱⋯(φ0,φn)(φ1,φn)⋮(φn,φn)⎦⎥⎥⎥⎤
需要明确的一点是:若以 x x x 为自变量,则 φ 0 ( x ) = x 0 = 1 , φ 1 ( x ) = x 1 , φ 2 ( x ) = x 2 , … , φ n ( x ) = x n φ_0 (x)=x^0=1,φ_1 (x)=x^1,φ_2 (x)=x^2,… ,φ_n (x)=x^n φ0(x)=x0=1,φ1(x)=x1,φ2(x)=x2,…,φn(x)=xn.
已知一组数据点 P = { p 1 ( x 1 , y 1 ) , p 2 ( x 2 , y 2 ) , … , p m ( x m , y m ) } P=\{p_1 (x_1,y_1 ),p_2 (x_2,y_2 ),…,p_m (x_m,y_m )\} P={p1(x1,y1),p2(x2,y2),…,pm(xm,ym)},每个点只观测一次,即 ω ( x ) ≡ 1 \omega(x)≡1 ω(x)≡1. 分别对其进行一阶、二阶、三阶多项式拟合,对应的拟合函数与法方程如下:
拟合阶数 | 拟合函数 | 法方程 |
---|---|---|
1 | S ( x ) = a 0 + a 1 x S(x)=a_0+a_1 x S(x)=a0+a1x | [ ∑ i = 1 m x i 0 ∑ i = 1 m x i 1 ∑ i = 1 m x i 1 ∑ i = 1 m x i 2 ] \begin{bmatrix}\sum_{i=1}^mx_i^0&\sum_{i=1}^mx_i^1&\\\sum_{i=1}^mx_i^1&\sum_{i=1}^mx_i^2&\\\end{bmatrix} [∑i=1mxi0∑i=1mxi1∑i=1mxi1∑i=1mxi2] [ a 0 a 1 ] \begin{bmatrix}a_0\\a_1\\\end{bmatrix} [a0a1]= [ ∑ i = 1 m x i 0 y i 1 ∑ i = 1 m x i 1 y i 1 ] \begin{bmatrix}\sum_{i=1}^mx_i^0y_i^1\\\sum_{i=1}^mx_i^1y_i^1\\\end{bmatrix} [∑i=1mxi0yi1∑i=1mxi1yi1] (详细) |
1 | S ( x ) = a 0 + a 1 x S(x)=a_0+a_1 x S(x)=a0+a1x | [ m ∑ x ∑ x ∑ x 2 ] \begin{bmatrix}m&\sum x\\\sum x&\sum x^2\\\end{bmatrix} [m∑x∑x∑x2] [ a 0 a 1 ] \begin{bmatrix}a_0\\a_1\\\end{bmatrix} [a0a1]= [ ∑ y ∑ x y ] \begin{bmatrix}\sum y\\\sum xy\\\end{bmatrix} [∑y∑xy] |
2 | S ( x ) = a 0 + a 1 x + a 2 x 2 S(x)=a_0+a_1 x+a_2 x^2 S(x)=a0+a1x+a2x2 | [ m ∑ x ∑ x 2 ∑ x ∑ x 2 ∑ x 3 ∑ x 2 ∑ x 3 ∑ x 4 ] \begin{bmatrix}m&\sum x&\sum x^2 \\\sum x&\sum x^2&\sum x^3\\\sum x^2&\sum x^3&\sum x^4\\\end{bmatrix} ⎣⎡m∑x∑x2∑x∑x2∑x3∑x2∑x3∑x4⎦⎤ [ a 0 a 1 a 2 ] \begin{bmatrix}a_0\\a_1\\a_2\\\end{bmatrix} ⎣⎡a0a1a2⎦⎤= [ ∑ y ∑ x y ∑ x 2 y ] \begin{bmatrix}\sum y\\\sum xy\\\sum x^2y\\\end{bmatrix} ⎣⎡∑y∑xy∑x2y⎦⎤ |
3 | S ( x ) = a 0 + a 1 x + a 2 x 2 + a 3 x 3 S(x)=a_0+a_1 x+a_2 x^2+a_3 x^3 S(x)=a0+a1x+a2x2+a3x3 | [ m ∑ x ∑ x 2 ∑ x 3 ∑ x ∑ x 2 ∑ x 3 ∑ x 4 ∑ x 2 ∑ x 3 ∑ x 4 ∑ x 5 ∑ x 3 ∑ x 4 ∑ x 5 ∑ x 6 ] \begin{bmatrix}m&\sum x&\sum x^2&\sum x^3 \\\sum x&\sum x^2&\sum x^3&\sum x^4\\\sum x^2&\sum x^3&\sum x^4&\sum x^5\\\sum x^3&\sum x^4&\sum x^5&\sum x^6\\\end{bmatrix} ⎣⎢⎢⎡m∑x∑x2∑x3∑x∑x2∑x3∑x4∑x2∑x3∑x4∑x5∑x3∑x4∑x5∑x6⎦⎥⎥⎤ [ a 0 a 1 a 2 a 3 ] \begin{bmatrix}a_0\\a_1\\a_2\\a_3\\\end{bmatrix} ⎣⎢⎢⎡a0a1a2a3⎦⎥⎥⎤= [ ∑ y ∑ x y ∑ x 2 y ∑ x 3 y ] \begin{bmatrix}\sum y\\\sum xy\\\sum x^2y\\\sum x^3y\\\end{bmatrix} ⎣⎢⎢⎡∑y∑xy∑x2y∑x3y⎦⎥⎥⎤ |
… | … | … |
为了便于理解,首先给出一阶多项式详细的法方程。
对比不同阶数的法方程可以看出, G G G 中的每个元素以一定规律对自变量的 k k k 次幂进行求和得到,低阶的法方程是高阶的法方程的一个 “子集”,且矩阵 G G G 为对称正定矩阵。
相关链接
基于Eigen库的线性方程组/矩阵方程求解(方法汇总)
最小二乘曲线拟合的C++实现
参考资料:
《数值分析》李庆扬