% eda03_13
% Illustrate calculating the variance of
% a joint distribution p(d1,d2)

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

% make a normal distribution
d1bar=15;
d2bar=25;
s1=3;
s2=5;
norm=1/(2*pi*s1*s2);
p1=exp(-((d1-d1bar).^2)/(2*s1*s1));
p2=exp(-((d2-d2bar).^2)/(2*s2*s2));
P=norm*p1*p2';

% sum along columns, which integrates P along d2 to get p1=p(d1)
p1 = Dd*sum(P,2);
% sum along rows, which integrates P along d1 to get p2=p(d2)
p2 = Dd*sum(P,1)';

% now compute one-dimensional means
d1mean = Dd*sum(d1 .* p1 );
d2mean = Dd*sum(d2 .* p2 );

% now compute one-dimensional variances
sigma12 = Dd*sum( ((d1-d1mean).^2) .* p1 );
sigma22 = Dd*sum( ((d2-d2mean).^2) .* p2 );

disp(sprintf('sigma1=%f sigma2=%f', sqrt(sigma12), sqrt(sigma22) ));

