Home > .. > 16-Profilereconstruction > Functions > estimate_profile_robust_flat.m

estimate_profile_robust_flat

PURPOSE ^

% estimate_profile_robust_flat

SYNOPSIS ^

function [xest,A,v,weights,Cov] = estimate_profile_robust_flat(N,ts,ys,sigma_e,sigma_n,Niter,type_out,type_r,print_type)

DESCRIPTION ^

% 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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         display(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 %     if iter == Niter+1
0052 %         rhs_w=[b(1:M)';weights(1:M)'];
0053 %     end
0054 
0055     k = M;
0056     for n = 1:N-1
0057         k = k+1;
0058         A(k,n:n+1)=[1 -1]/sigma_e;
0059     end
0060     
0061     % normal system
0062     Nmatrix = A'*A;
0063     nvector = A'*b;
0064     if N <= 1000 
0065         Cov = inv(Nmatrix);
0066     end
0067     dx = Nmatrix\nvector;
0068     xa = xa+dx;
0069     
0070     if iter == Niter+1 && print_type > 0
0071         xa_dx = [xa';dx'];
0072         display(xa_dx)
0073     end
0074     v = A*dx-b;
0075     
0076     %% adapt weights
0077     if iter < Niter
0078         
0079         if type_out == 0
0080             
0081             k=1;
0082             if iter < 6
0083                 weights(1:M) = min(1,1./abs(v(1:M))/k+0.0001);
0084             else
0085                 weights(1:M) = exp(-abs(v(1:M))/k);
0086             end
0087             
0088         else  % Kraus/Pfeifer
0089             
0090             k = 1;
0091             g = g_factor*k;
0092             w = w_factor*k;
0093             
0094             if type_robust == 0
0095                 
0096                 weights(1:M) = ...
0097                     max( min(1,1./abs(v(1:M))/k+0.0001),...
0098                     ((sign(v(1:M))/k)+1)/2 );
0099             else
0100                 
0101                 for m = 1:M
0102                     
0103                     if v(m) < g-w
0104                         weights(m) = 0;
0105                         
0106                     elseif v(m) > k
0107                         weights(m)=1;
0108                         
0109                     else
0110                         weights(m)=1/(1+(v(m)-g)^4);
0111                         
0112                     end
0113                 end
0114             end            
0115         end
0116     else
0117         weights(1:M) = weights(1:M)>0.5;
0118         if iter == Niter
0119             xa  = zeros(N,1);
0120         end
0121     end
0122     
0123 end
0124 
0125 xest = xa;
0126 
0127

Generated on Sat 17-Sep-2016 20:58:55 by m2html © 2005