diagnostics_GHM_multi_d

% Diagnostics Gauss-Helmert model multidimensional

 using sparse inverse

 params see GaussHelmertModel

 Wolfgang Förstner
 wfoerstn@uni-bonn.de
 
 See also GaussHelmertModel

0001 %% Diagnostics Gauss-Helmert model multidimensional
0002 %
0003 % using sparse inverse
0004 %
0005 % params see GaussHelmertModel
0006 %
0007 % Wolfgang Förstner
0008 % wfoerstn@uni-bonn.de
0009 %
0010 % See also GaussHelmertModel
0011 
0012 function [vi,Xv,Rim,nabla_lv,muv,muv1,U1mq,U2m] = diagnostics_GHM_multi_d...
0013     (r_U,I,d_I,d_G,Am,Bm,Cov_ll,W_ll,Cov_xx,W_xx,vv)
0014 
0015 [~,U] = size(Am);
0016 
0017 delta_0=4.13;  % all sizes referring to d_I.
0018 
0019 if ~isempty(r_U)
0020     % ranges
0021     r_C = r_U;
0022     r_D = setdiff(1:U,r_U);
0023     % partitioned design matrix A=[C,D]
0024     Cm  = Am(:,r_C);
0025     Dm  = Am(:,r_D);
0026     % reduced design matrix
0027     Cmr    = Cm - Dm / W_xx(r_D,r_D) * W_xx(r_D,r_C);  % C_reduced
0028     Cov_11 = Cov_xx(r_C,r_C);
0029 end
0030 %% observational groups
0031 Rim  = zeros(I,d_I^2);
0032 U1mq = zeros(I,d_I^2);
0033 U2m  = zeros(I,d_I^2);
0034 vi   = zeros(I,d_I);
0035 Xv   = zeros(I,1);
0036 nabla_lv = zeros(I,1);
0037 muv      = zeros(I,1);
0038 muv1     = zeros(I,1);
0039 for i=1:I
0040     % detectability, test statistics, sensitivity wrt all parameters
0041     i_range     = d_I*i-(d_I-1):d_I*i;
0042     g_range     = d_G*i-(d_G-1):d_G*i;
0043     Am_i        = Am(g_range,:)';
0044     Bm_i        = Bm(i_range,g_range);
0045     Cov_ll_ii   = reshape(Cov_ll(i,:),d_I,d_I);  % CovM of group
0046     W_ll_ii     = W_ll(i_range,i_range);         % Weight matrix of group
0047     W_gg_ii     = Bm_i' * Cov_ll_ii * Bm_i;      % Weight matrix of group
0048     % Cov_ll*Bm*W_gg*Bm'-Cov_ll*Bm*W_gg*Am*Cov_xx*Am'*W_gg*Bm'
0049     Cov_vv_ii   = Cov_ll_ii * Bm_i *...
0050         (eye(d_G) - W_gg_ii *Am_i'*Cov_xx*Am_i)*...
0051         W_gg_ii*Bm_i'*Cov_ll_ii;  % covariance matrix of resiudal
0052     W_vv_ii     = inv(Cov_vv_ii+eps);       % Weight matrix of residuals
0053     Rii         = Cov_vv_ii*W_ll_ii;        % redundancy matrix
0054     Rim(i,:)    = Rii(:)';                  % vectorized Rii
0055     Rii_inv     = full(inv(Rii));           % inverse of redundancy matrix
0056     vi(i,:)     = vv(i_range)';             % residual of group
0057     Xv(i)       = vv(i_range)'*W_vv_ii* vv(i_range);      %#ok<MINV>
0058     % Xhi^2-test statistic
0059     nabla_lv(i) = delta_0 * ...
0060         sqrt(real(max(eig(Rii_inv*Cov_ll_ii))));
0061     % minimum detectable outlier
0062     muv(i)      = sqrt(real(max(eig(Rii_inv-eye(d_I)))));
0063     % sensitivity factors
0064     
0065     % sensitivity wrt selected parameters range rU
0066     if ~isempty(r_U)
0067         % U's for constraints
0068         U1qc       = Cmr(g_range,:) * Cov_11  * Cmr(g_range,:)' * W_gg_ii;
0069         U2c        = Dm(g_range,:)  / (W_xx(r_D,r_D))  * ...
0070             Dm(g_range,:)' * W_gg_ii;
0071         % check      = U1qc + U2c + (eye(G)-Am_i'*Cov_xx*Am_i*inv(Cov_gg_ii));
0072         % I=U_C_bar+U_D+R, check-eye(d_I)
0073         
0074         % U's for observations
0075         U1q       = Cov_ll_ii*Bm_i*W_gg_ii * U1qc * Bm_i'; %idempotent, rank |r_U|
0076         U1mq(i,:) = U1q(:)';
0077         U2        = Cov_ll_ii*Bm_i*W_gg_ii * U2c * Bm_i'; %idempotent, rank U-|r_U|
0078         U2m(i,:)  = U2(:)';
0079         % check     = U1q + U2 + Rii;   % idempotent, rank G
0080         
0081         % sensitivity factors for seleced parameters
0082         muv1(i)   = sqrt(real(max(eig(U1q * Rii_inv))));
0083     end
0084 end