estimate_profile_robust_flat

% estimate_profile_robust_flat

 xest = estimate_profile_robust(N,ts,ys,sigma_e,sigma_n,
                    Niter,type_out,type_r,print_type)

 N        = scalar, length of profile
 ts,ys    = Mx1 vectors, observed profile
 sigma_e  = std of curvature
 sigma_n  = std of observations
 Niter    = number of iterations
 type_out = type of outlier (0=symm, 1=asymm)
 type_r   (1) > 0 = L1, 1 = Kraus 
          (2) = g for Kraus
          (3) = w for Kraus
 print_type   = 0,1,2 (no, low, high)

 Wolfgang Förstner 2014-10-05
 last changes: Susanne Wenzel 09/16
 wfoerstn@uni-bonn.de, wenzel@igg.uni-bonn.de

0001 %% estimate_profile_robust_flat
0002 %
0003 % xest = estimate_profile_robust(N,ts,ys,sigma_e,sigma_n,
0004 %                    Niter,type_out,type_r,print_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 % type_out = type of outlier (0=symm, 1=asymm)
0012 % type_r   (1) > 0 = L1, 1 = Kraus
0013 %          (2) = g for Kraus
0014 %          (3) = w for Kraus
0015 % print_type   = 0,1,2 (no, low, high)
0016 %
0017 % Wolfgang Förstner 2014-10-05
0018 % last changes: Susanne Wenzel 09/16
0019 % wfoerstn@uni-bonn.de, wenzel@igg.uni-bonn.de
0020 
0021 function [xest,A,v,weights,Cov] = estimate_profile_robust_flat...
0022     (N,ts,ys,sigma_e,sigma_n,Niter,type_out,type_r,print_type)
0023 
0024 % robust
0025 type_robust=type_r(1);
0026 g_factor   =type_r(2);
0027 w_factor   =type_r(3);
0028 
0029 % Number of observations
0030 M = length(ts);
0031 A = zeros(M+N-2,N);
0032 
0033 % initialize
0034 xa = zeros(N,1);
0035 weights = ones(M+N-2,1);
0036 Cov = 0;
0037 
0038 %% robust estimation
0039 for iter = 1:Niter+1
0040     if print_type > 0
0041         disp(iter)
0042     end
0043     
0044     b = zeros(M+N-1,1);
0045     for m = 1:M
0046         w = 1/sigma_n*sqrt(weights(m));
0047         A(m,ts(m)) = w;
0048         dl      = ys(m)-xa(ts(m));
0049         b(m)    = dl*w;
0050     end
0051 
0052     k = M;
0053     for n = 1:N-1
0054         k = k+1;
0055         A(k,n:n+1)=[1 -1]/sigma_e;
0056     end
0057     
0058     % normal system
0059     Nmatrix = A'*A;
0060     nvector = A'*b;
0061     if N <= 1000 
0062         Cov = inv(Nmatrix);
0063     end
0064     dx = Nmatrix\nvector;
0065     xa = xa+dx;
0066     
0067     if iter == Niter+1 && print_type > 0
0068         xa_dx = [xa';dx'];
0069         disp(xa_dx)
0070     end
0071     v = A*dx-b;
0072     
0073     %% adapt weights
0074     if iter < Niter
0075         
0076         if type_out == 0
0077             
0078             k=1;
0079             if iter < 6
0080                 weights(1:M) = min(1,1./abs(v(1:M))/k+0.0001);
0081             else
0082                 weights(1:M) = exp(-abs(v(1:M))/k);
0083             end
0084             
0085         else  % Kraus/Pfeifer
0086             
0087             k = 1;
0088             g = g_factor*k;
0089             w = w_factor*k;
0090             
0091             if type_robust == 0
0092                 
0093                 weights(1:M) = ...
0094                     max( min(1,1./abs(v(1:M))/k+0.0001),...
0095                     ((sign(v(1:M))/k)+1)/2 );
0096             else
0097                 
0098                 for m = 1:M
0099                     
0100                     if v(m) < g-w
0101                         weights(m) = 0;
0102                         
0103                     elseif v(m) > k
0104                         weights(m)=1;
0105                         
0106                     else
0107                         weights(m)=1/(1+(v(m)-g)^4);
0108                         
0109                     end
0110                 end
0111             end            
0112         end
0113     else
0114         weights(1:M) = weights(1:M)>0.5;
0115         if iter == Niter
0116             xa  = zeros(N,1);
0117         end
0118     end
0119     
0120 end
0121 
0122 xest = xa;
0123 
0124