+
−% 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