Open pboorm opened 1 year ago
I don't quite understand how the numbers work.
Here is the current code: https://github.com/JohannesBuchner/BXA/blob/master/bxa/xspec/solver.py#L64
Previous issues: https://github.com/JohannesBuchner/BXA/issues/12 https://github.com/JohannesBuchner/BXA/issues/10
There is probably a way to get the number out from the model and parameter in the transformation list.
Thanks! I had a go at tweaking the store_chain() function. A version I wrote that appears to work for models with one data group is:
def store_chain_universal(chainfilename, transformations, posterior, fit_statistic):
"""Writes a MCMC chain file in the same format as the Xspec chain command."""
import astropy.io.fits as pyfits
group_index = 1
old_model = transformations[0]['model']
indices = [t["index"] for t in transformations]
names = []
for t in transformations:
if t['model'] != old_model:
group_index += 1
old_model = t['model']
original_parname = AllModels(1)(t["index"]).name
names.append('%s__%d' % (original_parname, t['index']))# + (group_index - 1) * old_model.nParameters))
columns = [pyfits.Column(
name=name, format='D', unit=AllModels(1)(transformations[i]["index"]).unit, array=t['aftertransform'](posterior[:, i]))
for i, name in enumerate(names)]
columns = list(np.array(columns)[np.argsort(indices)])
columns.append(pyfits.Column(name='FIT_STATISTIC', format='D', array=fit_statistic))
table = pyfits.ColDefs(columns)
header = pyfits.Header()
header.add_comment("""Created by BXA (Bayesian X-ray spectal Analysis) for Xspec""")
header.add_comment("""refer to https://github.com/JohannesBuchner/""")
header['TEMPR001'] = 1.
header['STROW001'] = 1
header['EXTNAME'] = 'CHAIN'
tbhdu = pyfits.BinTableHDU.from_columns(table, header=header)
tbhdu.writeto(chainfilename, overwrite=True)
This saves the column names as the original parameter names (e.g., before adding extensions around the parameter name that depend on the prior) as well as the units stored in the model parameters. However, the code is using "AllModels(1)", so will not work for models with multiple datagroups/models. I've left the corresponding code commented that I think was being used to account for this in the original version (i.e. here: https://github.com/JohannesBuchner/BXA/blob/631e9f40f101aba781f89f1d1e53ffcfe9262519/bxa/xspec/solver.py#L68).
When I load the chain generated with this function, the error on flux can be computed from the chain.
Description
When using the BXA-generated chain files in Xspec after running a BXA fit, Xspec cannot estimate uncertainties from the chain with e.g., "flux 2. 10. err 1000". The issue seems to be that the BXA chain file column headings are incompatible with the model used.
What I Did
The model fit was:
And the free parameter numbers in the BXA fit were: 3, 6, 8, 9, 10, 12. The corresponding chain.fits file generated by BXA had column headings (note the nH prior was added out of order of the other parameters):
The problem was fixed by manually changing the BXA chain file column headings, column orders and parameter units to match an equivalent Xspec-generated MCMC chain file: