算法基础 - 最小二乘法（线性拟合）

1、最小二乘法的历史

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

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

2、最小二乘法的原理

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

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

使误差（所谓误差，当然是观察值与实际真实值的差量）平方和达到最小以寻求估计值的方法，就叫做最小二乘法，用最小二乘法得到的估计，叫做最小二乘估计。当然，取平方和作为目标函数只是众多可取的方法之一。这里说的实际真实值是指在运算过程中使用的**最小二乘估计(即：估算出的参数值)**计算出来的值。

3、典型应用场景 - 线性拟合

从一个简单的例子开始，已知平面上有三个点（1，2），（0，2），（2，3），我们想用一条直线去拟合它，像高中时一样，设这条直线的方程为 Y=kx+b (一次函数)，我们希望这条直线可以同时通过这三个点。

线性拟合是一种常见的数据分析方法，用于找到一条最佳拟合直线来描述数据点的趋势。

最小二乘法是一种通过最小化残差平方和来拟合数据的方法。

最小二乘法是一种常用的拟合方法，它通过最小化实际观测值与拟合值之间的残差平方和来确定拟合直线的参数。在线性拟合中，我们假设拟合直线的公式为y = kx + b，其中k是斜率，b是截距。

模型如下：

假设从总体中获取了n组观察值（X1，Y1），（X2，Y2）， …，（Xn，Yn）。对于平面中的这n个点，可以使用无数条曲线来拟合。综合起来看，这条直线处于样本数据的中心位置最合理。 选择最佳拟合曲线的标准可以确定为：使总的拟合误差（即总残差）达到最小。有以下三个标准可以选择：

（1）用“残差和最小”确定直线位置是一个途径。但很快发现计算“残差和”存在相互抵消的问题。

（2）用“残差绝对值和最小”确定直线位置也是一个途径。但绝对值的计算比较麻烦。

（3）最小二乘法的原则是以“残差平方和最小”确定直线位置。用最小二乘法除了计算比较方便外，得到的估计量还具有优良特性。这种方法对异常值非常敏感。

最常用的是普通最小二乘法（ Ordinary Least Square，OLS）：所选择的回归模型应该使所有观察值的残差平方和达到最小。（Q为残差平方和）

样本回归模型：

残差平方和：

则通过Q最小确定这条直线，即确定

以为变量，把它们看作是Q的函数，就变成了一个求极值的问题，可以通过求导数得到。求Q对两个待估参数的偏导数：

解得：

4、线性拟合 C++ 示例1 ： 求解 y = kx + b

#include<iostream>
#include<fstream>
#include<vector>
using namespace std;

class LeastSquare
{

private:
    double k, b;

public:
    LeastSquare(const vector<double>& x, const vector<double>& y)
    {
      double sumX = 0, sumY = 0, sumX2 = 0, sumXY = 0;

      for(int i=0; i< x.size(); ++i)
      {
        sumX += x[i];
        sumY += y[i];

        sumX2 += x[i] * x[i];
        sumXY += x[i] * y[i];
      }

      //依据推导公式: (n * sumXY - sumX * sumY) / (n * sumX2 - sumX * sumX)
      k = (sumXY * x.size() - sumX * sumY) / (sumX2 * x.size() - sumX * sumX);

      //依据推导公式: (sumX2 * sumY - sumX * sumXY) / (n * sumX2 - sumX * sumX)
      b = (sumX2 * sumY - sumX * sumXY) / (sumX2 * x.size() - sumX * sumX);

    }

    double getY(const double x) const
    {
      return k * x + b;
    }

    void print() const
    {
      std::cout << "y = " << k << "x + " <<b << "\n";
    }

};

int LeastSquare002()
{
    std::vector<double> x = {10, 15, 20, 25, 30, 35, 40, 45};
    std::vector<double> y = {2000.36, 2000.50, 2000.72, 2000.80, 2001.07, 2001.25, 2001.48, 2001.60};

    int sz = x.size();
    x.resize(sz / 2);
    LeastSquare ls(x, y);
    ls.print();

    cout<<"Input x:\n";
    double x0;
    while(cin>>x0)
    {
      cout<<"y = "<<ls.getY(x0)<<endl;
      cout<<"Input x:\n";
    }

  }

  return 0;

}

5、线性拟合 C++ 示例2 ： 求解 y = kx + b

#include <iostream>
#include <vector>
#include <numeric>

using Parameter = struct {
    double k; // 斜率
    double b; // 截距
};

// 最小二乘法计算过程
bool LeastSquares(std::vector<double>& X, std::vector<double>& Y, Parameter& param)
{
  if (X.empty() || Y.empty())
    return false;

  int n = X.size();

  double sumX = std::accumulate(X.begin(), X.end(), 0.0);
  double sumY = std::accumulate(Y.begin(), Y.end(), 0.0);
  std::cout << "sumX = " << sumX << ", sumY = " << sumY << std::endl;

  double meanX = sumX / n;
  double meanY = sumY / n;
  std::cout << "meanX = " << meanX << ", meanY = " << meanY << std::endl;


  double sumXY = 0.0;
  double sumX2 = 0.0;
  for (int i = 0; i < n; i++)
  {
    sumXY += X[i] * Y[i]; //对x与y的乘积求和
    sumX2 += X[i] * X[i]; //对x的平方求和
  }
  std::cout << "sumXY = " << sumXY << ", sumX2 = " << sumX2 << std::endl;

  //依据推导公式: (n * sumXY - sumX * sumY) / (n * sumX2 - sumX * sumX)
  param.k = (sumXY - n * meanX * meanY) / (sumX2 - n * meanX * meanX);

  //依据推导公式: y = kx + b ---> b = y - kx
  param.b = meanY - param.k * meanX;

  return true;
}

int Test001LeastSquares()
{
  std::vector<double> X = {1, 2, 3, 4, 5, 6, 7, 8, 9};
  std::vector<double> Y = {3, 5, 7, 9, 11, 13, 15, 17, 19};
  Parameter param;


  std::cout << " =============================== " << std::endl;

  if (LeastSquares(X, Y, param))
  {
    std::cout << "拟合直线的方程为: y = " << param.k << "x + " << param.b << std::endl;
  }
  else
  {
    std::cout << "拟合失败" << std::endl;
  }

  std::cout << " =============================== \n\n\n " << std::endl;

  return 0;
}