MCMC toolbox » Examples » Banana

Contents

Banana example

This techncal example constructs a non Gaussian target distribution by twisting two first dimensions of Gaussian distribution. The Jacobian of the transformation is 1, so it is easy to calculate the right probability regions for the banana and study different adaptive methods.

clear model data params options

'banana' sum-of-squares

bananafun = @(x,a,b) [a.*x(:,1),x(:,2)./a-b.*((a.*x(:,1)).^2+a^2),x(:,3:end)];
bananainv = @(x,a,b) [x(:,1)./a,x(:,2).*a+a.*b.*(x(:,1).^2+a^2),x(:,3:end)];
bananass  = @(x,d) bananainv(x-d.mu,d.a,d.b)*d.lam*bananainv(x-d.mu,d.a,d.b)';
a = 1; b = 1;         % banana parameters

npar = 12;             % number of unknowns
rho  = 0.9;            % target correlation
sig  = eye(npar); sig(1,2) = rho; sig(2,1) = rho;
lam  = inv(sig);       % target precision
mu   = zeros(1,npar);  % center

the data structure and parameters

data = struct('mu',mu,'a',a,'b',b,'lam',lam);
for i=1:npar
  params{i} = {sprintf('x_%d',i),0};
end

model.ssfun     = bananass;
model.N         = 1;

options.method  = 'dram';
options.nsimu   = 20000;
options.qcov    = eye(npar)*5; % [initial] proposal covariaance
[results,chain] = mcmcrun(model,data,params,options);
Setting nbatch to 1
Sampling these parameters:
name   start [min,max] N(mu,s^2)
x_1: 0 [-Inf,Inf] N(0,Inf)
x_2: 0 [-Inf,Inf] N(0,Inf)
x_3: 0 [-Inf,Inf] N(0,Inf)
x_4: 0 [-Inf,Inf] N(0,Inf)
x_5: 0 [-Inf,Inf] N(0,Inf)
x_6: 0 [-Inf,Inf] N(0,Inf)
x_7: 0 [-Inf,Inf] N(0,Inf)
x_8: 0 [-Inf,Inf] N(0,Inf)
x_9: 0 [-Inf,Inf] N(0,Inf)
x_10: 0 [-Inf,Inf] N(0,Inf)
x_11: 0 [-Inf,Inf] N(0,Inf)
x_12: 0 [-Inf,Inf] N(0,Inf)
figure(1); clf
mcmcplot(chain,[],results.names,'chainpanel')
figure(2); clf
mcmcplot(chain,[1,2],results.names,'pairs',0)
c50=1.3863; % critical values from chisq(2) distribution
c95=5.9915;
hold on
[xe,ye]=ellipse(mu,c50*sig(1:2,1:2));
xyplot(bananafun([xe,ye],a,b),'k-','LineWidth',2)
[xe,ye]=ellipse(mu,c95*sig(1:2,1:2));
xyplot(bananafun([xe,ye],a,b),'k-','LineWidth',2)
axis equal
hold off
title('2 first dimensions of the chain with 50% and 95% target contours')