Probing Parameter Space using Sampler in Sherpa

There is no fully Bayesian sampler inside Sherpa yet, however, there is a way to sample parameter space after the fit converged. The main assumption is that sampling starts from the best fit parameter values, the free parameters have Normal distributions with the scales  based on the covariance matrix with the off-diagonal elements to account for correlations.  The function "sample_energy_flux()" returns an array with sampled parameters and calculated flux for these parameters.

The thread   http://cxc.harvard.edu/sherpa/threads/flux_dist/index.py.html

describes how to use the functions to calculate the flux errors.

For the parameter errors it is easy to check the distribution of parameters returned by the sample_energy_flux()

Example:

sample100=sample_energy_flux(0.5,2.,num=1000)

print sample100

[...

[  2.90602963e-10   1.10100135e+00   8.39761377e-01   7.01873798e-01
   2.35626905e+00   1.03667988e+00]...

]

The first column contains the flux value for all the parameters in the other columns.

We can easily calculate 95% bounds for a given parameter:

sherpa-100> print sample100[:,1]

----------> print(sample100[:,1])

[ 1.10121073  1.1011715   1.10100135  1.10191317  1.10233657  1.10103386

  1.10360142  1.09774949  1.0984214   1.09927707  1.10041681  1.09919926

sherpa-102> min(sample100[:,1])

            1.0972469907219351

sherpa-103> max(sample100[:,1])

            1.1065345259787218

We can  use numpy to obtain the properties of the parameter distribution:

sherpa-104> numpy.median(sample100[:,1])

            1.1015496488552339

sherpa-105> numpy.mean(sample100[:,1])

            1.1015685708445886

sherpa-106> numpy.std(sample100[:,1])

            0.0015635753706038143

We can also sort and make a histogram and plot the distribution:

sort_par=numpy.sort(sample100[:,1])
sort_par[0.99*len(sample100[:,1])-1]

# make a histogram of parameter value and plot:
xspace=numpy.linspace(min(sample100[:,1]),max(sample100[:,1]),30)
xlo=xspace[:-1]
xhi=xspace[1:]
parhist=histogram1d(sample100[:,1],xlo,xhi)
parhist100=parhist/1000.
add_histogram(xlo,parhist100)

print_window("hist.png")

The histogram can be saved, but I could not place the figure into the Blogger - a known issue with the Google Blogger editor.