Live Script edama_03
This Live Script supports Chapter 3 of Environmental Data Analysis with MATLAB or PYTHON, Third Edition, by WIlliam Menke
The script naming system has the form edama_CC_SS, where CC is the chapter number and SS is the script number.
% edama_03_01: mode calculated from p(d)
clear all; % clears previously-defined variables
% performing a clear at the top of a LiveSript is a good idea
% Illustrative distribution p(d)
Dd = 1.0; % sampling interval
d = Dd*[0:14]' + 0.5*Dd; % vector of evenly-spaced d's
p = [5.00, 15.00, 10.00, 10.00, 10.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00]';
% normalize fistribution so its integral is unity
% note element-wise division
p = p./(Dd*sum(p));
% mode is d for which p is maximum
% max returns the maximum of p
% and the index k at which p is maximum
[pmax, k] = max(p);
% the mode is the d with index k
dmode = d(k);
% output to command window
fprintf('mode = %.2f\n', dmode);
% plot the p.d.f. with the mode shown with a dotted line
figure(1);
clf;
set(gca,'LineWidth',2);
set(gca,'Fontsize',14);
hold on;
axis( [min(d), max(d), 0, max(p)]);
plot(d,p,'k-','LineWidth',2);
plot( [dmode, dmode], [0, max(p)/10], 'k:', 'LineWidth', 2);
xlabel('d');
ylabel('p(d)');
% edama_03_02: median calculated from p(d)
% Illustrative distribution p(d)
Dd = 1.0; % sampling interval
d = Dd*[0:14]' + 0.5*Dd; % vector of evenly-spaced d's
p = [5.00, 15.00, 10.00, 10.00, 10.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00]';
N = length(p);
% normalize fistribution so its integral is unity
p = p./(Dd*sum(p));
% cumulative probability distribution
P = Dd*cumsum(p);
% median is d at the 50% area point of P
% find first element of P that's equal to or greatet than 0.5
% so search for first occurence of P>=0.5
for i=[1:N]
if P(i) >= 0.5
dmedian = d(i);
break;
end
end
% print results
fprintf('median = %.2f', dmedian);
% plot the p.d.f. with the median shown with a dotted line
figure(1);
clf;
set(gca,'LineWidth',2);
set(gca,'Fontsize',14);
hold on;
axis( [min(d), max(d), 0, max(p)]);
plot(d,p,'k-','LineWidth',2);
plot( [dmedian, dmedian], [0, max(p)/10], 'k:', 'LineWidth', 2);
xlabel('d');
ylabel('p(d)');
% edama_03_03: mean calculated from p(d)
% Illustrative distribution p(d)
Dd = 1.0;
d = Dd*[0:14]' + 0.5*Dd;
p = [5.00, 15.00, 10.00, 10.00, 10.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00]';
% normalize fistribution so its integral is unity
% note element-wise division
p = p./(Dd*sum(p));
% perform the interal for the mean
dmean = Dd*sum(d.*p);
% print results
fprintf('mean = %.2f', dmean);
% plot the p.d.f. with the mean shown with a dotted line
figure(1);
clf;
set(gca,'LineWidth',2);
set(gca,'Fontsize',14);
hold on;
axis( [min(d), max(d), 0, max(p)]);
plot(d,p,'k-','LineWidth',2);
plot( [dmean, dmean], [0, max(p)/10], 'k:', 'LineWidth', 2);
xlabel('d');
ylabel('p(d)');
% eda03_04: variance calculated from p(d)
% Illustrative distribution p(d)
Dd = 1.0; % sampling interval
d = Dd*[0:14]' + 0.5*Dd; % vector of evenly spaced d's
p = [5.00, 15.00, 10.00, 10.00, 10.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00]';
% normalize distribution so its integral is unity
% note element-wise division
p = p./(Dd*sum(p));
% perform the integral for the mean
dmean = Dd*sum(d.*p);
% perform the integral for the variance
q = (d-dmean).^2;
sigma2 = Dd*sum(q.*p);
sigma = sqrt(sigma2);
% print results
fprintf('sigma2 = %.2f sigma = %.2f\n', sigma2, sigma);
% plot the p.d.f. with the mean and a solid line and +/- 1 sigma
% as dotted lines
figure(1);
clf;
set(gca,'LineWidth',2);
set(gca,'Fontsize',14);
hold on;
axis( [min(d), max(d), 0, max(p)]);
plot(d,p,'k-','LineWidth',2);
plot( [dmean, dmean], [0, max(p)/10], 'k-', 'LineWidth', 2);
plot( [dmean-sigma, dmean-sigma], [0, max(p)/10], 'k:', 'LineWidth', 2);
plot( [dmean+sigma, dmean+sigma], [0, max(p)/10], 'k:', 'LineWidth', 2);
xlabel('d');
ylabel('p(d)');
% edama_03_05: visualation of process of computing variance
% The plot had three vectors, p(d), q(d) and p(d)q(d)
% depicted as grey-shaded coun vectors, side-by-side
% with mean and +/- one sigma levels marked
% all vectors have 40 rows
N=40;
% data, d
Dd=1.0;
d = Dd*[0:N]';
% distribution, p(d)
p = sin(pi*d/d(N))./(d+10); % cooked up to be skewed
p = p/(Dd*sum(p)); % normalize
% mean
dbar = Dd*sum(d.*p);
% quadratic
q = (d-dbar).^2;
% variance
qp = q.*p;
sigma2 = Dd*sum(qp);
sigma = sqrt(sigma2)
% Set up graphics
figure(1);
clf;
set(gca,'LineWidth',2);
hold on;
axis([-N/4, 7*N/4, 0, 2*N]);
axis([-N/4, 7*N/4, 0, 1.2*N]);
axis equal;
axis off;
set(gca,'XTick',[]); % turn off horizontal axis
set(gca,'YTick',[]); % turn off vertical axis
% black and white color map
bw=0.9*(256-linspace(0,255,256)')/256;
colormap([bw,bw,bw]);
hold on;
% plot the vecror p
range=max(p)-min(p);
if( range == 0.0 )
range=1;
end
imagesc( [1, N/16], [N-1, 0], (p-min(p))/range);
text(N/16,-N/16,'p');
hold on;
% plot the vector q
range=max(q)-min(q);
if( range == 0.0 )
range=1;
end
imagesc( [1+N/4, N/16+N/4], [N-1, 0], (q-min(q))/range);
text(N/32+N/4,-N/16,'q');
hold on;
% plot the vecror qp
range=max(qp)-min(qp);
if( range == 0.0 )
range=1;
end
imagesc( [1+2*N/4, N/16+2*N/4], [N-1, 0], (qp-min(qp))/range);
text(N/32+2*N/4,-N/16,'pq');
hold on;
% plot mean
y = dbar;
plot( [-N/8, 6*N/8], [N-y/Dd, N-y/Dd], 'k-', 'LineWidth', 2 );
text( 6*N/8, N-y/Dd, ' xbar');
y=dbar-sigma;
plot( [-N/8, 6*N/8], [N-y/Dd, N-y/Dd], 'k-', 'LineWidth', 2 );
text( 6*N/8, N-y/Dd, ' xbar-sigma');
y=dbar+sigma;
plot( [-N/8, 6*N/8], [N-y/Dd, N-y/Dd], 'k-', 'LineWidth', 2 );
text( 6*N/8, N-y/Dd, ' xbar+sigma');
% edama_03_06: Normal p.d.f.'s with different means
% the all these column-vectors have 40 rows
L=40;
sigma=5;
d=[1:L]';
% First normal curve
dbar=10;
norm = 1/(sqrt(2*pi)*sigma);
a=norm*exp( -(d-dbar).^2 / (2*sigma^2) ) / (sigma*sqrt(2*pi));
% Second normal curve
dbar=15;
norm = 1/(sqrt(2*pi)*sigma);
b=norm*exp( -(d-dbar).^2 / (2*sigma^2) ) / (sigma*sqrt(2*pi));
% Third normal curve
dbar=20;
norm = 1/(sqrt(2*pi)*sigma);
c=norm*exp( -(d-dbar).^2 / (2*sigma^2) ) / (sigma*sqrt(2*pi));
% Fourth normal curve
dbar=25;
norm = 1/(sqrt(2*pi)*sigma);
dd=norm*exp( -(d-dbar).^2 / (2*sigma^2) ) / (sigma*sqrt(2*pi));
% Fifth normal curve
dbar=30;
norm = 1/(sqrt(2*pi)*sigma);
e=norm*exp( -(d-dbar).^2 / (2*sigma^2) ) / (sigma*sqrt(2*pi));
% use simple drawing function that encapsulates all the graphics
eda_draw(' ',a,'caption a',b,'caption b',c,'caption c',dd,'caption d',e,'caption e');
% edama_03_07: Normal p.s.f.s with different variances
% the all these column-vectors have 40 rows
L=40;
dbar=20;
d=[1:L]';
% First normal curve
sigma=2.5;
norm = 1/(sqrt(2*pi)*sigma);
a=norm*exp( -(d-dbar).^2 / (2*sigma^2) ) / (sigma*sqrt(2*pi));
% Second normal curve
sigma=5;
norm = 1/(sqrt(2*pi)*sigma);
b=norm*exp( -(d-dbar).^2 / (2*sigma^2) ) / (sigma*sqrt(2*pi));
% Third normal curve
sigma=10;
norm = 1/(sqrt(2*pi)*sigma);
c=norm*exp( -(d-dbar).^2 / (2*sigma^2) ) / (sigma*sqrt(2*pi));
% Fourth normal curve
sigma=20;
norm = 1/(sqrt(2*pi)*sigma);
dd=norm*exp( -(d-dbar).^2 / (2*sigma^2) ) / (sigma*sqrt(2*pi));
% Fifth normal curve
sigma=40;
norm = 1/(sqrt(2*pi)*sigma);
e=norm*exp( -(d-dbar).^2 / (2*sigma^2) ) / (sigma*sqrt(2*pi));
% use simple drawing function that encapsulates all the graphics
eda_draw(' ', a,'caption a',b,'caption b',c,'caption c',dd,'caption d',e,'caption e');
% edama_03_08: transformation of variables m = d^2
% compare original p.d.f. is uniform; tat is, p(d)=1 on 0<d<1
% the transformation of variables is m = d^2
% with transformed p.d.f. is p(m)=0.5*m^(-0.5) on 0<m<1
% setup a vector d
L=40;
Dd=1/L;
d=Dd*[1:L]';
% the distribution p(d) is constant
pd = ones(L,1);
% setup a vector m
L=40;
Dm=1/L;
m=Dm*[1:L]';
% m is a function of d
pm = 0.5*m.^(-0.5);
% use simple drawing function that encapsulates all the graphics
eda_draw(' ',pd,'caption p(d)',' ',pm,'caption p(m)');
% edama_03_09: The joint p.d.f. p(d1,d2) and the two univariate
% p.d.f.'s p(d1) and p(d2) formed from it
% set up the d1 and d2 axes
L=40;
Dd = 1.0;
d1 = Dd*[0:L-1]';
d2 = Dd*[0:L-1]';
% make a Normal distribution centered at d1bar, d2bar
d1bar=10;
d2bar=20;
s1=5;
s2=5;
P=exp(-((d1-d1bar).^2)/(2*s1*s1)) * exp(-((d2-d2bar).^2)/(2*s1*s1))';
P = P / (Dd*Dd*sum(sum(P))); % normalize
% sum along columna, which integrates P along d2 to get p(d1)
p1 = Dd*sum(P,2);
% sum along rows, which integrates p along d2 to get p(d2)
p2 = Dd*sum(P,1)';
% use simple drawing function that encapsulates all the graphics
eda_draw(P,'caption p(d1,d2)',' ',p1,'caption p(d1)',' ',p2,'caption p(d2)');
% edama_03_10: The uniform p.d.f. p(d1,d2)=constant
% set up the d1 and d2 axes
Nd1=40;
Nd2=40;
Dd1 = 1.0;
Dd2 = 1.0;
d1 = Dd1*[0:Nd1-1]';
d2 = Dd2*[0:Nd2-1]';
% make a uniform p.d.f. and normalize it
P = ones(Nd1,Nd2);
% calculate the area under the p.d.f.
norm = (Dd1*Dd2)*sum(sum(P));
% normalize the p.d.f. by dividing by the area
P = P/norm;
% use simple drawing function that encapsulates all the graphics
eda_draw(P,'caption p(d1,d2)');
% eda03_11: construct a joint Normal p.f.d.
% and compute its mean and variance
% set up d1 axis
Nd1 = 101;
d1 = 5.0*[0:Nd1-1]'/(Nd1-1);
Dd1 = d1(2)-d1(1);
% set up d2 axis
Nd2 = 101;
d2 = 5.0*[0:Nd2-1]'/(Nd2-1);
Dd2 = d2(2) - d2(1);
% 2D normal p.f.d. setup
dbar1 = 1.5;
dbar2 = 2.5;
sigma1 = 0.25;
sigma2 = 0.25;
% since d1 and d2 are uncorrelated, the 2-d Normal
% p.d.f. p(d1,d2) is product of p(d1) and p(d2)
norm1 = (1/(sqrt(2*pi)*sigma1)); % normalization of p(d1)
norm2 = (1/(sqrt(2*pi)*sigma2)); % normalization of p(d2)
p1 = norm1 * exp(-(d1-dbar1).^2/(2*sigma1*sigma1)); % p(d1)
p2 = norm2 * exp(-(d2-dbar2).^2 /(2*sigma2*sigma2)); % p(d2)
P = p1 * p2'; % p(d1,d2) = p(d1) p(d2)
% integrate along d2-axis to get p1
p_1 = Dd2*sum(P,2);
% integrate along d1-axis to get p2
p_2 = Dd1*sum(P,1)';
% areas (should all be 1) (approximate integrals as sums)
A = Dd1*Dd2*sum(sum(p));
A1 = Dd1*sum(p_1);
A2 = Dd2*sum(p_2);
% vectors of ones
w1 = ones(Nd1,1);
w2 = ones(Nd2,1);
% mean along 1 axis, computed from p(d1,d2
D1 = d1 * w2';
dbar_1 = Dd1*Dd2*sum(sum(D1.*P));
% mean along 1 axis, computed from p_1
dbar_1a = Dd1*sum( d1 .* p_1);
% mean along 2 axis, computed from p(d1,d2
D2 = w1 * d2';
dbar_2 = Dd1*Dd2*sum(sum(D2.*P));
% mean along 2 axis, computed from p_2
dbar_2a = Dd2*sum(d2 .*p_2);
% variance along 1 axis, computed from p(d1,d2)
% use a tensor product to create a 2D version of the quantity (d1-d1bar)**2
D1s = ((d1-dbar1).^2) * w2';
% compute the variance with the usual integral, approximated as a sum
sigma_1_squared = Dd1*Dd2*sum(sum(D1s.*P));
sigma_1 = sqrt(sigma_1_squared);
% variance along 1 axis, computed from p(d1) via the usual integral, approximated as a sum
sigma_1a_squared = Dd1*sum(((d1-dbar_1).^2) .* p_1);
sigma_1a = sqrt(sigma_1_squared);
% variance along 2 axis, computed from p(d1,d2)
% use a tensor product to create a 2D version of the quantity (d2-d2bar)**2
D2s = w1 * ((d2-dbar2).^2)';
% compute the variance with the usual integral, approximated as a sum
sigma2_squared = Dd1*Dd2*sum(sum(D2s.*P));
sigma_2 = sqrt(sigma2_squared);
% variance along 2 axis, computed from p(d2) via the usual integral, approximated as a sum
sigma_2a_squared = Dd2*sum( ((d2-dbar_2).^2) .* p_2);
sigma_2a = sqrt(sigma_2a_squared);
eda_draw(p,'caption p(d1,d2)');
fprintf('area %.4f %.4f %.4f\n', A, A1, A2);
fprintf('dbar1 %.4f %.4f %.4f\n', dbar1, dbar_1, dbar_1a);
fprintf('dbar2 %.4f %.4f %.4f\n', dbar2, dbar_2, dbar_2a);
fprintf('sigma1 %.4f %.4f %.4f\n', sigma1, sigma_1, sigma_1a);
fprintf('sigma2 %.4f %.4f %.4f\n', sigma2, sigma_2, sigma_2a);
% edama_03_12: example of a conditional p.d.f.
% set up d1 axis
Nd1 = 101;
d1 = 5.0*[0:Nd1-1]'/(Nd1-1);
Dd1 = d1(2)-d1(1);
% set up d2 axis
Nd2 = 101;
d2 = 5.0*[0:Nd2-1]'/(Nd2-1);
Dd2 = d2(2) - d2(1);
% 2D normal p.f.d. setup
dbar1 = 1.5;
dbar2 = 3.5;
sigma1 = 0.75;
sigma2 = 0.75;
% since d1 and d2 are uncorrelated, the 2-d Normal
% p.d.f. p(d1,d2) is product of p(d1) and p(d2)
norm1 = (1/(sqrt(2*pi)*sigma1)); % normalization of p(d1)
norm2 = (1/(sqrt(2*pi)*sigma2)); % normalization of p(d2)
p1 = norm1 * exp( -(d1-dbar1).^2 / (2*sigma1*sigma1) ); % p(d1)
p2 = norm2 * exp( -(d2-dbar2).^2 / (2*sigma2*sigma2) ); % p(d2)
P = p1 * p2'; % p(d1,d2) i= p(d1) p(d2)
% integrate p(d1,d2) along d2-axis to get p1
p_1 = Dd1*sum(P,2);
% integrate p(d1,d2) along d1-axis to get p2
p_2 = Dd2*sum(P,1)';
% conditional p.d.f. P1g2 = P(d1|d2) = P(d1,d2)/p2
p1g2 = P ./ (ones(Nd1,1)*p_2');
% conditional p.d.f. P2g1 = P(d2|d1) = P(d1,d2)/p1
p2g1 = P ./ (p_1*ones(Nd2,1)');
% use simple drawing function that encapsulates all the graphics
eda_draw(P,'caption p(d1,d2)',p1g2,'caption p(d1|d2)',p2g1,'caption p(d2|d1)');
% edama_03_13: covariance calculated from p.d.f. p(d1,d2)
% set up d1 axis
Nd1 = 101;
d1 = -2.5 + 5.0*[0:Nd1-1]'/(Nd1-1);
Dd1 = d1(2)-d1(1);
% set up d2 axis
Nd2 = 101;
d2 = -2.5 + 5.0*[0:Nd2-1]'/(Nd2-1);
Dd2 = d2(2) - d2(1);
% create a highly and positively correlated p.d.f.
P = zeros(Nd1,Nd2);
for i=[1:Nd1]
for j=[1:Nd2]
d_lo = -1.0+d1(i);
d_hi = 1.0+d1(i);
r2 = 2.0 + d1(i)^2 + d2(j)^2;
P(i,j) = (d2(j)>=d_lo)*(d2(j)<=d_hi)/r2;
end
end
norm = Dd1*Dd2*sum(sum(P));
P = P/norm;
% compute means
d1bar = Dd1*Dd2*sum(sum((d1*ones(1,Nd2)).*P));
d2bar = Dd1*Dd2*sum(sum((ones(Nd2,1)*d2').*P));
% make the alternating sign function
S = (d1-d1bar) * (d2-d2bar)';
% perform covariance integral
sigma12 = (Dd1*Dd2)*sum(sum(S .* P));
fprintf('d1bar %.3f d2bar %.3f covariance %.3f\n', d1bar, d2bar, sigma12 );
eda_draw(S,'caption s(d1,d2)', P,'caption p(d1,d2)', S.*P,'caption s(d1,d2) p(d1,d2)');
% edama_03_14: Constructing Normal p.d.f. without for loop
% set up d1 axis
Nd1 = 101;
d1 = 5.0*[0:Nd1-1]'/(Nd1-1);
Dd1 = d1(2)-d1(1);
% set up d2 axis
Nd2 = 101;
d2 = 5.0*[0:Nd2-1]'/(Nd2-1);
Dd2 = d2(2) - d2(1);
% make a positively correlated normal distribution
d1bar=1.5;
d2bar=2.5;
sigma1 = 0.25;
sigma2 = 0.25;
covar = 0.05;
C=[[sigma1^2,covar]',[covar,sigma2^2]'];
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(Nd2,1)' ...
-0.5*CI(2,2)*ones(Nd1,1)*(dd2.^2)' ...
-CI(1,2)*dd1*dd2');
pp=P;
% make a negatively correlated normal distribution
d1bar=1.5;
d2bar=2.5;
sigma1 = 0.25;
sigma2 = 0.25;
covar = -0.05;
C=[[sigma1^2,covar]',[covar,sigma2^2]'];
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
pn=norm*exp(-0.5*CI(1,1)*(dd1.^2)*ones(Nd2,1)' ...
-0.5*CI(2,2)*ones(Nd1,1)*(dd2.^2)' ...
-CI(1,2)*dd1*dd2');
% make an uncorrelated normal distribution
d1bar=1.5;
d2bar=2.5;
sigma1 = 0.25;
sigma2 = 0.25;
covar = 0.0;
C=[[sigma1^2,covar]',[covar,sigma2^2]'];
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
pu=norm*exp(-0.5*CI(1,1)*(dd1.^2)*ones(Nd2,1)' ...
-0.5*CI(2,2)*ones(Nd1,1)*(dd2.^2)' ...
-CI(1,2)*dd1*dd2');
% use simple drawing function that encapsulates all the graphics
eda_draw(pp, 'caption +cor', pn, 'caption -cor', pu, 'caption uncor');
% edama_03_15: Constructing Normal p.d.f. with for loops
% set up d1 axis
Nd1 = 101;
d1 = 5.0*[0:Nd1-1]'/(Nd1-1);
Dd1 = d1(2)-d1(1);
% set up d2 axis
Nd2 = 101;
d2 = 5.0*[0:Nd2-1]'/(Nd2-1);
Dd2 = d2(2) - d2(1);
% make a positively correlated normal distribution
d1bar=1.5;
d2bar=2.5;
sigma1 = 0.25;
sigma2 = 0.25;
covar = 0.05;
covpos = covar;
C=[[sigma1^2,covar]',[covar,sigma2^2]'];
CI=inv(C);
norm=1/(2*pi*sqrt(det(C)));
P=zeros(Nd1,Nd2);
for i = [1:Nd1]
for j = [1:Nd2]
dd = [d1(i)-d1bar, d2(j)-d2bar]';
P(i,j)=norm*exp( -0.5 * dd' * CI * dd );
end
end
pp=P;
% make a negatively correlated normal distribution
d1bar=1.5;
d2bar=2.5;
sigma1 = 0.25;
sigma2 = 0.25;
covar = -0.05;
covneg = covar;
C=[[sigma1^2,covar]',[covar,sigma2^2]'];
CI=inv(C);
norm=1/(2*pi*sqrt(det(C)));
pn=zeros(Nd1,Nd2);
for i = [1:Nd1]
for j = [1:Nd2]
dd = [d1(i)-d1bar, d2(j)-d2bar]';
pn(i,j)=norm*exp( -0.5 * dd' * CI * dd );
end
end
% make a uncorrelated normal distribution
d1bar=1.5;
d2bar=2.5;
sigma1 = 0.25;
sigma2 = 0.25;
covar = 0.0;
covzero = covar;
C=[[sigma1^2,covar]',[covar,sigma2^2]'];
CI=inv(C);
norm=1/(2*pi*sqrt(det(C)));
pu=zeros(Nd1,Nd2);
for i = [1:Nd1]
for j = [1:Nd2]
dd = [d1(i)-d1bar, d2(j)-d2bar]';
pu(i,j)=norm*exp( -0.5 * dd' * CI * dd );
end
end
% use simple drawing function that encapsulates all the graphics
eda_draw(pp,'caption +cor',pn,'caption -cor',pu,'caption uncor');
% edama_03_16: compute covariance of p(d1,d2)
% warning: you must have run the previous section before this one will work
% make the alternating sign function
S = (d1-d1bar) * (d2-d2bar)';
% integrate
covpos_est = (Dd1*Dd2)*sum(sum(S .* pp));
covneg_est = (Dd1*Dd2)*sum(sum(S .* pn));
covzero_est = (Dd1*Dd2)*sum(sum(S .* pu));
eda_draw(pp,'caption P',S,'caption S',S.*P,'caption SP');
fprintf( 'cov +: true %.4f est %.4f\n', covpos, covpos_est);
fprintf( 'cov -: true %.4f est %.4f\n', covneg, covneg_est);
fprintf( 'cov 0: true %.4f est %.4f\n', covzero, covzero_est);
% edama_03_17: error propagation from data to model parameters
% set up d1 axis
Nd1 = 101;
d1 = 5.0*[0:Nd1-1]'/(Nd1-1);
Dd1 = d1(2)-d1(1);
% set up d2 axis
Nd2 = 101;
d2 = 5.0*[0:Nd2-1]'/(Nd2-1);
Dd2 = d2(2) - d2(1);
% 2D normal p.f.d. setup, uncorrelated
dbar1 = 1.0;
dbar2 = 4.0;
sigma1 = 0.25;
sigma2 = 0.25;
covar = 0.0;
Cd = [sigma1^2, covar; covar, sigma2^2]; % covariance matrix
CdI = inv(Cd); % its inverse
sqrtdet = sqrt(det(Cd)); % square root of its determinant
norm=(1/(2*pi*sqrtdet)); % normalization
dd1 = d1-dbar1; % d1 minus its mean
dd2 = d2-dbar2; % d2 minus its mean
dd11 = dd1.^2; % (d1 minus its mean) squared
dd22 = dd2.^2; % (d2 minus its mean) squared
pd=norm*exp(-0.5*CdI(1,1)*(dd1.^2)*ones(Nd2,1)' ...
-0.5*CdI(2,2)*ones(Nd1,1)*(dd2.^2)' ...
-CdI(1,2)*dd1*dd2');
% model is related to data by m=Md
M = [ -1, 1; -2, 1 ];
% error propagation
mbar = M*[dbar1; dbar2];
m1bar = mbar(1);
m2bar = mbar(2);
Cm = M*Cd*M';
% print results
fprintf('mbar1=%.2f mbar2=%.2f\n',m1bar,m2bar);
fprintf('varm1=%.2f varm2=%.2f covm=%.2f\n', Cm(1,1), Cm(2,2), Cm(1,2) );
% set up m1 axis
Nm1 = 101;
m1 = 5.0*[0:Nm1-1]'/(Nm1-1);
Dm1 = m1(2)-m1(1);
% set up m2 axis
Nm2 = 101;
m2 = 5.0*[0:Nm2-1]'/(Nm2-1);
Dm2 = m2(2) - m2(1);
CmI=inv(Cm);
norm=1/(2*pi*sqrt(det(Cm)));
dm1=m1-m1bar;
dm2=m2-m2bar;
pm=norm*exp(-0.5*CmI(1,1)*(dm1.^2)*ones(Nm2,1)' ...
-0.5*CmI(2,2)*ones(Nm1,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)');
