Home > 04-Estimation > GMM > Demo-GMM > fig_4_11_test_sensitivity_factors_GMM_similarity.m

fig_4_11_test_sensitivity_factors_GMM_similarity

PURPOSE ^

% Fig. 4.11 test sensitivity factors in GMM with similarity

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

% Fig. 4.11 test sensitivity factors in GMM with similarity

 see GMauss-Helmert model for similarity

 Wolfgang Förstner
 wfoerstn@uni-bonn.de

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 %% Fig. 4.11 test sensitivity factors in GMM with similarity
0002 %
0003 % see GMauss-Helmert model for similarity
0004 %
0005 % Wolfgang Förstner
0006 % wfoerstn@uni-bonn.de
0007 
0008 close all
0009 clearvars
0010 % clear all                                                                       %#ok<CLALL>
0011 
0012 addpath(genpath('../../../General-Functions/'));
0013 addpath(genpath('../Functions-GMM/'));
0014 
0015 disp(' ')
0016 disp(' ')
0017 disp('===========r========================================================')
0018 disp('---- Demo sensitivity factors in GMM with similarity Fig. 4.11 ----')
0019 disp('-------------------------------------------------------------------')
0020 
0021 %% Initiate parameters ####################################################
0022 
0023 % init_rand = 0;
0024 init_rand = 63;            % standard = 63
0025 init_rand_seed(init_rand);
0026 
0027 %type of generated data (1,2 or random) -----------------------------------
0028 
0029 %data_type = -10;        % type of generated data
0030 data_type =   2;         % standard = 2
0031 if data_type == 0
0032     data_type = -10;
0033 end
0034 
0035 % statistical parameters --------------------------------------------------
0036 
0037 % Significance level
0038 S = 0.99;
0039 
0040 % Tolerance
0041 tol = chi2inv(S,2);
0042 
0043 % selected parameters for sensititvity analysis ---------------------------
0044 %
0045 % 1,2 = translation (alternative choice)
0046 % 3,4 = rotation, scale
0047 
0048 rU = [3,4];
0049 
0050 %% ########################################################################
0051 disp('--------------------------------------------')
0052 disp(' planar similarity with Gauss-Markov Model  ')
0053 disp('--------------------------------------------')
0054 
0055 
0056 %% generate data ##########################################################
0057 % xis = a xi - b yi + c
0058 % yis = b xi + a yi + d
0059 switch data_type
0060     case 1
0061         tm  = [  1, 1; ...
0062             2, 2; ...
0063             -1,-2; ...
0064             -2,-1 ...
0065             ];
0066     case 2
0067         tm  = [  -7, 7 ;...
0068             1, 1; ...
0069             2, 2; ...
0070             -1,-2; ...
0071             -2,-1; ...
0072             ];
0073     otherwise
0074         tm=rand(-data_type,2)*5;
0075         
0076 end
0077 bb   = [min(tm(:,1)),max(tm(:,1)),min(tm(:,2)),max(tm(:,2))];
0078 maxd = max(max(tm(:,1))-min(tm(:,1)),max(tm(:,2))-min(tm(:,2)));
0079 
0080 % true parameters
0081 xt  = [2,0.5,3,-2]';
0082 disp('Parameters')
0083 disp('   [ x(1)     -x(2)      x(3) ]');
0084 disp('   [ x(2)      x(1)      x(4) ]');
0085 True_transformation = [xt(1) -xt(2) xt(3);xt(2) xt(1) xt(4)] %#ok<*NOPTS>
0086 
0087 %tm=rand(8,2)*2;
0088 N   = length(tm(:));
0089 d_I = 2;
0090 I   = N/d_I;
0091 Number_of_points = I
0092 U   = 4;
0093 sigma  = 0.01;
0094 Standard_deviation_observations_m = sigma
0095 
0096 
0097 %% setup elements of estimation
0098 % observed coordinates
0099 lm = zeros(I,d_I);
0100 % design matrix A =[A_i']
0101 Am = spalloc(N,U,N*U);
0102 Cov_ll_m = zeros(I,I-1);
0103 true_error = zeros(I,d_I);
0104 for i =1:I
0105     true_error(i,:) = randn(2,1)*sigma;
0106     lm(i,:)         = [xt(1) -xt(2) xt(3);xt(2) xt(1) xt(4)]*[tm(i,:)';1] ...
0107         + true_error(i,:)';
0108     Am(2*i-1:2*i,:) = [tm(i,1) -tm(i,2) 1 0;...
0109         tm(i,2)  tm(i,1) 0 1];                               %#ok<*SPRIX>
0110     Cov_ll_m(i,:)   = [1 0 0 1]*sigma^2;
0111 end
0112 true_errors=true_error
0113 true_and_observed_coordinates=[tm , lm]
0114 av  = zeros(N,1);
0115 
0116 %% factors for plotting
0117 factor_v  = 0.5*maxd/sqrt(N)/sigma;
0118 factor_r  = 0.5*maxd/sqrt(N);
0119 factor_X  = maxd/sqrt(N)/6;
0120 factor_mu = maxd/sqrt(N)/4;
0121 factor_nv = factor_v/10;
0122 
0123 %% estimate parameters
0124 
0125 [est_x,Cov_xx,sigma_0q,R,vm,Xv,Rim,nabla_lv,muv,muv1,Um1q,Um2]=...
0126     GaussMarkovModelLinear_groups(lm,Cov_ll_m,Am,av,rU);
0127 
0128 
0129 
0130 %% plot results
0131 screensize = plot_init;
0132 figure('Color','w','Position',[100 100 screensize+[-300 -300]]);
0133 
0134 subplot(2,3,1);
0135 hold on;
0136 plot(tm(:,1),tm(:,2),'.r','Markersize',12);       % given data
0137 for i=1:I
0138     plot([tm(i,1),tm(i,1)+factor_v*vm(i,1)],[tm(i,2),tm(i,2)+factor_v*vm(i,2)],'-k');
0139 end
0140 title({'residuals [m]', strcat('max=',num2str(max(norm(vm(:,1:2)))))});
0141 xlim([min(tm(:,1))-1,max(tm(:,1))+1]);
0142 ylim([min(tm(:,2))-1,max(tm(:,2))+1]);
0143 axis equal;
0144 
0145 
0146 subplot(2,3,2);
0147 hold on;
0148 for i=1:I
0149     plot_circle(tm(i,1),tm(i,2),factor_r*Rim(i,1),'-r');
0150 end
0151 title({'redundancy number', strcat('min=',num2str(min(Rim(:,1))))});
0152 xlim([min(tm(:,1))-factor_r*1,max(tm(:,1))+factor_r*1]);
0153 ylim([min(tm(:,2))-factor_r*1,max(tm(:,2))+factor_r*1]);
0154 axis equal;
0155 
0156 
0157 subplot(2,3,4);
0158 hold on;
0159 for i=1:I
0160     plot_circle(tm(i,1),tm(i,2),factor_X*sqrt(Xv(i)/d_I),'-k','LineWidth',2);
0161     plot_circle(tm(i,1),tm(i,2),factor_X*sqrt(tol/d_I),'--r');
0162 end
0163 title({'test statistics', strcat('max=',num2str(max(sqrt(Xv/d_I)))),'- - - critical values'});
0164 xlim([min(tm(:,1))-factor_X*5,max(tm(:,1))+factor_X*5]);
0165 ylim([min(tm(:,2))-factor_X*5,max(tm(:,2))+factor_X*5]);
0166 axis equal;
0167 
0168 subplot(2,3,3);
0169 hold on;
0170 for i=1:I
0171     plot_circle(tm(i,1),tm(i,2),factor_mu*muv(i),'-k');
0172 end
0173 title({'sensitivity factors', strcat('max=',num2str(max(muv)))});
0174 xlim([min(tm(:,1))-factor_mu*5,max(tm(:,1))+factor_mu*5]);
0175 ylim([min(tm(:,2))-factor_mu*5,max(tm(:,2))+factor_mu*5]);
0176 axis equal;
0177 
0178 
0179 subplot(2,3,5);
0180 hold on;
0181 for i=1:I
0182     plot_circle(tm(i,1),tm(i,2),factor_nv*nabla_lv(i),'-r');
0183 end
0184 title({'min. detectable outliers [m]';strcat('max=',num2str(max(nabla_lv)))});
0185 xlim([min(tm(:,1))-1,max(tm(:,1))+1]);
0186 ylim([min(tm(:,2))-1,max(tm(:,2))+1]);
0187 axis equal;
0188 
0189 if ~isempty(rU)
0190     subplot(2,3,6);
0191     hold on;
0192     for i=1:I
0193         plot_circle(tm(i,1),tm(i,2),factor_mu*muv1(i),'-k');
0194     end
0195     title({'sens. factors partial';strcat('max=',num2str(max(muv1)),', param=',num2str(rU))});
0196     xlim([min(tm(:,1))-factor_mu*5,max(tm(:,1))+factor_mu*5]);
0197     ylim([min(tm(:,2))-factor_mu*5,max(tm(:,2))+factor_mu*5]);
0198     axis equal;
0199 end
0200 
0201 
0202 %%
0203 
0204 disp('.............................')
0205 disp('        diagnostics          ')
0206 disp('.............................')
0207 Estimated_transformation = ...
0208     [est_x(1) -est_x(2) est_x(3);est_x(2) est_x(1) est_x(4)]
0209 Theoretical_precision_of_paramters   = full(sqrt(diag(Cov_xx)))'
0210 Redundancy              =R
0211 Estimated_sigma_0       = sqrt(sigma_0q)
0212 
0213 a________I________vx________vy_________z=...
0214     [(1:I)',vm,sqrt(Xv/d_I)]
0215 
0216 [m,i]=max(sqrt(diag(true_error*true_error')));
0217 disp(horzcat('Maximal true error                 = ',num2str(m,'% 12.5f'),...
0218     ' at: ', num2str(i,'% 3.0f')));
0219 [m,i]=max(sqrt(diag(vm*vm')));
0220 disp(horzcat('Maximal residual                   = ',num2str(m,'% 12.5f'),...
0221     ' at: ', num2str(i,'% 3.0f')));
0222 [m,i]=min(Rim(:,1));
0223 disp(horzcat('Minimal redundancy number          = ',num2str(m,'% 12.5f'),...
0224     ' at: ', num2str(i,'% 3.0f')));
0225 
0226 [m,i]=max(sqrt(Xv/d_I));
0227 if max(Xv) > tol
0228     disp(horzcat('Maximal test statistic             = ',num2str(m,'% 12.5f'),...
0229         ' at: ', num2str(i,'% 3.0f')), ' ***');
0230 else
0231     disp(horzcat('Maximal test statistic             = ',num2str(m,'% 12.5f'),...
0232         ' at: ', num2str(i,'% 3.0f')));
0233 end
0234 [m,i]=max(nabla_lv);
0235 disp(horzcat('Max. of minimal detectable outlier = ',num2str(m,'% 12.5f'),...
0236     ' at: ', num2str(i,'% 3.0f')));
0237 
0238 [m,i]=max(muv);
0239 disp(horzcat('Max. sensitivity factor            = ',num2str(m,'% 12.5f'),...
0240     ' at: ', num2str(i,'% 3.0f')));
0241 if ~isempty(rU)
0242     [m,i]=max(muv1);
0243     disp(horzcat('Max. sensitivity factor [',num2str(rU),']     = ',num2str(m,'% 12.5f'),...
0244         ' at: ', num2str(i,'% 3.0f')));
0245 end
0246 
0247 
0248

Generated on Sat 21-Jul-2018 20:56:10 by m2html © 2005