c++椭圆最小二乘法原理_算法知识：最小二乘法理论及其编程实现

【最小二乘法的历史】

1801年，意大利天文学家朱赛普·皮亚齐发现了第一颗小行星谷神星。经过40天的跟踪观测后，由于谷神星运行至太阳背后，使得皮亚齐失去了谷神星的位置。随后全世界的科学家利用皮亚齐的观测数据开始寻找谷神星，但是根据大多数人计算的结果来寻找谷神星都没有结果。

时年24岁的高斯也计算了谷神星的轨道。奥地利天文学家海因里希·奥尔伯斯根据高斯计算出来的轨道重新发现了谷神星。高斯使用的最小二乘法的方法发表于1809年他的著作《天体运动论》中。法国科学家勒让德于1806年独立发明“最小二乘法”，但因不为世人所知而默默无闻。勒让德曾与高斯为谁最早创立最小二乘法原理发生争执。1829年，高斯提供了最小二乘法的优化效果强于其他方法的证明，因此被称为高斯-马尔可夫定理。

【最小二乘法的原理】

我们口头中经常说：一般来说，平均来说。如平均来说，不吸烟的健康优于吸烟者，之所以要加“平均”二字，是因为凡事皆有例外，总存在某个特别的人，他吸烟但由于经常锻炼所以他的健康状况可能会优于他身边不吸烟的朋友。而最小二乘法的一个最简单的例子便是算术平均。

最小二乘法(又称最小平方法)是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。

用函数表示为：

使误差(所谓误差，当然是观察值与实际真实值的差量)平方和达到最小以寻求估计值的方法，就叫做最小二乘法，用最小二乘法得到的估计，叫做最小二乘估计。当然，取平方和作为目标函数只是众多可取的方法之一。

最小二乘法的一般形式可表示为：

有效的最小二乘法是勒让德在 1805

年发表的，基本思想就是认为测量中有误差，所以所有方程的累积误差为：

我们求解出导致累积误差最小的参数即可：

勒让德在论文中对最小二乘法的优良性做了几点说明：•最小二乘使得误差平方和最小，并在各个方程的误差之间建立了一种平衡，从而防止某一个极端误差取得支配地位；•计算中只要求偏导后求解线性方程组，计算过程明确便捷；•最小二乘可以导出算术平均值作为估计值。对于最后一点，从统计学的角度来看是很重要的一个性质。推理如下：假设真值为

θ, x1,?,xn 为 n 次测量值, 每次测量的误差为ei=xi?θ，按最小二乘法，误差累积为

求解 θ 使达到最小 L(θ)，正好是算术平均 。

由于算术平均是一个历经考验的方法，而以上的推理说明，算术平均是最小二乘的一个特例，所以从另一个角度说明了最小二乘方法的优良性，使我们对最小二乘法更加有信心。

【最小二乘法的几何解释】

从一个简单的例子开始，已知平面上有三个点(1，2)，(0，2)，(2，3)，我们想用一条直线去拟合它，像高中时一样，设这条直线的方程为

Y=kx+b (一次函数)，我们希望这条直线可以同时通过这三个点，也就是这条直线的参数要满足：

学过初中解方程组的同学知道，这个方程组是无解的。怎么解一个无解的方程组呢？

为了解释的方便，我们用 X1 表示 k，用 X2 表示

b，则：

写出矩阵的形式：

从列的角度看：

一旦写出列的形式，我们很自然的想到把向量 a1、a2 和 b

画到图上。

要找到解，就要找到 a1 和 a2 的一个线性组合，使得组合的向量刚好等于

b，可惜任何的 a1 和 a2 的线性组合，只能出现在 a1 和 a2

所在的平面s上(高中必修二第二章平面基本性质二的推论二两条相交直线确定一个平面，a1 和 a2 可以认为是两条相交直线)，但是向量 b

不在平面 s 上，不可能找到解，怎么办呢？

找不到完美的解，就只能找到一个最接近的解，所以我们想在平面 S 上找到一个最接近向量

b 的向量来代替向量 b，记这个替代品位向量 p，我们知道最接近的肯定是它的投影，即过直线 b 的终点做平面 S

的垂线，垂足就是替代向量 p 的终点，p 和 b 之间的误差 e=b-p。

原来的方程是无解的，我们用 p 代替 b 后，p 在 a1 和 a2

的所在的面上，所以现在方程就一定有解了。

【最小二乘法的求解流程】

要求的解就是式子中间带帽子的 X(是真实值的近似解或者说最接近的那个值)，要使 p 和 b 之间的差距(误差

e)最小，那么 e 一定是垂直于平面 S 的，也就是要垂直于 a1 和 a2

的，想想我们高中是怎么表示两个向量垂直的？只要他们的点乘等于 0 就行了。也就是：

用矩阵表示就是：

即：

把上式代入值前面的e的方程得到最终的结果：

化简一下，反解出 x，最终最佳的近似解就是：

【最小二乘法的 C++ 语言实现】

#include

#include

#include

using namespace std;

class LeastSquare{

double a,

b;

public:

LeastSquare(const vector& x, const vector& y)

{

double t1=0, t2=0, t3=0, t4=0;

for(int i=0; i

{

t1 += x[i]*x[i];

t2 += x[i];

t3 += x[i]*y[i];

t4 += y[i];

}

a = (t3*x.size() - t2*t4) / (t1*x.size() -

t2*t2); // 求得β1

b = (t1*t4 - t2*t3) / (t1*x.size() -

t2*t2); // 求得β2

}

double getY(const double x) const

{

return a*x + b;

}

void print() const

{

cout<

"<

}

};

int main(int argc, char *argv[])

{

if(argc !=

2)

{

cout<

return -1;

}

else

{

vector x;

ifstream in(argv[1]);

for(double d; in>>d; )

x.push_back(d);

int sz = x.size();

vector y(x.begin()+sz/2, x.end());

x.resize(sz/2);

LeastSquare ls(x, y);

ls.print();

cout<

double x0;

while(cin>>x0)

{

cout<

cout<

}

}

}