最小二乘法

文章目录
原文地址:https://ngunauj.github.io

最小二乘法(又称最小平方法)是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据,并使得这些求得的数据与实际数据之间误差的平方和为最小。
“最小二乘法”是对过度确定系统,即其中存在比未知数更多的方程组,以回归分析求得近似解的标准方法。在这整个解决方案中,最小二乘法演算为每一方程式的结果中,将残差平方和的总和最小化。
以上来源维基百科。


从200多年前,最小二乘法一直是科技工作者整理实验数据、建立数学模型的有力工具。
最小二乘法主要应用在曲线拟合方面。
1.直线拟合
已知x-y平面上有一些点,怎样作出一条直线\(y=a+bx\)去拟合这些数据点,反映这组数据的总趋势?
从观察可以看出直线不可能通过所有点,在\(x_i\)处有误差:
$$e_i=y_i-(a+bx_i)(i=1,2…,n).$$
称\(e_i\)为直线\(y=a+bx\)在\(x_i\)处的残差,它反映了这条拟合直线在点\(x_i\)处拟合的好坏。怎样衡量总体上拟合的好坏?常用的方法有三种:
(1)残差的最大绝对值\(\max \limits_{1\le i \le n} |e_i| \);
(2) 残差的绝对值之和\(\sum _{i=1}^n|e_i|\)
(3) 残差的平方和\(\sum _{i=1}^ne^2_i\)
以上三种准则1,2两种由于含有绝对值计算比较复杂,一般都使用3准则。
那么原问题等价于求$$Q=\sum _{i=1}^n(y_i-(a+bx_i))^2$$的最小值。
由于Q与一次式有关,也就是与一次式的系数a,b有关,所以Q是a,b的二元函数。求一次式使Q最小,就是求二元函数Q的最小点。用微分法球Q的最小值点。令\(\frac{\partial Q}{\partial a}=0,\frac{\partial Q}{\partial b}=0\)
最后计算得到如下方程组。
$$
\begin{cases}
na+(\sum _{i=1}^nx_i)b=\sum _{i=1}^ny_i,\\
(\sum _{i=1}^nx_i)a+(\sum _{i=1}^nx_i^2)b=\sum _{i=1}^ny_ix_i\end{cases}
$$
直接计算求和公式的值,然后解方程组就得到a,b的值了。
哈哈,记得初高中那时候不知道偏微分,拟合直线直接记公式:
112
2.多项式拟合
多项式拟合的最小二乘法:
求一个m次多项式\(P_m(x)= \sum _{j=0}^ma_jx^j\)使残差平方和
$$Q=\sum _{i=1}^n(y_i- \sum _{j=0}^ma_jx^j_i)^2$$为最小,称这样作出的\(P_m(x)\)为多项式拟合的最小二乘法解。
求\(P_m(x)\)的过程同直线拟合一样(直线拟合不就是一次多项式拟合嘛),求m+1元函数的最小值点(使m+1元函数取得最小值的\(a_j,j=0,1,…m\))。令\(\frac{\partial Q}{\partial a_k}=0(k=0,1,…,m)\),计算偏导数后得:
$$\sum _{i=1}^n(y_i- \sum _{j=0}^ma_jx^j_i)x_i^k=0(k=0,1,…m)$$移项得
$$\sum _{i=1}^n(\sum _{j=0}^ma_jx^j_i)x_i^k=\sum _{i=1}^ny_ix_i^k(k=0,1,…m)$$交换顺序得$$\sum _{j=0}^m(\sum _{i=1}^nx_i^{j+k})a_j=\sum _{i=1}^ny_ix_i^k(k=0,1,…m)$$上面等式可以改写成如下矩阵相乘的形式
113
这样看着能舒服点,然后就是解方程组。分别计算\(\sum _{i=1}^nx_i^j(j=0,1,2,…2m)\)和\(\sum _{i=1}^nx_i^jy_i(j=0,1,2…m)\)带入方程组求解。(原式子可以写成矩阵形式XA=Y 那么A=\((X^TX)^{-1}X^TY\,如果不是方阵逆矩阵用伪逆矩阵代替,这个公式可以大胆去用),上面的方程组也可以用高斯消元去解.
事实上python numpy包中已经提供了用于解线性方程组XA=Y的方法numpy.linalg.solve(),或者matlab中直接X\Y求得。当然matlab中也提供了曲线拟合的方法lsqcurvefit()。
3.指数类型曲线拟合
已知一组数据点(采样点)(xi,yi)(i=1,2,3…n),且已知这组数据满足某一函数原型(经验公式)\(y=F(a,x)\),其中向量a待定.曲线拟合的最小二乘法就是要确定待定向量a,使残差平方和$$Q=\sum _{i=1}^n(y_i-F(a,x_i))^2$$最小。称使总误差最小的F(a,x)为曲线拟合的最小二乘解。
这里假设\(F(a,x)=a_1e^{a_2/x}\)。即\(y=a_1e^{a_2/x}\)
由于F(a,x)是指数函数,所以可以先把它化为线性模型,再按残差平方和最小规则确定向量a。具体做法如下:
对于函数y两边取对数,得:
$$lny=lna_1+\frac{a_2}{x}等价于c_0+c_1x^{-1}$$这是一个线性模型,按直线拟合的最小二乘法,求\(c_0,c_1\),使总误差$$Q=\sum _{i=1}^n(c_0+c_1x_i^{-1}-lny_i)^2$$为最小。
$$
\begin{cases}
nc_0+(\sum _{i=1}^nx_i^{-1})c_1=\sum _{i=1}^nlny_i,\\
(\sum _{i=1}^nx_i^{-1})c_0+(\sum _{i=1}^n(x_i^{-1})^2)c_1=\sum _{i=1}^nlny_ix_i^{-1}\end{cases}
$$把数据带入,求\(c_0,c_1\)然后再求\(a_1,a_2\)最后求得原函数。

4.平面拟合
平面方程的一般表达式为:
$$Ax+By+Cz+d=0(C\neq0)$$
$$z=-\frac{A}{C}-\frac{B}{C}-\frac{D}{C}$$
记:\(a0=-\frac{A}{C},a1=-\frac{B}{C},a2=-\frac{D}{C}\)则\(z=a0x+a1y+a2\)
平面方程拟合:
对于一系列的n个点\((n \ge3)\):
\((x_i,y_i,z_i),i=0,1,…n-1\)
使残差的平方和
$$S=\sum _{i=1}^n(a_0x+a_1y+a_2-z)^2$$最小,应满足:\(\frac{\partial S}{\partial a_k}=0(k=0,1,2)\)
推导后可以得到如下方程组
\(\left[ \begin{matrix} \sum x_i^2 \ \sum x_iy_i \ \sum x_i \\ \sum x_iy_i \ \sum y_i^2 \ \sum y_i \\ \sum x_i \ \sum y_i \ n
\end{matrix}\right]\)\(\left[ \begin{matrix} a_0 \\ a_1 \\
a_2 \end{matrix}\right]\)=\(\left[ \begin{matrix} \sum x_iz_i \\ \sum y_iz_i \\ \sum z_i \end{matrix}
\right]\)
可以看成X*A=Y已知XY求A,由于X存在逆矩阵那么代码就好写多了。下面是之前写的代码纯c++.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
/* ***********************************************
Author :guanjun
Created Time :2017/4/1 18:09:46
File Name :LeastSquare.cpp
************************************************ */
#include <bits/stdc++.h>
struct node{
double x,y,z;
};
using namespace std;
vector<node>point;
class LeastSquare{
double A[3][3];
double B[3][1];
public:
double a,b,c,d;
// z=ax+by+c
void getParameter(vector<node>&point){
memset(A,0,sizeof A);
memset(B,0,sizeof B);
int n=point.size();
for(int i=0;i<n;i++){
A[0][0]+=point[i].x*point[i].x;
A[0][1]+=point[i].x*point[i].y;
A[0][2]+=point[i].x;
A[1][1]+=point[i].y*point[i].y;
A[1][2]+=point[i].y;
}
A[1][0]=A[0][1];
A[2][0]=A[0][2];
A[2][1]=A[1][2];
A[2][2]=n*1.0;
for(int i=0;i<n;i++){
B[0][0]+=point[i].x*point[i].z;
B[1][0]+=point[i].y*point[i].z;
B[2][0]+=point[i].z;
}
double des[3][3],a0,a1,a2;
if(!getMatrixInverse(A,3,des)){
cout<<"奇异矩阵!"<<endl;return ;
}
double arr[3]={0.0};
for(int i=0;i<3;i++){
for(int j=0;j<3;j++){
arr[i]+=des[i][j]*B[j][0];
}
}
cout<<"inv(A): \n";
for(int i=0;i<3;i++){
for(int j=0;j<3;j++){
cout<<des[i][j]<<" ";
}
cout<<endl;
}
a=-arr[0];b=-arr[1];d=-arr[2];c=1.0;
cout<<"a0 "<<arr[0]<<" "<<arr[1]<<" "<<arr[2]<<endl;
}
double getdetA(double a[3][3],int n){
if(n==1){
return a[0][0];
}
double ans=0;
double tmp[3][3]={0.0};
for(int k=0;k<n;k++){
for(int i=0;i<n-1;i++){
for(int j=0;j<n-1;j++){
tmp[i][j]=a[i+1][(j<k)?j:j+1];
}
}
double t=getdetA(tmp,n-1);
if(k&1){
ans-=a[0][k]*t;
}
else ans+=a[0][k]*t;
}
return ans;
}
void getAstar(double a[3][3],int n,double ans[3][3]){
// ans 是a的伴随矩阵
if(n==1){
ans[0][0]=1;
return ;
}
double tmp[3][3];
for(int i=0;i<n;i++){
for(int j=0;j<n;j++){
for(int k=0;k<n-1;k++){
for(int t=0;t<n-1;t++){
tmp[k][t]=a[(k<i)?k:k+1][(t<j)?t:t+1];
}
}
ans[j][i]=getdetA(tmp,n-1);
if((i+j)&1){
ans[j][i]=-ans[j][i];
}
}
}
cout<<"A* "<<endl;
for(int i=0;i<3;i++){
for(int j=0;j<3;j++){
cout<<ans[i][j]<<" ";
}
cout<<endl;
}
}
bool getMatrixInverse(double A[3][3],int n,double des[3][3]){
double detA=getdetA(A,n);
double t[3][3];
if(detA==0){
return false;
}
else{
getAstar(A,n,t);
for(int i=0;i<n;i++){
for(int j=0;j<n;j++){
des[i][j]=t[i][j]/detA;
}
}
return true;
}
}
};
int main()
{
ifstream in("21out.txt");
if(!in){
cout<<"error,can't open\n";
}
else{
node tmp;
point.clear();
while(in>>tmp.x>>tmp.y>>tmp.z){
point.push_back(tmp);
}
LeastSquare ls;
ls.getParameter(point);
cout<<"line Parameter:\n";
cout<<ls.a<<" "<<ls.b<<" "<<ls.c<<" "<<ls.d<<endl;
}
in.close();
return 0;
}

×

纯属好玩

扫码支持
扫码打赏,你说多少就多少

打开支付宝扫一扫,即可进行扫码打赏哦

文章目录
,