generate_observed_AR2

% generate_observed_AR2 with noise and ouliers

 [x,y,out_in,select,xs,ys] = generate_observed_AR2(N, sigma_e, sigma_n, Pout, Mout, type_out, dens)

 N         = number of points
 sigma_e   = process noise
 sigma_n   = observation noise
 Pout      = probability for outliers
 Mout      = maximal outlier
 type_out  = type otliers (0=symm., 1=asymm.)
 dens      = percentage of observed points in (0,1]

 x         = true profile with slope 0
 y         = noisy profile with slope 0 (no outliers)
 out_in    = Nx1 vector, 0/1: o inlier, 1 outlier
 select    = index of observed points
 xs        = observed true profile
 ys        = observed noisy profile with outliers

 Wolfgang Förstner
 wfoerstn@uni-bonn.de

0001 %% generate_observed_AR2 with noise and ouliers
0002 %
0003 % [x,y,out_in,select,xs,ys] = generate_observed_AR2(N, sigma_e, sigma_n, Pout, Mout, type_out, dens)
0004 %
0005 % N         = number of points
0006 % sigma_e   = process noise
0007 % sigma_n   = observation noise
0008 % Pout      = probability for outliers
0009 % Mout      = maximal outlier
0010 % type_out  = type otliers (0=symm., 1=asymm.)
0011 % dens      = percentage of observed points in (0,1]
0012 %
0013 % x         = true profile with slope 0
0014 % y         = noisy profile with slope 0 (no outliers)
0015 % out_in    = Nx1 vector, 0/1: o inlier, 1 outlier
0016 % select    = index of observed points
0017 % xs        = observed true profile
0018 % ys        = observed noisy profile with outliers
0019 %
0020 % Wolfgang Förstner
0021 % wfoerstn@uni-bonn.de
0022 
0023 function [x, y, out_in, select, xs, ys] = ...
0024     generate_observed_AR2(N, sigma_e, sigma_n, Pout, Mout, type_out, dens)
0025 
0026 % initiate true profiles
0027 x = zeros(N, 1);
0028 y = zeros(N, 1);
0029 
0030 % generate true and noisy profile
0031 for n = 3:N
0032     x(n) = 1.9998 * x(n - 1) - 0.9999 * x(n - 2) + randn(1) * sigma_e;
0033     y(n) = x(n) + randn(1) * sigma_n;
0034 end
0035 
0036 % enforce slope 0
0037 for i = 1:N
0038     y(i) = y(i) - (i - 1) * (x(N) - x(1)) / (N - 1);
0039     x(i) = x(i) - (i - 1) * (x(N) - x(1)) / (N - 1);
0040 end
0041 % add ouliers and select observations
0042 out_in = zeros(N, 1);
0043 m = 0;
0044 for i = 1:N
0045     % add outliers
0046     if rand(1) < Pout
0047         out_in(i) = 1;
0048         if type_out == 0 % symmetric in [-Mout,+Mout]
0049             y(i) = y(i) + 2 * (rand(1) - 0.5) * Mout;
0050         else % asymmetric in [0,Mout]
0051             y(i) = y(i) + rand(1) * Mout;
0052         end
0053     end
0054     
0055     % select observations
0056     if rand(1) < dens
0057         m = m + 1;
0058         select(m) = i;                                                      %#ok<*AGROW>
0059         xs(m) = x(i);
0060         ys(m) = y(i);
0061     end
0062 end
0063 
0064 
0065