diff --git a/docs/api.rst b/docs/api.rst index 6118f6095a..048e73728e 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -140,6 +140,7 @@ Calculators :nosignatures: calculators.generate_asimov_data + calculators.HypoTestFitResults calculators.AsymptoticTestStatDistribution calculators.EmpiricalDistribution calculators.AsymptoticCalculator diff --git a/docs/contributors.rst b/docs/contributors.rst index 264c4405fa..51d05c1363 100644 --- a/docs/contributors.rst +++ b/docs/contributors.rst @@ -25,4 +25,5 @@ Contributors include: - Saransh Chopra - Sviatoslav Sydorenko - Mason Proffitt +- Lars Henkelmann - Aryan Roy diff --git a/src/pyhf/infer/__init__.py b/src/pyhf/infer/__init__.py index 851e19adb0..501c55a588 100644 --- a/src/pyhf/infer/__init__.py +++ b/src/pyhf/infer/__init__.py @@ -28,6 +28,7 @@ def hypotest( return_tail_probs=False, return_expected=False, return_expected_set=False, + return_calculator=False, **kwargs, ): r""" @@ -66,9 +67,12 @@ def hypotest( return_tail_probs (:obj:`bool`): Bool for returning :math:`\mathrm{CL}_{s+b}` and :math:`\mathrm{CL}_{b}` return_expected (:obj:`bool`): Bool for returning :math:`\mathrm{CL}_{\mathrm{exp}}` return_expected_set (:obj:`bool`): Bool for returning the :math:`(-2,-1,0,1,2)\sigma` :math:`\mathrm{CL}_{\mathrm{exp}}` --- the "Brazil band" + return_calculator (:obj:`bool`): Bool for returning calculator. Returns: - Tuple of Floats and lists of Floats: + Tuple of Floats and lists of Floats and + a :py:class:`~pyhf.infer.calculators.AsymptoticCalculator` + or :py:class:`~pyhf.infer.calculators.ToyCalculator` instance: - :math:`\mathrm{CL}_{s}`: The modified :math:`p`-value compared to the given threshold :math:`\alpha`, typically taken to be :math:`0.05`, @@ -139,6 +143,12 @@ def hypotest( referred to as the "Brazil band". Only returned when ``return_expected_set`` is ``True``. + - a calculator: The calculator instance used in the computation of the :math:`p`-values. + Either an instance of :py:class:`~pyhf.infer.calculators.AsymptoticCalculator` + or :py:class:`~pyhf.infer.calculators.ToyCalculator`, + depending on the value of ``calctype``. + Only returned when ``return_calculator`` is ``True``. + """ init_pars = init_pars or pdf.config.suggested_init() par_bounds = par_bounds or pdf.config.suggested_bounds() @@ -188,6 +198,8 @@ def hypotest( _returns.append(pvalues_exp_band) elif return_expected: _returns.append(tb.astensor(pvalues_exp_band[2])) + if return_calculator: + _returns.append(calc) # Enforce a consistent return type of the observed CLs return tuple(_returns) if len(_returns) > 1 else _returns[0] diff --git a/src/pyhf/infer/calculators.py b/src/pyhf/infer/calculators.py index 15d6ec24f3..1988f9872c 100644 --- a/src/pyhf/infer/calculators.py +++ b/src/pyhf/infer/calculators.py @@ -12,6 +12,7 @@ from pyhf.infer import utils import tqdm +from dataclasses import dataclass import logging log = logging.getLogger(__name__) @@ -29,7 +30,9 @@ def __dir__(): return __all__ -def generate_asimov_data(asimov_mu, data, pdf, init_pars, par_bounds, fixed_params): +def generate_asimov_data( + asimov_mu, data, pdf, init_pars, par_bounds, fixed_params, return_fitted_pars=False +): """ Compute Asimov Dataset (expected yields at best-fit values) for a given POI value. @@ -46,6 +49,14 @@ def generate_asimov_data(asimov_mu, data, pdf, init_pars, par_bounds, fixed_para >>> pyhf.infer.calculators.generate_asimov_data(mu_test, data, model, None, None, None) array([ 60.61229858, 56.52802479, 270.06832542, 48.31545488]) + It is possible to access the Asimov parameters as well: + + >>> pyhf.infer.calculators.generate_asimov_data( + ... mu_test, data, model, None, None, None, + ... return_fitted_pars = True + ... ) + (array([ 60.61229858, 56.52802479, 270.06832542, 48.31545488]), array([1. , 0.97224597, 0.87553894])) + Args: asimov_mu (:obj:`float`): The value for the parameter of interest to be used. data (:obj:`tensor`): The observed data. @@ -56,15 +67,23 @@ def generate_asimov_data(asimov_mu, data, pdf, init_pars, par_bounds, fixed_para The shape should be ``(n, 2)`` for ``n`` model parameters. fixed_params (:obj:`tensor` of :obj:`bool`): The flag to set a parameter constant to its starting value during minimization. + return_fitted_pars (:obj:`bool`): Return the best-fit parameter values for the given ``asimov_mu``. + Returns: - Tensor: The Asimov dataset. + A Tensor or a Tuple of two Tensors: + + - The Asimov dataset. + - The Asimov parameters. Only returned if ``return_fitted_pars`` is ``True``. """ bestfit_nuisance_asimov = fixed_poi_fit( asimov_mu, data, pdf, init_pars, par_bounds, fixed_params ) - return pdf.expected_data(bestfit_nuisance_asimov) + asimov_data = pdf.expected_data(bestfit_nuisance_asimov) + if return_fitted_pars: + return asimov_data, bestfit_nuisance_asimov + return asimov_data class AsymptoticTestStatDistribution: @@ -188,6 +207,21 @@ def expected_value(self, nsigma): ) +@dataclass(frozen=True) +class HypoTestFitResults: + """ + Fitted model parameters of the fits in + :py:meth:`AsymptoticCalculator.teststatistic ` + """ + + # ignore "F821 undefined name 'Tensor'" so as to avoid typing.Any + asimov_pars: 'Tensor' # noqa: F821 + free_fit_to_data: 'Tensor' # noqa: F821 + free_fit_to_asimov: 'Tensor' # noqa: F821 + fixed_poi_fit_to_data: 'Tensor' # noqa: F821 + fixed_poi_fit_to_asimov: 'Tensor' # noqa: F821 + + class AsymptoticCalculator: """The Asymptotic Calculator.""" @@ -251,6 +285,7 @@ def __init__( self.test_stat = test_stat self.calc_base_dist = calc_base_dist self.sqrtqmuA_v = None + self.fitted_pars = None def distributions(self, poi_test): r""" @@ -297,9 +332,13 @@ def distributions(self, poi_test): return sb_dist, b_dist def teststatistic(self, poi_test): - """ + r""" Compute the test statistic for the observed data under the studied model. + The fitted parameters of the five fits that are implicitly ran at every call + of this method are afterwards accessible through ``self.fitted_pars``, + which is a :py:class:`~pyhf.infer.calculators.HypoTestFitResults` instance. + Example: >>> import pyhf @@ -314,6 +353,16 @@ def teststatistic(self, poi_test): >>> asymptotic_calculator.teststatistic(mu_test) array(0.14043184) + Access the best-fit parameters afterwards: + + >>> asymptotic_calculator.fitted_pars + HypoTestFitResults(asimov_pars=array([0. , 1.0030482 , 0.96264534]), free_fit_to_data=array([0. , 1.0030512 , 0.96266961]), free_fit_to_asimov=array([0. , 1.00304893, 0.96263365]), fixed_poi_fit_to_data=array([1. , 0.97224597, 0.87553894]), fixed_poi_fit_to_asimov=array([1. , 0.97276864, 0.87142047])) + + E.g. the :math:`\hat{\mu}` and :math:`\hat{\theta}` fitted to the asimov dataset: + + >>> asymptotic_calculator.fitted_pars.free_fit_to_asimov + array([0. , 1.00304893, 0.96263365]) + Args: poi_test (:obj:`float` or :obj:`tensor`): The value for the parameter of interest. @@ -325,35 +374,45 @@ def teststatistic(self, poi_test): teststat_func = utils.get_test_stat(self.test_stat) - qmu_v = teststat_func( + qmu_v, (mubhathat, muhatbhat) = teststat_func( poi_test, self.data, self.pdf, self.init_pars, self.par_bounds, self.fixed_params, + return_fitted_pars=True, ) sqrtqmu_v = tensorlib.sqrt(qmu_v) asimov_mu = 1.0 if self.test_stat == 'q0' else 0.0 - asimov_data = generate_asimov_data( + asimov_data, asimov_mubhathat = generate_asimov_data( asimov_mu, self.data, self.pdf, self.init_pars, self.par_bounds, self.fixed_params, + return_fitted_pars=True, ) - qmuA_v = teststat_func( + qmuA_v, (mubhathat_A, muhatbhat_A) = teststat_func( poi_test, asimov_data, self.pdf, self.init_pars, self.par_bounds, self.fixed_params, + return_fitted_pars=True, ) self.sqrtqmuA_v = tensorlib.sqrt(qmuA_v) + self.fitted_pars = HypoTestFitResults( + asimov_pars=asimov_mubhathat, + free_fit_to_data=muhatbhat, + free_fit_to_asimov=muhatbhat_A, + fixed_poi_fit_to_data=mubhathat, + fixed_poi_fit_to_asimov=mubhathat_A, + ) if self.test_stat in ["q", "q0"]: # qmu or q0 teststat = sqrtqmu_v - self.sqrtqmuA_v diff --git a/src/pyhf/infer/test_statistics.py b/src/pyhf/infer/test_statistics.py index c0e4c51212..18d75d1738 100644 --- a/src/pyhf/infer/test_statistics.py +++ b/src/pyhf/infer/test_statistics.py @@ -13,7 +13,9 @@ def __dir__(): return __all__ -def _qmu_like(mu, data, pdf, init_pars, par_bounds, fixed_params): +def _qmu_like( + mu, data, pdf, init_pars, par_bounds, fixed_params, return_fitted_pars=False +): """ Clipped version of _tmu_like where the returned test statistic is 0 if muhat > 0 else tmu_like_stat. @@ -22,12 +24,14 @@ def _qmu_like(mu, data, pdf, init_pars, par_bounds, fixed_params): qmu_tilde. Otherwise this is qmu (no tilde). """ tensorlib, optimizer = get_backend() - tmu_like_stat, (_, muhatbhat) = _tmu_like( + tmu_like_stat, (mubhathat, muhatbhat) = _tmu_like( mu, data, pdf, init_pars, par_bounds, fixed_params, return_fitted_pars=True ) qmu_like_stat = tensorlib.where( muhatbhat[pdf.config.poi_index] > mu, tensorlib.astensor(0.0), tmu_like_stat ) + if return_fitted_pars: + return qmu_like_stat, (mubhathat, muhatbhat) return qmu_like_stat @@ -56,7 +60,7 @@ def _tmu_like( return tmu_like_stat -def qmu(mu, data, pdf, init_pars, par_bounds, fixed_params): +def qmu(mu, data, pdf, init_pars, par_bounds, fixed_params, return_fitted_pars=False): r""" The test statistic, :math:`q_{\mu}`, for establishing an upper limit on the strength parameter, :math:`\mu`, as defiend in @@ -94,6 +98,11 @@ def qmu(mu, data, pdf, init_pars, par_bounds, fixed_params): >>> pyhf.infer.test_statistics.qmu(test_mu, data, model, init_pars, par_bounds, fixed_params) array(3.9549891) + Access the best-fit parameter tensors: + + >>> pyhf.infer.test_statistics.qmu(test_mu, data, model, init_pars, par_bounds, fixed_params, return_fitted_pars = True) + (array(3.9549891), (array([1. , 0.97224597, 0.87553894]), array([-0.06679525, 1.00555369, 0.96930896]))) + Args: mu (Number or Tensor): The signal strength parameter data (Tensor): The data to be considered @@ -104,9 +113,18 @@ def qmu(mu, data, pdf, init_pars, par_bounds, fixed_params): The shape should be ``(n, 2)`` for ``n`` model parameters. fixed_params (:obj:`list` of :obj:`bool`): The flag to set a parameter constant to its starting value during minimization. + return_fitted_pars (:obj:`bool`): Return the best-fit parameter tensors + the fixed-POI and unconstrained fits have converged on + (i.e. :math:`\mu, \hat{\hat{\theta}}` and :math:`\hat{\mu}, \hat{\theta}`) Returns: - Float: The calculated test statistic, :math:`q_{\mu}` + Tuple of a Float and a Tuple of Tensors: + + - The calculated test statistic, :math:`q_{\mu}` + + - The parameter tensors corresponding to the constrained and unconstrained best fit, + :math:`\mu, \hat{\hat{\theta}}` and :math:`\hat{\mu}, \hat{\theta}`. + Only returned if ``return_fitted_pars`` is ``True``. """ if pdf.config.poi_index is None: raise UnspecifiedPOI( @@ -117,10 +135,20 @@ def qmu(mu, data, pdf, init_pars, par_bounds, fixed_params): 'qmu test statistic used for fit configuration with POI bounded at zero.\n' + 'Use the qmu_tilde test statistic (pyhf.infer.test_statistics.qmu_tilde) instead.' ) - return _qmu_like(mu, data, pdf, init_pars, par_bounds, fixed_params) + return _qmu_like( + mu, + data, + pdf, + init_pars, + par_bounds, + fixed_params, + return_fitted_pars=return_fitted_pars, + ) -def qmu_tilde(mu, data, pdf, init_pars, par_bounds, fixed_params): +def qmu_tilde( + mu, data, pdf, init_pars, par_bounds, fixed_params, return_fitted_pars=False +): r""" The "alternative" test statistic, :math:`\tilde{q}_{\mu}`, for establishing an upper limit on the strength parameter, :math:`\mu`, for models with @@ -163,6 +191,11 @@ def qmu_tilde(mu, data, pdf, init_pars, par_bounds, fixed_params): >>> pyhf.infer.test_statistics.qmu_tilde(test_mu, data, model, init_pars, par_bounds, fixed_params) array(3.93824492) + Access the best-fit parameter tensors: + + >>> pyhf.infer.test_statistics.qmu_tilde(test_mu, data, model, init_pars, par_bounds, fixed_params, return_fitted_pars = True) + (array(3.93824492), (array([1. , 0.97224597, 0.87553894]), array([0. , 1.0030512 , 0.96266961]))) + Args: mu (Number or Tensor): The signal strength parameter data (:obj:`tensor`): The data to be considered @@ -173,9 +206,18 @@ def qmu_tilde(mu, data, pdf, init_pars, par_bounds, fixed_params): The shape should be ``(n, 2)`` for ``n`` model parameters. fixed_params (:obj:`list` of :obj:`bool`): The flag to set a parameter constant to its starting value during minimization. + return_fitted_pars (:obj:`bool`): Return the best-fit parameter tensors + the fixed-POI and unconstrained fits have converged on + (i.e. :math:`\mu, \hat{\hat{\theta}}` and :math:`\hat{\mu}, \hat{\theta}`) Returns: - Float: The calculated test statistic, :math:`\tilde{q}_{\mu}` + Tuple of a Float and a Tuple of Tensors: + + - The calculated test statistic, :math:`\tilde{q}_{\mu}` + + - The parameter tensors corresponding to the constrained and unconstrained best fit, + :math:`\mu, \hat{\hat{\theta}}` and :math:`\hat{\mu}, \hat{\theta}`. + Only returned if ``return_fitted_pars`` is ``True``. """ if pdf.config.poi_index is None: raise UnspecifiedPOI( @@ -186,10 +228,18 @@ def qmu_tilde(mu, data, pdf, init_pars, par_bounds, fixed_params): 'qmu_tilde test statistic used for fit configuration with POI not bounded at zero.\n' + 'Use the qmu test statistic (pyhf.infer.test_statistics.qmu) instead.' ) - return _qmu_like(mu, data, pdf, init_pars, par_bounds, fixed_params) + return _qmu_like( + mu, + data, + pdf, + init_pars, + par_bounds, + fixed_params, + return_fitted_pars=return_fitted_pars, + ) -def tmu(mu, data, pdf, init_pars, par_bounds, fixed_params): +def tmu(mu, data, pdf, init_pars, par_bounds, fixed_params, return_fitted_pars=False): r""" The test statistic, :math:`t_{\mu}`, for establishing a two-sided interval on the strength parameter, :math:`\mu`, as defiend in Equation (8) @@ -221,6 +271,11 @@ def tmu(mu, data, pdf, init_pars, par_bounds, fixed_params): >>> pyhf.infer.test_statistics.tmu(test_mu, data, model, init_pars, par_bounds, fixed_params) array(3.9549891) + Access the best-fit parameter tensors: + + >>> pyhf.infer.test_statistics.tmu(test_mu, data, model, init_pars, par_bounds, fixed_params, return_fitted_pars = True) + (array(3.9549891), (array([1. , 0.97224597, 0.87553894]), array([-0.06679525, 1.00555369, 0.96930896]))) + Args: mu (Number or Tensor): The signal strength parameter data (Tensor): The data to be considered @@ -231,9 +286,18 @@ def tmu(mu, data, pdf, init_pars, par_bounds, fixed_params): The shape should be ``(n, 2)`` for ``n`` model parameters. fixed_params (:obj:`list` of :obj:`bool`): The flag to set a parameter constant to its starting value during minimization. + return_fitted_pars (:obj:`bool`): Return the best-fit parameter tensors + the fixed-POI and unconstrained fits have converged on + (i.e. :math:`\mu, \hat{\hat{\theta}}` and :math:`\hat{\mu}, \hat{\theta}`) Returns: - Float: The calculated test statistic, :math:`t_{\mu}` + Tuple of a Float and a Tuple of Tensors: + + - The calculated test statistic, :math:`t_{\mu}` + + - The parameter tensors corresponding to the constrained and unconstrained best fit, + :math:`\mu, \hat{\hat{\theta}}` and :math:`\hat{\mu}, \hat{\theta}`. + Only returned if ``return_fitted_pars`` is ``True``. """ if pdf.config.poi_index is None: raise UnspecifiedPOI( @@ -244,10 +308,20 @@ def tmu(mu, data, pdf, init_pars, par_bounds, fixed_params): 'tmu test statistic used for fit configuration with POI bounded at zero.\n' + 'Use the tmu_tilde test statistic (pyhf.infer.test_statistics.tmu_tilde) instead.' ) - return _tmu_like(mu, data, pdf, init_pars, par_bounds, fixed_params) + return _tmu_like( + mu, + data, + pdf, + init_pars, + par_bounds, + fixed_params, + return_fitted_pars=return_fitted_pars, + ) -def tmu_tilde(mu, data, pdf, init_pars, par_bounds, fixed_params): +def tmu_tilde( + mu, data, pdf, init_pars, par_bounds, fixed_params, return_fitted_pars=False +): r""" The test statistic, :math:`\tilde{t}_{\mu}`, for establishing a two-sided interval on the strength parameter, :math:`\mu`, for models with @@ -270,6 +344,7 @@ def tmu_tilde(mu, data, pdf, init_pars, par_bounds, fixed_params): \end{equation} Example: + >>> import pyhf >>> pyhf.set_backend("numpy") >>> model = pyhf.simplemodels.uncorrelated_background( @@ -284,6 +359,11 @@ def tmu_tilde(mu, data, pdf, init_pars, par_bounds, fixed_params): >>> pyhf.infer.test_statistics.tmu_tilde(test_mu, data, model, init_pars, par_bounds, fixed_params) array(3.93824492) + Access the best-fit parameter tensors: + + >>> pyhf.infer.test_statistics.tmu_tilde(test_mu, data, model, init_pars, par_bounds, fixed_params, return_fitted_pars = True) + (array(3.93824492), (array([1. , 0.97224597, 0.87553894]), array([0. , 1.0030512 , 0.96266961]))) + Args: mu (Number or Tensor): The signal strength parameter data (:obj:`tensor`): The data to be considered @@ -294,9 +374,18 @@ def tmu_tilde(mu, data, pdf, init_pars, par_bounds, fixed_params): The shape should be ``(n, 2)`` for ``n`` model parameters. fixed_params (:obj:`list` of :obj:`bool`): The flag to set a parameter constant to its starting value during minimization. + return_fitted_pars (:obj:`bool`): Return the best-fit parameter tensors + the fixed-POI and unconstrained fits have converged on + (i.e. :math:`\mu, \hat{\hat{\theta}}` and :math:`\hat{\mu}, \hat{\theta}`) Returns: - Float: The calculated test statistic, :math:`\tilde{t}_{\mu}` + Tuple of a Float and a Tuple of Tensors: + + - The calculated test statistic, :math:`\tilde{t}_{\mu}` + + - The parameter tensors corresponding to the constrained and unconstrained best fit, + :math:`\mu, \hat{\hat{\theta}}` and :math:`\hat{\mu}, \hat{\theta}`. + Only returned if ``return_fitted_pars`` is ``True``. """ if pdf.config.poi_index is None: raise UnspecifiedPOI( @@ -307,10 +396,18 @@ def tmu_tilde(mu, data, pdf, init_pars, par_bounds, fixed_params): 'tmu_tilde test statistic used for fit configuration with POI not bounded at zero.\n' + 'Use the tmu test statistic (pyhf.infer.test_statistics.tmu) instead.' ) - return _tmu_like(mu, data, pdf, init_pars, par_bounds, fixed_params) + return _tmu_like( + mu, + data, + pdf, + init_pars, + par_bounds, + fixed_params, + return_fitted_pars=return_fitted_pars, + ) -def q0(mu, data, pdf, init_pars, par_bounds, fixed_params): +def q0(mu, data, pdf, init_pars, par_bounds, fixed_params, return_fitted_pars=False): r""" The test statistic, :math:`q_{0}`, for discovery of a positive signal as defined in Equation (12) in :xref:`arXiv:1007.1727`, for :math:`\mu=0`. @@ -340,6 +437,11 @@ def q0(mu, data, pdf, init_pars, par_bounds, fixed_params): >>> pyhf.infer.test_statistics.q0(test_mu, data, model, init_pars, par_bounds, fixed_params) array(2.98339447) + Access the best-fit parameter tensors: + + >>> pyhf.infer.test_statistics.q0(test_mu, data, model, init_pars, par_bounds, fixed_params, return_fitted_pars = True) + (array(2.98339447), (array([0. , 1.03050845, 1.12128752]), array([0.95260667, 0.99635345, 1.02140172]))) + Args: mu (Number or Tensor): The signal strength parameter (must be set to zero) data (Tensor): The data to be considered @@ -350,9 +452,18 @@ def q0(mu, data, pdf, init_pars, par_bounds, fixed_params): The shape should be ``(n, 2)`` for ``n`` model parameters. fixed_params (:obj:`list` of :obj:`bool`): The flag to set a parameter constant to its starting value during minimization. + return_fitted_pars (:obj:`bool`): Return the best-fit parameter tensors + the fixed-POI and unconstrained fits have converged on + (i.e. :math:`\mu, \hat{\hat{\theta}}` and :math:`\hat{\mu}, \hat{\theta}`) Returns: - Float: The calculated test statistic, :math:`q_{0}` + Tuple of a Float and a Tuple of Tensors: + + - The calculated test statistic, :math:`q_{0}` + + - The parameter tensors corresponding to the constrained and unconstrained best fit, + :math:`\mu, \hat{\hat{\theta}}` and :math:`\hat{\mu}, \hat{\theta}`. + Only returned if ``return_fitted_pars`` is ``True``. """ if pdf.config.poi_index is None: @@ -367,10 +478,12 @@ def q0(mu, data, pdf, init_pars, par_bounds, fixed_params): tensorlib, optimizer = get_backend() - tmu_like_stat, (_, muhatbhat) = _tmu_like( + tmu_like_stat, (mubhathat, muhatbhat) = _tmu_like( mu, data, pdf, init_pars, par_bounds, fixed_params, return_fitted_pars=True ) q0_stat = tensorlib.where( muhatbhat[pdf.config.poi_index] < 0, tensorlib.astensor(0.0), tmu_like_stat ) + if return_fitted_pars: + return q0_stat, (mubhathat, muhatbhat) return q0_stat diff --git a/tests/test_calculator.py b/tests/test_calculator.py index 65021a89dc..97b5fa4e3d 100644 --- a/tests/test_calculator.py +++ b/tests/test_calculator.py @@ -1,3 +1,5 @@ +import pytest + import pyhf import pyhf.infer.calculators @@ -5,3 +7,83 @@ def test_calc_dist(): asymptotic_dist = pyhf.infer.calculators.AsymptoticTestStatDistribution(0.0) assert asymptotic_dist.pvalue(-1) == 1 - asymptotic_dist.cdf(-1) + + +@pytest.mark.parametrize("return_fitted_pars", [False, True]) +def test_generate_asimov_can_return_fitted_pars(return_fitted_pars): + model = pyhf.simplemodels.uncorrelated_background([1, 1], [1, 1], [1, 1]) + data = [2, 2, 1, 1] # [main x 2, aux x 2] + init_pars = model.config.suggested_init() + par_bounds = model.config.suggested_bounds() + fixed_params = model.config.suggested_fixed() + + result = pyhf.infer.calculators.generate_asimov_data( + 1.0, + data, + model, + init_pars, + par_bounds, + fixed_params, + return_fitted_pars=return_fitted_pars, + ) + + if return_fitted_pars: + assert len(result) == 2 + result, asimov_pars = result + assert pytest.approx([1.0, 1.0, 1.0]) == pyhf.tensorlib.tolist(asimov_pars) + assert pytest.approx([2.0, 2.0, 1.0, 1.0]) == pyhf.tensorlib.tolist(result) + + +# test different test stats because those affect the control flow +# in AsymptotiCalculator.teststatistic, where the fit results should be set +# the other kwargs don't impact the logic of that method, +# so leave them at the default so as not to put a burden on future changes +@pytest.mark.parametrize('test_stat', ['qtilde', 'q', 'q0']) +def test_asymptotic_calculator_has_fitted_pars(test_stat): + model = pyhf.simplemodels.uncorrelated_background([1], [1], [1]) + data = [2, 1] # [main, aux] + + calc = pyhf.infer.calculators.AsymptoticCalculator(data, model, test_stat=test_stat) + calc.teststatistic(0 if test_stat == 'q0' else 1) + + assert hasattr(calc, 'fitted_pars') + fitted_pars = calc.fitted_pars + assert hasattr(fitted_pars, 'asimov_pars') + assert hasattr(fitted_pars, 'fixed_poi_fit_to_data') + assert hasattr(fitted_pars, 'fixed_poi_fit_to_asimov') + assert hasattr(fitted_pars, 'free_fit_to_data') + assert hasattr(fitted_pars, 'free_fit_to_asimov') + + rtol = 1e-5 + if test_stat == 'q0': + assert pytest.approx([1.0, 1.0], rel=rtol) == pyhf.tensorlib.tolist( + fitted_pars.asimov_pars + ) + assert pytest.approx([0.0, 1.5], rel=rtol) == pyhf.tensorlib.tolist( + fitted_pars.fixed_poi_fit_to_data + ) + assert pytest.approx([0.0, 1.5], rel=rtol) == pyhf.tensorlib.tolist( + fitted_pars.fixed_poi_fit_to_asimov + ) + assert pytest.approx([1.0, 1.0], rel=rtol) == pyhf.tensorlib.tolist( + fitted_pars.free_fit_to_data + ) + assert pytest.approx([1.0, 1.0], rel=rtol) == pyhf.tensorlib.tolist( + fitted_pars.free_fit_to_asimov + ) + else: + assert pytest.approx([0.0, 1.5], rel=rtol) == pyhf.tensorlib.tolist( + fitted_pars.asimov_pars + ) + assert pytest.approx([1.0, 1.0], rel=rtol) == pyhf.tensorlib.tolist( + fitted_pars.fixed_poi_fit_to_data + ) + assert pytest.approx([1.0, 1.1513553], rel=rtol) == pyhf.tensorlib.tolist( + fitted_pars.fixed_poi_fit_to_asimov + ) + assert pytest.approx([1.0, 1.0], rel=rtol) == pyhf.tensorlib.tolist( + fitted_pars.free_fit_to_data + ) + assert pytest.approx( + [7.6470499e-05, 1.4997178], rel=rtol + ) == pyhf.tensorlib.tolist(fitted_pars.free_fit_to_asimov) diff --git a/tests/test_infer.py b/tests/test_infer.py index 32f10277c4..520c4a9356 100644 --- a/tests/test_infer.py +++ b/tests/test_infer.py @@ -173,6 +173,59 @@ def test_hypotest_return_expected_set(tmpdir, hypotest_args, test_stat): assert check_uniform_type(result[3]) +@pytest.mark.parametrize( + 'calctype,kwargs,expected_type', + [ + ('asymptotics', {}, pyhf.infer.calculators.AsymptoticCalculator), + ('toybased', dict(ntoys=1), pyhf.infer.calculators.ToyCalculator), + ], +) +@pytest.mark.parametrize('return_tail_probs', [True, False]) +@pytest.mark.parametrize('return_expected', [True, False]) +@pytest.mark.parametrize('return_expected_set', [True, False]) +def test_hypotest_return_calculator( + tmpdir, + hypotest_args, + calctype, + kwargs, + expected_type, + return_tail_probs, + return_expected, + return_expected_set, +): + """ + Check that the return structure of pyhf.infer.hypotest with the + additon of the return_calculator keyword arg is as expected + """ + *_, model = hypotest_args + + # only those return flags where the toggled return value + # is placed in front of the calculator in the returned tuple + extra_returns = sum( + int(return_flag) + for return_flag in ( + return_tail_probs, + return_expected, + return_expected_set, + ) + ) + + result = pyhf.infer.hypotest( + *hypotest_args, + return_calculator=True, + return_tail_probs=return_tail_probs, + return_expected=return_expected, + return_expected_set=return_expected_set, + calctype=calctype, + **kwargs, + ) + + assert len(list(result)) == 2 + extra_returns + # not *_, calc = result b.c. in future, there could be additional optional returns + calc = result[1 + extra_returns] + assert isinstance(calc, expected_type) + + @pytest.mark.parametrize( "kwargs", [{'calctype': 'asymptotics'}, {'calctype': 'toybased', 'ntoys': 5}], diff --git a/tests/test_teststats.py b/tests/test_teststats.py index 640e8a8794..f3681b2fd8 100644 --- a/tests/test_teststats.py +++ b/tests/test_teststats.py @@ -172,3 +172,40 @@ def test_get_teststat_by_name(test_stat): def test_get_teststat_error(): with pytest.raises(pyhf.exceptions.InvalidTestStatistic): pyhf.infer.utils.get_test_stat("look at me i'm not real") + + +@pytest.mark.parametrize("return_fitted_pars", [False, True]) +@pytest.mark.parametrize( + "test_stat", + [ + pyhf.infer.test_statistics.q0, + pyhf.infer.test_statistics.qmu, + pyhf.infer.test_statistics.qmu_tilde, + pyhf.infer.test_statistics.tmu, + pyhf.infer.test_statistics.tmu_tilde, + ], +) +def test_return_fitted_pars(test_stat, return_fitted_pars): + mu = 0.0 if test_stat is pyhf.infer.test_statistics.q0 else 1.0 + model = pyhf.simplemodels.uncorrelated_background([6], [9], [3]) + data = [9] + model.config.auxdata + init_pars = model.config.suggested_init() + par_bounds = model.config.suggested_bounds() + fixed_params = model.config.suggested_fixed() + + result = test_stat( + mu, + data, + model, + init_pars, + par_bounds, + fixed_params, + return_fitted_pars=return_fitted_pars, + ) + if return_fitted_pars: + assert len(result) == 2 + assert len(result[1]) == 2 + result, (pars_bestfit, pars_constrained_fit) = result + assert len(pars_bestfit) == len(init_pars) + assert len(pars_constrained_fit) == len(init_pars) + assert result > -1e4 # >= 0 but with generous tolerance