This page describes how to use the chain command to run Markov Chain Monte Carlo in XSPEC. You might also like to see our MCMC example.
Basics
The basic idea of MCMC is to produce a chain of parameter values whose density gives the probability distribution for that parameter. The XSPEC chain command is used to create and manipulate these chains, which are stored as FITS or ascii files. The fundamental parameter for a chain is its length which we set by eg :
XSPEC12> chain length 50000
It is generally useful to erase the effects of the starting parameters of the chain by discarding an initial run. So eg :
XSPEC12> chain burn 1000
tells XSPEC to perform 1000 steps and discard the results then perform 50000 steps for the chain proper. The command to run the chain is :
XSPEC12> chain run mychain.fits
which will place the output chain in the FITS file mychain.fits. In general it is a good idea to run a number of chains instead of one long one to provide a consistency check. This also allows a crude parallelization if multiple machines are available.
The proposal distribution
The critical practical issue with MCMC is to choose the proposal distribution from which trial parameters are drawn. It is a remarkable fact that any distribution will work in the sense that it will eventually generate a valid chain, however the correct choice makes a huge difference in how long this takes.
The basic command to set up the proposal distribution is chain proposal The two distribution types available are gaussian or cauchy (alias Lorentzian - a gaussian with fat tails). Both distributions are multivariate and require an input covariance matrix. If a fit has been performed then the covariance at best-fit can be selected by using the fit option e.g. :
XSPEC12> chain proposal gaussian fit
Alternatively, a matrix can be entered on the command line or read in from a file. If an MCMC chain has already been run then a new covariance matrix can be calculated from it. This last option is specified by :
XSPEC12> chain proposal gaussian chain
Our experience is that a good strategy is to use the fit option and then divide the entire covariance matrix by a factor of 8. A script (rescalecov.tcl) to do this is included with the current release and can be run by :
XSPEC12> rescalecov 0.125
The standard rule-of-thumb for a good proposal distribution is that the chain should have a repeat fraction of 0.75. The chain stat command writes out the repeat fraction.
Tempering
The Metropolis-Hastings algorithm is used to generate the next element in the chain. The basic algorithm can be modified to include tempering by setting a temperature value using chain temperature. The algorithm is then as follows. Suppose that our current set of parameters give a fit statistic value of S_0. We generate a new proposal set of parameters and calculate a new fit statistic of S_1. If S_1 <= S_0 we accept the new set of parameters. If not, then we generate random number R distributed uniformly between 0 and 1. If the temperature is T then we accept the new set of parameters if R < exp(0.5*(S_0-S_1)/T).
Using the MCMC chain results
There are a number of ways to use the MCMC chains inside XSPEC. If any chains are loaded the error, flux error, lumin error, and eqwidth error commands will all use the chains instead of performing their own independent calculations. This should be a more accurate estimate and in particular will eliminate the problem reported by some users of the flux error range not including the best-fit value.
The chains can be turned into a probability distribution using the margin command. This takes the same arguments as steppar and calculates probabilities on a multi-dimensional grid. Parameters not specified in the arguments to margin are automatically integrated over (marginalized).
The results of margin can be plotted using plot margin. If only one parameter was specified in margin the result is a simple line plot, if two parameters then a contour plot. This contour plot is of the probability density not of the integrated probability within the contour and hence is not directly comparable to the results of plot contour after steppar.