最⼩⼆乘法拟合直线
最⼩⼆乘法拟合直线
在科学实验和⽣产实践中，经常需要从⼀组实验数据出发寻求函数y=f(x)的⼀个近似表达式，也称为经验公式。从⼏何上看，就是希望根据给定的m个点，求曲线y=f(x)的⼀条近似曲线。因此这是个曲线拟合问题。
当我们要求近似曲线严格通过给定的每个点时，这是插值算法。对于本⽂所述的直线拟合来说，如果⽤插值算法，则只需要两个点就够了。实际直线拟合数据可能满⾜不了这个条件，为了便于计算，分析与应⽤，我们较多地根据“使测量点到直线距离的平⽅和最⼩”的原则来拟合。按最⼩⼆乘原则选择拟合曲线的⽅法，称为最⼩⼆乘法（Method of Least Squares）。
利⽤最⼩⼆乘法拟合曲线的⼀般步骤是：
将实验数据显⽰出来，分析曲线的形式；
确定拟合曲线的形式。对于本⽂来说，曲线形式是直线，y=a+bx;
建⽴法⽅程组并对其进⾏求解；
因为OpenCASCADE中数据结构及算法丰富，所以⽤OpenCASCADE可以快速实现直线的最⼩⼆乘法拟合算法。下⾯给出具体实现代码：
#include <iostream>
#include <gp_Lin2d.hxx>
#include <gp_Pnt2d.hxx>
#include <TColgp_Array1OfPnt2d.hxx>
#include <math_Vector.hxx>
#include <math_SVD.hxx>
#include <math_Gauss.hxx>
#include <math_GaussLeastSquare.hxx>
#include <OSD_Chronometer.hxx>
void fitLine(const TColgp_Array1OfPnt2d& thePoints,
const std::string& theFileName,
gp_Lin2d& theLine)
{
math_Vector aB(1, 2, 0.0);
math_Vector aX(1, 2, 0.0);
math_Matrix aM(1, 2, 1, 2);
Standard_Real aSxi = 0.0;
Standard_Real aSyi = 0.0;
Standard_Real aSxx = 0.0;
Standard_Real aSxy = 0.0;
std::ofstream aDrawFile(theFileName);
for (Standard_Integer i = thePoints.Lower(); i <= thePoints.Upper(); ++i)
{
missingconst gp_Pnt2d& aPoint = thePoints.Value(i);
aSxi += aPoint.X();
aSyi += aPoint.Y();
aSxx += aPoint.X() * aPoint.X();
aSxy += aPoint.X() * aPoint.Y();
aDrawFile << "vpoint p" << i << "" <<
aPoint.X() << "" << aPoint.Y() << " 0" << std::endl;
}
aM(1, 1) = thePoints.Size();
aM(1, 2) = aSxi;
发放英语aM(2, 1) = aSxi;
aM(2, 2) = aSxx;
aB(1) = aSyi;
aB(2) = aSxy;
OSD_Chronometer aChronometer;
wrongaChronometer.Start();
math_Gauss aSolver(aM);
//math_GaussLeastSquare aSolver(aM);
大学英语四级试题//math_SVD aSolver(aM);
aSolver.Solve(aB, aX);
记录英文if (aSolver.IsDone())
{
domestic是什么意思
Standard_Real aA = aX(1);
Standard_Real aB = aX(2);
gp_Pnt2d aP1(0.0, aA);
gp_Pnt2d aP2(-aA/aB, 0.0);
汗流浃背的读音theLine.SetLocation(aP1);
theLine.SetDirection(gp_Vec2d(aP1, aP2).XY());
aDrawFile << "vaxis l "
<< aP1.X() << "" << aP1.Y() << " 0 "
<< aP2.X() << "" << aP2.Y() << " 0 " << std::endl;
std::cout << "===================" << std::endl;
aX.Dump(std::cout);
}
aChronometer.Stop();
辅音字母大学英语四级查询aChronometer.Show();
}
int main()
{
gp_Lin2d aLine;
总统日// Test data 1
TColgp_Array1OfPnt2d aPoints1(1, 6);
aPoints1.SetValue(1, gp_Pnt2d(36.9, 181.0));
aPoints1.SetValue(2, gp_Pnt2d(46.7, 197.0));
aPoints1.SetValue(3, gp_Pnt2d(63.7, 235.0));
aPoints1.SetValue(4, gp_Pnt2d(77.8, 270.0));
aPoints1.SetValue(5, gp_Pnt2d(84.0, 283.0));
aPoints1.SetValue(6, gp_Pnt2d(87.5, 292.0));
fitLine(aPoints1, "l", aLine);
// Test data 2
TColgp_Array1OfPnt2d aPoints2(0, 7);
aPoints2.SetValue(0, gp_Pnt2d(0.0, 27.0));
aPoints2.SetValue(1, gp_Pnt2d(1.0, 26.8));
aPoints2.SetValue(2, gp_Pnt2d(2.0, 26.5));
aPoints2.SetValue(3, gp_Pnt2d(3.0, 26.3));
aPoints2.SetValue(4, gp_Pnt2d(4.0, 26.1));
aPoints2.SetValue(5, gp_Pnt2d(5.0, 25.7));
aPoints2.SetValue(6, gp_Pnt2d(6.0, 25.3));
aPoints2.SetValue(7, gp_Pnt2d(7.0, 24.8));
fitLine(aPoints2, "l", aLine);
return0;
}
在函数fitLine()中，根据拟合点建⽴法⽅程组，并使⽤math_Gauss来对法⽅程组进⾏求解。其实也可以使⽤math_GaussLeastSquare或者math_SVD等求解法⽅程组。在主函数main()中测试了两组数据。测试数据1来⾃易⼤义等《计算⽅法》，测试数据2来⾃《⾼等数学》。程序运⾏结果如下图所⽰：
与书中计算结果吻合。
由于需要将计算结果显⽰出来，所以在fitLine()函数中增加了输出Draw脚本⽂件的代码，实际运⽤时可将这部分代码去掉。将程序⽣成的脚本⽂件加载到Draw中，即可得到下⾯两个图：
测试数据1拟合直线
测试数据2拟合直线
将计算结果导出Draw脚本可视化，可以⽅便直观地查看拟合结果。如果熟悉其他脚本库如Python的matplotlib，也可以类似处理来将结果可视化。
