This page shows a simple example of MCMC analysis in XSPEC. We use for data file1.pha in spectral/session and a simple absorbed power-law for the model :
XSPEC12> data file1 XSPEC12> model phabs(pow)
start by doing a fit
XSPEC12> fit
to give the result
======================================================================== Model phabs<1>*powerlaw<2> Source No.: 1 Active/On Model Model Component Parameter Unit Value par comp 1 1 phabs nH 10^22 9.5831062 +/- 0.29171265 2 2 powerlaw PhoIndex 1.3142755 +/- 3.3257374E-02 3 2 powerlaw norm 0.12037190 +/- 8.2219031E-03 ________________________________________________________________________ Chi-Squared = 78.09 using 64 PHA bins. Reduced chi-squared = 1.280 for 61 degrees of freedom Null hypothesis probability = 6.917881e-02
We will start by estimating a good proposal distribution. Start by using a multivariate cauchy using the covariance matrix from the fit. Set up the chain length and burn-in for a short test chain and run it :
XSPEC12> chain proposal cauchy fit XSPEC12> chain length 1000 XSPEC12> chain burn 0 XSPEC12> chain run test1.fits
To find out how well the chain worked we use chain stat on the first parameter :
Means in chains : 9.6619207 Mean over all chains : 9.6619207 and variance of chain means : 0. Variance over all chains : 0.075465715 Rubin-Gelman convergence measure : 0.99900000 Fraction of repeated values: 0.823000 (rule of thumb target is 0.75)
This looks pretty good: the fraction of repeated values is close to the target of 0.75. Let's do another test by halving all the covariance matrix values and see what effect that has :
XSPEC12> rescalecov 0.5 XSPEC12> chain run test2.fits XSPEC12> chain stat 1 Means in chains : 9.6619207 9.5525420 Mean over all chains : 9.6072314 and variance of chain means : 0.0059818499 Variance over all chains : 0.089958227 Rubin-Gelman convergence measure : 1.0655624 Fraction of repeated values: 0.823000 0.695000 (rule of thumb target is 0.75)
This looks OK so let's try a longer run :
XSPEC12> chain length 50000 XSPEC12> chain burn 1000 XSPEC12> chain run run1.fits