Figure 4.10:

The multivariate normal, marginals, marginal box, and bounding box.

Code for Figure 4.10

Text of the GNU GPL.

main.m

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
m = [1;2];
n = length(m);
P = [2, 0.75; 0.75, 0.5];
%P = eye(2);
nsam = 1000;
randn('seed', 0)
% generate the samples
x = repmat(m, 1, nsam) + sqrtm(P)*randn(n,nsam);

% compute the 95% confidence interval ellipse
alpha = 0.95;
A = inv(P);
b = chi2inv(alpha, n);
[xe, ye] = ellipse(A, b, 100, m);

% count how many samples are outside the  ellipse
e = x - repmat(m, 1, nsam);
elpts = 1 - sum( diag(e'*A*e)> b )/nsam

% 95% bounding box
xd = sqrt(b*diag(P))';
bbox = [-xd;  [xd(1), -xd(2)]; xd; [-xd(1), xd(2)]; -xd];
bbox = bbox + kron(ones(2*n+1,1), m');
%count samples
z = x-m;
bbpts = 1 - sum(sum(abs(z)<xd')<2)/nsam

% 95% marginal intervals
b2 = chi2inv(alpha, 1);
xd2 = sqrt(b2.*diag(P))';
mbox = [-xd2;  [xd2(1), -xd2(2)]; xd2; [-xd2(1), xd2(2)]; -xd2];
mbox = mbox + kron(ones(2*n+1,1), m');
% count samples
mbpts = 1 - sum(sum(abs(z)<xd2')<2)/nsam

% make the two 1-d plots of the marginals
xrange = round([min(x'); max(x')]');
nplot = 100;
% workaround; matlab's linspace is braindead
% octave: xplot = linspace(xrange(:,1),xrange(:,2), nplot);
tmp = linspace(0, 1, nplot);
xplot = xrange(:,1) + tmp.*(xrange(:,2)-xrange(:,1));

z = xplot-m;
fx = exp(-(1/2)*(z./diag(P).*z))./sqrt(2*pi*diag(P));

% compute closed curves for using filledcurve in gnuplot
ind = (abs(z) <= xd2');
xp1 = xplot(1,:);
xp2 = xplot(2,:);
fx1 = fx(1,:);
fx2 = fx(2,:);
ind1 = ind(1,:);
ind2 = ind(2,:);

fill1 = [xp1(ind1); fx1(ind1)];
fill1 = [ [fill1(1,1);0], fill1, [fill1(1,end);0] ]';
fill2 = [xp2(ind2); fx2(ind2)];
fill2 = [ [fill2(1,1);0], fill2, [fill2(1,end);0] ]';


% plot ellipse, bounding box, and marginal box
figure()
plot(x(1,:), x(2,:), 'o', xe, ye, bbox(:,1),bbox(:,2), mbox(:,1),mbox(:,2))

% histograms for the marginals of the samples
nbin = 20;
figure()
hist(x(1,:),nbin)
figure()
hist(x(2,:),nbin)

figure()


plot(xplot', fx')
mtab = [xplot', fx'];
%
samtab = x;
eltab = [xe ye];
samtab = x';
save confmarg.dat samtab eltab mtab bbox mbox fill1 fill2