#acl All:read

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)
}}}