0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026 function [xest,A,v,weights,Cov] = estimate_profile_robust...
0027 (N,ts,ys,sigma_e,sigma_n,Niter,type_out,type_r,...
0028 print_type,plot_type)
0029
0030
0031 k = 2.0;
0032
0033 type_robust = type_r(1);
0034 g_factor = type_r(2);
0035 w_factor = type_r(3);
0036 rob_p_d = type_r(4);
0037
0038
0039 M = length(ts);
0040 A = zeros(M+N-2,N);
0041
0042
0043 xa = zeros(N,1);
0044 weights = ones(M+N-2,1);
0045 Cov = 0;
0046 sigma_e0 = sigma_e;
0047
0048
0049
0050
0051
0052 gain=1;
0053
0054
0055 for iter = 1:Niter+1
0056 if plot_type> 0
0057 figure
0058 end
0059
0060 sigma_e = sigma_e0*gain^(Niter-iter);
0061
0062 if print_type > 0
0063 iter__sigmae_sigma_n_ratio = [iter,sigma_e,sigma_n,sigma_e/sigma_n]
0064 weight_before_iteration = weights(1:20)'
0065 end
0066
0067 b=zeros(M+N-2,1);
0068
0069 for m = 1:M
0070 w = 1/sigma_n*sqrt(weights(m));
0071 A(m,ts(m)) = w;
0072 dl = ys(m)-xa(ts(m));
0073 b(m) = dl*w;
0074 end
0075
0076
0077
0078
0079 kk = M;
0080
0081 for n = 2:N-1
0082 kk = kk+1;
0083 w = 1/sigma_e*sqrt(weights(kk));
0084 A(kk,n-1:n+1)= [1 -2 1]*w;
0085 dl = -(xa(n+1)-2*xa(n)+xa(n-1));
0086 b(kk) = dl*w;
0087 end
0088
0089
0090 Nmatrix = A'*A;
0091 nvector = A'*b;
0092 if N <= 1000
0093 Cov = inv(Nmatrix);
0094 end
0095
0096
0097 dx = Nmatrix\nvector;
0098
0099 xa = xa+dx;
0100
0101 if iter == Niter+1 && print_type > 0
0102 xa_dx = [xa';dx']
0103 end
0104 v = A*dx-b;
0105
0106 eso = norm(v)/sqrt(M-2);
0107
0108 vori = (v*sigma_e./sqrt(weights));
0109
0110 vs = vori;
0111
0112 eso_p = median(abs(vs(1:M)))*1.48;
0113
0114 eso_k = median(abs(vs(M+1:end)))*1.48;
0115
0116 if print_type > 0
0117 estimated_s0_s0p_s0k=[eso,eso_p,eso_k]
0118 end
0119
0120
0121 kp = k;
0122 kk = k*eso_k;
0123
0124
0125 if iter < Niter
0126 if type_out == 0
0127 if print_type > 0
0128 display('rob_p_d')
0129 end
0130 if iter < 4
0131 if rob_p_d == 1 || rob_p_d == 3
0132 weights(M+1:M+N-2) = ...
0133 min(1,1./(abs(vs(M+1:M+N-2))/kk+0.0001));
0134 if print_type > 1
0135 iter_kk_max_v_min_w = ...
0136 [iter,kk,max(abs(vs(M+1:M+N-2))),...
0137 min(weights(M+1:M+N-2))];
0138 display('iter_kk_max_v_min_w')
0139 display('modify dem weights')
0140 residuals_new_weights_of_dem=...
0141 [vs(M+1:M+N-2)';weights(M+1:M+N-2)'];
0142 display('residuals_new_weights_of_dem')
0143 end
0144 if plot_type > 0
0145 subplot(3,1,1)
0146 plot(2:N-1,xa(2:N-1),'-b')
0147 title(num2str(sigma_e/sigma_n));
0148 subplot(3,1,2);
0149 hold on
0150 plot(2:N-1,vs(M+1:M+N-2),'-b')
0151 plot(2:N-1,vori(M+1:M+N-2),'-r')
0152 plot(2:N-1,ones(N-2,1)*kk,'-b')
0153 plot(2:N-1,-ones(N-2,1)*kk,'-b')
0154 title(num2str(kk))
0155 subplot(3,1,3);
0156 plot(2:N-1,weights(M+1:M+N-2),'-b')
0157 title(num2str(iter))
0158 end
0159 end
0160 if rob_p_d == 2 || rob_p_d == 3
0161 weights(1:M) = min(1,1./(abs(vs(1:M))/kp+0.0001));
0162 if print_type > 1
0163 residuals_new_weights_of_points_L1 = ...
0164 [(1:20);v(1:20)';vs(1:20)';weights(1:20)']
0165 end
0166 end
0167 else
0168 if rob_p_d == 1 || rob_p_d == 3
0169 weights(M+1:M+N-2) = exp(-abs(vs(M+1:M+N-2)).^2./kk^2/2)+0.0001;
0170 if print_type > 1
0171 iter_kk_max_v_min_w = ...
0172 [iter,kk,max(abs(vs(M+1:M+N-2))),min(weights(M+1:M+N-2))]
0173 display('modify dem weights')
0174 residuals_new_weights_of_dem = ...
0175 [vs(M+1:M+N-2)';weights(M+1:M+N-2)']
0176 end
0177
0178 if plot_type > 0
0179 subplot(3,1,1)
0180 plot(2:N-1,xa(2:N-1),'-b')
0181 title(num2str(sigma_e/sigma_n));
0182 subplot(3,1,2);
0183 hold on
0184 plot(2:N-1,vs(M+1:M+N-2),'-b')
0185 plot(2:N-1,vori(M+1:M+N-2),'-r')
0186 plot(2:N-1,ones(N-2,1)*kk,'-b')
0187 plot(2:N-1,-ones(N-2,1)*kk,'-b')
0188 title(num2str(kk))
0189 subplot(3,1,3);
0190 plot(2:N-1,weights(M+1:M+N-2),'-b')
0191 title(num2str(iter))
0192 end
0193 end
0194 if rob_p_d == 2 || rob_p_d == 3
0195 weights(1:M) = exp(-abs(vs(1:M)).^2/kp^2/2)+0.0001;
0196 if print_type > 1
0197 residuals_new_weights_of_points_exp = ...
0198 [(1:20);v(1:20)';vs(1:20)'/kp;weights(1:20)']
0199 end
0200 end
0201 end
0202 else
0203 k = 1;
0204 g = g_factor*k;
0205 w = w_factor*k;
0206 if type_robust==0
0207 weights(1:M) = max(min(1,1./abs(v(1:M))/k+0.0001),((sign(v(1:M))/k)+1)/2);
0208 else
0209 for m=1:M
0210 if v(m) < g-w
0211 weights(m) = 0;
0212 elseif v(m) > g
0213 weights(m)=1;
0214 else
0215 weights(m)=1/(1+(v(m)-g)^4);
0216 end
0217 end
0218 end
0219 end
0220 else
0221
0222 weights = 0.9999*ceil(weights-0.1)+0.0001;
0223
0224 if iter == Niter
0225 xa = zeros(N,1);
0226 end
0227 end
0228 end
0229
0230 xest = xa;
0231
0232
0233