% eda03_18
% Illustrate error propagation from data
% to model parameters

% set up vectors d1 and d2
L=40;
Dd = 1.0;
d1 = Dd*[0:L-1]';
d2 = Dd*[0:L-1]';

% make an uncorrelated normal distribution p(d1,d2)
% with equal variance
d1bar=15;
d2bar=25;
dbar = [d1bar, d2bar]';
Cd=[[5,0]',[0,5]'];
CdI=inv(Cd);
norm=1/(2*pi*sqrt(det(Cd)));
dd1=d1-d1bar;
dd2=d2-d2bar;
% note that a quadratic form v'Mv can be written
% M(1,1)*v(1)^2 + M(2,2)*v2(2)^2 + 2*M(1,2)*v(1)*v(2)
% M is a 2x2 symmetric matrix and v is a column-vector
Pd=norm*exp(-0.5*CdI(1,1)*(dd1.^2)*ones(L,1)' ...
    -0.5*CdI(2,2)*ones(L,1)*(dd2.^2)' ...
    -CdI(1,2)*dd1*dd2');

% model is related to data by m=Md
M = [ [-1, -2]', [1, 1]' ];

% error propagation
mbar = M*dbar;
m1bar = mbar(1);
m2bar = mbar(2);
Cm = M*Cd*M';

% print results
disp(sprintf('mbar1=%f mbar2=%f',m1bar,m2bar));
disp(sprintf('varm1=%f varm2=%f covm=%f', Cm(1,1), Cm(2,2), Cm(1,2) ));

% evaluate normal distribution p(m1,m2)
% with equal varian
Dm = 1.0;
m1 = Dd*[-L/2:L/2-1]';  % note change in range of m1 and m2
m2 = Dd*[-L/2:L/2-1]';
CmI=inv(Cm);
norm=1/(2*pi*sqrt(det(Cm)));
dm1=m1-m1bar;
dm2=m2-m2bar;
% note that a quadratic form v'Mv can be written
% M(1,1)*v(1)^2 + M(2,2)*v2(2)^2 + 2*M(1,2)*v(1)*v(2)
% M is a 2x2 symmetric matrix and v is a column-vector
Pm=norm*exp(-0.5*CmI(1,1)*(dm1.^2)*ones(L,1)' ...
    -0.5*CmI(2,2)*ones(L,1)*(dm2.^2)' ...
    -CmI(1,2)*dm1*dm2');

% use simple drawing function that encapsulates all the graphics
eda_draw(Pd, 'caption p(d1,d2)', Pm, 'caption p(m1,m2)');
