对于给定的数据点(xi, yi)(i = 1, 2, …, N)求作一次式 y = a + b x y=a+bx y=a+bx,使总误差
Q = ∑ i = 1 N [ y i − ( a + b x i ) ] 2 Q=\sum_{i=1}^{N} \left [y_{i}-(a+bx_{i})\right ] ^{2} Q=i=1∑N[yi−(a+bxi)]2
为最小。
这个问题是不难求解的。由微积分中求极值的方法知,使 Q 达到极值的参数 a,b 应满足
∂ Q ∂ a = 0 , ∂ Q ∂ b = 0 \frac{\partial Q}{\partial a} =0,\frac{\partial Q}{\partial b} =0 ∂a∂Q=0,∂b∂Q=0
即成立
{ a N + b ∑ x i = ∑ y i a ∑ x i + b ∑ x i 2 = ∑ x i y i (1) \begin{cases} aN+b {\textstyle \sum x_{i}} =\sum y_{i}\\ a\sum x_{i}+b\sum x_{i}^{2}=\sum x_{i}y_{i} \end{cases}\tag1 {aN+b∑xi=∑yia∑xi+b∑xi2=∑xiyi(1)
式中 ∑ \sum ∑ 表示关于下表 i 从 1 到 N 的求和。
例:炼钢是个氧化脱碳的过程,钢液含碳量的多少直接影响冶炼时间的长短。下面是某平炉的生产记录,表中 i 为试验次数,xi 为全部炉料熔化完毕时钢液的含碳量,yi 为熔毕至出钢的冶炼时间(以分为单位)。
i | 1 | 2 | 3 | 4 | 5 |
---|---|---|---|---|---|
xi | 165 | 123 | 150 | 123 | 141 |
yi | 187 | 126 | 172 | 125 | 148 |
把表中所给数据画在坐标纸上,将会看到,数据点的分布可以用一条直线来近似描述。设所求的拟合直线为 y = a + b x y=a+bx y=a+bx,则方程组(1)的具体形式是
{ 5 a + 702 b = 758 702 a + 99864 b = 108396 \begin{cases} 5a+702b=758\\ 702a+99864b=108396 \end{cases} {5a+702b=758702a+99864b=108396
解出 a,b,即得拟合直线 y = − 60.0392 + 1.5138 x y=-60.0392+1.5138x 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;
}