0001
0002
0003
0004
0005
0006
0007
0008 close all
0009 clearvars
0010
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
0022
0023
0024 init_rand = 63;
0025 init_rand_seed(init_rand);
0026
0027
0028
0029
0030 data_type = 2;
0031 if data_type == 0
0032 data_type = -10;
0033 end
0034
0035
0036
0037
0038 S = 0.99;
0039
0040
0041 tol = chi2inv(S,2);
0042
0043
0044
0045
0046
0047
0048 rU = [3,4];
0049
0050
0051 disp('--------------------------------------------')
0052 disp(' planar similarity with Gauss-Markov Model ')
0053 disp('--------------------------------------------')
0054
0055
0056
0057
0058
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
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)]
0086
0087
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
0098
0099 lm = zeros(I,d_I);
0100
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];
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
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
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
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);
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