% diagnostics GMM 1d r_U = range of paraemters of interest Am = design matrix Cov_ll = CovM of observations Cov_xx = CovM of parameters W_xx = weight matrix of parameters vv = vector of residuals riv = redundancy numbers of observations zv = standardized resiuals nabla_lv = minimum size of detectable outliers muv = sensitivity factors wrt to all parameters uv1q = effect onty parameters of interest uv2 = effect onto nuisance parameters Wolfgang Förstner 6/6/2015 wfoerstn@uni-bonn.de
0001 %% diagnostics GMM 1d 0002 % 0003 % r_U = range of paraemters of interest 0004 % Am = design matrix 0005 % Cov_ll = CovM of observations 0006 % Cov_xx = CovM of parameters 0007 % W_xx = weight matrix of parameters 0008 % vv = vector of residuals 0009 % 0010 % riv = redundancy numbers of observations 0011 % zv = standardized resiuals 0012 % nabla_lv = minimum size of detectable outliers 0013 % muv = sensitivity factors wrt to all parameters 0014 % uv1q = effect onty parameters of interest 0015 % uv2 = effect onto nuisance parameters 0016 % 0017 % Wolfgang Förstner 6/6/2015 0018 % wfoerstn@uni-bonn.de 0019 0020 function [riv,zv,nabla_lv,muv,muv1,uv1q,uv2] = ... 0021 diagnostics_GHM_1d(r_U,Am,Bm,Cov_ll,Cov_xx,W_xx,W_gg,vv) 0022 0023 delta_0 = 4.13; 0024 0025 U = size(Cov_xx,1); 0026 0027 %% observations: detectatbility, statistical test 0028 Rm = Cov_ll*Bm*W_gg*Bm'-Cov_ll*Bm*W_gg*Am*Cov_xx*Am'*W_gg*Bm'; % redundancy matrix 0029 riv = diag(Rm); % redundancy numbers 0030 zv = vv./sqrt(riv.*diag(Cov_ll)+eps); % normalized residuals 0031 nabla_lv = delta_0 * sqrt(1./riv); 0032 0033 %% sensitivity wrt all parameters 0034 muv = sqrt((1-riv)./(riv+eps)); % sensitivity factors 0035 0036 %% sensitivity wrt selected parameters 0037 if ~isempty(r_U) 0038 % ranges 0039 r_C= r_U; 0040 r_D= setdiff(1:U,r_U); 0041 % partitioned Jacobian matrix A=[C,D], G x [U_C,U_D] 0042 Cm = Am(:,r_C); % G * U_C 0043 Dm = Am(:,r_D); % G * U_D 0044 % reduced design matrix 0045 Cmr = Cm - Dm * inv(W_xx(r_D,r_D)) * W_xx(r_D,r_C); %#ok<MINV> % C_reduced 0046 Cov_11 = Cov_xx(r_C,r_C); 0047 % effect of constraints onto parameters 0048 U1qc = Cmr * Cov_11 * Cmr' * W_gg; % U_C_bar 0049 U2c = Dm * inv(W_xx(r_D,r_D)) * Dm' * W_gg; %#ok<MINV> % U_D 0050 0051 % --- for checking 0052 % uv1qc = diag(U1qc); % u_C_bar 0053 % uv2c = diag(U2c) ; % u_D 0054 % Rmc = eye(G) - Am*Cov_xx*Am'*W_gg; 0055 % rivc = diag(Rmc); 0056 % check =uv1qc+uv2c+rivc ; % should be ones 0057 0058 % sensitivity factors for constraints 0059 % wrt observations 0060 U1q = Cov_ll*Bm*W_gg * U1qc *Bm'; 0061 uv1q = diag(U1q); 0062 U2 = Cov_ll*Bm*W_gg * U2c * Bm'; 0063 uv2 = diag(U2); 0064 0065 % --- for checking 0066 % muv1c = sqrt(uv1qc ./ (rivc+eps)); 0067 % check =uv1q+uv2+riv; % should be ones 0068 0069 % sensitivity factors for observations 0070 muv1 = sqrt(uv1q ./ (riv+eps)); 0071 end