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