var_prop_classical

% variance propagation by numerical differentiation

 Usage: 
   [my,Syy,J]=var_prop_classical(@f,mx,Sxx,p)%

 Input:
    f   = nonlinear function y=f(x,p) as handle of m-file
    mx  = mean of x
    Sxx = covaraince of x
    p   = parameters for function

 Output 
    my  = mean of y
    Syy = covariance of y
    J = Jacobian matrix

 condition: Sxx must not have zero diagonals, otherwise numerically unstable

0001 %% variance propagation by numerical differentiation
0002 %
0003 % Usage:
0004 %   [my,Syy,J]=var_prop_classical(@f,mx,Sxx,p)%
0005 %
0006 % Input:
0007 %    f   = nonlinear function y=f(x,p) as handle of m-file
0008 %    mx  = mean of x
0009 %    Sxx = covaraince of x
0010 %    p   = parameters for function
0011 %
0012 % Output
0013 %    my  = mean of y
0014 %    Syy = covariance of y
0015 %    J = Jacobian matrix
0016 %
0017 % condition: Sxx must not have zero diagonals, otherwise numerically unstable
0018 %
0019 function [my,Syy,J] = var_prop_classical(f,mx,Sxx,p)
0020 
0021 % determine
0022 % mean of output
0023 % dimension of output (must be coded in function)
0024 if nargin == 3
0025     my = f(mx);
0026 else
0027     my = f(mx,p);
0028 end
0029 nf = size(my,1);
0030 
0031 % determine dimension of input
0032 n = size(mx,1);
0033 
0034 % factor for deviations from mean for Jacobian
0035 k = 0.1;
0036  
0037 % Jacobian
0038 J = zeros(nf,n);
0039 for ii = 1:n
0040     % finite difference of input
0041     s = sqrt(Sxx(ii,ii));
0042     % vector of effect: in direction if inout
0043     eff = zeros(n,1);
0044     eff(ii) = k;
0045     % symmetric finite difference
0046     if nargin == 3
0047         t = (f(mx+s*eff)-f(mx-s*eff))/(2*k*s+eps);
0048     else
0049         t = (f(mx+s*eff,p)-f(mx-s*eff,p))/(2*k*s+eps);
0050     end
0051     % Jacobian element
0052     J(:,ii) = t(:);
0053 end;
0054 
0055 % Variance
0056 Syy = J*Sxx*J';