Wstat (joint likelihood) crashes with sinning.py assertion error

[kchoX14]

It seems like it is impossible to have the posterior plots when the joint likelihood of the source+background and the background are used. BXAsolver() performs fine, but the binning expects integral counts in the plot background subtraction. This fails each time because it is likely that one or more source channels have counts lower than the corresponding background channel. Or the source may have 0 counts in a channel.

Could the background having 0 counts lead to the mess? Should a tiny bit of rebinning of the original data (to have at least 1 count per bin, like XSPEC recommends) be allowed to take this into account? Ignoring artificial features or smearing of information, what other problems may this lead to?

[JohannesBuchner]

I think what is happening is that bxa.binning tries to get the plot from xspec, but xspec typically plots a background-subtracted spectrum (even if the analysis does not subtract). This leads to problems, because the bins are not integer. Rebinning will not help. bxa.binning requires the counts to be integers for meaningful (Poisson) analyses. I don't know if there is a pyxspec switch to disable the subtraction.

Maybe the BXA plotting capabilities should be removed, it is quite limited.

This comment https://github.com/JohannesBuchner/BXA/issues/28#issuecomment-871572545 explains how you can do your own plotting.

[JohannesBuchner]

To expand, posterior_predictions_convolved uses posterior_predictions_plot, which extracts from xspec with:

xspec.Plot('counts')

[JohannesBuchner]

I believe the ideal solution would be if we were able to generate counts from the model, and then plot both predicted counts and observed counts on an area-corrected y-axis. But it is very difficult to do in xspec, especially without a background model.

[kchoX14]

I think what is happening is that bxa.binning tries to get the plot from xspec, but xspec typically plots a background-subtracted spectrum (even if the analysis does not subtract). This leads to problems, because the bins are not integer. Rebinning will not help. bxa.binning requires the counts to be integers for meaningful (Poisson) analyses. I don't know if there is a pyxspec switch to disable the subtraction.

Maybe the BXA plotting capabilities should be removed, it is quite limited.

This comment https://github.com/JohannesBuchner/BXA/issues/28#issuecomment-871572545 explains how you can do your own plotting.

Yes, either (maybe) the Xspec-imported plotting needs to go or the BXA analysis should be modified to temporarily store the counts information for its own posterior binning. The current mess does not allow anything with Wstat beyond optimized fit retrieval, corner plots, trace plots, and the plot for the evolution of the likelihood.

In addition, the Cstat variant (background modeled independent of the source) treatment makes a mess of the QQ visual inspection with multiple instruments/data over the same energy band. Even a PPC plot may not be that visually helpful here (not to forget AIC has a frequentist base and should not be used for model selection/comparison off Bayesian data treatment).

[JohannesBuchner]

BXA does not have access to the counts information during the fit, the statistic is calculated by xspec.

[JohannesBuchner]

If you have a idea or proof of concept how we could make a PPC-like plot work, it would be greatly appreciated.

Reading the counts from the input file is probably a brittle mess too.

PS: AIC is based on information theory, not frequentism, perhaps you are thinking of likelihood ratio or F tests.

[kchoX14]

PS: AIC is based on information theory, not frequentism, perhaps you are thinking of likelihood ratio or F tests.

"The AIC should not be used with bayesian fits."

This is an excerpt from one of the comments on the threeML repository by Michael. It is, as also mentioned later, suitable for MLE based fits and not Bayesian. What do you have to say on that?

[JohannesBuchner]

Well, AIC is based on the maximum likelihood, and it is not Bayesian.

I find it useful for empirical (ad-hoc) models, when model prior odds and parameter priors are difficult to justify. But visualisations and predictive quality are also powerful model comparison tools, especially for M-open (true model is not part of the specified set of models).

[kchoX14]

Agreed that it is MLE-driven. But visualizations (as it turns out with the PPC way over threeML) can get messy and crowded with multiple instruments spectra, like the QQ plots do. So, it is a compromizing trade off. Better optimization with Bayesian but relatively poorer model selection/comparison checks. Poor conversion with the MLE methods but reliable checks using toold like AIC. The only dependable saviour (in terms of numbers) may be parametric bootstrapping with MCMC. ThreeML allows for that, but BXA does not seem to do it since XSPEC has its own way of doing it. The multinormal assumption of MCMC/emcee is probably an issue here.

[JohannesBuchner]

Lets refocus on the issue: how to fix visualisations for wstat/cstat.

I don't know a way with wstat, because xspec does not have allow access to the estimated background count rate (I think?). Perhaps there is some secret option. Or what if we temporarily set background to none, plot counts, and then bring the background back. The difference should be the wstat contribution?

But for multi-spectrum Cstat with no background, it must be fixable. Perhaps we should have another issue for that. I think the bins can just be sorted by energy, and then QQ and posterior predictions could be useful.

[kchoX14]

Hiding the background to get the bins right would effectively render Wstat as Cstat, but post-analysis. I wonder what difference it'd make. A workable example will be needed, and I may give it a shot.

On QQ with Cstat, that's possibly another thing to look at. Maybe a workable example on your end would be nice for multiple spectra/instruments? It'd create a new BXA branch than the existing, single spectrum analyses. Cumulative of counts across spectra or instruments just dilutes all information with current modifications.