[计算机数值分析]直线拟合的最小二乘法

Spring-_-Bear 的 CSDN 博客导航

---

对于给定的数据点（x_i, y_i）（i = 1, 2, …, N）求作一次式 $y=a+bx$，使总误差

$$Q=\sum _{i=1}^N{\left[y_i-(a+b{x}_i)\right]}^2$$

为最小。

这个问题是不难求解的。由微积分中求极值的方法知，使 Q 达到极值的参数 a，b 应满足

$$\frac{\partial Q}{\partial a}=0，\frac{\partial Q}{\partial b}=0$$

即成立

$$\{\begin{array}{l}aN+b\sum x_i=\sum y_i\\ a\sum x_i+b\sum x_i^2=\sum x_i{y}_i\end{array}\text{\hspace{1em}(1)}$$

式中 $\sum$ 表示关于下表 i 从 1 到 N 的求和。

例：炼钢是个氧化脱碳的过程，钢液含碳量的多少直接影响冶炼时间的长短。下面是某平炉的生产记录，表中 i 为试验次数，x_i 为全部炉料熔化完毕时钢液的含碳量，y_i 为熔毕至出钢的冶炼时间（以分为单位）。

i	1	2	3	4	5
xi	165	123	150	123	141
yi	187	126	172	125	148

把表中所给数据画在坐标纸上，将会看到，数据点的分布可以用一条直线来近似描述。设所求的拟合直线为 $y=a+bx$，则方程组（1）的具体形式是

$$\{\begin{array}{l}5a+702b=758\\ 702a+99864b=108396\end{array}$$

解出 a，b，即得拟合直线 $y=-60.0392+1.5138x$。

运行示例：

程序源码：

#include <iostream>

using namespace std;

#define MAX 10

typedef struct Point
{
    double x;
    double y;
} point;

int main(void)
{
    int N;
    cout << "请输入点的个数:";
    cin >> N;

    point p[MAX];
    for (int i = 1; i <= N; i++)
    {
        cout << "请输入第 " << i << " 个点的坐标:";
        cin >> p[i].x;
        cin >> p[i].y;
    }

    // ∑X, ∑Y, ∑X², ∑XY
    double sumX = 0;
    double sumY = 0;
    double sumX2 = 0;
    double sumXY = 0;
    for (int i = 1; i <= N; i++)
    {
        sumX += p[i].x;
        sumY += p[i].y;
        sumX2 += p[i].x * p[i].x;
        sumXY += p[i].x * p[i].y;
    }

    cout << "\n一元二次方程方程如下：" << endl;
    // 根据 sumX 的正负性控制输出 + 或 -
    if (sumX > 0)
    {
        cout << N << " * a + " << sumX << " = " << sumY << endl;
    }
    else
    {
        cout << N << " * a " << sumX << " = " << sumY << endl;
    }
    cout << sumX << " * a + " << sumX2 << " = " << sumXY << endl;

    /*
     * 根据以上定义，原一元二次方程变为以下形式
     * - N * a + sumX * b = sumY
     * - sumX * a + sumX2 * b =sumXY
     * 解出 a, b, 得到拟合直线的方程
     */
    double b = (sumXY * N - sumX * sumY) / (sumX2 * N - sumX * sumX);
    double a = (sumY / N) - (sumX * sumXY * N - sumX * sumX * sumY) / (sumX2 * N * N - sumX * sumX * N);
    cout << "\n拟合直线 y = " << b << "x " << a << endl;

    return 0;
}