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