0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
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
0025 type_robust=type_r(1);
0026 g_factor =type_r(2);
0027 w_factor =type_r(3);
0028
0029
0030 M = length(ts);
0031 A = zeros(M+N-2,N);
0032
0033
0034 xa = zeros(N,1);
0035 weights = ones(M+N-2,1);
0036 Cov = 0;
0037
0038
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
0052
0053
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
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
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
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