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)
However, the chain provides us with far more information than just a single uncertainty. We can make the entire posterior probability function by histogramming the chain using the margin command e.g.
XSPEC12>margin 2 1.2 1.5 20
Probability Mod param
2
0.000000 1.207500
1.000000e-04 1.222500
2.600000e-04 1.237500
3.620000e-03 1.252500
9.800000e-03 1.267500
0.030480 1.282500
0.056680 1.297500
0.100820 1.312500
0.150400 1.327500
0.165320 1.342500
0.161040 1.357500
0.132880 1.372500
0.095400 1.387500
0.052060 1.402500
0.026580 1.417500
9.800000e-03 1.432500
2.580000e-03 1.447500
1.340000e-03 1.462500
4.600000e-04 1.477500
1.800000e-04 1.492500
The margin command has the same arguments as steppar. We can also plot the result using
XSPEC12>plot margin
In this case the probability distribution looks gaussian so the error command does indeed provide an adequate measure of the uncertainty.