c++ - 如何用 C++ 的线条拟合 2D 散点图数据

Question

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

score 11 · Accepted Answer

此页面描述的算法比维基百科更容易，无需额外的步骤来计算平均值等：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;
  }
};

请注意，如果点的“最佳”描述是垂直线，则此算法和来自 Wikipedia ( 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检查分母是否等于 0 的“then”块中，而不是// 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。

score 7 · Accepted Answer

我建议从头开始编码。这是一个非常简单的 C++ 实现。polyfit您可以直接从此处的公式中从数据中编写最小二乘拟合的截距和梯度（与相同的方法）

http://en.wikipedia.org/wiki/Simple_linear_regression#Fitting_the_regression_line

这些是封闭形式的公式，您可以使用循环轻松评估自己。如果您使用更高程度的拟合，那么我建议您使用矩阵库或更复杂的算法，但是对于您上面描述的简单线性回归，这就是您所需要的。对于这样的问题（在我看来），矩阵和线性代数例程将是多余的。

score 5 · Accepted Answer

直线方程为 A x + B y + C=0。

所以它可以很容易（当 B 不是那么接近零时）转换为 y = (-A/B)*x + (-C/B)

typedef double scalar_type;
typedef std::array< scalar_type, 2 > point_type;
typedef std::vector< point_type > cloud_type;

bool fit( scalar_type & A, scalar_type & B, scalar_type & C, cloud_type const& cloud )
{
    if( cloud.size() < 2 ){ return false; }

    scalar_type X=0, Y=0, XY=0, X2=0, Y2=0;

    for( auto const& point: cloud )
    { // Do all calculation symmetric regarding X and Y
        X  += point[0];
        Y  += point[1];
        XY += point[0] * point[1];
        X2 += point[0] * point[0];
        Y2 += point[1] * point[1];
    }

    X  /= cloud.size();
    Y  /= cloud.size();
    XY /= cloud.size();
    X2 /= cloud.size();
    Y2 /= cloud.size();

    A = - ( XY - X * Y ); //!< Common for both solution

    scalar_type Bx = X2 - X * X;
    scalar_type By = Y2 - Y * Y;

    if( fabs( Bx ) < fabs( By ) ) //!< Test verticality/horizontality
    { // Line is more Vertical.
        B = By;
        std::swap(A,B);
    }
    else
    {   // Line is more Horizontal.
        // Classical solution, when we expect more horizontal-like line
        B = Bx;
    }
    C = - ( A * X + B * Y );

    //Optional normalization:
    // scalar_type  D = sqrt( A*A + B*B );
    // A /= D;
    // B /= D;
    // C /= D;
    return true;
}

score 1 · Accepted Answer

1

您还可以使用或查看此实现，这里也有文档。

于 2012-07-12T10:31:14.830 回答

score 1 · Accepted Answer

拟合线可以通过不同的方式完成。最小二乘意味着最小化距离的平方和。但是您可以将另一个成本函数作为（非平方）距离的示例。但通常你使用平方距离（最小二乘）。也有可能以不同的方式定义距离。通常，您只需使用“y”轴作为距离。但您也可以使用总/正交距离。那里的距离是在 x 和 y 方向计算的。如果您在 x 方向上也有错误（假设它是测量时间）并且您没有在保存在数据中的确切时间开始测量，这可能会更合适。对于最小二乘法和总最小二乘法线拟合存在封闭形式的算法。因此，如果您安装其中之一，您将获得与数据点的距离平方和最小的线。就你的辩护而言，你再合适不过了。您可以将定义更改为采用另一个成本函数或以另一种方式定义距离的示例。

有很多关于将模型拟合到您能想到的数据中的东西，但通常它们都使用“最小二乘法拟合”，而且您在大多数情况下应该没问题。但是，如果您有特殊情况，则可能有必要考虑您在做什么。可能在几分钟内完成最小二乘。考虑哪种方法最适合您解决问题涉及理解数学，这可能需要无限时间:-)。

score 1 · Accepted Answer

注意：这个答案不是这个问题的答案，而是这个“最接近一组点的线”被标记为这个“重复”的“线”（我认为不正确），无法添加新的答案.

问题要求：

找到与所有点的距离最小的线？我所说的距离是指点和线之间的最短距离。

“点和线之间”距离的最常见解释是欧几里德距离，“从所有点”最常见的解释是距离之和（绝对值或平方值）。

当目标是最小化欧式距离平方和时，线性回归 (LST) 不是要使用的算法。此外，线性回归不能产生垂直线。要使用的算法是“总最小二乘法”。有关问题描述的示例维基百科和数学堆栈交换中的此答案有关公式的详细信息，请参见示例。

score 0 · Accepted Answer

要适合一条线，y=param[0]x+param[1]只需执行以下操作：

// loop over data:
{               
sum_x += x[i];
sum_y += y[i];
sum_xy += x[i] * y[i];
sum_x2 += x[i] * x[i];
}

// means
double mean_x = sum_x / ninliers;
double mean_y = sum_y / ninliers;

float varx = sum_x2 - sum_x * mean_x;
float cov = sum_xy - sum_x * mean_y;

// 检查零 varx

param[0] = cov / varx;
param[1] = mean_y - param[0] * mean_x;

有关该主题的更多信息http://easycalculation.com/statistics/learn-regression.php （公式相同，它们只是乘以并除以 N，样本 sz。）。如果您想将平面拟合到 3D 数据，请使用类似的方法 - http://www.mymathforum.com/viewtopic.php?f=13&t=8793

免责声明：从某种意义上说，所有二次拟合都是线性的并且是最优的，因为它们可以减少参数中的噪声。但是，您可能对减少数据中的噪声感兴趣。您可能还想忽略异常值，因为它们会极大地影响您的解决方案。这两个问题都可以通过 RANSAC 解决。请参阅我的帖子：

c++ - 如何用 C++ 的线条拟合 2D 散点图数据

7 回答 7

Related

Reference