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