一.说明
本文主要讲解利用最小二乘法解决曲线拟合的问题,并利用C#代码实现。重点在代码实现,原理简单介绍下,本文所有解释都以二次曲线举例。
二.简单介绍最小二乘法原理
最小二乘法(又称最小平方法)是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据,并使得这些求得的数据与实际数据之间误差的平方和为最小。
在曲线拟合的应用场景中,相当于拟合的曲线后的方程利用给定x坐标计算出结果y1与给的y坐标的所有平方和最小即:∑i=1n[y(xi)−yi]2\sum_{i=1}^{n}[y(x_i)-y_i]^2∑i=1n[y(xi)−yi]2最小。
最终转化对系数a1 a2 a3求偏导,多次迭代后算出估算值。
偏导公式:∂Q/∂aj=2∑i=1n[y(xi)−yi]xj∂Q/∂a_j = 2\sum_{i=1}^{n}[y(x_i)-y_i]x^j∂Q/∂aj=2∑i=1n[y(xi)−yi]xj
1.对常数a3偏导:∂Q/∂a3=2∑i=1n[(a1xi2+a2xi+a3)−yi]=0∂Q/∂a_3 = 2\sum_{i=1}^{n}[(a_1x_i^2+a_2x_i+a_3)-y_i]=0∂Q/∂a3=2∑i=1n[(a1xi2+a2xi+a3)−yi]=0
2.对一次系数a2偏导:∂Q/∂a2=2∑i=1n[(a1xi2+a2xi+a3)−yi]x=0∂Q/∂a_2 = 2\sum_{i=1}^{n}[(a_1x_i^2+a_2x_i+a_3)-y_i]x=0∂Q/∂a2=2∑i=1n[(a1xi2+a2xi+a3)−yi]x=0
3.对二次系数a3偏导:∂Q/∂a1=2∑i=1n[(a1xi2+a2xi+a3)−yi]x2=0∂Q/∂a_1 = 2\sum_{i=1}^{n}[(a_1x_i^2+a_2x_i+a_3)-y_i]x^2=0∂Q/∂a1=2∑i=1n[(a1xi2+a2xi+a3)−yi]x2=0
由此推导出
a3=(∑i=1nyi−a1∑i=1nxi2−a2∑i=1nxi)/na3 = ( \sum_{i=1}^{n}y_i - a_1 \sum_{i=1}^{n}x_i^2-a_2 \sum_{i=1}^{n}x_i)/na3=(∑i=1nyi−a1∑i=1nxi2−a2∑i=1nxi)/n
a2=(∑i=1nyixi−a3∑i=1nxi−a1∑i=1nxi3)/∑i=1nxi2a2 = ( \sum_{i=1}^{n}y_ix_i - a_3\sum_{i=1}^{n}x_i-a_1 \sum_{i=1}^{n}x_i^3)/ \sum_{i=1}^{n}x_i^2a2=(∑i=1nyixi−a3∑i=1nxi−a1∑i=1nxi3)/∑i=1nxi2
a1=(∑i=1nyixi2−a3∑i=1nxi2−a2∑i=1nxi3)/∑i=1nxi4a1 = ( \sum_{i=1}^{n}y_ix_i^2 - a_3 \sum_{i=1}^{n}x_i^2-a_2 \sum_{i=1}^{n}x_i^3)/\sum_{i=1}^{n}x_i^4a1=(∑i=1nyixi2−a3∑i=1nxi2−a2∑i=1nxi3)/∑i=1nxi4
三.代码实现
代码中的a代表二次系数,b代表一次系数,c代表常数。
这里构造了x参数xarray[],y参数yarray[]进行演示,最终结果也与我们构造的参数基本一致。
这里的z1,z2,z3用来评估该次运算精度是否达到我们的预期值,值越小计算时间越久。
var xarray = new double[100];
var yarray = new double[100];
for (int i = 0; i < xarray.Length; i++)
{
xarray[i] = i;
yarray[i] = 1.356*i*i+7.969 * i - 0.26;
}
double a, b, c, m1, m2, m3, z1, z2, z3;
a = b = c = 0;
double sumx = 0, sumx2 = 0, sumx3 = 0, sumx4 = 0, sumy = 0, sumxy = 0, sumx2y = 0;
for (var i = 0; i < xarray.Length; i++)
{
sumx += xarray[i];
sumy += yarray[i];
sumx2 += Math.Pow(xarray[i], 2);
sumxy += xarray[i] * yarray[i];
sumx3 += Math.Pow(xarray[i], 3);
sumx2y += Math.Pow(xarray[i], 2) * yarray[i];
sumx4 += Math.Pow(xarray[i], 4);
}
do
{
m1 = a;
a = (sumx2y - sumx3 * b - sumx2 * c) / sumx4;
z1 = (a - m1) * (a - m1);
m2 = b;
b = (sumxy - sumx * c - sumx3 * a) / sumx2;
z2 = (b - m2) * (b - m2);
m3 = c;
c = (sumy - sumx2 * a - sumx * b) / xarray.Length;
z3 = (c - m3) * (c - m3);
} while (z1 > 1e-13 || z2 > 1e-13 || z3 > 1e-13);
Debug.WriteLine($"a:{a} b;{b} c:{c}");
}
四.其他方法
其他次数曲线同理,除了利用最小二乘法,还可以利用高斯牛顿方程求解,下一节介绍利用高斯牛顿方程及场景应用,各位大佬有其他更好的方法也可告知。