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 This chain has a different length than all other loaded Chains. Do you still want to load it? (y/n): y New chain run1.fits is now loaded.
Note that you get a warning because you are loading a chain whose length differs from those already in XSPEC. Let's tidy up by unloading the two unwanted chains.
XSPEC12>chain info
Output file format : fits
Current chain length setting : 50000
Burn-in length : 1000
Chain width : 4
Current temperature : 1.0000e+00
(Model parameters 1 2 3)
Randomization off
Current chain proposal distribution setting: cauchy <cmdline>
Loaded chains: length:
1. run1.fits 50000
2. test1.fits 1000
3. test2.fits 1000
XSPEC12>chain unload 2-3
Chains unloaded:
test1.fits
test2.fits
Now check the statistics on the full-length chain :
XSPEC12>chain stat 1
Means in chains : 9.6170305
Mean over all chains : 9.6170305 and variance of chain means : 0.
Variance over all chains : 0.090427714
Rubin-Gelman convergence measure : 0.99998000
Fraction of repeated values: 0.716380
(rule of thumb target is 0.75)
We can use this chain within XSPEC now to estimate confidence ranges eg
XSPEC12>error 1
Errors calculated from chains
Parameter Confidence Range (2.706)
1 9.1270669 10.129351 (-0.45603932,0.54624449)