### 39.12 Computing the covariance matrix of best fit parameters

Function: int gsl_multifit_fdfsolver_jac (gsl_multifit_fdfsolver * s, gsl_matrix * J)

This function stores the n-by-p Jacobian matrix for the current iteration of the solver s into the output J.

Function: int gsl_multifit_covar (const gsl_matrix * J, const double epsrel, gsl_matrix * covar)

This function computes the covariance matrix of best-fit parameters using the Jacobian matrix J and stores it in covar. The parameter epsrel is used to remove linear-dependent columns when J is rank deficient.

The covariance matrix is given by,

covar = (J^T J)^{-1}


or in the weighted case,

covar = (J^T W J)^{-1}


and is computed by QR decomposition of J with column-pivoting. Any columns of R which satisfy

|R_{kk}| <= epsrel |R_{11}|


are considered linearly-dependent and are excluded from the covariance matrix (the corresponding rows and columns of the covariance matrix are set to zero).

If the minimisation uses the weighted least-squares function f_i = (Y(x, t_i) - y_i) / \sigma_i then the covariance matrix above gives the statistical error on the best-fit parameters resulting from the Gaussian errors \sigma_i on the underlying data y_i. This can be verified from the relation \delta f = J \delta c and the fact that the fluctuations in f from the data y_i are normalised by \sigma_i and so satisfy <\delta f \delta f^T> = I.

For an unweighted least-squares function f_i = (Y(x, t_i) - y_i) the covariance matrix above should be multiplied by the variance of the residuals about the best-fit \sigma^2 = \sum (y_i - Y(x,t_i))^2 / (n-p) to give the variance-covariance matrix \sigma^2 C. This estimates the statistical error on the best-fit parameters from the scatter of the underlying data.