% eda03_17
% Illustrate the computation of covariance

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

% make a positively correlated normal distribution
d1bar=15;
d2bar=25;
C=[[6,4]',[4,10]'];
CI=inv(C);
norm=1/(2*pi*sqrt(det(C)));
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
P=norm*exp(-0.5*CI(1,1)*(dd1.^2)*ones(L,1)' ...
    -0.5*CI(2,2)*ones(L,1)*(dd2.^2)' ...
    -CI(1,2)*dd1*dd2');

% make the alternating sign function
S = (d1-d1bar) * (d2-d2bar)';

% form the product
SP = S .* P;

% integrate
cov = (Dd^2)*sum(sum(SP));
disp(sprintf('covariance=%f', cov));

% use simple drawing function that encapsulates all the graphics
eda_draw(P, 'caption p', S, 'caption s', SP, 'caption product');


