by $$\rmn{Pr}(\bmath{\theta}|\bmath{d},\mathcal{M}_\rmn{j}) \propto
  \rmn{Pr}(\bmath{d}|\bmath{\theta},\mathcal{M}_\rmn{j})^{\lambda}\rmn{Pr}(\bmath{\theta}|\mathcal{M}_\rmn{j}),$$ where the cooling factor ($\lambda$) is slowly increased from 0 to 1. The rate of cooling (equal to the change in $\lambda$ over one step) should be much lower than unity in order to ensure that BAYESYS is robust to finding regions of high posterior value in large prior volumes. Therefore the rate of cooling is set to a value of 0.1 and multiple runs are performed with different seed values to check for good convergence. The robustness of the process is improved by using an ensemble of 10 simultaneous MCMC chains that communicate after each step and so ensure that they do not become trapped in local maxima of likelihood.

After burn-in, the program draws samples from the posterior distribution at $\lambda = 1$, where the total number of samples drawn is equal to the product of the number of MCMC chains and steps. The number density of the model samples should then be a good estimate of the true posterior probability distribution. The resulting output includes a list of the samples which can be marginalised over to obtain posterior distributions for the cluster properties, as well as an estimate of the evidence.

Testing the 