fig_16_13_profile_reconstruction_demo_precision

% Fig 16.13 page 750
 profile recondstruction with precision

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

0001 %% Fig 16.13 page 750
0002 % profile recondstruction with precision
0003 %
0004 % Wolfgang Förstner 2014-08-06
0005 % last changes: Susanne Wenzel 09/16
0006 % wfoerstn@uni-bonn.de, wenzel@igg.uni-bonn.de
0007 
0008 function fig_16_13_profile_reconstruction_demo_precision()
0009 
0010 addpath(genpath('../../General-Functions/'));
0011 addpath('Functions')
0012 
0013 close all
0014 
0015 %% initialize random number generation by fixed seed
0016 init_rand_seed(23);
0017 
0018 %% plot settings
0019 ss = plot_init;
0020 
0021 %% set parameters
0022 N = 200;                    % number of grid points
0023 sigma_e=0.5;                % process noise
0024 sigma_n=0.5;                % observation noise
0025 factor_sigma = 1;           % factor for sigma_n
0026 type_outlier = 0;           % symmetric
0027 type_robust  = [0,0,0,0];   % not robust
0028 Niter = 0;                  % not robust
0029 print_type = 0;
0030 plot_type  = 0;
0031 
0032 %% generate profile
0033 [x,y,select,xs,ys] = ...
0034     generate_observed_AR2_demo(N,sigma_e,sigma_n);
0035 
0036 %% reconstruct profile
0037 [xest,aa,ab,ac,Cov] = estimate_profile_robust...
0038     (N,select,ys,sigma_e/factor_sigma,sigma_n,Niter,...
0039     type_outlier,type_robust,...
0040     plot_type,print_type);
0041 
0042 %% show precision
0043 
0044 factor = 1;               % blow-up factor for error band
0045 
0046 figure('name','Fig. 16.13: Profile Reconstruction and precision','color','w',...
0047     'Position',[0.2*ss(1),0.2*ss(2),0.5*ss(1),0.5*ss(2)]);
0048 hold on;
0049 % for n=1:N-1
0050 %     plotrect([n,x(n)-3*factor*sqrt((Cov(n,n)))],...
0051 %              [n+1,x(n+1)+3*factor*sqrt((Cov(n,n)))],'g' );
0052 % end
0053 xbandtop = x+3*factor*sqrt(diag(Cov));
0054 xbandbottom = x-3*factor*sqrt(diag(Cov));
0055 fill([(1:N)';(N:-1:1)'],[xbandtop;xbandbottom(end:-1:1)],[0.3,0.97,0.02],'FaceAlpha',0.5,'EdgeAlpha',0);
0056 plot(1:N,xest,'-k','LineWidth',2)
0057 plot(1:N,x,'--r','LineWidth',2)
0058 plot(1:N,xbandtop,'-r','LineWidth',1)
0059 plot(1:N,xbandbottom,'-r','LineWidth',1)
0060 plot(select,ys,'.b','MarkerSize',15)
0061 plot(1:N,xest,'.k','LineWidth',2)
0062 title(['Fig 16.13: $N = ',num2str(N),'$, $s_e = ',num2str(sigma_e),'$, $s_n =',num2str(sigma_e),'$'])
0063 
0064 
0065 return
0066 
0067 
0068 %% generate_observed_AR2
0069 %
0070 % N         = number of points
0071 % sigma_e   = process noise
0072 % sigma_n   = observation noise
0073 %
0074 % x,y       = true signal
0075 % select    = indices for observed points
0076 % xy,ys     = observed points
0077 
0078 function [x,y,select,xs,ys] = ...
0079     generate_observed_AR2_demo(N,sigma_e,sigma_n)
0080 
0081 x=zeros(N,1);
0082 y=zeros(N,1);
0083 
0084 for n = 3:N
0085     x(n) = 1.9998*x(n-1)-0.9999*x(n-2)+randn(1)*sigma_e;
0086     y(n) = x(n)+randn(1)*sigma_n;
0087 end
0088 for i=1:N
0089     y(i) = y(i) - (i-1)*(x(N)-x(1))/(N-1);
0090     x(i) = x(i) - (i-1)*(x(N)-x(1))/(N-1);
0091 end
0092 
0093 M = [6,12,24,81,95,124,138,176];
0094 select = zeros(8,1);
0095 xs = zeros(8,1);
0096 ys = zeros(8,1);
0097 for m = 1:8
0098     select(m) = M(m);
0099     xs(m) = x(M(m));
0100     ys(m) = y(M(m));
0101 end
0102 
0103 return