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
m = [1;2];
n = length(m);
P = [2, 0.75; 0.75, 0.5];
%P = eye(2);
nsam = 1000;
rng(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 += 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 += 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;
xplot = linspace(xrange(:,1),xrange(:,2), nplot);
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');
fill1 = [xplot(1,:)(ind(1,:)); fx(1,:)(ind(1,:))];
fill1 = [ [fill1(1,1);0], fill1, [fill1(1,end);0] ]';
fill2 = [xplot(2,:)(ind(2,:)); fx(2,:)(ind(2,:))];
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