Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/jr2cov.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 % JR2COV Calculates coefficients covariance matrix from Jacobian and Residuals.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION
5 %
6 % Calculates coefficients covariance matrix from Jacobian and
7 % Residuals. The algorithm uses the QR factorization of J to perform
8 % the calculation inv(J'J)*s wehre J is the Jacobian matrix and s is
9 % the mean squared error.
10 %
11 % CALL:
12 %
13 % covmat = jr2cov(J,resid)
14 %
15 % INPUT:
16 %
17 % J - Jacobian of the function with respect to the given coefficients
18 % resid - Fit residuals
19 %
20 %
21 % OUTPUT:
22 %
23 % covmat - covariance matrix of the fit coefficients
24 %
25 % Note: Resid should be a column vector, Number of rows of J must be equal
26 % to the number of rows of Resid. The number of columns of J defines the
27 % number of corresponding fit parameters for which we want to calculate the
28 % covariance matrix
29 % Estimate covariance from J and residuals. Look at the nlparci.m function
30 % of the stats toolbox for further details
31 %
32 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
33 % VERSION: $Id: jr2cov.m,v 1.2 2009/03/27 14:48:54 luigi Exp $
34 %
35 % HISTORY: 26-03-2009 L Ferraioli
36 % Creation
37 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
38
39 function covmat = jr2cov(J,resid)
40
41 missing = isnan(resid);
42 if ~isempty(missing)
43 resid(missing) = [];
44 end
45
46 n = length(resid);
47
48 J(missing,:) = [];
49
50 if size(J,1)~=n
51 error('The number of rows of J does not match the number of rows of RESID.');
52 end
53 [n,p] = size(J);
54 v = n-p; % degrees of freedom of the parameters estimation
55
56 % Approximation when a column is zero vector
57 temp = find(max(abs(J)) == 0);
58 if ~isempty(temp)
59 J(temp,:) = J(temp,:) + sqrt(eps(class(J)));
60 end
61
62 % Calculate covariance matrix
63 [Q,R] = qr(J,0);
64 Rinv = R\eye(size(R));
65
66 mse = norm(resid)^2 / v; % mean square error
67 covmat = mse * Rinv*Rinv';
68
69 end