#

MCMC toolbox » Examples » Normal distribution

## MCMC toolbox example

In this example, we generate Gaussian target with known covariance matrix. The target distribution has a known form and could be calculated analytically, so this simulation is mainly for testing of the algorithms. We study the correctness of the chain by calculating points that lie inside 50% and 95% probability contours.

    clear model data params options

    nsimu = 10000;  % number of simulations
    npar  = 100;     % dimension of the unknown

    data.x0 = zeros(1,npar); % mean vector

We create the needed prameter structure in a loop

    for i=1:npar
      %            'name'                     initial
      params{i} = {sprintf('\\theta_{%d}',i), data.x0(i)};
    end

Create covariance and precision matrises. Function `covcond` creates a covariance matrix with given condition number and a direction of the first principal axis. It returns the covariance matrix and its inverse, the precision matrix.

    [Sig,Lam] = covcond(100,ones(npar,1));

Store the precision matrix in `data` structure so we can use it in the `ssfun`.

    data.Lam = Lam;

The function `ssfun` for `mcmcrun` is the quadratic form in the Gaussian distribution.

    model.ssfun     = @(x,data) (x-data.x0)*data.Lam*(x-data.x0)';
    model.N         = 1;

For `mcmcrun` we use scaled version of the known target covariance as the proposal covariance. This scaling is known to produce an optimal proposal.

    options.nsimu   = nsimu;
    options.qcov    = 2.4^2/npar*Sig;
    options.method  = 'ram';  % use the (default) DRAM method
    options.verbosity = 1;

Generate the MCMC chain.

    [results,chain] = mcmcrun(model,data,params,options);

    Setting nbatch to 1
    Sampling these parameters:
    name   start [min,max] N(mu,s^2)
    \theta_{1}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{2}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{3}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{4}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{5}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{6}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{7}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{8}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{9}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{10}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{11}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{12}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{13}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{14}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{15}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{16}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{17}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{18}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{19}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{20}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{21}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{22}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{23}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{24}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{25}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{26}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{27}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{28}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{29}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{30}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{31}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{32}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{33}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{34}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{35}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{36}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{37}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{38}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{39}: 0 [-Inf,Inf] N(0,Inf)
    \theta_{40}: 0 [-Inf,Inf] N(0,Inf)
    ...

    figure(1); clf
    mcmcplot(chain,[1:4],results.names,'chainpanel')

From the generated chain we calculate the relative distances of the chain points from the origin and count the points that are inside given probability limits. Then we plot the first two dimensions of the chain together with the correct probability contours.

The title of the 2d plot shows the rejection rate and the propotion of points inside the ellipsoids. Number tau*t in the title tells how many seconds it takes to generate 1000 independent samples according to the integrated autocorrelation time (`iact`) estimate.

    d    = mahalanobis(chain(:,1:npar),data.x0,Lam,1);
    c50  = chiqf(0.50,npar);
    c95  = chiqf(0.95,npar);
    cc50 = sum(d<c50)./nsimu;
    cc95 = sum(d<c95)./nsimu;

    figure(2); clf
    mcmcplot(chain,[1,2],results.names,'pairs',0)

    title(sprintf('Rejected = %.1f%%, c50 = %.1f%%, c95 = %.1f%%, \\tau*t=%.1f', ...
                  results.rejected*100, cc50*100, cc95*100, ...
                  results.simutime/results.nsimu*1000*mean(iact(chain))))

    c50  = chiqf(0.50,2);
    c95  = chiqf(0.95,2);

    hold on
    ellipse(data.x0(1:2),c50*Sig(1:2,1:2),'r--','LineWidth',2);
    ellipse(data.x0(1:2),c95*Sig(1:2,1:2),'r-','LineWidth',2);
    axis equal
    hold off