GaussMarkovModelLinear_groups

% Gauss-Markov model linear, groups of uncorrel. observations

 all observational groups have the same dimension
 all observational groups are mutually uncorrelated

 input matrix Am may be sparse
 covariance matrix of estimates is assumed to be full

 lm       = I x d_I matrix of observational groups
 Cov_ll_m = I x d_i^2 matrix of vectorized covariance matrices
 Am       = NxU Jacobian (may be sparse)
 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=0)
 R        = redundancy
 vm       = I x d_I matrix of estimated residuals
 Xv       = Ix1 vector of test statistics  (Xhi^2_d_I using sigma_0=1)
 Rim      = I x d_i^2 matrix of redundacy matrices R_ii
 nabla_lv = size of minimal detectable errors
            (from max eigenvalue of C_ll W_vv C_ll * delta_0
 muv      = Ix1 vector of senstivity factors
               (mu_i:=mu_ix, cf. (4.307))
 muv1     = Ix1 vector of senstivity factors for selected parameters
               (mu_i:=mu_ik, cf. (4.315))
 Um1q     = I x d_i^2 matrix of reduced U_ii-matrices
              for selected parameters (cf. bar(U)_k(4.127))
 Um2      = I x d_i^2 matrix of U_ii-matrices
              for not selected parameters (cf. U_t (4.127))

 [est_x,Cov_xx,sigma_0q,R,vm,Xv,Rim,muv]=...
            GaussMarkovModelLinear_groups(lm,Cov_ll,Am,av)

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

0001 %% Gauss-Markov model linear, groups of uncorrel. observations
0002 %
0003 % all observational groups have the same dimension
0004 % all observational groups are mutually uncorrelated
0005 %
0006 % input matrix Am may be sparse
0007 % covariance matrix of estimates is assumed to be full
0008 %
0009 % lm       = I x d_I matrix of observational groups
0010 % Cov_ll_m = I x d_i^2 matrix of vectorized covariance matrices
0011 % Am       = NxU Jacobian (may be sparse)
0012 % av       = Nx1 additive vector
0013 % rU       = range of unknown parameters of interest
0014 %
0015 % est_x    = Ux1 estimated parameter
0016 % Cov_xx   = UxU theoretical covariance matrix of estimated parameters
0017 % sigma_0q = estimated variance factor (=1 if R=0)
0018 % R        = redundancy
0019 % vm       = I x d_I matrix of estimated residuals
0020 % Xv       = Ix1 vector of test statistics  (Xhi^2_d_I using sigma_0=1)
0021 % Rim      = I x d_i^2 matrix of redundacy matrices R_ii
0022 % nabla_lv = size of minimal detectable errors
0023 %            (from max eigenvalue of C_ll W_vv C_ll * delta_0
0024 % muv      = Ix1 vector of senstivity factors
0025 %               (mu_i:=mu_ix, cf. (4.307))
0026 % muv1     = Ix1 vector of senstivity factors for selected parameters
0027 %               (mu_i:=mu_ik, cf. (4.315))
0028 % Um1q     = I x d_i^2 matrix of reduced U_ii-matrices
0029 %              for selected parameters (cf. bar(U)_k(4.127))
0030 % Um2      = I x d_i^2 matrix of U_ii-matrices
0031 %              for not selected parameters (cf. U_t (4.127))
0032 %
0033 % [est_x,Cov_xx,sigma_0q,R,vm,Xv,Rim,muv]=...
0034 %            GaussMarkovModelLinear_groups(lm,Cov_ll,Am,av)
0035 %
0036 % Wolfgang Förstner 2015-06-04
0037 % wfoerstn@uni-bonn.de
0038 
0039 
0040 function [est_x,Cov_xx,sigma_0q,R,vi,Xv,Rim,nabla_lv,muv,muv1,Um1q,Um2] = ...
0041     GaussMarkovModelLinear_groups(lm,Cov_ll_m,Am,av,rU)
0042 
0043 %% initialization
0044 
0045 % numbers I, d_I
0046 [I,d_I]= size(lm);
0047 
0048 % number N and U
0049 [N,U] = size(Am);
0050 
0051 % redundancy
0052 R = N-U;
0053 if R < 0
0054     disp('not enough observations')
0055     return;
0056 end
0057 
0058 % reshape observations and covariance matrix
0059 lv     = zeros(N,1);
0060 Cov_ll = spalloc(N,N,I*d_I^2);
0061 for i = 1:I
0062     lv(d_I*i-(d_I-1):d_I*i)=lm(i,:)';
0063     Cov_ll(d_I*i-(d_I-1):d_I*i,d_I*i-(d_I-1):d_I*i) = ...
0064         reshape(Cov_ll_m(i,:),d_I,d_I);
0065 end
0066 %% estimation
0067 W_ll   = inv(Cov_ll);                  % weight matrix of lv
0068 Bm     = W_ll*Am;                      %#ok<*MINV> % ancillary matrix
0069 Nm     = Bm'*Am;                       % normal equation matrix
0070 mv     = Bm'*(lv-av);                  % right hand sides
0071 
0072 % use sparseinv in Matlab, otherwise inv
0073 %Cov_xx = sparseinv(Nm);               % covariance matrix of est. parameters
0074 Cov_xx = inv(Nm);                      % covariance matrix of est. parameters
0075 est_x  = Nm\mv;                        % estimated parameters
0076 vv     = Am*est_x+av-lv;               % estimated residuals
0077 if R > 0
0078     sigma_0q = vv'*W_ll*vv/R;          % estimated variance factor
0079 else
0080     sigma_0q = 1;
0081 end
0082 %% diagnostics
0083 [vi,Xv,Rim,nabla_lv,muv,muv1,Um1q,Um2]=...
0084     diagnostics_GMM_multi_d(rU,I,d_I,Am,Cov_ll,W_ll,Cov_xx,Nm,vv);