A way to lift best-fit-params up from the depths of the pyhf.infer.hypotest implementation?

[lhenkelm]

Dear pyhf experts,

Is there a way to access the best-fit POI and NP values resulting from the fits that pyhf runs when one calls pyhf.infer.hypotest (or the underlying test statistics/calculators)? I have some funny/unexpected p-values coming out of hypotest, and I hope being able to see where those fits place mu-hat-hat (etc.) would help me understand what is going on.

The next best thing to a solution I could find is to use cabinetry.model_utils.get_asimov_parameters, andpyhf.infer.mle.fit and pyhf.infer.mle.fixed_poi_fit to manually duplicate the fits happening inside the pyhf implementation of the test statistics.
But then I would still have to trust that there are no subtle differences in implementation.

Thanks, Lars

[kratsg]

There are multiple fits happening in hypotest -- rather than a single one -- and the comparison of these fit is what allows you to compute the corresponding test statistic. Which one are you looking for?

[lhenkelm]

Ideally, I want both the NPs maximised for the fixed POI value, and the POI and NPs for the unconstrained fit, for whatever fixed mu values would be appropriate for the teststatistic I am looking at. (so mu = 0 for a discovery p0, mu = 0 and 1 for a CLs using qmu, etc.). That way I could feed them back into the model to get an idea of how the physical distributions are different under different hypotheses and the constrained vs, the free optimisation.

One reason I asked is that I was wondering if there was some way to avoid running each of the fits twice (once implicitly in hypotest, then by hand to access the parameters). If possible (not to imply it needs to be) that would be convenient, because I would not need to worry so much that the way I call pyhf.infer.mle.fit and pyhf.infer.mle.fixed_poi_fit actually shadows the test statistics implementations correctly.

[matthewfeickert]

Thanks for clarifying! At the moment there isn't any API to do this.

All of the inference functions try to take advantage of the pyhf.infer.calculators which themselves use the test statistic distributions. To get down to where fit and fixed_poi_fit are being called you need to go into the test statistic functions themselves:

pyhf/src/pyhf/infer/test_statistics.py

Lines 44 to 49 in 4b66cd2

	mubhathat, fixed_poi_fit_lhood_val = fixed_poi_fit(
	mu, data, pdf, init_pars, par_bounds, fixed_params, return_fitted_val=True
	)
	muhatbhat, unconstrained_fit_lhood_val = fit(
	data, pdf, init_pars, par_bounds, fixed_params, return_fitted_val=True
	)

[lukasheinrich]

options on the test stat functions to return thoose could be an interesting PR.. but also @lhenkelm - I think running the fits by had should be OK too, the chance to have a diveging cofig should be low

[lhenkelm]

Thanks @matthewfeickert for the answer, and pointing out the fit calls in the source! I think that is what I'll do for now.

About @lukasheinrich's suggestion to change the API, I (obviously) agree this might be useful, but I don't want to rush in and impose my preferences. So if you guys think this might be worthwhile adding, I am happy to put some work into a PR (or open an issue first to discuss details).

[matthewfeickert]

Thanks for the answer, and pointing out the fit calls in the source! I think that is what I'll do for now.

Sure. I hope that's helpful, but I'll admit that you need to do some call tracing to get from hypotest all the way down to the test stats, so apologies in advance if anything is a bit confusing there.

About @lukasheinrich's suggestion to change the API, I (obviously) agree this might be useful, but I don't want to rush in and impose my preferences. So if you guys think this might be worthwhile adding, I am happy to put some work into a PR (or open an issue first to discuss details).

Yeah, I have nothing against this as a feature PR as it wouldn't slow things down and could just be an additional kwarg and so wouldn't be disruptive to the API. As I'm already very behind on things for release v0.6.3 I wouldn't add this to the core dev project boards until at least v0.7.0, but @lhenkelm if you're interested in taking this on then please open an Issue like you suggested and I think we'd welcome a PR whenever from it. 👍

(Aside: We're trying to get better about figuring out how to accept new feature PRs while staying close to SemVer which will require a change in release branch strategy probably, but for this I think it could go in whenever. Example of where we've failed to do this well is PR #1274 which I feel bad about).