% Diagnostics Gauss-Helmert model multidimensional
referring to constraints
assuming each group of observations belongs to one group of constraints
number = I, with d_I observations and d_G constraints
[vgi,Xv,Rgm,nabla_cg,muv,muv1,U1mqc,U2mc]=
diagnostics_GHM_constraints_multi_d...
(r_U,I,d_I,d_G,Am,Bm,Cov_xx,W_xx,W_gg,vv)
r_U = index set for parameters to be estimated (others are nuisance)
I = number of observational groups
d_I = dimension of observatioanl groups (must be the same for all)
Am = G x U Jacobian
Bm = G x N Jacobain
Cov_xx = U x U CovM of all estimated parameters
W_xx = U x U WeightM of all estimated parameters
Wgg =
vv = N x 1 residuals
vgi = I x dI matrix of residuals
Xv = I x 1 vector standardized test statistics (chi^2)
Rim = I * d_I^2 redundancy matrices as row vectors
nabla_lv = I x 1 lower bound for detectable outliers (maximum eigenvalue)
muv = I x 1 sensitivity factor for all parameters
muv1 = I x 1 sensitivity factor for parameters specified by r_U
Um1q = I * d_I^2 diagonal blocks of \bar U_1 (numerator of muv1^2)
Um2 = I * d_I^2 diagonal blocks of U2
Wolfgang Förstner
wfoerstn@uni-bonn.de
See also GaussHelmertModel0001 %% Diagnostics Gauss-Helmert model multidimensional 0002 % 0003 % referring to constraints 0004 % assuming each group of observations belongs to one group of constraints 0005 % number = I, with d_I observations and d_G constraints 0006 % 0007 % [vgi,Xv,Rgm,nabla_cg,muv,muv1,U1mqc,U2mc]= 0008 % diagnostics_GHM_constraints_multi_d... 0009 % (r_U,I,d_I,d_G,Am,Bm,Cov_xx,W_xx,W_gg,vv) 0010 % 0011 % r_U = index set for parameters to be estimated (others are nuisance) 0012 % I = number of observational groups 0013 % d_I = dimension of observatioanl groups (must be the same for all) 0014 % Am = G x U Jacobian 0015 % Bm = G x N Jacobain 0016 % Cov_xx = U x U CovM of all estimated parameters 0017 % W_xx = U x U WeightM of all estimated parameters 0018 % Wgg = 0019 % vv = N x 1 residuals 0020 % 0021 % vgi = I x dI matrix of residuals 0022 % Xv = I x 1 vector standardized test statistics (chi^2) 0023 % Rim = I * d_I^2 redundancy matrices as row vectors 0024 % nabla_lv = I x 1 lower bound for detectable outliers (maximum eigenvalue) 0025 % muv = I x 1 sensitivity factor for all parameters 0026 % muv1 = I x 1 sensitivity factor for parameters specified by r_U 0027 % Um1q = I * d_I^2 diagonal blocks of \bar U_1 (numerator of muv1^2) 0028 % Um2 = I * d_I^2 diagonal blocks of U2 0029 % 0030 % 0031 % Wolfgang Förstner 0032 % wfoerstn@uni-bonn.de 0033 % 0034 % See also GaussHelmertModel 0035 0036 function [vgi,Xv,Rgm,nabla_cg,muv,muv1,U1mqc,U2mc]=diagnostics_GHM_constraints_multi_d... 0037 (r_U,I,d_I,d_G,Am,Bm,Cov_xx,W_xx,W_gg,vv) 0038 0039 [~,U] = size(Am); 0040 0041 delta_0=4.13; % all sizes referring to d_J. 0042 0043 if ~isempty(r_U) 0044 % ranges 0045 r_C= r_U; 0046 r_D= setdiff(1:U,r_U); 0047 % partitioned design matrix A=[C,D] 0048 Cm = Am(:,r_C); 0049 Dm = Am(:,r_D); 0050 % reduced design matrix 0051 Cmr = Cm - Dm / W_xx(r_D,r_D) * W_xx(r_D,r_C); % C_reduced 0052 Cov_11 = Cov_xx(r_C,r_C); 0053 end 0054 0055 %% observational groups 0056 Rgm = zeros(I,d_G^2); 0057 U1mqc = zeros(I,d_G^2); 0058 U2mc = zeros(I,d_G^2); 0059 Xv = zeros(I,1); 0060 vgi = zeros(I,d_G); 0061 nabla_cg = zeros(I,1); 0062 muv = zeros(I,1); 0063 muv1 = zeros(I,1); 0064 0065 for i=1:I 0066 % detectability, test statistics, sensitivity wrt all parameters 0067 i_range = d_I*i-(d_I-1):d_I*i; 0068 g_range = d_G*i-(d_G-1):d_G*i; 0069 Amt_i = Am(g_range,:); 0070 Bmt_i = Bm(i_range,g_range)'; 0071 %W_gg_ii = inv(Bmt_i * Cov_ll_ii * Bmt_i'); % Weight matrix of group 0072 W_gg_ii = W_gg(g_range,g_range); % Weight matrix of group 0073 Cov_gg_ii = inv(W_gg_ii); 0074 Cov_vgvg_ii = Cov_gg_ii - Amt_i * Cov_xx * Amt_i'; 0075 W_vgvg_ii = inv(Cov_vgvg_ii+eps); 0076 % redundancy matrix of constraints 0077 R_gg_ii = full(Cov_vgvg_ii * W_gg_ii); 0078 R_gg_ii_inv = inv(R_gg_ii); 0079 Rgm(i,:) = R_gg_ii(:)'; % vectorized Rgg 0080 vgi(i,:) = Bmt_i*vv(i_range); % residual of group 0081 Xv(i) = vgi(i,:) * W_vgvg_ii * vgi(i,:)'; %#ok<MINV> 0082 % Xhi^2-test statistic 0083 nabla_cg(i) = delta_0 * ... 0084 sqrt(real(max(eig( R_gg_ii \ inv(W_gg_ii) )))); 0085 % minimum detectable outlier 0086 muv(i) = sqrt(real(max(eig(R_gg_ii_inv - eye(d_G) )))); 0087 % sensitivity factors 0088 0089 % sensitivity wrt selected parameters range rU 0090 0091 if ~isempty(r_U) 0092 % U's for constraints 0093 U1qc = Cmr(g_range,:) * Cov_11 * Cmr(g_range,:)' * W_gg_ii; 0094 U1mqc(i,:) = U1qc(:)'; 0095 U2c = Dm(g_range,:) / W_xx(r_D,r_D) * ... 0096 Dm(g_range,:)' * W_gg_ii; 0097 U2mc(i,:) = U2c(:)'; 0098 % check = U1qc + U2c + R_gg_ii - eye(d_G) 0099 % I=U_C_bar+U_D+R, check-eye(d_I) 0100 0101 % sensitivity factors 0102 muv1(i) = sqrt(real(max(eig(U1qc / R_gg_ii )))); 0103 end 0104 end 0105