#acl All:read

The release of XSPEC includes greatly enhanced support for Markov Chain Monte Carlo with more flexibility and improved integration with the rest of the XSPEC capabilities.

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 Metropolis-Hastings algorithm is used to generate the next element in the chain.

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. 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.

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 new release. 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.

The number of links in the chain is set using chain length. The beginning of the chain can be discarded using chain burn; this is useful to erase any dependency on the starting parameters and is referred to in the literature as the burn-in. 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.

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.

Alternatively, the analysis can continue outside XSPEC, using fv or other ftools. The fv histogram option can be used to create plots of the probability density as in the illustrations below.

attachment:fv_2d_histo.jpg  attachment:fv_1d_histo.jpg