Chris@49: // Tutorial for linear least square fitting
Chris@49: // Author: Pierre Moulon
Chris@49: // Date: 8 December 2009
Chris@49: // Objective:
Chris@49: // Fit a 2D line with to a set of points
Chris@49: // Specifically, find a and b in the model y = ax + b
Chris@49: // 
Chris@49: // Direct application of the example in:
Chris@49: // http://en.wikipedia.org/wiki/Linear_least_squares#Motivational_example
Chris@49: 
Chris@49: 
Chris@49: #include <iostream>
Chris@49: 
Chris@49: #include "armadillo"
Chris@49: 
Chris@49: using namespace arma;
Chris@49: using namespace std;
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: int main(int argc, char** argv)
Chris@49:   {
Chris@49:   // points to which we will fit the line
Chris@49:   mat data = "1 6; 2 5; 3 7; 4 10";
Chris@49: 
Chris@49:   cout << "Points used for the estimation:" << endl;
Chris@49:   cout << data << endl;
Chris@49: 
Chris@49:   // Build matrices to solve Ax = b problem:
Chris@49:   vec b(data.n_rows);
Chris@49:   mat C(data.n_rows, 2);
Chris@49: 
Chris@49:   for(u32 i=0; i<data.n_rows; ++i)
Chris@49:     {
Chris@49:     b(i)   = data(i,1);
Chris@49:     
Chris@49:     C(i,0) = 1;
Chris@49:     C(i,1) = data(i,0);
Chris@49:     }
Chris@49:   
Chris@49:   cout << "b:" << endl;
Chris@49:   cout << b << endl;
Chris@49:   
Chris@49:   cout << "Constraint matrix:" << endl;
Chris@49:   cout << C << endl;
Chris@49:   
Chris@49:   // Compute least-squares solution:
Chris@49:   vec solution = solve(C,b);
Chris@49:   
Chris@49:   // solution should be "3.5; 1.4"
Chris@49:   cout << "solution:" << endl;
Chris@49:   cout << solution << endl;
Chris@49:   
Chris@49: 
Chris@49:   cout << "Reprojection error:" << endl;
Chris@49:   
Chris@49:   for(u32 i=0; i<data.n_rows; ++i)
Chris@49:     {
Chris@49:     cout << "  residual: " << ( data(i,1) - (solution(0) + solution(1) * data(i,0)) ) << endl;
Chris@49:     }
Chris@49:     
Chris@49:   return 0;
Chris@49:   }
Chris@49: