Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/getjacobian.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 % GETJACOBIAN Calculate Jacobian of a given model function.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION
5 %
6 % Calculate Jacobian of a given model function for the given set of
7 % coefficients. Jacobian is approximated with finite difference method.
8 %
9 % CALL:
10 %
11 % J = getjacobian(coeff,model,X)
12 %
13 % INPUT:
14 %
15 % J - fit coefficients
16 % model - model function
17 % X - x vactor (abscissa)
18 %
19 %
20 % OUTPUT:
21 %
22 % J - Jacobian
23 %
24 % Note: Look at nlinfit.m of the stats toolbox. Model should be a matlab
25 % function calculating Model values as a function of X and of coefficients
26 % NOTE: The function prefer to work with column objects. Therefore it is
27 % good practise to directly input coeff and X as column objects
28 %
29 % Examples:
30 %
31 % Use @ to specify MODELFUN:
32 % load reaction;
33 % J = getjacobian(coeff,@mymodel,X);
34 %
35 % where MYMODEL is a MATLAB function such as:
36 % function yhat = mymodel(beta, X)
37 % yhat = (beta(1)*X(:,2) - X(:,3)/beta(5)) ./ ...
38 % (1+beta(2)*X(:,1)+beta(3)*X(:,2)+beta(4)*X(:,3));
39 % or
40 %
41 % mymodel = @(beta,X)(beta(1).*X(:,1)+beta(2).*X(:,2)+...)
42 %
43 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
44 % VERSION: $Id: getjacobian.m,v 1.1 2009/03/26 16:26:05 luigi Exp $
45 %
46 % HISTORY: 26-03-2009 L Ferraioli
47 % Creation
48 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
49
50 function J = getjacobian(coeff,model,X)
51 % checking size, willing to work with columns
52 [a,b] = size(coeff);
53 if a<b
54 coeff = coeff.';
55 end
56 [a,b] = size(X);
57 if a<b
58 X = X.';
59 end
60
61 % finite difference relative step
62 fdiffstep = eps^(1/3);
63 % evaluate model on the input coefficients
64 yfit = model(coeff,X);
65 % check for NaNs
66 nans = isnan(yfit(:));
67
68 % initialization
69 p = numel(coeff);
70 delta = zeros(size(coeff));
71 J = zeros(numel(yfit),p);
72 for k = 1:p
73 if (coeff(k) == 0)
74 nb = sqrt(norm(coeff));
75 delta(k) = fdiffstep * (nb + (nb==0));
76 else
77 delta(k) = fdiffstep*coeff(k);
78 end
79 yplus = model(coeff+delta,X);
80 dy = yplus(:) - yfit(:);
81 dy(nans) = [];
82 J(:,k) = dy/delta(k);
83 delta(k) = 0;
84 end
85 end