--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/+utils/@math/getjacobian.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,85 @@
+% GETJACOBIAN Calculate Jacobian of a given model function.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% DESCRIPTION
+% 
+%     Calculate Jacobian of a given model function for the given set of
+%     coefficients. Jacobian is approximated with finite difference method.
+% 
+% CALL:
+% 
+%     J = getjacobian(coeff,model,X)
+% 
+% INPUT:
+% 
+%     J - fit coefficients
+%     model - model function
+%     X - x vactor (abscissa)
+%       
+% 
+% OUTPUT:
+% 
+%     J - Jacobian
+% 
+% Note: Look at nlinfit.m of the stats toolbox. Model should be a matlab
+% function calculating Model values as a function of X and of coefficients
+% NOTE: The function prefer to work with column objects. Therefore it is
+% good practise to directly input coeff and X as column objects
+%
+%   Examples:
+%
+%      Use @ to specify MODELFUN:
+%         load reaction;
+%         J = getjacobian(coeff,@mymodel,X);
+%
+%      where MYMODEL is a MATLAB function such as:
+%         function yhat = mymodel(beta, X)
+%         yhat = (beta(1)*X(:,2) - X(:,3)/beta(5)) ./ ...
+%                        (1+beta(2)*X(:,1)+beta(3)*X(:,2)+beta(4)*X(:,3));
+%      or
+% 
+%         mymodel = @(beta,X)(beta(1).*X(:,1)+beta(2).*X(:,2)+...)
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% VERSION: $Id: getjacobian.m,v 1.1 2009/03/26 16:26:05 luigi Exp $
+%
+% HISTORY:     26-03-2009 L Ferraioli
+%                 Creation
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function J = getjacobian(coeff,model,X)
+  % checking size, willing to work with columns
+  [a,b] = size(coeff);
+  if a<b
+    coeff = coeff.';
+  end
+  [a,b] = size(X);
+  if a<b
+    X = X.';
+  end
+  
+  % finite difference relative step
+  fdiffstep = eps^(1/3);
+  % evaluate model on the input coefficients
+  yfit = model(coeff,X);
+  % check for NaNs
+  nans = isnan(yfit(:));
+  
+  % initialization
+  p = numel(coeff);
+  delta = zeros(size(coeff));
+  J = zeros(numel(yfit),p);
+  for k = 1:p
+    if (coeff(k) == 0)
+      nb = sqrt(norm(coeff));
+      delta(k) = fdiffstep * (nb + (nb==0));
+    else
+      delta(k) = fdiffstep*coeff(k);
+    end
+    yplus = model(coeff+delta,X);
+    dy = yplus(:) - yfit(:);
+    dy(nans) = [];
+    J(:,k) = dy/delta(k);
+    delta(k) = 0;
+  end
+end
\ No newline at end of file