Ordinary least-squares fit

1. (6 points)

   Implement a routine which makes a least-squares fit of a given data-set, {x_i, y_i, δy_i | i=1...n}, with a linear combination of given functions, F(x) =∑ _k=1..m c_k f_k(x).

   The subroutine must calculate the vector of the coefficients, c_k, and the covariance matrix. The fitting functions {f_k(x)} can be programmed, e.g., as

   • a vector of functions (C++11,JavaScript,Python),
     std::vector<std::function<double(double)>> fitfunctions = {f1,...,fm};

   • a function of two parameters (C,Fortran,C++,Python),
     double fitfunctions(int i, double x){return fi(x);}
     For example:

     #include<math.h> // NAN is here
     double funs(int i, double x){
        switch(i){
        case 0: return 1.0; break;
        case 1: return x;   break;
        case 2: return x*x; break;
        default: {fprintf(stderr,"funs: wrong i:%d",i); return NAN;}
        }
     }
     

   • In the worst case you can simply precalculate the matrix A_ik=f_k(x_i) and pass it as a parameter.
   Make up some interesting fits and check that the reported errors of the fitting coefficients are reasonable. For example,

   

2. (3 points) Linear least-squares solution using thin singular-value decomposition.

   Here under the thin singular value decomposition we shall understand a representation of a tall n×m (n>m) matrix A in the form A=UΣV^T where U is an orthogonal n×m matrix (U^TU=1), Σ is a square m×m diagonal matrix with non-negative real numbers on the diagonal (called singular values of matrix A), and V is a square m×m orthoginal matrix (V^TV=1).

   Singular value decomposition can be found by diagonalising the A^TA matrix, A^TA=VDV^T, where D is a diagonal matrix with eigenvalues of A^TA on the diagonal and V is the matrix of the corresponding eigenvectors. Indeed it is easy to check (do it) that the sought decomposition can the be constructed as A=UΣV^T where Σ=D^1/2, U=AVD^-1/2.

   Singular value decomposition can be used to solve the linear least squares problem Ax=b. Indeed inserting the decomposition into the equation gives UΣV^Tx=b. Multiplying from the left with U^T and using the orthogonality of U one gets the projected equation ΣV^Tx=U^Tb. This is a square system which can be easily solved first by solving the diagonal system Σy=U^Tb for y and then obtaining x as x=Vy.

   1. Implement a function that makes a singular value decomposition of a (tall) matrix A by diagonalising A^TA (with your implementation of the Jacobi eigenvalue algorithm), A^TA=VDV^T (where D is a diagonal matrix with eigenvalues of A^TA and V is the matrix of the corresponding eigenvectors) and then using A=UΣV where U=AVD^-1/2 and Σ=D^1/2.
   2. Implement a function that solves the linear least squares problem Ax=b using you implementation of singular-value decomposition.
   3. Implement a function that makes ordinary least-squares fit using your implementation of singular-value decomposition.

3. (1 point)

   Make up two experimental points with error bars, {x_i, y_i ± Δ y_i | i=1,2}. Assume that the theoretical curve is a straight line. Make an ordinary least-squares fit of your experimental data with a linear function y(x)=a+bx. Make a plot of your results. On the plot illustrate the probable distribution of the variable y as function of x (in the assumption of normal distribution of y_i, a, b if you need that).

   Hint: the variation of y with respect to variations of a and b as function of x is given as
   δy = δa + δb x.
   Therefore
   δy² = δa² + δb² x² +2δa δb x.
   For normally distributed a and b this translates into
   Δy²(x) = Δa² + Δb²x² + 2<δa δb>x
   where <δa δb> is the covariance.

4. (0 points)

   Make up an interesting exercise here.