Hessian of POI around its maximum

[WolfgangWaltenberger]

Hey guys,

in SModelS we have a use case where we fit the signal strength multiplier mu that maximizes the likelhood via
pyhf.infer.mle.fit. Is there a way to obtain an estimate of the hessian, or its inverse, at that point?
So d^2 NLL / d mu^2 | _{mu=mu_hat} ? I would want to employ the Cramer-Rao bound and use this Hessian to get a rough estimate of sigma_mu around mu_hat. Or is there another way to get the standard deviation for mu, so sigma_mu ( mu = mu_hat ), even if just crudely approximated?
Cheers

Wolfgang (SModelS)

[kratsg]

You should be able to get it through minuit.

>>> import pyhf
>>> model = pyhf.simplemodels.uncorrelated_background([5.0], [10.0], [2.5])
>>> data = [12.5] + model.config.auxdata
>>> pyhf.set_backend(pyhf.tensorlib, 'minuit')
>>> pyhf.infer.mle.fit(data, model, return_correlations=True)
(array([0.5015113 , 0.99947511]), array([[ 1.        , -0.57697422],
       [-0.57697422,  1.        ]]))

Keyword arguments for most/all fit functions are listed in the Optimizer API (https://scikit-hep.org/pyhf/_generated/pyhf.optimize.mixins.OptimizerMixin.html#pyhf.optimize.mixins.OptimizerMixin.minimize)

[kratsg]

You can also return the minuit object instead and interact with it as you expect with return_result_obj=True.

[WolfgangWaltenberger]

Gocha, thanks @kratsg . That settles my question.

Wolfgang