How to fit the 2D scatter data with a line with C++

Question

I used to work with MATLAB, and for the question I raised I can use p = polyfit(x,y,1) to estimate the best fit line for the scatter data in a plate. I was wondering which resources I can rely on to implement the line fitting algorithm with C++. I understand there are a lot of algorithms for this subject, and for me I expect the algorithm should be fast and meantime it can obtain the comparable accuracy of polyfit function in MATLAB.

Solution 1

I would suggest coding it from scratch. It is a very simple implementation in C++. You can code up both the intercept and gradient for least-squares fit (the same method as polyfit) from your data directly from the formulas here

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

These are closed form formulas that you can easily evaluate yourself using loops. If you were using higher degree fits then I would suggest a matrix library or more sophisticated algorithms but for simple linear regression as you describe above this is all you need. Matrices and linear algebra routines would be overkill for such a problem (in my opinion).

Solution 2

This page describes the algorithm easier than Wikipedia, without extra steps to calculate the means etc. : http://faculty.cs.niu.edu/~hutchins/csci230/best-fit.htm . Almost quoted from there, in C++ it's:

#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;
  }
};

Please, be aware that both this algorithm and the algorithm from Wikipedia ( http://en.wikipedia.org/wiki/Simple_linear_regression#Fitting_the_regression_line ) fail in case the "best" description of points is a vertical line. They fail because they use

y = k*x + b 

line equation which intrinsically is not capable to describe vertical lines. If you want to cover also the cases when data points are "best" described by vertical lines, you need a line fitting algorithm which uses

A*x + B*y + C = 0

line equation. You can still modify the current algorithm to produce that equation:

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

In terms of the above code:

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

And in "then" block of the if checking for denominator equal to 0, instead of // Fail: it seems a vertical line, produce the following line equation:

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

I've just noticed that the original article I was referring to has been deleted. And this web page proposes a little different formula for line fitting: 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;

The formulas are identical because nPoints*xMean == sumX and nPoints*xMean*yMean == sumX * yMean == sumY * xMean.

Solution 3

Equation of line is Ax + By + C=0.

So it can be easily( when B is not so close to zero ) convert to 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;
}

Solution 4

You can also use or go over this implementation there is also documentation here.

Solution 5

Fitting a Line can be acomplished in different ways. Least Square means minimizing the sum of the squared distance. But you could take another cost function as example the (not squared) distance. But normaly you use the squred distance (Least Square). There is also a possibility to define the distance in different ways. Normaly you just use the "y"-axis for the distance. But you could also use the total/orthogonal distance. There the distance is calculated in x- and y-direction. This can be a better fit if you have also errors in x direction (let it be the time of measurment) and you didn't start the measurment on the exact time you saved in the data. For Least Square and Total Least Square Line fit exist algorithms in closed form. So if you fitted with one of those you will get the line with the minimal sum of the squared distance to the datapoints. You can't fit a better line in the sence of your defenition. You could just change the definition as examples taking another cost function or defining distance in another way.

There is a lot of stuff about fitting models into data you could think of, but normaly they all use the "Least Square Line Fit" and you should be fine most times. But if you have a special case it can be necessary to think about what your doing. Taking Least Square done in maybe a few minutes. Thinking about what Method fits you best to the problem envolves understanding the math, which can take indefinit time :-).

View more solutions