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