#

MCMC toolbox » Examples » Banana

## 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   = 10000;
    options.qcov    = eye(npar)*0.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')