% JR2COV Calculates coefficients covariance matrix from Jacobian and Residuals.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION
% 
%     Calculates coefficients covariance matrix from Jacobian and
%     Residuals. The algorithm uses the QR factorization of J to perform
%     the calculation inv(J'J)*s wehre J is the Jacobian matrix and s is
%     the mean squared error.
% 
% CALL:
% 
%     covmat = jr2cov(J,resid)
% 
% INPUT:
% 
%     J - Jacobian of the function with respect to the given coefficients
%     resid - Fit residuals
%       
% 
% OUTPUT:
% 
%     covmat - covariance matrix of the fit coefficients
% 
% Note: Resid should be a column vector, Number of rows of J must be equal
% to the number of rows of Resid. The number of columns of J defines the
% number of corresponding fit parameters for which we want to calculate the
% covariance matrix
% Estimate covariance from J and residuals. Look at the nlparci.m function
% of the stats toolbox for further details
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% VERSION: $Id: jr2cov.m,v 1.2 2009/03/27 14:48:54 luigi Exp $
%
% HISTORY:     26-03-2009 L Ferraioli
%                 Creation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function covmat = jr2cov(J,resid)
  
  missing = isnan(resid);
  if ~isempty(missing)
      resid(missing) = [];
  end
  
  n = length(resid);
  
  J(missing,:) = [];
  
  if size(J,1)~=n
    error('The number of rows of J does not match the number of rows of RESID.');
  end
  [n,p] = size(J);
  v = n-p; % degrees of freedom of the parameters estimation
  
  % Approximation when a column is zero vector
  temp = find(max(abs(J)) == 0);
  if ~isempty(temp)
    J(temp,:) = J(temp,:) + sqrt(eps(class(J)));
  end
  
  % Calculate covariance matrix
  [Q,R] = qr(J,0);
  Rinv = R\eye(size(R));
  
  mse = norm(resid)^2 / v; % mean square error
  covmat = mse * Rinv*Rinv';
  
end