GaussMarkovModelLinear

% Gauss-Markov model linear, individual observations

 lv       = Nx1 vector of observations
 Cov_ll   = NxN covariance matrix of observations (Cov_ll^a)
 Am       = NxU Jacobian
 av       = Nx1 additive vector
 rU       = range of unknown parameters of interest

 est_x    = Ux1 estimated parameter
 Cov_xx   = UxU theoretical covariance matrix of estimated parameters
 sigma_0q = estimated variance factor (=1 if R=1)
 R        = redundancy
 vv       = Nx1 vector of estimated corrections
 zv       = Nx1 vector of standardized residuals (using sigma_0=1)
 riv      = Nx1 vector of redundacy numbers
 mu       = Nx1 vector of senstivity factors

 [est_x,Cov_xx,sigma_0q,R,vv,zv,riv,muv]=...
            GaussMarkovModelLinear(l,Cov_ll,A,a)

 Wolfgang Förstner 2015-06-04
 wfoerstn@uni-bonn.de

0001 %% Gauss-Markov model linear, individual observations
0002 %
0003 % lv       = Nx1 vector of observations
0004 % Cov_ll   = NxN covariance matrix of observations (Cov_ll^a)
0005 % Am       = NxU Jacobian
0006 % av       = Nx1 additive vector
0007 % rU       = range of unknown parameters of interest
0008 %
0009 % est_x    = Ux1 estimated parameter
0010 % Cov_xx   = UxU theoretical covariance matrix of estimated parameters
0011 % sigma_0q = estimated variance factor (=1 if R=1)
0012 % R        = redundancy
0013 % vv       = Nx1 vector of estimated corrections
0014 % zv       = Nx1 vector of standardized residuals (using sigma_0=1)
0015 % riv      = Nx1 vector of redundacy numbers
0016 % mu       = Nx1 vector of senstivity factors
0017 %
0018 % [est_x,Cov_xx,sigma_0q,R,vv,zv,riv,muv]=...
0019 %            GaussMarkovModelLinear(l,Cov_ll,A,a)
0020 %
0021 % Wolfgang Förstner 2015-06-04
0022 % wfoerstn@uni-bonn.de
0023 
0024 
0025 function [est_x,Cov_xx,sigma_0q,R,vv,zv,riv,nabla_lv,muv,muv1,uv1q,uv2]=...
0026     GaussMarkovModelLinear(lv,Cov_ll,Am,av,rU)
0027 
0028 %% initialization
0029 
0030 % number N and U
0031 [N,U] = size(Am);
0032 
0033 % redundancy
0034 R = N-U;
0035 if R < 0
0036     disp('not enough observations')
0037     return;
0038 end
0039 
0040 %% estimation
0041 W_ll   = inv(Cov_ll);                  % weight matrix of lv
0042 Bm     = W_ll*Am;                         %#ok<*MINV> % ancillary matrix
0043 Nm     = Bm'*Am;                       % normal equation matrix
0044 mv     = Bm'*(lv-av);                  % right hand sides
0045 Cov_xx =inv(Nm);                       % covariance matrix of est. parameters
0046 est_x  = Cov_xx*mv;                    % estimated parameters
0047 vv     = Am*est_x+av-lv;               % estimated residuals
0048 if R > 0
0049     sigma_0q = vv'*W_ll*vv/R;          % estimated variance factor
0050 else
0051     sigma_0q = 1;
0052 end
0053 
0054 %% diagnostics (useful if observations are uncorrelated)
0055 [riv,zv,nabla_lv,muv,muv1,uv1q,uv2] = diagnostics_GMM_1d(rU,Am,Cov_ll,W_ll,Cov_xx,Nm,vv);