demo_GHM_2D_line

% Demo nonlinear GHM without constraints

 Wolfgang Förstner 6/6/2016
 wfoerstn@uni-bonn.de

 g = ^l2 - a ^l1 -b = 0 straight line  y_i - (a x_i + b) = 0

0001 %% Demo nonlinear GHM without constraints
0002 %
0003 % Wolfgang Förstner 6/6/2016
0004 % wfoerstn@uni-bonn.de
0005 %
0006 % g = ^l2 - a ^l1 -b = 0 straight line  y_i - (a x_i + b) = 0
0007 %
0008 clearvars
0009 close all
0010 
0011 addpath(genpath('../../../General-Functions/'));
0012 addpath(genpath('../Functions-GHM/'));
0013 
0014 % initialization of random numbers
0015 % = 0 CPU-time dependent
0016 % > 0 fixed
0017 %init_rand  = 0;
0018 init_rand   = 15;     % standard = 15
0019 init_rand = init_rand_seed(init_rand);
0020 
0021 disp(' =============================================')
0022 disp(' ----------  Demo GHM fitting line -----------')
0023 disp(' ----------- y_i - (a x_i + b) = 0 -----------')
0024 disp(' ---------------------------------------------')
0025 disp(['Seed for random numbers                              = ' ,num2str(init_rand)]);
0026 
0027 % choose type of covariance matrix
0028 % cov_type = 0; % homogeneous, isotropic, uncorrelated
0029 % cov_type = 1; % homogeneous, anisotropic, uncorrelated
0030 cov_type = 2; % inhomogeneous, anisotropic, uncorrelated
0031 %cov_type = 3; % random, correlated
0032 
0033 % choose type of approximate values
0034 xa = [0.8,1.0]';
0035 % appr_type=0;
0036 appr_type = 1;
0037 
0038 %% generate data
0039 
0040 % line parameters
0041 a = 1;
0042 b = 2;
0043 disp(['true parameters  [a, b]                              = [', num2str(a),',',num2str(b),']'])
0044 
0045 % number of observations
0046 N = 4; .... cannot be changed
0047     
0048 % true observations
0049 l1_t  = randn(N,1);
0050 l2_t  = a*l1_t + b;
0051 % covariance matrix
0052 
0053 switch cov_type
0054     case 0 % l1 and l2 independent, isotropic
0055         sigma_1 = 0.03;
0056         Cov_ll = [sigma_1^2*eye(N), zeros(N);zeros(N),sigma_1^2*eye(N)];
0057         % standard deviations
0058         disp(['Standard deviations for x_i and y_i                  =  ', num2str(sigma_1)])
0059         
0060     case 1
0061         % l1 and l2 independent, isotropic
0062         % standard deviations
0063         sigma_1 = 0.02;
0064         sigma_2 = 0.05;
0065         Cov_ll = [sigma_1^2*eye(N), zeros(N);zeros(N),sigma_2^2*eye(N)];
0066         disp(['Standard deviations for x_i and y_i                  =  ', num2str(sigma_1),',',num2str(sigma_2)])
0067         
0068     case 2
0069         % l1 and l2 independend, non isotropic
0070         var=[1,2,3,4,5,6,7,8];
0071         sigma_1=0.02;
0072         Cov_ll = diag(var)*sigma_1^2;
0073         disp(['Standard deviations for observations                 =  ', num2str(sqrt(var))])
0074     case 3
0075         % fully correlated
0076         sigma_1=0.01;
0077         A = randn(8);
0078         Cov_ll=(3*eye(8)+A*A')*sigma_1^2;
0079         disp('Standard deviations for observations                 =  random')
0080 end
0081 
0082 % all observations
0083 lvt = [l1_t;l2_t];
0084 % noisy observations
0085 lv = rand_gauss(lvt,Cov_ll,1);
0086 l1 = lv(1:N);
0087 l2 = lv(N+1:2*N);
0088 
0089 
0090 ScrS = plot_init;
0091 figure('Color','w','Position',[100 100  ScrS+[ -600 -300]])
0092 hold on
0093 plot(l1,l2,'.r','MarkerSize',5)
0094 for n=1:N
0095     ell.cov  = Cov_ll([n,n+4],[n,n+4]);
0096     ell.mean = [l1(n);l2(n)];
0097     plot_ellipse(ell,'-b');
0098 end
0099 title('Given points with standard ellipses, fitted line')
0100 axis equal
0101 
0102 r_U     = 1;                         % we are interested in the slope only
0103 %% determine approximate values
0104 if appr_type ==1 % estimate approximate values
0105     x0    = mean([l1,l2])';
0106     C0    = cov([l1,l2]);
0107     phi   = 0.5*atan2(2*C0(1,2),C0(2,2)-C0(1,1));
0108     xa(1) = tan(phi);
0109     xa(2) = x0(2) - x0(1) * xa(1); % xa(1) = tan(phi) = (x0(2)-xa(2))/(x0(1)-0)
0110 end
0111 disp(['approximate values                                   = [', num2str(xa(1)),',',num2str(xa(2)),']'])
0112 
0113 %% estimate
0114 
0115 sx      = [sigma_1,sigma_1]/sqrt(N); % approxmate standard deviations
0116 Tx      = 0.02;                     % theshold for convergence
0117 maxiter = 10;                        % maximum number of iterations
0118 
0119 [est_x,Cov_xx,sigma_0q,R,vv,zv,riv,nabla_lv,muv,muv1,uv1q,uv2]...
0120     = GaussHelmertModel(lv,Cov_ll,@cg_2D_line,xa,sx,Tx,maxiter,r_U);
0121 
0122 %% provide result
0123 disp(' ')
0124 disp(['estimated parameters  [a, b]   = [', num2str(est_x(1)),',',num2str(est_x(2)),']'])
0125 disp(['estimated sigms_0              =  ', num2str(sqrt(sigma_0q))])
0126 disp(['Covariance matrix              = [',num2str(Cov_xx(1,:)),']'])
0127 disp(['                               = [',num2str(Cov_xx(2,:)),']'])
0128 disp(['Standard deviations for [a, b] =  ', num2str(sqrt(Cov_xx(1,1))),', ',num2str(sqrt(Cov_xx(2,2)))])
0129 
0130 xmin=min(l1)-0.2;
0131 xmax=max(l1)+0.2;
0132 ymin=min(l2);
0133 ymax=max(l2);
0134 
0135 plot([xmin,xmax],[est_x(1)*xmin+est_x(2),est_x(1)*xmax+est_x(2)],'-r','LineWidth',2);
0136 xlim([xmin-0.2,xmax+0.2])
0137 ylim([ymin-0.2,ymax+0.2])