% Jyvaskylä Summer School 2014
% Aku Seppanen
% University of Eastern Finland, Kuopio
% Department of Applied Physics

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example 2.1: pair of linear equations %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

a = 1e-6;
f_true = [1 1]'; % True parameter values
K = [0 1; a 1];
g_true = K*f_true;
var_n = 1e-6; % variance of the observation noise

% 100 noise realizations
for j = 1:100
   n = sqrt(var_n)*randn(2,1);
   g_noisy = g_true + n;
   g_noisy_all(:,j) = g_noisy;
   f_est = inv(K)*g_noisy;
   f_est_all(:,j) = f_est;
end

figure(1), clf, plot(f_true(1),f_true(2),'bo'), hold on
plot(f_est_all(1,:),f_est_all(2,:),'r+')
xlabel('f_1'); ylabel('f_2');
legend('f_{true}','estimates')
title('Estimates for f corresponding to 100 noise realizations')

axis equal

figure(2), clf, plot(g_true(1),g_true(2),'bo'), hold on
plot(g_noisy_all(1,:),g_noisy_all(2,:),'r+')
xlabel('g_1'); ylabel('g_2');
legend('g_{true}','g^{\delta}')
title('Noisy data, 100 realizations')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example 2.2: modeling error %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

a = 1e-6;
f_true = [1 1]'; % True parameter values
K_accurate = [0 .999; a 1];
K_approximate = [0 1; a 1];
g_true = K_accurate*f_true;
var_n = 1e-10; % variance of the observation noise
% Note: the noise level is extremely low

% 100 noise realizations
for j = 1:100
   n = sqrt(var_n)*randn(2,1);
   g_noisy = g_true + n;
   g_noisy_all(:,j) = g_noisy;
   f_est = inv(K_approximate)*g_noisy;
   f_est_all(:,j) = f_est;
end

figure(1), clf, plot(f_true(1),f_true(2),'bo'), hold on
plot(f_est_all(1,:),f_est_all(2,:),'r+')
xlabel('f_1'); ylabel('f_2');
legend('f_{true}','estimates')
title('Estimates for f corresponding to 100 noise realizations')

axis equal

figure(2), clf, plot(g_true(1),g_true(2),'bo'), hold on
plot(g_noisy_all(1,:),g_noisy_all(2,:),'r+')
xlabel('g_1'); ylabel('g_2');
legend('g_{true}','g^{\delta}')
title('Noisy data, 100 realizations')






