% Diagnostics Gauss-Markov model multidimensional
given all parameters of a linearized estimation problem
determine the essential diagnostic values
assuming the number N of observations is small enough
(requires several NxN matrices)
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 = Jacobian
Cov_ll = CovM of observations
W_ll = WeightM of observations
Cov_xx = CovM of all estimated parameters
W_xx = WeightM of all estimated parameters
vv = residuals
Xv = standardized test statistics (chi^2)
Rim = redundancy matrices as row vectors (I * d_I^2)
nabla_lv = lower bound for detectable outliers (maximum eigenvalue)
muv = sensitivity factor for all parameters
muv1 = sensitivity factor for parameters specified by r_U
Um1q = diagonal blocks of \bar U_1 (numerator of muv1^2)
Um2 = diagonal blocks of U2
Wolfgang Förstner 10/2016
wfoerstn@uni-bonn.de0001 %% Diagnostics Gauss-Markov model multidimensional 0002 % 0003 % given all parameters of a linearized estimation problem 0004 % determine the essential diagnostic values 0005 % assuming the number N of observations is small enough 0006 % (requires several NxN matrices) 0007 % 0008 % r_U = index set for parameters to be estimated (others are nuisance) 0009 % I = number of observational groups 0010 % d_I = dimension of observatioanl groups (must be the same for all) 0011 % Am = Jacobian 0012 % Cov_ll = CovM of observations 0013 % W_ll = WeightM of observations 0014 % Cov_xx = CovM of all estimated parameters 0015 % W_xx = WeightM of all estimated parameters 0016 % vv = residuals 0017 % 0018 % 0019 % Xv = standardized test statistics (chi^2) 0020 % Rim = redundancy matrices as row vectors (I * d_I^2) 0021 % nabla_lv = lower bound for detectable outliers (maximum eigenvalue) 0022 % muv = sensitivity factor for all parameters 0023 % muv1 = sensitivity factor for parameters specified by r_U 0024 % Um1q = diagonal blocks of \bar U_1 (numerator of muv1^2) 0025 % Um2 = diagonal blocks of U2 0026 % 0027 % 0028 % Wolfgang Förstner 10/2016 0029 % wfoerstn@uni-bonn.de 0030 0031 0032 function [vi,Xv,Rim,nabla_lv,muv,muv1,U1mq,U2m] = diagnostics_GMM_multi_d... 0033 (r_U,I,d_I,Am,Cov_ll,W_ll,Cov_xx,W_xx,vv) 0034 0035 [~,U] = size(Am); 0036 0037 delta_0 = 4.13; % all sizes referring to d=1. 0038 0039 if ~isempty(r_U) 0040 % ranges 0041 r_C = r_U; 0042 r_D = setdiff(1:U,r_U); 0043 % partitioned design matrix A=[C,D] 0044 Cm = Am(:,r_C); 0045 Dm = Am(:,r_D); 0046 % reduced design matrix 0047 Cmr = Cm - Dm * inv(W_xx(r_D,r_D)) * W_xx(r_D,r_C); % C_reduced 0048 Cov_11 = Cov_xx(r_C,r_C); 0049 end 0050 %% observational groups 0051 Rim = zeros(I,d_I^2); 0052 U1mq = zeros(I,d_I^2); 0053 U2m = zeros(I,d_I^2); 0054 vi = zeros(I,d_I); 0055 Xv = zeros(I,1); 0056 nabla_lv = zeros(I,1); 0057 muv = zeros(I,1); 0058 muv1 = zeros(I,1); 0059 for i=1:I 0060 % detectability, test statistics, sensitivity wrt all parameters 0061 i_range = d_I*i-(d_I-1):d_I*i; 0062 Ai = Am(i_range,:)'; 0063 Cov_ll_ii = Cov_ll(i_range,i_range); % CovM of group 0064 W_ll_ii = W_ll(i_range,i_range); % Weight matrix of group 0065 Cov_vv_ii = Cov_ll_ii-Ai'*Cov_xx*Ai; % covariance matrix of resiudal 0066 W_vv_ii = inv(Cov_vv_ii+eps); % Weihgt matrix of residuals 0067 Rii = Cov_vv_ii*W_ll_ii; % redundancy matrix 0068 Rim(i,:) = Rii(:)'; % vectorized Rii for storing 0069 Rii_inv = full(inv(Rii)); % inverse of redundancy matrix 0070 vi(i,:) = vv(i_range)'; % residual of group 0071 Xv(i) = vv(i_range)'*W_vv_ii* vv(i_range); %#ok<*MINV> 0072 % Xhi^2-test statistic 0073 nabla_lv(i) = delta_0 * ... 0074 sqrt(real(max(eig(Rii_inv*Cov_ll_ii)))); 0075 % minimum detectable outlier 0076 muv(i) = sqrt(real(max(eig(Rii_inv-eye(d_I))))); 0077 % sensitivity factors 0078 0079 % sensitivity wrt selected parameters range rU 0080 if ~isempty(r_U) 0081 U1q = Cmr(i_range,:) * Cov_11 * Cmr(i_range,:)' * W_ll_ii; 0082 U1mq(i,:) = U1q(:)'; 0083 U2 = Dm(i_range,:) / W_xx(r_D,r_D) * ... 0084 Dm(i_range,:)' * W_ll_ii; 0085 U2m(i,:) = U2(:)'; 0086 %check = U1q + U2 + Rii-eye(d_I) % I=U_C_bar+U_D+R 0087 %i,check-eye(d_I) 0088 muv1(i) = sqrt(real(max(eig(U1q * Rii_inv)))); 0089 % sensitivity factors 0090 else 0091 U1mq=0; 0092 U2m =0; 0093 muv1=0; 0094 end 0095 end