XSPEC v12.8 - Markov Chain Monte Carlo

The Metropolis-Hastings algorithm used by XSPEC for Markov Chain Monte Carlo requires a choice of proposal distribution. Finding the best distribution can be difficult and makes MCMC harder to use. A new algorithm Goodman-Weare is included in v12.8. This algorithm does not require a choice of proposal distribution. It works by holding multiple sets of parameters, called walkers, for each step in the chain and generating walkers for the next step using those from the current step. This appears to work very well in many cases. The same algorithm is the basis of the emcee python package.

To use the new option in XSPEC do the following

XSPEC12> chain type gw
XSPEC12> chain walkers 8
XSPEC12> chain length 10000

The actual number of steps performed will be length/walkers so length should be chosen as a multiple of the number of walkers. The output chain combines all the walkers so in this case the first eight rows will be the walkers for the first step, the next eight for the second step and so on. A similar rule holds for the burn-in length set by chain burn.

Note that the output from running a chain is fully integrated into XSPEC. So, if a chain is loaded then tclout simpars will draw from the distribution defined by the chain, the error command will be based on the chain, as will the error options on flux and eqwidth.

Another new feature related to MCMC is the thin option on plot chain. It can be difficult to see anything on a plot with many thousands of points so this option applies a thinning factor to only show every Nth step in the chain. For instance:

XSPEC12> plot chain thin 100 1

will plot the values of parameter 1 every 100th step.