No dynamic range in posteriors with MultiNest (BXA v2.9)

[pboorm]

• BXA version: 2.9
• Python version: 3.8.12
• Xspec version: 12.12.1
• Operating System: macOS Monterey v12.6

Description

I am fitting 9 spectra simultaneously with PyXspec (one XMM/PN, one XMM/MOS1, one XMM/MOS2 for three different epochs). Since the spectra are all low signal to noise (source is <~50% of total counts), I am fitting a simple model but have a cross-calibration constant for each dataset that are free to vary relative to the first one (i.e. 8 free parameters in addition to the free parameters of the main model). For these cross-calibration constants I am using a custom log-Gaussian prior (code attached) in the BXA fit.

After fitting, the posteriors for a lot of the constants and some of the model parameters are all a single value which then gives an error when trying to produce the corner plot with corner after the fit since the parameter posteriors have no dynamic range.

Are there some suggested practices to overcome this? E.g., would increasing the number of live points help for the fit? This issue does not occur when I only fit 3 spectra simultaneously (i.e. with two cross-calibration constants), so I am wondering if the custom log-Gaussian prior may be too restrictive for fitting with a lot of spectra?

Code

Here is the code I am using for the custom log-Gaussian prior

def create_loggaussian_prior_for(model, par, mean, std):
    """
    This combines the Gaussian and loguniform BXA priors to give
    a Gaussian distribution of the log-parameter.

    NOTE: mean and std are also in log units
    """
    import scipy.stats
    pval, pdelta, pmin, pbottom, ptop, pmax = par.values

    if pmin == 0:
        raise Exception('You forgot to set reasonable parameter limits on %s' % par.name)
    low = np.log10(pmin)
    spread = np.log10(pmax) - np.log10(pmin)
    hi = np.log10(pmax)
    if spread > 10:
        print('   note: this parameter spans *many* dex. Double-check the limits are reasonable.')

    print('  log-gaussian prior for %s of %f +- %f' % (par.name, mean, std))
    rv = scipy.stats.norm(mean, std)
    def loggauss_transform(x): 
        return max(low, min(hi, rv.ppf(x * spread + low)))
    def loggauss_after_transform(x): return 10**x
    return dict(model=model, index=par._Parameter__index, name='log(%s)' % par.name, 
        transform=loggauss_transform, aftertransform=loggauss_after_transform)

[JohannesBuchner]

Have you tried plotting the best fit?

[pboorm]

Thanks! After plotting it, the best-fit for all datasets seems to be below the data systematically.

[JohannesBuchner]

Probably it is fighting the data or the prior range, and your likelihood is very low. Try to by hand find a not-horrible fit and make sure those parameters are allowed in your prior. That means making your model more flexible.

[pboorm]

Thanks a lot for your help, Johannes. Increasing the custom log-Gaussian prior standard deviation allowed sufficient dynamic range for all parameters in the fit to generate a corner plot.