% Gauss Helmert Model
g(^l,^x) = 0
lv = Nx1 vector of observations
Cov_ll = NxN covariance matrix of observations (Cov_ll^a)
cg_f = Gx1 constraint function -> cg, A, B
xa = Ux1 vector of approximate parameters
sx = Ux1 vector of standard deviations of parameters
Tx = thershold for convergence
maxiter = maximum number of iterations
rU = range of unknown parameters of interest
est_x = Ux1 estimated parameter
Cov_xx = UxU theoretical covariance matrix of estimated parameters
may be full
sigma_0q = estimated variance factor (=1 if R=0)
R = redundancy
vv = Nx1 vector of estimated corrections
zv = Nx1 vector of standardized residuals (using sigma_0=1)
riv = Nx1 vector of redundancy numbers
nabla_lv = minimum detectabel outlier
muv = Nx1 vector of senstivity factors wrt to all parameters
muv1 = Nx1 vector of senstivity factors wrt to selected group
uv1q = Nx1 vector of contributions to selected group
uv2 = Nx1 vector of contribution to remaining parameters
[est_x,Cov_xx,sigma_0q,R,vv,zv,riv,nabla_lv,muv,muv1,uv1q,uv2]=...
GaussHelmertModelNonlinearNotrobust(lv,Cov_ll,cg_f,xa,sx,Tx,maxiter,rU)
Wolfgang Förstner 2016-06-06
wfoerstn@uni-bonn.de0001 %% Gauss Helmert Model 0002 % 0003 % g(^l,^x) = 0 0004 % 0005 % lv = Nx1 vector of observations 0006 % Cov_ll = NxN covariance matrix of observations (Cov_ll^a) 0007 % cg_f = Gx1 constraint function -> cg, A, B 0008 % xa = Ux1 vector of approximate parameters 0009 % sx = Ux1 vector of standard deviations of parameters 0010 % Tx = thershold for convergence 0011 % maxiter = maximum number of iterations 0012 % rU = range of unknown parameters of interest 0013 % 0014 % est_x = Ux1 estimated parameter 0015 % Cov_xx = UxU theoretical covariance matrix of estimated parameters 0016 % may be full 0017 % sigma_0q = estimated variance factor (=1 if R=0) 0018 % R = redundancy 0019 % vv = Nx1 vector of estimated corrections 0020 % zv = Nx1 vector of standardized residuals (using sigma_0=1) 0021 % riv = Nx1 vector of redundancy numbers 0022 % nabla_lv = minimum detectabel outlier 0023 % muv = Nx1 vector of senstivity factors wrt to all parameters 0024 % muv1 = Nx1 vector of senstivity factors wrt to selected group 0025 % uv1q = Nx1 vector of contributions to selected group 0026 % uv2 = Nx1 vector of contribution to remaining parameters 0027 % 0028 % [est_x,Cov_xx,sigma_0q,R,vv,zv,riv,nabla_lv,muv,muv1,uv1q,uv2]=... 0029 % GaussHelmertModelNonlinearNotrobust(lv,Cov_ll,cg_f,xa,sx,Tx,maxiter,rU) 0030 % 0031 % Wolfgang Förstner 2016-06-06 0032 % wfoerstn@uni-bonn.de 0033 0034 0035 function [est_x,Cov_xx,sigma_0q,R,vv,zv,riv,nabla_lv,muv,muv1,uv1q,uv2] = ... 0036 GaussHelmertModel(lv,Cov_ll,cg_f,xa,sx,Tx,maxiter,rU) 0037 0038 %% initialization of estimation 0039 0040 % numbers G, and U 0041 [cg,Am,Bm] = cg_f(lv,lv,xa); 0042 G = size(cg,1); 0043 U = size(xa,1); 0044 0045 % redundancy 0046 R = G-U; 0047 if R < 0 0048 disp('not enough observations') 0049 return; 0050 end 0051 %% initialization of iterations 0052 0053 nu = 0; % number of iterations 0054 la = lv; 0055 lnu = la; 0056 xnu = xa; 0057 s = 0; 0058 0059 %% iterations 0060 0061 for iter = 1: maxiter 0062 [cg,Am,Bm] = cg_f(lv,lnu,xnu); % residual if constraints, Jacobians 0063 % W_ll = inv(Cov_ll); % weight matrix of lv 0064 W_gg = inv(Bm' * Cov_ll * Bm); 0065 ABWB = Am' * W_gg; %#ok<MINV> 0066 Nm = ABWB * Am; % normal equation matrix 0067 nv = ABWB * cg; % right hand sides 0068 lambda_N = eigs(Nm); 0069 if abs(log(lambda_N(1)/lambda_N(U))) > 10 0070 disp('normal equation matrix nearly singular') 0071 return; 0072 end 0073 Cov_xx = inv(Nm); % CovM of parameters 0074 delta_x = Cov_xx * nv; %#ok<MINV> % correction of parameters 0075 0076 nu = nu+1; 0077 0078 % check for final iteration 0079 if max(abs(delta_x(:)./sx(:))) < Tx || nu == maxiter 0080 s = 2; 0081 end 0082 % correction of parameters 0083 xnu = xnu + delta_x; 0084 0085 disp([num2str(nu),'. Iteration: delta_x = [',num2str(delta_x'),... 0086 '], est_x = [',num2str(xnu'),']']) 0087 0088 % correction to fitted observations 0089 vv = lnu - lv; % residuals 0090 delta_l = Cov_ll * Bm * W_gg * (cg - Am * delta_x) - vv; %#ok<MINV> 0091 lnu = lnu + delta_l; % fitted observations 0092 0093 if s==2 0094 break; 0095 end 0096 end 0097 if R > 0 0098 sigma_0q = (cg' * W_gg *cg)/R; %#ok<MINV> 0099 else 0100 sigma_0q = 1; 0101 end 0102 est_x = xnu; 0103 0104 %% diagnostics (useful if observations are uncorrelated) 0105 [riv,zv,nabla_lv,muv,muv1,uv1q,uv2] = ... 0106 diagnostics_GHM_1d(rU,Am,Bm,Cov_ll,Cov_xx,Nm,W_gg,vv);