c++ - 如何用 C++ 拟合二维散点数据

Question

我以前使用 MATLAB，对于我提出的问题，我可以使用 p = polyfit(x,y,1) 来估计板中散点数据的最佳拟合线。我想知道我可以依靠哪些资源来用 C++ 实现线拟合算法。我知道这个主题有很多算法，对我来说，我希望算法应该很快，同时它可以在 MATLAB 中获得与 polyfit 函数相当的精度。

Answer 1

这个页面描述的算法比维基百科更简单，没有额外的步骤来计算均值等:http://faculty.cs.niu.edu/~hutchins/csci230/best-fit.htm .几乎从那里引用，在 C++ 中是:

#include <vector>
#include <cmath>

struct Point {
  double _x, _y;
};
struct Line {
  double _slope, _yInt;
  double getYforX(double x) {
    return _slope*x + _yInt;
  }
  // Construct line from points
  bool fitPoints(const std::vector<Point> &pts) {
    int nPoints = pts.size();
    if( nPoints < 2 ) {
      // Fail: infinitely many lines passing through this single point
      return false;
    }
    double sumX=0, sumY=0, sumXY=0, sumX2=0;
    for(int i=0; i<nPoints; i++) {
      sumX += pts[i]._x;
      sumY += pts[i]._y;
      sumXY += pts[i]._x * pts[i]._y;
      sumX2 += pts[i]._x * pts[i]._x;
    }
    double xMean = sumX / nPoints;
    double yMean = sumY / nPoints;
    double denominator = sumX2 - sumX * xMean;
    // You can tune the eps (1e-7) below for your specific task
    if( std::fabs(denominator) < 1e-7 ) {
      // Fail: it seems a vertical line
      return false;
    }
    _slope = (sumXY - sumX * yMean) / denominator;
    _yInt = yMean - _slope * xMean;
    return true;
  }
};

请注意，如果对点的“最佳”描述是垂直线，则此算法和维基百科 (http://en.wikipedia.org/wiki/Simple_linear_regression#Fitting_the_regression_line) 中的算法都会失败。他们失败是因为他们使用了

y = k*x + b 

本质上不能描述垂直线的线方程。如果你还想涵盖数据点被垂直线“最佳”描述的情况，你需要一个使用的线拟合算法

A*x + B*y + C = 0

线方程。您仍然可以修改当前算法以生成该等式:

y = k*x + b <=>
y - k*x - b = 0 <=>
B=1, A=-k, C=-b

就上面的代码而言:

B=1, A=-_slope, C=-_yInt

然后在 if 的“then” block 中检查分母是否等于 0，而不是 //Fail: it seems a vertical line，生成以下直线方程:

x = xMean <=>
x - xMean = 0 <=>
A=1, B=0, C=-xMean

我刚刚注意到我所指的原始文章已被删除。这个网页提出了一些不同的线拟合公式:http://hotmath.com/hotmath_help/topics/line-of-best-fit.html

double denominator = sumX2 - 2 * sumX * xMean + nPoints * xMean * xMean;
...
_slope = (sumXY - sumY*xMean - sumX * yMean + nPoints * xMean * yMean) / denominator;

公式是相同的，因为 nPoints*xMean == sumX 和 nPoints*xMean*yMean == sumX * yMean == sumY * xMean。