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