next up previous contents
Next: INTEGRAL/SPI: A Walk Through Up: Walks through XSPEC Previous: Producing Plots: Modifying the

Markov Chain Monte Carlo Example

To illustrate MCMC methods we will use the same data as the first walkthrough.

XSPEC12> data s54405
... 
XSPEC12> model phabs(pow)
... 
XSPEC12> fit
...
XSPEC12> chain type gw
XSPEC12> chain walkers 8
XSPEC12> chain length 10000

We use the Goodman-Weare algorithm with 8 walkers and a total length of 10,000. For the G-W algorithm the actual number of steps are 10,000/8 but we combine the results from all 8 walkers into a single output file. We start the chain at the default model parameters except that we use the renorm command to make sure that the model and the data have the same normalization. If we had multiple additive parameters with their own norms then a good starting place would be to use the fit 1 command to initially set the normalizations to something sensible.

XSPEC12> chain run test1.fits

The first thing to check is what happened to the fit statistic during the run.

XSPEC12> plot chain 0

Figure 4.11: The statistic from an MCMC run showing the initial burn-in phase.
Image figk

The result is shown in Figure 4.11, which plots the statistic value against the chain step. It is clear that after about 2000 steps the chain reached a steady state. We would usually have told XSPEC to discard the first few thousand steps but included them for illustrative purposes. Let us do this again but specifying a burn-in phase that will not be stored.

XSPEC12> chain burn 5000
XSPEC12> chain run test1.fits

The output chain now comprises 10,000 steady-state samples of the parameter probability distribution. Repeating plot chain 0 will confirm that the chain is in a steady state. The other parameter values can be plotted either singly using eg plot chain 2 for the power-law index or in pairs eg plot chain 1 2 giving a scatter plot as shown in Figure 4.12.

Using the error command at this point will generate errors based on the chain values.

XSPEC12>error 1 2 3
Errors calculated from chains

 Parameter   Confidence Range (2.706)
     1     0.264971     0.919546
     2       2.1134      2.41307
     3    0.0107304    0.0171814

The 90% confidence ranges are determined by ordering the parameter values in the chain then finding the center 90%.

Figure 4.12: The scatter plot from a 10,000 step MCMC run.
Image figl

To make confidence regions in two dimensions the margin command takes the same arguments as steppar.

XSPEC12>margin 1 0.0 1.5 25 2 1.5 3.0 25

Then use plot integprob to produce a plot of confidence regions.


next up previous contents
Next: INTEGRAL/SPI: A Walk Through Up: Walks through XSPEC Previous: Producing Plots: Modifying the