Description of estimate_profile

0001 %% estimate_profile_robust
0002 %
0003 % xest = estimate_profile_robust(N,ts,ys,sigma_e,sigma_n,
0004 %                    Niter,type_out,type_r,print_type,plot_type)
0005 %
0006 % N        = scalar, length of profile
0007 % ts,ys    = Mx1 vectors, observed profile
0008 % sigma_e  = std of curvature
0009 % sigma_n  = std of observations
0010 % Niter    = number of iterations
0011 %          = 1 -> no robust estimation
0012 % type_out = type of outlier (0=symm, 1=asymm)
0013 % type_r   (1) > 0 = L1, 1 = Kraus
0014 %          (2) = g for Kraus
0015 %          (3) = w for Kraus
0016 %          (4) = 1 only dem
0017 %                2 only points
0018 % plot_type    = 0,1,2 (no, low, high)
0019 % print_type   = 0,1,2 (no, low, high)
0020 %
0021 % Wolfgang Förstner 2014-10-05
0022 % last changes: Susanne Wenzel 09/16
0023 % wfoerstn@uni-bonn.de, wenzel@igg.uni-bonn.de
0024 
0025 function [xest,A,v,weights,Cov] = estimate_profile_robust...
0026     (N,ts,ys,sigma_e,sigma_n,Niter,type_out,type_r,...
0027     print_type,plot_type)
0028 
0029 % robust - explore sensistivity of result wrt k! Choose k in [1:3]
0030 k = 2; 
0031 
0032 type_robust = type_r(1);
0033 g_factor    = type_r(2);
0034 w_factor    = type_r(3);
0035 rob_p_d     = type_r(4);
0036 
0037 % Number of observations
0038 M = length(ts);
0039 A = zeros(M+N-2,N);
0040 
0041 % initialize
0042 xa = zeros(N,1);
0043 weights = ones(M+N-2,1);
0044 Cov = 0;
0045 sigma_e0 = sigma_e;
0046 gain = 1;
0047 
0048 %% robust estimation
0049 for iter = 1:Niter+1
0050     if plot_type> 0
0051         figure
0052     end
0053     
0054     sigma_e = sigma_e0*gain^(Niter-iter);
0055     
0056     if print_type > 0
0057         iter__sigmae_sigma_n_ratio = [iter,sigma_e,sigma_n,sigma_e/sigma_n] %#ok<NOPRT,NASGU>
0058         weight_before_iteration = weights(1:20)' %#ok<NOPRT,NASGU>
0059     end
0060     
0061     b=zeros(M+N-2,1);
0062     % set A and b for observed points
0063     for m = 1:M
0064         w = 1/sigma_n*sqrt(weights(m));
0065         A(m,ts(m)) = w;
0066         dl      = ys(m)-xa(ts(m));
0067         b(m)    = dl*w;
0068     end
0069     
0070     kk = M;
0071     % set A and b for regularizing observations
0072     for n = 2:N-1
0073         kk = kk+1;
0074         w = 1/sigma_e*sqrt(weights(kk));
0075         A(kk,n-1:n+1)= [1 -2 1]*w;
0076         dl           = -(xa(n+1)-2*xa(n)+xa(n-1));
0077         b(kk)        = dl*w;
0078     end
0079     
0080     % normal equation system
0081     Nmatrix = A'*A;
0082     nvector = A'*b;
0083     if N <= 1000
0084         Cov = inv(Nmatrix);
0085     end
0086     
0087     % solution of correction
0088     dx = Nmatrix\nvector;
0089     % corrected estimate
0090     xa = xa+dx;
0091     
0092 %     % ----------------
0093 %     % apply conjugate gradients
0094 %     reorth=0;
0095 %     k=2*N;
0096 %     start_time=cputime
0097 %     [dx_cg,rho,eta,F] = cgls_mod(A,b,k,reorth);
0098 %
0099 %     check_conjugate_gradients = norm(dx-dx_cg)
0100 %         time_for_solution=cputime-start_time
0101 %     % -----------------------
0102 
0103     if iter == Niter+1 && print_type > 0
0104         xa_dx = [xa';dx'] %#ok<NOPRT,NASGU>
0105     end
0106     v = A*dx-b;                             % estimated standardized residuals
0107     
0108     %eso        = norm(v)/sqrt(M-2);              % estimated sigma0
0109     eso_robust = median(abs(v()))*1.48;
0110     
0111     vori=v;
0112     vori(1:M)      = v(1:M)*sigma_n./sqrt(weights(1:M));
0113     vori(M+1:end)  = v(M+1:end)*sigma_e./sqrt(weights(M+1:end));
0114     vs    = vori;
0115     
0116     %eso_n = median(abs(v(1:M)))*1.48;      % estimated RMS of point residuals
0117     eso_n = median(abs(vs(1:M)))*1.48;      % estimated RMS of point residuals
0118     
0119     %eso_e = median(abs(v(M+1:end)))*1.48;  % estimated RMS curvature residuals
0120     eso_e = median(abs(vs(M+1:end)))*1.48;  % estimated RMS curvature residuals
0121     
0122     if print_type > 0
0123         estimated_s0_s0p_s0k=[eso,eso_p,eso_k] %#ok<NOPRT,NASGU>
0124     end
0125 
0126     % critical values
0127     kn = k*eso_n;
0128     ke = k*eso_e;
0129 
0130 %% adapt weights
0131     if iter < Niter
0132         if type_out == 0                            % symmetric
0133             if print_type > 0
0134                 disp('rob_p_d')
0135             end
0136             if iter < 4                             % first iterations
0137                 if rob_p_d == 1 || rob_p_d == 3     % dem
0138                     weights(M+1:M+N-2) = ...
0139                         min(1,1./(abs(vs(M+1:M+N-2))/ke+0.0001));
0140 %                     weights(M+1:M+N-2) = ...
0141 %                         min(1,1./(abs(v(M+1:M+N-2))/ke+0.0001));
0142                    
0143                     if print_type > 1
0144                         iter_ke_max_v_min_w = ...
0145                             [iter,ke,max(abs(vs(M+1:M+N-2))),...
0146                             min(weights(M+1:M+N-2))]; %#ok<NASGU>
0147                         disp('iter_ke_max_v_min_w')
0148                         disp('modify dem weights')
0149                         residuals_new_weights_of_dem=...
0150                             [vs(M+1:M+N-2)';weights(M+1:M+N-2)']; %#ok<NASGU>
0151                         disp('residuals_new_weights_of_dem')
0152                     end
0153                     if plot_type > 0
0154                         subplot(3,1,1)
0155                         plot(2:N-1,xa(2:N-1),'-b')
0156                         title(num2str(sigma_e/sigma_n));
0157                         subplot(3,1,2);
0158                         hold on
0159                         plot(2:N-1,vs(M+1:M+N-2),'-b')
0160                         plot(2:N-1,vori(M+1:M+N-2),'-r')
0161                         plot(2:N-1,ones(N-2,1)*ke,'-b')
0162                         plot(2:N-1,-ones(N-2,1)*ke,'-b')
0163                         title(num2str(ke))
0164                         subplot(3,1,3);
0165                         plot(2:N-1,weights(M+1:M+N-2),'-b')
0166                         title(num2str(iter))
0167                     end
0168                 end
0169                 if rob_p_d == 2 || rob_p_d == 3     % points
0170                     weights(1:M) = min(1,1./(abs(vs(1:M))/kn+0.0001));
0171                     %weights(1:M) = min(1,1./(abs(v(1:M))/kn+0.0001));
0172                     
0173                     
0174                     if print_type > 1
0175                         residuals_new_weights_of_points_L1 = ...
0176                             [(1:20);v(1:20)';vs(1:20)';weights(1:20)'] %#ok<NOPRT,NASGU>
0177                     end
0178                 end               
0179             else                                    % iteration 4 following
0180                 if rob_p_d == 1 || rob_p_d == 3     % dem
0181                     weights(M+1:M+N-2) = exp(-abs(vs(M+1:M+N-2)).^2./ke^2/2)+0.0001;
0182                     %weights(M+1:M+N-2) = exp(-abs(v(M+1:M+N-2)).^2./ke^2/2)+0.0001;
0183                     
0184                     
0185                     if print_type > 1
0186                         iter_kk_max_v_min_w = ...
0187                             [iter,ke,max(abs(vs(M+1:M+N-2))),min(weights(M+1:M+N-2))] %#ok<NOPRT,NASGU>
0188                         disp('modify dem weights')
0189                         residuals_new_weights_of_dem = ...
0190                             [vs(M+1:M+N-2)';weights(M+1:M+N-2)'] %#ok<NOPRT,NASGU>
0191                     end
0192                     %
0193                     if plot_type > 0
0194                         subplot(3,1,1)
0195                         plot(2:N-1,xa(2:N-1),'-b')
0196                         title(num2str(sigma_e/sigma_n));
0197                         subplot(3,1,2);
0198                         hold on
0199                         plot(2:N-1,vs(M+1:M+N-2),'-b')
0200                         plot(2:N-1,vori(M+1:M+N-2),'-r')
0201                         plot(2:N-1,ones(N-2,1)*ke,'-b')
0202                         plot(2:N-1,-ones(N-2,1)*ke,'-b')
0203                         title(num2str(ke))
0204                         subplot(3,1,3);
0205                         plot(2:N-1,weights(M+1:M+N-2),'-b')
0206                         title(num2str(iter))
0207                     end
0208                 end
0209                 if rob_p_d == 2 || rob_p_d == 3         %points
0210                     weights(1:M) = exp(-abs(vs(1:M)).^2/kn^2/2)+0.0001;
0211                     %weights(1:M) = exp(-abs(v(1:M)).^2/kn^2/2)+0.0001;
0212                     
0213                     
0214                     if print_type > 1
0215                         residuals_new_weights_of_points_exp = ...
0216                             [(1:20);v(1:20)';vs(1:20)'/kn;weights(1:20)'] %#ok<NOPRT,NASGU>
0217                     end
0218                 end                
0219             end
0220         else  %  aymmetric
0221             k = 1;
0222             g = g_factor*k;
0223             w = w_factor*k;
0224             if type_robust==0   % asymmetric L1
0225                 weights(1:M) = max(min(1,1./abs(vs(1:M))/k+0.0001),((sign(vs(1:M))/k)+1)/2);
0226             else                % Kraus/Pfeifer
0227                 for m=1:M
0228                     if vs(m) < g-w
0229                         weights(m) = 0;
0230                     elseif vs(m) > g
0231                         weights(m)=1;
0232                     else
0233                         weights(m)=1/(1+(vs(m)-g)^4);
0234                     end
0235                 end
0236             end
0237         end
0238     else
0239 
0240         weights = 0.9999*ceil(weights-0.5)+0.0001;
0241         % initiate last iteration (Niter+1) with approximate values 0
0242         if iter == Niter
0243             xa = zeros(N,1);
0244         end
0245     end
0246 end
0247 
0248 xest = xa;
0249 
0250 
0251
estimate_profile_robust

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
