Home > General-Functions > Statistics > var_prop_classical.m

var_prop_classical

PURPOSE ^

% variance propagation by numerical differentiation

SYNOPSIS ^

function [my,Syy,J] = var_prop_classical(f,mx,Sxx,p)

DESCRIPTION ^

% 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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';

Generated on Mon 08-Jan-2018 17:21:49 by m2html © 2005