Ordinary least-squares fit

Fit a linear combination of given functions, ∑ c_kf_k(x)|_k=1..m , to a (measured) data set {x_i, y_i, δy_i}_i=1..n , where δy_i are the uncertainties of the measurement (error-bars).

You can pass around a set of functions {fi(x)}k=1..m in a variable of the type Func<double,double>[] (that is, as an array of functions). For example,

var fs = new Func<double,double>[] { z => 1.0, z => z, z => z*z };

1. (6 points) Ordinary least-squares fit by QR-decomposition

   • Make sure that your QR-decomposition routines work for tall matrices.

   • Implement a routine that makes a least-squares fit—using your QR-decomposition routines—of a given data-set,

   {xi, yi, δyi}i=1...n ,
   

   with a linear combination

   Fc(x) ≐ ∑k=1..m ck fk(x)
   

   of given functions f_k(x)|_k=1..m . The interface could be something like

   vector lsfit(Func<double,double>[] fs, vector x, vector y, vector dy) 

   The routine takes as arguments the data to fit, {x_i, y_i, δy_i}, and the set of functions, {f_k}, the linear combination of which should fit the data. The routine must calculate and return the vector of the best fit coefficients, {c_k}.

   • The law of radioactive decay. In 1902 [Rutherford and Soddy] measured the radioactivity of the (then not well explored) element, called ThX at the time, and obtained the following results,

   Time t (days)                     : 1,  2,  3, 4, 6, 9,   10,  13,  15
   Activity y of ThX (relative units): 117,100,88,72,53,29.5,25.2,15.2,11.1
   

   From this data they correctly deduced that radioactive decay follows exponentil law, y(t)=ae^-λt (equation (1) in the article).

   Now, assume that the incertainty δy of the measurement (which is actually nowhere to find in the article – not good at all) was determined by the last-but-one digit of the measurement,

   δy: 5,5,5,4,4,3,3,2,2 

   and fit the data with exponential function in the usual logarithmic way, ln(y)=ln(a)-λt. The uncertainty of the logarithm should be probably taken as δln(y)=δy/y.

   Plot the experimental data (with error-bars) and your best fit.

   From your fit find out the half-life time, T_½ = ^ln(2)/_λ, of ThX.

   ThX is today known as ²²⁴Ra – compare your result with the modern value.

2. (3 points) Uncertainties of the fitting coefficients
   • Modify you least-squares fitting function such that it also calculates the covariance matrix and the uncertainties of the fitting coefficients.

   vector lsfit(Func<double,double>[] fs, vector x, vector y, vector dy) 

   
   • Calculate the uncertainty of the half-life value for ThX from the given data.
   • Does it agree with the modern value within the estimated uncertainty?

3. (1 points) Evaluation of the quality of the uncertainties on the fit coefficients
   • Plot your best fit,
   ∑ _k=1..m c_k f_k(x),
   together with the fits where you change the fit coefficients by the estimated δc, that is,
   ∑ _k=1..m (c_k±δc_k) f_k(x).
   You can either add(subtract) δc_k to all coefficients at once, or individually.