From c31919ba240647a7cb23d2c313f8091fbf9297ad Mon Sep 17 00:00:00 2001 From: Carlo Gottardo Date: Sun, 26 Apr 2020 14:24:03 +0200 Subject: [PATCH 01/24] fixed_vals add as argument to infer functions --- src/pyhf/infer/mle.py | 17 +++++++++++++---- src/pyhf/infer/test_statistics.py | 30 ++++++++++++++++++------------ 2 files changed, 31 insertions(+), 16 deletions(-) diff --git a/src/pyhf/infer/mle.py b/src/pyhf/infer/mle.py index 1391c67e02..82448d5958 100644 --- a/src/pyhf/infer/mle.py +++ b/src/pyhf/infer/mle.py @@ -47,7 +47,7 @@ def twice_nll(pars, data, pdf): return -2 * pdf.logpdf(pars, data) -def fit(data, pdf, init_pars=None, par_bounds=None, **kwargs): +def fit(data, pdf, init_pars=None, par_bounds=None, fixed_vals=None, **kwargs): r""" Run a unconstrained maximum likelihood fit. This is done by minimizing the objective function :func:`~pyhf.infer.mle.twice_nll` @@ -96,10 +96,14 @@ def fit(data, pdf, init_pars=None, par_bounds=None, **kwargs): _, opt = get_backend() init_pars = init_pars or pdf.config.suggested_init() par_bounds = par_bounds or pdf.config.suggested_bounds() - return opt.minimize(twice_nll, data, pdf, init_pars, par_bounds, **kwargs) + return opt.minimize( + twice_nll, data, pdf, init_pars, par_bounds, fixed_vals, **kwargs + ) -def fixed_poi_fit(poi_val, data, pdf, init_pars=None, par_bounds=None, **kwargs): +def fixed_poi_fit( + poi_val, data, pdf, init_pars=None, par_bounds=None, fixed_vals=None, **kwargs +): r""" Run a maximum likelihood fit with the POI value fixed. This is done by minimizing the objective function of :func:`~pyhf.infer.mle.twice_nll` @@ -152,6 +156,11 @@ def fixed_poi_fit(poi_val, data, pdf, init_pars=None, par_bounds=None, **kwargs) raise UnspecifiedPOI( 'No POI is defined. A POI is required to fit with a fixed POI.' ) + + fixed_params = [(pdf.config.poi_index, poi_val)] + if fixed_vals: + fixed_params += fixed_vals + _, opt = get_backend() init_pars = init_pars or pdf.config.suggested_init() par_bounds = par_bounds or pdf.config.suggested_bounds() @@ -161,6 +170,6 @@ def fixed_poi_fit(poi_val, data, pdf, init_pars=None, par_bounds=None, **kwargs) pdf, init_pars, par_bounds, - [(pdf.config.poi_index, poi_val)], + fixed_params, **kwargs, ) diff --git a/src/pyhf/infer/test_statistics.py b/src/pyhf/infer/test_statistics.py index c636d04ab6..d1cd663ecd 100644 --- a/src/pyhf/infer/test_statistics.py +++ b/src/pyhf/infer/test_statistics.py @@ -7,7 +7,7 @@ log = logging.getLogger(__name__) -def _qmu_like(mu, data, pdf, init_pars, par_bounds): +def _qmu_like(mu, data, pdf, init_pars, par_bounds, fixed_vals): """ Clipped version of _tmu_like where the returned test statistic is 0 if muhat > 0 else tmu_like_stat. @@ -25,7 +25,9 @@ def _qmu_like(mu, data, pdf, init_pars, par_bounds): return qmu_like_stat -def _tmu_like(mu, data, pdf, init_pars, par_bounds, return_fitted_pars=False): +def _tmu_like( + mu, data, pdf, init_pars, par_bounds, fixed_vals, return_fitted_pars=False +): """ Basic Profile Likelihood test statistic. @@ -34,10 +36,10 @@ def _tmu_like(mu, data, pdf, init_pars, par_bounds, return_fitted_pars=False): """ tensorlib, optimizer = get_backend() mubhathat, fixed_poi_fit_lhood_val = fixed_poi_fit( - mu, data, pdf, init_pars, par_bounds, return_fitted_val=True + mu, data, pdf, init_pars, par_bounds, fixed_vals, return_fitted_val=True ) muhatbhat, unconstrained_fit_lhood_val = fit( - data, pdf, init_pars, par_bounds, return_fitted_val=True + data, pdf, init_pars, par_bounds, fixed_vals, return_fitted_val=True ) log_likelihood_ratio = fixed_poi_fit_lhood_val - unconstrained_fit_lhood_val tmu_like_stat = tensorlib.astensor( @@ -48,7 +50,7 @@ def _tmu_like(mu, data, pdf, init_pars, par_bounds, return_fitted_pars=False): return tmu_like_stat -def qmu(mu, data, pdf, init_pars, par_bounds): +def qmu(mu, data, pdf, init_pars, par_bounds, fixed_vals): r""" The test statistic, :math:`q_{\mu}`, for establishing an upper limit on the strength parameter, :math:`\mu`, as defiend in @@ -91,6 +93,7 @@ def qmu(mu, data, pdf, init_pars, par_bounds): pdf (~pyhf.pdf.Model): The HistFactory statistical model used in the likelihood ratio calculation init_pars (`list`): Values to initialize the model parameters at for the fit par_bounds (`list` of `list`\s or `tuple`\s): The extrema of values the model parameters are allowed to reach in the fit + fixed_vals(`list`): Parameters held constant in the fit Returns: Float: The calculated test statistic, :math:`q_{\mu}` @@ -104,10 +107,10 @@ def qmu(mu, data, pdf, init_pars, par_bounds): '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) + return _qmu_like(mu, data, pdf, init_pars, par_bounds, fixed_vals) -def qmu_tilde(mu, data, pdf, init_pars, par_bounds): +def qmu_tilde(mu, data, pdf, init_pars, par_bounds, fixed_vals): r""" The test statistic, :math:`\tilde{q}_{\mu}`, for establishing an upper limit on the strength parameter, :math:`\mu`, for models with @@ -155,6 +158,7 @@ def qmu_tilde(mu, data, pdf, init_pars, par_bounds): pdf (~pyhf.pdf.Model): The statistical model adhering to the schema model.json init_pars (`list`): Values to initialize the model parameters at for the fit par_bounds (`list` of `list`\s or `tuple`\s): The extrema of values the model parameters are allowed to reach in the fit + fixed_vals(`list`): Parameters held constant in the fit Returns: Float: The calculated test statistic, :math:`\tilde{q}_{\mu}` @@ -168,10 +172,10 @@ def qmu_tilde(mu, data, pdf, init_pars, par_bounds): '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) + return _qmu_like(mu, data, pdf, init_pars, par_bounds, fixed_vals) -def tmu(mu, data, pdf, init_pars, par_bounds): +def tmu(mu, data, pdf, init_pars, par_bounds, fixed_vals): r""" The test statistic, :math:`t_{\mu}`, for establishing a two-sided interval on the strength parameter, :math:`\mu`, as defiend in Equation (8) @@ -208,6 +212,7 @@ def tmu(mu, data, pdf, init_pars, par_bounds): pdf (~pyhf.pdf.Model): The statistical model adhering to the schema model.json init_pars (`list`): Values to initialize the model parameters at for the fit par_bounds (`list` of `list`\s or `tuple`\s): The extrema of values the model parameters are allowed to reach in the fit + fixed_vals(`list`): Parameters held constant in the fit Returns: Float: The calculated test statistic, :math:`t_{\mu}` @@ -221,10 +226,10 @@ def tmu(mu, data, pdf, init_pars, par_bounds): '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) + return _tmu_like(mu, data, pdf, init_pars, par_bounds, fixed_vals) -def tmu_tilde(mu, data, pdf, init_pars, par_bounds): +def tmu_tilde(mu, data, pdf, init_pars, par_bounds, fixed_vals): r""" The test statistic, :math:`\tilde{t}_{\mu}`, for establishing a two-sided interval on the strength parameter, :math:`\mu`, for models with @@ -266,6 +271,7 @@ def tmu_tilde(mu, data, pdf, init_pars, par_bounds): pdf (~pyhf.pdf.Model): The statistical model adhering to the schema model.json init_pars (`list`): Values to initialize the model parameters at for the fit par_bounds (`list` of `list`\s or `tuple`\s): The extrema of values the model parameters are allowed to reach in the fit + fixed_vals(`list`): Parameters held constant in the fit Returns: Float: The calculated test statistic, :math:`\tilde{t}_{\mu}` @@ -279,4 +285,4 @@ def tmu_tilde(mu, data, pdf, init_pars, par_bounds): '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) + return _tmu_like(mu, data, pdf, init_pars, par_bounds, fixed_vals) From 80dafdb320a6a0ed6a3a324fe253f1093b5d6d91 Mon Sep 17 00:00:00 2001 From: Carlo Gottardo Date: Sun, 26 Apr 2020 18:54:53 +0200 Subject: [PATCH 02/24] fixed_vals add as argument to infer functions --- src/pyhf/infer/calculators.py | 27 ++++++++++++++++++++++----- 1 file changed, 22 insertions(+), 5 deletions(-) diff --git a/src/pyhf/infer/calculators.py b/src/pyhf/infer/calculators.py index 62f90329ab..528a3ba023 100644 --- a/src/pyhf/infer/calculators.py +++ b/src/pyhf/infer/calculators.py @@ -12,7 +12,7 @@ from .test_statistics import qmu, qmu_tilde -def generate_asimov_data(asimov_mu, data, pdf, init_pars, par_bounds): +def generate_asimov_data(asimov_mu, data, pdf, init_pars, par_bounds, fixed_vals): """ Compute Asimov Dataset (expected yields at best-fit values) for a given POI value. @@ -22,12 +22,15 @@ def generate_asimov_data(asimov_mu, data, pdf, init_pars, par_bounds): pdf (~pyhf.pdf.Model): The statistical model adhering to the schema ``model.json``. init_pars (`tensor`): The initial parameter values to be used for fitting. par_bounds (`tensor`): The parameter value bounds to be used for fitting. + fixed_vals (`tensor`): Parameters to be held constant in the fit. Returns: Tensor: The Asimov dataset. """ - bestfit_nuisance_asimov = fixed_poi_fit(asimov_mu, data, pdf, init_pars, par_bounds) + bestfit_nuisance_asimov = fixed_poi_fit( + asimov_mu, data, pdf, init_pars, par_bounds, fixed_vals + ) return pdf.expected_data(bestfit_nuisance_asimov) @@ -118,7 +121,9 @@ def expected_value(self, nsigma): class AsymptoticCalculator(object): """The Asymptotic Calculator.""" - def __init__(self, data, pdf, init_pars=None, par_bounds=None, qtilde=False): + def __init__( + self, data, pdf, init_pars=None, par_bounds=None, fixed_vals=None, qtilde=False + ): """ Asymptotic Calculator. @@ -127,6 +132,7 @@ def __init__(self, data, pdf, init_pars=None, par_bounds=None, qtilde=False): pdf (~pyhf.pdf.Model): The statistical model adhering to the schema ``model.json``. init_pars (`tensor`): The initial parameter values to be used for fitting. par_bounds (`tensor`): The parameter value bounds to be used for fitting. + fixed_vals (`tensor`): Parameters to be held constant in the fit. qtilde (`bool`): Whether to use qtilde as the test statistic. Returns: @@ -137,6 +143,7 @@ def __init__(self, data, pdf, init_pars=None, par_bounds=None, qtilde=False): self.pdf = pdf self.init_pars = init_pars or pdf.config.suggested_init() self.par_bounds = par_bounds or pdf.config.suggested_bounds() + self.fixed_vals = fixed_vals self.qtilde = qtilde self.sqrtqmuA_v = None @@ -173,7 +180,12 @@ def teststatistic(self, poi_test): teststat_func = qmu_tilde if self.qtilde else qmu qmu_v = teststat_func( - poi_test, self.data, self.pdf, self.init_pars, self.par_bounds + poi_test, + self.data, + self.pdf, + self.init_pars, + self.par_bounds, + self.fixed_vals, ) sqrtqmu_v = tensorlib.sqrt(qmu_v) @@ -182,7 +194,12 @@ def teststatistic(self, poi_test): asimov_mu, self.data, self.pdf, self.init_pars, self.par_bounds ) qmuA_v = teststat_func( - poi_test, asimov_data, self.pdf, self.init_pars, self.par_bounds + poi_test, + asimov_data, + self.pdf, + self.init_pars, + self.par_bounds, + self.fixed_vals, ) self.sqrtqmuA_v = tensorlib.sqrt(qmuA_v) From 32a7955f0abbe2717c10d8f1de01b00a45e10ab9 Mon Sep 17 00:00:00 2001 From: Carlo Gottardo Date: Sun, 26 Apr 2020 18:56:40 +0200 Subject: [PATCH 03/24] fixed_vals add as argument to infer functions --- src/pyhf/infer/__init__.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/pyhf/infer/__init__.py b/src/pyhf/infer/__init__.py index 246d2df3b0..61700b9d84 100644 --- a/src/pyhf/infer/__init__.py +++ b/src/pyhf/infer/__init__.py @@ -6,7 +6,7 @@ def hypotest( - poi_test, data, pdf, init_pars=None, par_bounds=None, qtilde=False, **kwargs + poi_test, data, pdf, init_pars=None, par_bounds=None, fixed_vals=None, qtilde=False, **kwargs ): r""" Compute :math:`p`-values and test statistics for a single value of the parameter of interest. @@ -34,6 +34,7 @@ def hypotest( pdf (~pyhf.pdf.Model): The statistical model adhering to the schema ``model.json`` init_pars (Array or Tensor): The initial parameter values to be used for minimization par_bounds (Array or Tensor): The parameter value bounds to be used for minimization + fixed_vals (list of tuples): Parameters to be held constant and their value qtilde (Bool): When ``True`` perform the calculation using the alternative test statistic, :math:`\tilde{q}_{\mu}`, as defined under the Wald approximation in Equation (62) of :xref:`arXiv:1007.1727`. @@ -120,7 +121,7 @@ def hypotest( par_bounds = par_bounds or pdf.config.suggested_bounds() tensorlib, _ = get_backend() - calc = AsymptoticCalculator(data, pdf, init_pars, par_bounds, qtilde=qtilde) + calc = AsymptoticCalculator(data, pdf, init_pars, par_bounds, fixed_vals, qtilde=qtilde) teststat = calc.teststatistic(poi_test) sig_plus_bkg_distribution, b_only_distribution = calc.distributions(poi_test) From eb90ed5885f9afc395ef607f4ebcfdbc9ff77c9d Mon Sep 17 00:00:00 2001 From: Carlo Gottardo Date: Sun, 26 Apr 2020 22:24:24 +0200 Subject: [PATCH 04/24] tested --- src/pyhf/infer/calculators.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pyhf/infer/calculators.py b/src/pyhf/infer/calculators.py index 528a3ba023..b6e1c1404b 100644 --- a/src/pyhf/infer/calculators.py +++ b/src/pyhf/infer/calculators.py @@ -191,7 +191,7 @@ def teststatistic(self, poi_test): asimov_mu = 0.0 asimov_data = generate_asimov_data( - asimov_mu, self.data, self.pdf, self.init_pars, self.par_bounds + asimov_mu, self.data, self.pdf, self.init_pars, self.par_bounds, self.fixed_vals ) qmuA_v = teststat_func( poi_test, From bf72f330b01afa97893a95ff99c306e8f0b4e919 Mon Sep 17 00:00:00 2001 From: Carlo Gottardo Date: Sun, 26 Apr 2020 23:28:31 +0200 Subject: [PATCH 05/24] black formatting --- src/pyhf/infer/__init__.py | 13 +++++++++++-- src/pyhf/infer/calculators.py | 7 ++++++- 2 files changed, 17 insertions(+), 3 deletions(-) diff --git a/src/pyhf/infer/__init__.py b/src/pyhf/infer/__init__.py index 61700b9d84..89c2d24773 100644 --- a/src/pyhf/infer/__init__.py +++ b/src/pyhf/infer/__init__.py @@ -6,7 +6,14 @@ def hypotest( - poi_test, data, pdf, init_pars=None, par_bounds=None, fixed_vals=None, qtilde=False, **kwargs + poi_test, + data, + pdf, + init_pars=None, + par_bounds=None, + fixed_vals=None, + qtilde=False, + **kwargs, ): r""" Compute :math:`p`-values and test statistics for a single value of the parameter of interest. @@ -121,7 +128,9 @@ def hypotest( par_bounds = par_bounds or pdf.config.suggested_bounds() tensorlib, _ = get_backend() - calc = AsymptoticCalculator(data, pdf, init_pars, par_bounds, fixed_vals, qtilde=qtilde) + calc = AsymptoticCalculator( + data, pdf, init_pars, par_bounds, fixed_vals, qtilde=qtilde + ) teststat = calc.teststatistic(poi_test) sig_plus_bkg_distribution, b_only_distribution = calc.distributions(poi_test) diff --git a/src/pyhf/infer/calculators.py b/src/pyhf/infer/calculators.py index b6e1c1404b..fb07d089e6 100644 --- a/src/pyhf/infer/calculators.py +++ b/src/pyhf/infer/calculators.py @@ -191,7 +191,12 @@ def teststatistic(self, poi_test): asimov_mu = 0.0 asimov_data = generate_asimov_data( - asimov_mu, self.data, self.pdf, self.init_pars, self.par_bounds, self.fixed_vals + asimov_mu, + self.data, + self.pdf, + self.init_pars, + self.par_bounds, + self.fixed_vals, ) qmuA_v = teststat_func( poi_test, From 5f8e16ff1b3173e2be606dae7adf4b5a2571f830 Mon Sep 17 00:00:00 2001 From: Carlo Gottardo Date: Sun, 26 Apr 2020 23:56:49 +0200 Subject: [PATCH 06/24] added fixed_vals positional argument to test backends --- tests/test_backend_consistency.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/tests/test_backend_consistency.py b/tests/test_backend_consistency.py index 61f9b8579c..86fd3578fe 100644 --- a/tests/test_backend_consistency.py +++ b/tests/test_backend_consistency.py @@ -115,7 +115,12 @@ def test_hypotest_qmu_tilde( pyhf.set_backend(backend) qmu_tilde = pyhf.infer.test_statistics.qmu_tilde( - 1.0, data, pdf, pdf.config.suggested_init(), pdf.config.suggested_bounds() + 1.0, + data, + pdf, + pdf.config.suggested_init(), + pdf.config.suggested_bounds(), + None, ) test_statistic.append(qmu_tilde) From 9437118b76125a905f695808f83b5c1abc64efb0 Mon Sep 17 00:00:00 2001 From: Carlo Gottardo Date: Mon, 27 Apr 2020 11:45:20 +0200 Subject: [PATCH 07/24] Adding NP fix and shift in test_backend_consistency --- tests/test_backend_consistency.py | 104 ++++++++++++++++++------------ 1 file changed, 64 insertions(+), 40 deletions(-) diff --git a/tests/test_backend_consistency.py b/tests/test_backend_consistency.py index 86fd3578fe..5c46f5ab24 100644 --- a/tests/test_backend_consistency.py +++ b/tests/test_backend_consistency.py @@ -75,6 +75,8 @@ def test_hypotest_qmu_tilde( """ source = generate_source_static(n_bins) + bkg_unc_10pc_up = [x + 0.10 * x for x in source['bindata']['bkg']] + bkg_unc_10pc_dn = [x - 0.10 * x for x in source['bindata']['bkg']] signal_sample = { 'name': 'signal', @@ -90,7 +92,13 @@ def test_hypotest_qmu_tilde( 'name': 'uncorr_bkguncrt', 'type': 'shapesys', 'data': source['bindata']['bkgerr'], - } + }, + { + 'name': 'norm_bkgunc', + 'type': 'histosys', + 'hi_data': bkg_unc_10pc_up, + 'lo_data': bkg_unc_10pc_dn, + }, ], } samples = ( @@ -101,6 +109,19 @@ def test_hypotest_qmu_tilde( spec = {'channels': [{'name': 'singlechannel', 'samples': samples}]} pdf = pyhf.Model(spec) + # get norm_bkgunc index to set it constatnt + norm_bkgunc_idx = -1 + for idx, par in enumerate(pdf.config.par_map): + if par == "norm_bkgunc": + norm_bkgunc_idx = idx + # Floating norm_bkgunc, fixed at nominal, plus/minus 1 sigma + param_tests = [ + None, + (norm_bkgunc_idx, 0.0), + (norm_bkgunc_idx, -1.0), + (norm_bkgunc_idx, 1.0), + ] + data = source['bindata']['data'] + pdf.config.auxdata backends = [ @@ -110,44 +131,47 @@ def test_hypotest_qmu_tilde( pyhf.tensor.jax_backend(precision='64b'), ] - test_statistic = [] - for backend in backends: - pyhf.set_backend(backend) - - qmu_tilde = pyhf.infer.test_statistics.qmu_tilde( - 1.0, - data, - pdf, - pdf.config.suggested_init(), - pdf.config.suggested_bounds(), - None, - ) - test_statistic.append(qmu_tilde) - - # compare to NumPy/SciPy - test_statistic = np.array(test_statistic) - numpy_ratio = np.divide(test_statistic, test_statistic[0]) - numpy_ratio_delta_unity = np.absolute(np.subtract(numpy_ratio, 1)) - - # compare tensor libraries to each other - tensors_ratio = np.divide(test_statistic[1], test_statistic[2]) - tensors_ratio_delta_unity = np.absolute(np.subtract(tensors_ratio, 1)) - - try: - assert (numpy_ratio_delta_unity < tolerance['numpy']).all() - except AssertionError: - print( - 'Ratio to NumPy+SciPy exceeded tolerance of {}: {}'.format( - tolerance['numpy'], numpy_ratio_delta_unity.tolist() + for p in param_tests: + fixed_vals = [p] + + test_statistic = [] + for backend in backends: + pyhf.set_backend(backend) + + qmu_tilde = pyhf.infer.test_statistics.qmu_tilde( + 1.0, + data, + pdf, + pdf.config.suggested_init(), + pdf.config.suggested_bounds(), + fixed_vals, + ) + test_statistic.append(qmu_tilde) + + # compare to NumPy/SciPy + test_statistic = np.array(test_statistic) + numpy_ratio = np.divide(test_statistic, test_statistic[0]) + numpy_ratio_delta_unity = np.absolute(np.subtract(numpy_ratio, 1)) + + # compare tensor libraries to each other + tensors_ratio = np.divide(test_statistic[1], test_statistic[2]) + tensors_ratio_delta_unity = np.absolute(np.subtract(tensors_ratio, 1)) + + try: + assert (numpy_ratio_delta_unity < tolerance['numpy']).all() + except AssertionError: + print( + 'Ratio to NumPy+SciPy exceeded tolerance of {}: {}'.format( + tolerance['numpy'], numpy_ratio_delta_unity.tolist() + ) ) - ) - assert False - try: - assert (tensors_ratio_delta_unity < tolerance['tensors']).all() - except AssertionError: - print( - 'Ratio between tensor backends exceeded tolerance of {}: {}'.format( - tolerance['tensors'], tensors_ratio_delta_unity.tolist() + assert False + try: + assert (tensors_ratio_delta_unity < tolerance['tensors']).all() + except AssertionError: + print( + 'Ratio between tensor backends exceeded tolerance of {}: {}'.format( + tolerance['tensors'], tensors_ratio_delta_unity.tolist() + ) ) - ) - assert False + assert False From 40b781b4f0f5c830848c253e15927c811301fca4 Mon Sep 17 00:00:00 2001 From: Carlo Gottardo Date: Mon, 27 Apr 2020 12:41:11 +0200 Subject: [PATCH 08/24] correct histosys data structure in test_backend_consistency --- tests/test_backend_consistency.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_backend_consistency.py b/tests/test_backend_consistency.py index 5c46f5ab24..bf4edecc85 100644 --- a/tests/test_backend_consistency.py +++ b/tests/test_backend_consistency.py @@ -96,11 +96,11 @@ def test_hypotest_qmu_tilde( { 'name': 'norm_bkgunc', 'type': 'histosys', - 'hi_data': bkg_unc_10pc_up, - 'lo_data': bkg_unc_10pc_dn, + 'data': {'hi_data': bkg_unc_10pc_up, 'lo_data': bkg_unc_10pc_dn,}, }, ], } + print(background_sample) samples = ( [background_sample, signal_sample] if invert_order From bcd477d7e3f16554766502db1f362feea208c3e8 Mon Sep 17 00:00:00 2001 From: Carlo Gottardo Date: Mon, 27 Apr 2020 13:08:44 +0200 Subject: [PATCH 09/24] fix test_backend_consistency --- tests/test_backend_consistency.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/test_backend_consistency.py b/tests/test_backend_consistency.py index bf4edecc85..32d8faae93 100644 --- a/tests/test_backend_consistency.py +++ b/tests/test_backend_consistency.py @@ -116,10 +116,10 @@ def test_hypotest_qmu_tilde( norm_bkgunc_idx = idx # Floating norm_bkgunc, fixed at nominal, plus/minus 1 sigma param_tests = [ - None, - (norm_bkgunc_idx, 0.0), - (norm_bkgunc_idx, -1.0), - (norm_bkgunc_idx, 1.0), + [], + [(norm_bkgunc_idx, 0.0)], + [(norm_bkgunc_idx, 1.0)], + [(norm_bkgunc_idx, -1.0)], ] data = source['bindata']['data'] + pdf.config.auxdata @@ -132,7 +132,7 @@ def test_hypotest_qmu_tilde( ] for p in param_tests: - fixed_vals = [p] + fixed_vals = p test_statistic = [] for backend in backends: From cc774e7ca1744c04fe510ef4995941bdb725ed9c Mon Sep 17 00:00:00 2001 From: Carlo Gottardo Date: Tue, 5 May 2020 15:26:25 +0200 Subject: [PATCH 10/24] NP index fix in test_backend_consistency.py --- tests/test_backend_consistency.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/test_backend_consistency.py b/tests/test_backend_consistency.py index 32d8faae93..736ac5ca00 100644 --- a/tests/test_backend_consistency.py +++ b/tests/test_backend_consistency.py @@ -96,7 +96,7 @@ def test_hypotest_qmu_tilde( { 'name': 'norm_bkgunc', 'type': 'histosys', - 'data': {'hi_data': bkg_unc_10pc_up, 'lo_data': bkg_unc_10pc_dn,}, + 'data': {'hi_data': bkg_unc_10pc_up, 'lo_data': bkg_unc_10pc_dn}, }, ], } @@ -114,6 +114,7 @@ def test_hypotest_qmu_tilde( for idx, par in enumerate(pdf.config.par_map): if par == "norm_bkgunc": norm_bkgunc_idx = idx + # Floating norm_bkgunc, fixed at nominal, plus/minus 1 sigma param_tests = [ [], From d9e197208037da567207ab3f1a44da82d4134ade Mon Sep 17 00:00:00 2001 From: Carlo Gottardo Date: Tue, 5 May 2020 15:42:05 +0200 Subject: [PATCH 11/24] Need 5 percent tolerane for 500 bin fit with TF --- tests/test_backend_consistency.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/tests/test_backend_consistency.py b/tests/test_backend_consistency.py index 736ac5ca00..7824bacf38 100644 --- a/tests/test_backend_consistency.py +++ b/tests/test_backend_consistency.py @@ -58,8 +58,11 @@ def generate_source_poisson(n_bins): @pytest.mark.parametrize('n_bins', bins, ids=bin_ids) @pytest.mark.parametrize('invert_order', [False, True], ids=['normal', 'inverted']) +@pytest.mark.parametrize( + 'fix_param', [None, 0.0, 1.0, -1.0], ids=['none', 'central', 'up', 'down'] +) def test_hypotest_qmu_tilde( - n_bins, invert_order, tolerance={'numpy': 1e-02, 'tensors': 5e-03} + n_bins, invert_order, fix_param, tolerance={'numpy': 1e-02, 'tensors': 5e-03} ): """ Check that the different backends all compute a test statistic From 9d45958650615785e326e072e530d37c5d3b6b2c Mon Sep 17 00:00:00 2001 From: Giordon Stark Date: Thu, 3 Sep 2020 18:54:07 -0400 Subject: [PATCH 12/24] pass in fixed_vals by converting to fixed_params --- src/pyhf/infer/mle.py | 38 +++++++++++++++++++++++++++++++++----- 1 file changed, 33 insertions(+), 5 deletions(-) diff --git a/src/pyhf/infer/mle.py b/src/pyhf/infer/mle.py index 82448d5958..f30dbe1e9b 100644 --- a/src/pyhf/infer/mle.py +++ b/src/pyhf/infer/mle.py @@ -96,8 +96,23 @@ def fit(data, pdf, init_pars=None, par_bounds=None, fixed_vals=None, **kwargs): _, opt = get_backend() init_pars = init_pars or pdf.config.suggested_init() par_bounds = par_bounds or pdf.config.suggested_bounds() + + # get fixed params from the model + model_fixed_params = [ + (index, init) + for index, (init, is_fixed) in enumerate( + zip(init_pars, pdf.config.suggested_fixed()) + ) + if is_fixed + ] + # add user-defined ones at the end + fixed_params = model_fixed_params + (fixed_vals or []) + + # de-dupe and use last-appended result for each index + fixed_params = list(dict(fixed_params)) + return opt.minimize( - twice_nll, data, pdf, init_pars, par_bounds, fixed_vals, **kwargs + twice_nll, data, pdf, init_pars, par_bounds, fixed_params, **kwargs ) @@ -157,13 +172,26 @@ def fixed_poi_fit( 'No POI is defined. A POI is required to fit with a fixed POI.' ) - fixed_params = [(pdf.config.poi_index, poi_val)] - if fixed_vals: - fixed_params += fixed_vals - _, opt = get_backend() init_pars = init_pars or pdf.config.suggested_init() par_bounds = par_bounds or pdf.config.suggested_bounds() + + # get fixed params from the model + model_fixed_params = [ + (index, init) + for index, (init, is_fixed) in enumerate( + zip(init_pars, pdf.config.suggested_fixed()) + ) + if is_fixed + ] + # add user-defined ones at the end + fixed_params = model_fixed_params + (fixed_vals or []) + # add the fixed POI + fixed_params += [(pdf.config.poi_index, poi_val)] + + # de-dupe and use last-appended result for each index + fixed_params = list(dict(fixed_params)) + return opt.minimize( twice_nll, data, From d3c3b47e046c354702e70e8a05b4d733fd9eb497 Mon Sep 17 00:00:00 2001 From: Giordon Stark Date: Thu, 3 Sep 2020 19:45:55 -0400 Subject: [PATCH 13/24] make sure we propagate fixed_vals correctly everywhere --- src/pyhf/infer/__init__.py | 16 +++++++++++++++- src/pyhf/infer/calculators.py | 15 +++++++++++++++ src/pyhf/infer/mle.py | 10 +++++----- src/pyhf/infer/test_statistics.py | 19 +++++++++++-------- 4 files changed, 46 insertions(+), 14 deletions(-) diff --git a/src/pyhf/infer/__init__.py b/src/pyhf/infer/__init__.py index 89c2d24773..3b5c739eeb 100644 --- a/src/pyhf/infer/__init__.py +++ b/src/pyhf/infer/__init__.py @@ -126,7 +126,19 @@ def hypotest( """ init_pars = init_pars or pdf.config.suggested_init() par_bounds = par_bounds or pdf.config.suggested_bounds() - tensorlib, _ = get_backend() + # get fixed vals from the model + model_fixed_vals = [ + (index, init) + for index, (init, is_fixed) in enumerate( + zip(init_pars, pdf.config.suggested_fixed()) + ) + if is_fixed + ] + # add user-defined ones at the end + fixed_vals = model_fixed_vals + (fixed_vals or []) + + # de-dupe and use last-appended result for each index + fixed_vals = list(dict(fixed_vals)) calc = AsymptoticCalculator( data, pdf, init_pars, par_bounds, fixed_vals, qtilde=qtilde @@ -137,6 +149,8 @@ def hypotest( CLsb = sig_plus_bkg_distribution.pvalue(teststat) CLb = b_only_distribution.pvalue(teststat) CLs = CLsb / CLb + + tensorlib, _ = get_backend() # Ensure that all CL values are 0-d tensors CLsb, CLb, CLs = ( tensorlib.astensor(CLsb), diff --git a/src/pyhf/infer/calculators.py b/src/pyhf/infer/calculators.py index fb07d089e6..aafd33c3f9 100644 --- a/src/pyhf/infer/calculators.py +++ b/src/pyhf/infer/calculators.py @@ -143,6 +143,21 @@ def __init__( self.pdf = pdf self.init_pars = init_pars or pdf.config.suggested_init() self.par_bounds = par_bounds or pdf.config.suggested_bounds() + + # get fixed vals from the model + model_fixed_vals = [ + (index, init) + for index, (init, is_fixed) in enumerate( + zip(init_pars, pdf.config.suggested_fixed()) + ) + if is_fixed + ] + # add user-defined ones at the end + fixed_vals = model_fixed_vals + (fixed_vals or []) + + # de-dupe and use last-appended result for each index + fixed_vals = list(dict(fixed_vals)) + self.fixed_vals = fixed_vals self.qtilde = qtilde self.sqrtqmuA_v = None diff --git a/src/pyhf/infer/mle.py b/src/pyhf/infer/mle.py index f30dbe1e9b..a6166398ce 100644 --- a/src/pyhf/infer/mle.py +++ b/src/pyhf/infer/mle.py @@ -97,8 +97,8 @@ def fit(data, pdf, init_pars=None, par_bounds=None, fixed_vals=None, **kwargs): init_pars = init_pars or pdf.config.suggested_init() par_bounds = par_bounds or pdf.config.suggested_bounds() - # get fixed params from the model - model_fixed_params = [ + # get fixed vals from the model + model_fixed_vals = [ (index, init) for index, (init, is_fixed) in enumerate( zip(init_pars, pdf.config.suggested_fixed()) @@ -106,13 +106,13 @@ def fit(data, pdf, init_pars=None, par_bounds=None, fixed_vals=None, **kwargs): if is_fixed ] # add user-defined ones at the end - fixed_params = model_fixed_params + (fixed_vals or []) + fixed_vals = model_fixed_vals + (fixed_vals or []) # de-dupe and use last-appended result for each index - fixed_params = list(dict(fixed_params)) + fixed_vals = list(dict(fixed_vals)) return opt.minimize( - twice_nll, data, pdf, init_pars, par_bounds, fixed_params, **kwargs + twice_nll, data, pdf, init_pars, par_bounds, fixed_vals, **kwargs ) diff --git a/src/pyhf/infer/test_statistics.py b/src/pyhf/infer/test_statistics.py index d1cd663ecd..7bd4c4fdb3 100644 --- a/src/pyhf/infer/test_statistics.py +++ b/src/pyhf/infer/test_statistics.py @@ -17,7 +17,7 @@ def _qmu_like(mu, data, pdf, init_pars, par_bounds, fixed_vals): """ tensorlib, optimizer = get_backend() tmu_like_stat, (_, muhatbhat) = _tmu_like( - mu, data, pdf, init_pars, par_bounds, return_fitted_pars=True + mu, data, pdf, init_pars, par_bounds, fixed_vals, return_fitted_pars=True ) qmu_like_stat = tensorlib.where( muhatbhat[pdf.config.poi_index] > mu, tensorlib.astensor(0.0), tmu_like_stat @@ -84,7 +84,8 @@ def qmu(mu, data, pdf, init_pars, par_bounds, fixed_vals): >>> init_pars = model.config.suggested_init() >>> par_bounds = model.config.suggested_bounds() >>> par_bounds[model.config.poi_index] = [-10.0, 10.0] - >>> pyhf.infer.test_statistics.qmu(test_mu, data, model, init_pars, par_bounds) + >>> fixed_vals = [] + >>> pyhf.infer.test_statistics.qmu(test_mu, data, model, init_pars, par_bounds, []) array(3.9549891) Args: @@ -93,7 +94,7 @@ def qmu(mu, data, pdf, init_pars, par_bounds, fixed_vals): pdf (~pyhf.pdf.Model): The HistFactory statistical model used in the likelihood ratio calculation init_pars (`list`): Values to initialize the model parameters at for the fit par_bounds (`list` of `list`\s or `tuple`\s): The extrema of values the model parameters are allowed to reach in the fit - fixed_vals(`list`): Parameters held constant in the fit + fixed_vals (`list`): Parameters held constant in the fit Returns: Float: The calculated test statistic, :math:`q_{\mu}` @@ -158,7 +159,7 @@ def qmu_tilde(mu, data, pdf, init_pars, par_bounds, fixed_vals): pdf (~pyhf.pdf.Model): The statistical model adhering to the schema model.json init_pars (`list`): Values to initialize the model parameters at for the fit par_bounds (`list` of `list`\s or `tuple`\s): The extrema of values the model parameters are allowed to reach in the fit - fixed_vals(`list`): Parameters held constant in the fit + fixed_vals (`list`): Parameters held constant in the fit Returns: Float: The calculated test statistic, :math:`\tilde{q}_{\mu}` @@ -203,7 +204,8 @@ def tmu(mu, data, pdf, init_pars, par_bounds, fixed_vals): >>> init_pars = model.config.suggested_init() >>> par_bounds = model.config.suggested_bounds() >>> par_bounds[model.config.poi_index] = [-10.0, 10.0] - >>> pyhf.infer.test_statistics.tmu(test_mu, data, model, init_pars, par_bounds) + >>> fixed_vals = [] + >>> pyhf.infer.test_statistics.tmu(test_mu, data, model, init_pars, par_bounds, []) array(3.9549891) Args: @@ -212,7 +214,7 @@ def tmu(mu, data, pdf, init_pars, par_bounds, fixed_vals): pdf (~pyhf.pdf.Model): The statistical model adhering to the schema model.json init_pars (`list`): Values to initialize the model parameters at for the fit par_bounds (`list` of `list`\s or `tuple`\s): The extrema of values the model parameters are allowed to reach in the fit - fixed_vals(`list`): Parameters held constant in the fit + fixed_vals (`list`): Parameters held constant in the fit Returns: Float: The calculated test statistic, :math:`t_{\mu}` @@ -262,7 +264,8 @@ def tmu_tilde(mu, data, pdf, init_pars, par_bounds, fixed_vals): >>> test_mu = 1.0 >>> init_pars = model.config.suggested_init() >>> par_bounds = model.config.suggested_bounds() - >>> pyhf.infer.test_statistics.tmu_tilde(test_mu, data, model, init_pars, par_bounds) + >>> fixed_vals = [] + >>> pyhf.infer.test_statistics.tmu_tilde(test_mu, data, model, init_pars, par_bounds, []) array(3.93824492) Args: @@ -271,7 +274,7 @@ def tmu_tilde(mu, data, pdf, init_pars, par_bounds, fixed_vals): pdf (~pyhf.pdf.Model): The statistical model adhering to the schema model.json init_pars (`list`): Values to initialize the model parameters at for the fit par_bounds (`list` of `list`\s or `tuple`\s): The extrema of values the model parameters are allowed to reach in the fit - fixed_vals(`list`): Parameters held constant in the fit + fixed_vals (`list`): Parameters held constant in the fit Returns: Float: The calculated test statistic, :math:`\tilde{t}_{\mu}` From e0c82d1e5f6433abdf9e80ef7e13d9296be78482 Mon Sep 17 00:00:00 2001 From: Giordon Stark Date: Thu, 3 Sep 2020 19:48:10 -0400 Subject: [PATCH 14/24] fix up doc a bit --- src/pyhf/infer/test_statistics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pyhf/infer/test_statistics.py b/src/pyhf/infer/test_statistics.py index 7bd4c4fdb3..b4d03182a2 100644 --- a/src/pyhf/infer/test_statistics.py +++ b/src/pyhf/infer/test_statistics.py @@ -150,7 +150,7 @@ def qmu_tilde(mu, data, pdf, init_pars, par_bounds, fixed_vals): >>> test_mu = 1.0 >>> init_pars = model.config.suggested_init() >>> par_bounds = model.config.suggested_bounds() - >>> pyhf.infer.test_statistics.qmu_tilde(test_mu, data, model, init_pars, par_bounds) + >>> pyhf.infer.test_statistics.qmu_tilde(test_mu, data, model, init_pars, par_bounds, []) array(3.93824492) Args: From 0c232a6ba51d4d854188c9779669dfe6c262f44c Mon Sep 17 00:00:00 2001 From: Giordon Stark Date: Thu, 3 Sep 2020 19:54:25 -0400 Subject: [PATCH 15/24] fix up notational difference --- src/pyhf/infer/mle.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/pyhf/infer/mle.py b/src/pyhf/infer/mle.py index a6166398ce..bf581eeb51 100644 --- a/src/pyhf/infer/mle.py +++ b/src/pyhf/infer/mle.py @@ -176,8 +176,8 @@ def fixed_poi_fit( init_pars = init_pars or pdf.config.suggested_init() par_bounds = par_bounds or pdf.config.suggested_bounds() - # get fixed params from the model - model_fixed_params = [ + # get fixed vals from the model + model_fixed_vals = [ (index, init) for index, (init, is_fixed) in enumerate( zip(init_pars, pdf.config.suggested_fixed()) @@ -185,12 +185,13 @@ def fixed_poi_fit( if is_fixed ] # add user-defined ones at the end - fixed_params = model_fixed_params + (fixed_vals or []) + fixed_vals = model_fixed_vals + (fixed_vals or []) + # add the fixed POI - fixed_params += [(pdf.config.poi_index, poi_val)] + fixed_vals = fixed_vals + [(pdf.config.poi_index, poi_val)] # de-dupe and use last-appended result for each index - fixed_params = list(dict(fixed_params)) + fixed_vals = list(dict(fixed_vals)) return opt.minimize( twice_nll, @@ -198,6 +199,6 @@ def fixed_poi_fit( pdf, init_pars, par_bounds, - fixed_params, + fixed_vals, **kwargs, ) From aada2e97fabbdb1128f6f30e5182532f3bc9f879 Mon Sep 17 00:00:00 2001 From: Giordon Stark Date: Thu, 3 Sep 2020 20:00:25 -0400 Subject: [PATCH 16/24] undo changes to tests for now --- tests/test_backend_consistency.py | 105 ++++++++++-------------------- 1 file changed, 36 insertions(+), 69 deletions(-) diff --git a/tests/test_backend_consistency.py b/tests/test_backend_consistency.py index 7824bacf38..61f9b8579c 100644 --- a/tests/test_backend_consistency.py +++ b/tests/test_backend_consistency.py @@ -58,11 +58,8 @@ def generate_source_poisson(n_bins): @pytest.mark.parametrize('n_bins', bins, ids=bin_ids) @pytest.mark.parametrize('invert_order', [False, True], ids=['normal', 'inverted']) -@pytest.mark.parametrize( - 'fix_param', [None, 0.0, 1.0, -1.0], ids=['none', 'central', 'up', 'down'] -) def test_hypotest_qmu_tilde( - n_bins, invert_order, fix_param, tolerance={'numpy': 1e-02, 'tensors': 5e-03} + n_bins, invert_order, tolerance={'numpy': 1e-02, 'tensors': 5e-03} ): """ Check that the different backends all compute a test statistic @@ -78,8 +75,6 @@ def test_hypotest_qmu_tilde( """ source = generate_source_static(n_bins) - bkg_unc_10pc_up = [x + 0.10 * x for x in source['bindata']['bkg']] - bkg_unc_10pc_dn = [x - 0.10 * x for x in source['bindata']['bkg']] signal_sample = { 'name': 'signal', @@ -95,15 +90,9 @@ def test_hypotest_qmu_tilde( 'name': 'uncorr_bkguncrt', 'type': 'shapesys', 'data': source['bindata']['bkgerr'], - }, - { - 'name': 'norm_bkgunc', - 'type': 'histosys', - 'data': {'hi_data': bkg_unc_10pc_up, 'lo_data': bkg_unc_10pc_dn}, - }, + } ], } - print(background_sample) samples = ( [background_sample, signal_sample] if invert_order @@ -112,20 +101,6 @@ def test_hypotest_qmu_tilde( spec = {'channels': [{'name': 'singlechannel', 'samples': samples}]} pdf = pyhf.Model(spec) - # get norm_bkgunc index to set it constatnt - norm_bkgunc_idx = -1 - for idx, par in enumerate(pdf.config.par_map): - if par == "norm_bkgunc": - norm_bkgunc_idx = idx - - # Floating norm_bkgunc, fixed at nominal, plus/minus 1 sigma - param_tests = [ - [], - [(norm_bkgunc_idx, 0.0)], - [(norm_bkgunc_idx, 1.0)], - [(norm_bkgunc_idx, -1.0)], - ] - data = source['bindata']['data'] + pdf.config.auxdata backends = [ @@ -135,47 +110,39 @@ def test_hypotest_qmu_tilde( pyhf.tensor.jax_backend(precision='64b'), ] - for p in param_tests: - fixed_vals = p - - test_statistic = [] - for backend in backends: - pyhf.set_backend(backend) - - qmu_tilde = pyhf.infer.test_statistics.qmu_tilde( - 1.0, - data, - pdf, - pdf.config.suggested_init(), - pdf.config.suggested_bounds(), - fixed_vals, - ) - test_statistic.append(qmu_tilde) - - # compare to NumPy/SciPy - test_statistic = np.array(test_statistic) - numpy_ratio = np.divide(test_statistic, test_statistic[0]) - numpy_ratio_delta_unity = np.absolute(np.subtract(numpy_ratio, 1)) - - # compare tensor libraries to each other - tensors_ratio = np.divide(test_statistic[1], test_statistic[2]) - tensors_ratio_delta_unity = np.absolute(np.subtract(tensors_ratio, 1)) - - try: - assert (numpy_ratio_delta_unity < tolerance['numpy']).all() - except AssertionError: - print( - 'Ratio to NumPy+SciPy exceeded tolerance of {}: {}'.format( - tolerance['numpy'], numpy_ratio_delta_unity.tolist() - ) + test_statistic = [] + for backend in backends: + pyhf.set_backend(backend) + + qmu_tilde = pyhf.infer.test_statistics.qmu_tilde( + 1.0, data, pdf, pdf.config.suggested_init(), pdf.config.suggested_bounds() + ) + test_statistic.append(qmu_tilde) + + # compare to NumPy/SciPy + test_statistic = np.array(test_statistic) + numpy_ratio = np.divide(test_statistic, test_statistic[0]) + numpy_ratio_delta_unity = np.absolute(np.subtract(numpy_ratio, 1)) + + # compare tensor libraries to each other + tensors_ratio = np.divide(test_statistic[1], test_statistic[2]) + tensors_ratio_delta_unity = np.absolute(np.subtract(tensors_ratio, 1)) + + try: + assert (numpy_ratio_delta_unity < tolerance['numpy']).all() + except AssertionError: + print( + 'Ratio to NumPy+SciPy exceeded tolerance of {}: {}'.format( + tolerance['numpy'], numpy_ratio_delta_unity.tolist() ) - assert False - try: - assert (tensors_ratio_delta_unity < tolerance['tensors']).all() - except AssertionError: - print( - 'Ratio between tensor backends exceeded tolerance of {}: {}'.format( - tolerance['tensors'], tensors_ratio_delta_unity.tolist() - ) + ) + assert False + try: + assert (tensors_ratio_delta_unity < tolerance['tensors']).all() + except AssertionError: + print( + 'Ratio between tensor backends exceeded tolerance of {}: {}'.format( + tolerance['tensors'], tensors_ratio_delta_unity.tolist() ) - assert False + ) + assert False From 61b97ec0741a0e0c534c15bca57ec7b574fea1d0 Mon Sep 17 00:00:00 2001 From: Giordon Stark Date: Thu, 3 Sep 2020 20:16:09 -0400 Subject: [PATCH 17/24] fixed it up --- src/pyhf/infer/__init__.py | 2 +- src/pyhf/infer/calculators.py | 2 +- src/pyhf/infer/mle.py | 4 ++-- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/pyhf/infer/__init__.py b/src/pyhf/infer/__init__.py index 3b5c739eeb..f9f6e2570a 100644 --- a/src/pyhf/infer/__init__.py +++ b/src/pyhf/infer/__init__.py @@ -138,7 +138,7 @@ def hypotest( fixed_vals = model_fixed_vals + (fixed_vals or []) # de-dupe and use last-appended result for each index - fixed_vals = list(dict(fixed_vals)) + fixed_vals = list(dict(fixed_vals).items()) calc = AsymptoticCalculator( data, pdf, init_pars, par_bounds, fixed_vals, qtilde=qtilde diff --git a/src/pyhf/infer/calculators.py b/src/pyhf/infer/calculators.py index aafd33c3f9..5f7e53be58 100644 --- a/src/pyhf/infer/calculators.py +++ b/src/pyhf/infer/calculators.py @@ -156,7 +156,7 @@ def __init__( fixed_vals = model_fixed_vals + (fixed_vals or []) # de-dupe and use last-appended result for each index - fixed_vals = list(dict(fixed_vals)) + fixed_vals = list(dict(fixed_vals).items()) self.fixed_vals = fixed_vals self.qtilde = qtilde diff --git a/src/pyhf/infer/mle.py b/src/pyhf/infer/mle.py index bf581eeb51..b5cdb8a803 100644 --- a/src/pyhf/infer/mle.py +++ b/src/pyhf/infer/mle.py @@ -109,7 +109,7 @@ def fit(data, pdf, init_pars=None, par_bounds=None, fixed_vals=None, **kwargs): fixed_vals = model_fixed_vals + (fixed_vals or []) # de-dupe and use last-appended result for each index - fixed_vals = list(dict(fixed_vals)) + fixed_vals = list(dict(fixed_vals).items()) return opt.minimize( twice_nll, data, pdf, init_pars, par_bounds, fixed_vals, **kwargs @@ -191,7 +191,7 @@ def fixed_poi_fit( fixed_vals = fixed_vals + [(pdf.config.poi_index, poi_val)] # de-dupe and use last-appended result for each index - fixed_vals = list(dict(fixed_vals)) + fixed_vals = list(dict(fixed_vals).items()) return opt.minimize( twice_nll, From 702b2171b73f04bb743c2a1e4a50a265cc6edf9d Mon Sep 17 00:00:00 2001 From: Giordon Stark Date: Thu, 3 Sep 2020 20:28:22 -0400 Subject: [PATCH 18/24] get tests passing --- tests/test_backend_consistency.py | 7 ++++++- tests/test_infer.py | 3 +++ tests/test_teststats.py | 16 ++++++++-------- 3 files changed, 17 insertions(+), 9 deletions(-) diff --git a/tests/test_backend_consistency.py b/tests/test_backend_consistency.py index 61f9b8579c..b42f7afaed 100644 --- a/tests/test_backend_consistency.py +++ b/tests/test_backend_consistency.py @@ -115,7 +115,12 @@ def test_hypotest_qmu_tilde( pyhf.set_backend(backend) qmu_tilde = pyhf.infer.test_statistics.qmu_tilde( - 1.0, data, pdf, pdf.config.suggested_init(), pdf.config.suggested_bounds() + 1.0, + data, + pdf, + pdf.config.suggested_init(), + pdf.config.suggested_bounds(), + [], ) test_statistic.append(qmu_tilde) diff --git a/tests/test_infer.py b/tests/test_infer.py index 527ec60f1c..2a2b4c9a2b 100644 --- a/tests/test_infer.py +++ b/tests/test_infer.py @@ -138,6 +138,9 @@ def suggested_init(self): def suggested_bounds(self): return [[0.0, 10.0], [0.0, 10.0]] + def suggested_fixed(self): + return [False, False] + class NonPyhfModel(object): def __init__(self, spec): self.sig, self.nominal, self.uncert = spec diff --git a/tests/test_teststats.py b/tests/test_teststats.py index bd942152e4..aca8f9f361 100644 --- a/tests/test_teststats.py +++ b/tests/test_teststats.py @@ -12,7 +12,7 @@ def test_qmu(caplog): par_bounds = model.config.suggested_bounds() with caplog.at_level(logging.WARNING, "pyhf.infer.test_statistics"): - pyhf.infer.test_statistics.qmu(mu, data, model, init_pars, par_bounds) + pyhf.infer.test_statistics.qmu(mu, data, model, init_pars, par_bounds, []) assert "WARNING qmu test statistic used for fit" in caplog.text caplog.clear() @@ -26,7 +26,7 @@ def test_qmu_tilde(caplog): par_bounds[model.config.poi_index] = [-10, 10] with caplog.at_level(logging.WARNING, "pyhf.infer.test_statistics"): - pyhf.infer.test_statistics.qmu_tilde(mu, data, model, init_pars, par_bounds) + pyhf.infer.test_statistics.qmu_tilde(mu, data, model, init_pars, par_bounds, []) assert "WARNING qmu_tilde test statistic used for fit" in caplog.text caplog.clear() @@ -38,7 +38,7 @@ def test_tmu(caplog): init_pars = model.config.suggested_init() par_bounds = model.config.suggested_bounds() with caplog.at_level(logging.WARNING, "pyhf.infer.test_statistics"): - pyhf.infer.test_statistics.tmu(mu, data, model, init_pars, par_bounds) + pyhf.infer.test_statistics.tmu(mu, data, model, init_pars, par_bounds, []) assert "WARNING tmu test statistic used for fit" in caplog.text caplog.clear() @@ -52,7 +52,7 @@ def test_tmu_tilde(caplog): par_bounds[model.config.poi_index] = [-10, 10] with caplog.at_level(logging.WARNING, "pyhf.infer.test_statistics"): - pyhf.infer.test_statistics.tmu_tilde(mu, data, model, init_pars, par_bounds) + pyhf.infer.test_statistics.tmu_tilde(mu, data, model, init_pars, par_bounds, []) assert "WARNING tmu_tilde test statistic used for fit" in caplog.text caplog.clear() @@ -86,7 +86,7 @@ def test_no_poi_test_stats(): par_bounds = model.config.suggested_bounds() with pytest.raises(pyhf.exceptions.UnspecifiedPOI) as excinfo: - pyhf.infer.test_statistics.qmu(test_poi, data, model, init_pars, par_bounds) + pyhf.infer.test_statistics.qmu(test_poi, data, model, init_pars, par_bounds, []) assert ( "No POI is defined. A POI is required for profile likelihood based test statistics." in str(excinfo.value) @@ -94,7 +94,7 @@ def test_no_poi_test_stats(): with pytest.raises(pyhf.exceptions.UnspecifiedPOI) as excinfo: pyhf.infer.test_statistics.qmu_tilde( - test_poi, data, model, init_pars, par_bounds + test_poi, data, model, init_pars, par_bounds, [] ) assert ( "No POI is defined. A POI is required for profile likelihood based test statistics." @@ -102,7 +102,7 @@ def test_no_poi_test_stats(): ) with pytest.raises(pyhf.exceptions.UnspecifiedPOI) as excinfo: - pyhf.infer.test_statistics.tmu(test_poi, data, model, init_pars, par_bounds) + pyhf.infer.test_statistics.tmu(test_poi, data, model, init_pars, par_bounds, []) assert ( "No POI is defined. A POI is required for profile likelihood based test statistics." in str(excinfo.value) @@ -110,7 +110,7 @@ def test_no_poi_test_stats(): with pytest.raises(pyhf.exceptions.UnspecifiedPOI) as excinfo: pyhf.infer.test_statistics.tmu_tilde( - test_poi, data, model, init_pars, par_bounds + test_poi, data, model, init_pars, par_bounds, [] ) assert ( "No POI is defined. A POI is required for profile likelihood based test statistics." From 680502b482ffe070dfe17b6910c74d753f088287 Mon Sep 17 00:00:00 2001 From: Giordon Stark Date: Thu, 3 Sep 2020 20:57:05 -0400 Subject: [PATCH 19/24] update a test and provide a fixed_params to deal with the pdf issue --- src/pyhf/infer/calculators.py | 14 +++++++++++--- tests/test_infer.py | 2 +- 2 files changed, 12 insertions(+), 4 deletions(-) diff --git a/src/pyhf/infer/calculators.py b/src/pyhf/infer/calculators.py index 5f7e53be58..7f9b0a8261 100644 --- a/src/pyhf/infer/calculators.py +++ b/src/pyhf/infer/calculators.py @@ -122,7 +122,14 @@ class AsymptoticCalculator(object): """The Asymptotic Calculator.""" def __init__( - self, data, pdf, init_pars=None, par_bounds=None, fixed_vals=None, qtilde=False + self, + data, + pdf, + init_pars=None, + par_bounds=None, + fixed_params=None, + fixed_vals=None, + qtilde=False, ): """ Asymptotic Calculator. @@ -132,7 +139,8 @@ def __init__( pdf (~pyhf.pdf.Model): The statistical model adhering to the schema ``model.json``. init_pars (`tensor`): The initial parameter values to be used for fitting. par_bounds (`tensor`): The parameter value bounds to be used for fitting. - fixed_vals (`tensor`): Parameters to be held constant in the fit. + fixed_params (`tensor`): A list of booleans for each parameter on whether it should be fixed in the fit or not. + fixed_vals (list of tuples): The fixed parameters and fixed values to be used for fitting. qtilde (`bool`): Whether to use qtilde as the test statistic. Returns: @@ -148,7 +156,7 @@ def __init__( model_fixed_vals = [ (index, init) for index, (init, is_fixed) in enumerate( - zip(init_pars, pdf.config.suggested_fixed()) + zip(init_pars, fixed_params or pdf.config.suggested_fixed()) ) if is_fixed ] diff --git a/tests/test_infer.py b/tests/test_infer.py index 2a2b4c9a2b..291e403702 100644 --- a/tests/test_infer.py +++ b/tests/test_infer.py @@ -187,7 +187,7 @@ def logpdf(self, pars, data): @pytest.mark.parametrize("qtilde", [True, False]) def test_calculator_distributions_without_teststatistic(qtilde): calc = pyhf.infer.AsymptoticCalculator( - [0.0], {}, [1.0], [(0.0, 10.0)], qtilde=qtilde + [0.0], {}, [1.0], [(0.0, 10.0)], [False], qtilde=qtilde ) with pytest.raises(RuntimeError): calc.distributions(1.0) From d6831b6e4e59b3cd3e3d1cf415964d2a5dfe41e1 Mon Sep 17 00:00:00 2001 From: Giordon Stark Date: Thu, 3 Sep 2020 21:03:13 -0400 Subject: [PATCH 20/24] scipy is able to minimize better so bump to 1 maxiter to force it to fail --- tests/test_scripts.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_scripts.py b/tests/test_scripts.py index ce74887d1b..3a2ff3fca4 100644 --- a/tests/test_scripts.py +++ b/tests/test_scripts.py @@ -238,7 +238,7 @@ def test_testpoi(tmpdir, script_runner): 'optimizer', ['scipy', 'minuit', 'scipy_optimizer', 'minuit_optimizer'] ) @pytest.mark.parametrize( - 'opts,success', [(['maxiter=1000'], True), (['maxiter=10'], False)] + 'opts,success', [(['maxiter=1000'], True), (['maxiter=1'], False)] ) def test_cls_optimizer(tmpdir, script_runner, optimizer, opts, success): temp = tmpdir.join("parsed_output.json") From 91847be9229d797e4e75fa4149ac9adbb4851505 Mon Sep 17 00:00:00 2001 From: Giordon Stark Date: Fri, 4 Sep 2020 10:55:50 -0400 Subject: [PATCH 21/24] switch to fixed_params for public-facing calls in pyhf.infer --- src/pyhf/infer/__init__.py | 24 ++++----------- src/pyhf/infer/calculators.py | 32 +++++--------------- src/pyhf/infer/mle.py | 27 +++++------------ src/pyhf/infer/test_statistics.py | 49 ++++++++++++++++--------------- 4 files changed, 47 insertions(+), 85 deletions(-) diff --git a/src/pyhf/infer/__init__.py b/src/pyhf/infer/__init__.py index f9f6e2570a..ce64eb819b 100644 --- a/src/pyhf/infer/__init__.py +++ b/src/pyhf/infer/__init__.py @@ -11,7 +11,7 @@ def hypotest( pdf, init_pars=None, par_bounds=None, - fixed_vals=None, + fixed_params=None, qtilde=False, **kwargs, ): @@ -39,9 +39,9 @@ def hypotest( poi_test (Number or Tensor): The value of the parameter of interest (POI) data (Number or Tensor): The data considered pdf (~pyhf.pdf.Model): The statistical model adhering to the schema ``model.json`` - init_pars (Array or Tensor): The initial parameter values to be used for minimization - par_bounds (Array or Tensor): The parameter value bounds to be used for minimization - fixed_vals (list of tuples): Parameters to be held constant and their value + init_pars (`tensor`): The initial parameter values to be used for minimization + par_bounds (`tensor`): The parameter value bounds to be used for minimization + fixed_params (`tensor`): Whether to fix the parameter to the init_pars value during minimization qtilde (Bool): When ``True`` perform the calculation using the alternative test statistic, :math:`\tilde{q}_{\mu}`, as defined under the Wald approximation in Equation (62) of :xref:`arXiv:1007.1727`. @@ -126,22 +126,10 @@ def hypotest( """ init_pars = init_pars or pdf.config.suggested_init() par_bounds = par_bounds or pdf.config.suggested_bounds() - # get fixed vals from the model - model_fixed_vals = [ - (index, init) - for index, (init, is_fixed) in enumerate( - zip(init_pars, pdf.config.suggested_fixed()) - ) - if is_fixed - ] - # add user-defined ones at the end - fixed_vals = model_fixed_vals + (fixed_vals or []) - - # de-dupe and use last-appended result for each index - fixed_vals = list(dict(fixed_vals).items()) + fixed_params = fixed_params or pdf.config.suggested_fixed() calc = AsymptoticCalculator( - data, pdf, init_pars, par_bounds, fixed_vals, qtilde=qtilde + data, pdf, init_pars, par_bounds, fixed_params, qtilde=qtilde ) teststat = calc.teststatistic(poi_test) sig_plus_bkg_distribution, b_only_distribution = calc.distributions(poi_test) diff --git a/src/pyhf/infer/calculators.py b/src/pyhf/infer/calculators.py index 7f9b0a8261..fd961e05f3 100644 --- a/src/pyhf/infer/calculators.py +++ b/src/pyhf/infer/calculators.py @@ -12,7 +12,7 @@ from .test_statistics import qmu, qmu_tilde -def generate_asimov_data(asimov_mu, data, pdf, init_pars, par_bounds, fixed_vals): +def generate_asimov_data(asimov_mu, data, pdf, init_pars, par_bounds, fixed_params): """ Compute Asimov Dataset (expected yields at best-fit values) for a given POI value. @@ -22,14 +22,14 @@ def generate_asimov_data(asimov_mu, data, pdf, init_pars, par_bounds, fixed_vals pdf (~pyhf.pdf.Model): The statistical model adhering to the schema ``model.json``. init_pars (`tensor`): The initial parameter values to be used for fitting. par_bounds (`tensor`): The parameter value bounds to be used for fitting. - fixed_vals (`tensor`): Parameters to be held constant in the fit. + fixed_params (`tensor`): Parameters to be held constant in the fit. Returns: Tensor: The Asimov dataset. """ bestfit_nuisance_asimov = fixed_poi_fit( - asimov_mu, data, pdf, init_pars, par_bounds, fixed_vals + asimov_mu, data, pdf, init_pars, par_bounds, fixed_params ) return pdf.expected_data(bestfit_nuisance_asimov) @@ -128,7 +128,6 @@ def __init__( init_pars=None, par_bounds=None, fixed_params=None, - fixed_vals=None, qtilde=False, ): """ @@ -139,8 +138,7 @@ def __init__( pdf (~pyhf.pdf.Model): The statistical model adhering to the schema ``model.json``. init_pars (`tensor`): The initial parameter values to be used for fitting. par_bounds (`tensor`): The parameter value bounds to be used for fitting. - fixed_params (`tensor`): A list of booleans for each parameter on whether it should be fixed in the fit or not. - fixed_vals (list of tuples): The fixed parameters and fixed values to be used for fitting. + fixed_params (`tensor`): Whether to fix the parameter to the init_pars value during minimization qtilde (`bool`): Whether to use qtilde as the test statistic. Returns: @@ -151,22 +149,8 @@ def __init__( self.pdf = pdf self.init_pars = init_pars or pdf.config.suggested_init() self.par_bounds = par_bounds or pdf.config.suggested_bounds() + self.fixed_params = fixed_params or pdf.config.suggested_fixed() - # get fixed vals from the model - model_fixed_vals = [ - (index, init) - for index, (init, is_fixed) in enumerate( - zip(init_pars, fixed_params or pdf.config.suggested_fixed()) - ) - if is_fixed - ] - # add user-defined ones at the end - fixed_vals = model_fixed_vals + (fixed_vals or []) - - # de-dupe and use last-appended result for each index - fixed_vals = list(dict(fixed_vals).items()) - - self.fixed_vals = fixed_vals self.qtilde = qtilde self.sqrtqmuA_v = None @@ -208,7 +192,7 @@ def teststatistic(self, poi_test): self.pdf, self.init_pars, self.par_bounds, - self.fixed_vals, + self.fixed_params, ) sqrtqmu_v = tensorlib.sqrt(qmu_v) @@ -219,7 +203,7 @@ def teststatistic(self, poi_test): self.pdf, self.init_pars, self.par_bounds, - self.fixed_vals, + self.fixed_params, ) qmuA_v = teststat_func( poi_test, @@ -227,7 +211,7 @@ def teststatistic(self, poi_test): self.pdf, self.init_pars, self.par_bounds, - self.fixed_vals, + self.fixed_params, ) self.sqrtqmuA_v = tensorlib.sqrt(qmuA_v) diff --git a/src/pyhf/infer/mle.py b/src/pyhf/infer/mle.py index b5cdb8a803..8a8cff69aa 100644 --- a/src/pyhf/infer/mle.py +++ b/src/pyhf/infer/mle.py @@ -47,7 +47,7 @@ def twice_nll(pars, data, pdf): return -2 * pdf.logpdf(pars, data) -def fit(data, pdf, init_pars=None, par_bounds=None, fixed_vals=None, **kwargs): +def fit(data, pdf, init_pars=None, par_bounds=None, fixed_params=None, **kwargs): r""" Run a unconstrained maximum likelihood fit. This is done by minimizing the objective function :func:`~pyhf.infer.mle.twice_nll` @@ -96,20 +96,14 @@ def fit(data, pdf, init_pars=None, par_bounds=None, fixed_vals=None, **kwargs): _, opt = get_backend() init_pars = init_pars or pdf.config.suggested_init() par_bounds = par_bounds or pdf.config.suggested_bounds() + fixed_params = fixed_params or pdf.config.suggested_fixed() # get fixed vals from the model - model_fixed_vals = [ + fixed_vals = [ (index, init) - for index, (init, is_fixed) in enumerate( - zip(init_pars, pdf.config.suggested_fixed()) - ) + for index, (init, is_fixed) in enumerate(zip(init_pars, fixed_params)) if is_fixed ] - # add user-defined ones at the end - fixed_vals = model_fixed_vals + (fixed_vals or []) - - # de-dupe and use last-appended result for each index - fixed_vals = list(dict(fixed_vals).items()) return opt.minimize( twice_nll, data, pdf, init_pars, par_bounds, fixed_vals, **kwargs @@ -117,7 +111,7 @@ def fit(data, pdf, init_pars=None, par_bounds=None, fixed_vals=None, **kwargs): def fixed_poi_fit( - poi_val, data, pdf, init_pars=None, par_bounds=None, fixed_vals=None, **kwargs + poi_val, data, pdf, init_pars=None, par_bounds=None, fixed_params=None, **kwargs ): r""" Run a maximum likelihood fit with the POI value fixed. @@ -175,21 +169,16 @@ def fixed_poi_fit( _, opt = get_backend() init_pars = init_pars or pdf.config.suggested_init() par_bounds = par_bounds or pdf.config.suggested_bounds() + fixed_params = fixed_params or pdf.config.suggested_fixed() # get fixed vals from the model - model_fixed_vals = [ + fixed_vals = [ (index, init) - for index, (init, is_fixed) in enumerate( - zip(init_pars, pdf.config.suggested_fixed()) - ) + for index, (init, is_fixed) in enumerate(zip(init_pars, fixed_params)) if is_fixed ] - # add user-defined ones at the end - fixed_vals = model_fixed_vals + (fixed_vals or []) - # add the fixed POI fixed_vals = fixed_vals + [(pdf.config.poi_index, poi_val)] - # de-dupe and use last-appended result for each index fixed_vals = list(dict(fixed_vals).items()) diff --git a/src/pyhf/infer/test_statistics.py b/src/pyhf/infer/test_statistics.py index b4d03182a2..0f122f152f 100644 --- a/src/pyhf/infer/test_statistics.py +++ b/src/pyhf/infer/test_statistics.py @@ -7,7 +7,7 @@ log = logging.getLogger(__name__) -def _qmu_like(mu, data, pdf, init_pars, par_bounds, fixed_vals): +def _qmu_like(mu, data, pdf, init_pars, par_bounds, fixed_params): """ Clipped version of _tmu_like where the returned test statistic is 0 if muhat > 0 else tmu_like_stat. @@ -17,7 +17,7 @@ def _qmu_like(mu, data, pdf, init_pars, par_bounds, fixed_vals): """ tensorlib, optimizer = get_backend() tmu_like_stat, (_, muhatbhat) = _tmu_like( - mu, data, pdf, init_pars, par_bounds, fixed_vals, return_fitted_pars=True + 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 @@ -26,7 +26,7 @@ def _qmu_like(mu, data, pdf, init_pars, par_bounds, fixed_vals): def _tmu_like( - mu, data, pdf, init_pars, par_bounds, fixed_vals, return_fitted_pars=False + mu, data, pdf, init_pars, par_bounds, fixed_params, return_fitted_pars=False ): """ Basic Profile Likelihood test statistic. @@ -36,10 +36,10 @@ def _tmu_like( """ tensorlib, optimizer = get_backend() mubhathat, fixed_poi_fit_lhood_val = fixed_poi_fit( - mu, data, pdf, init_pars, par_bounds, fixed_vals, return_fitted_val=True + 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_vals, return_fitted_val=True + data, pdf, init_pars, par_bounds, fixed_params, return_fitted_val=True ) log_likelihood_ratio = fixed_poi_fit_lhood_val - unconstrained_fit_lhood_val tmu_like_stat = tensorlib.astensor( @@ -50,7 +50,7 @@ def _tmu_like( return tmu_like_stat -def qmu(mu, data, pdf, init_pars, par_bounds, fixed_vals): +def qmu(mu, data, pdf, init_pars, par_bounds, fixed_params): r""" The test statistic, :math:`q_{\mu}`, for establishing an upper limit on the strength parameter, :math:`\mu`, as defiend in @@ -84,8 +84,8 @@ def qmu(mu, data, pdf, init_pars, par_bounds, fixed_vals): >>> init_pars = model.config.suggested_init() >>> par_bounds = model.config.suggested_bounds() >>> par_bounds[model.config.poi_index] = [-10.0, 10.0] - >>> fixed_vals = [] - >>> pyhf.infer.test_statistics.qmu(test_mu, data, model, init_pars, par_bounds, []) + >>> fixed_params = model.config.suggested_fixed() + >>> pyhf.infer.test_statistics.qmu(test_mu, data, model, init_pars, par_bounds, fixed_params) array(3.9549891) Args: @@ -94,7 +94,7 @@ def qmu(mu, data, pdf, init_pars, par_bounds, fixed_vals): pdf (~pyhf.pdf.Model): The HistFactory statistical model used in the likelihood ratio calculation init_pars (`list`): Values to initialize the model parameters at for the fit par_bounds (`list` of `list`\s or `tuple`\s): The extrema of values the model parameters are allowed to reach in the fit - fixed_vals (`list`): Parameters held constant in the fit + fixed_params (`list`): Parameters held constant in the fit Returns: Float: The calculated test statistic, :math:`q_{\mu}` @@ -108,10 +108,10 @@ def qmu(mu, data, pdf, init_pars, par_bounds, fixed_vals): '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_vals) + return _qmu_like(mu, data, pdf, init_pars, par_bounds, fixed_params) -def qmu_tilde(mu, data, pdf, init_pars, par_bounds, fixed_vals): +def qmu_tilde(mu, data, pdf, init_pars, par_bounds, fixed_params): r""" The test statistic, :math:`\tilde{q}_{\mu}`, for establishing an upper limit on the strength parameter, :math:`\mu`, for models with @@ -150,7 +150,8 @@ def qmu_tilde(mu, data, pdf, init_pars, par_bounds, fixed_vals): >>> test_mu = 1.0 >>> init_pars = model.config.suggested_init() >>> par_bounds = model.config.suggested_bounds() - >>> pyhf.infer.test_statistics.qmu_tilde(test_mu, data, model, init_pars, par_bounds, []) + >>> fixed_params = model.config.suggested_fixed() + >>> pyhf.infer.test_statistics.qmu_tilde(test_mu, data, model, init_pars, par_bounds, fixed_params) array(3.93824492) Args: @@ -159,7 +160,7 @@ def qmu_tilde(mu, data, pdf, init_pars, par_bounds, fixed_vals): pdf (~pyhf.pdf.Model): The statistical model adhering to the schema model.json init_pars (`list`): Values to initialize the model parameters at for the fit par_bounds (`list` of `list`\s or `tuple`\s): The extrema of values the model parameters are allowed to reach in the fit - fixed_vals (`list`): Parameters held constant in the fit + fixed_params (`list`): Parameters held constant in the fit Returns: Float: The calculated test statistic, :math:`\tilde{q}_{\mu}` @@ -173,10 +174,10 @@ def qmu_tilde(mu, data, pdf, init_pars, par_bounds, fixed_vals): '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_vals) + return _qmu_like(mu, data, pdf, init_pars, par_bounds, fixed_params) -def tmu(mu, data, pdf, init_pars, par_bounds, fixed_vals): +def tmu(mu, data, pdf, init_pars, par_bounds, fixed_params): r""" The test statistic, :math:`t_{\mu}`, for establishing a two-sided interval on the strength parameter, :math:`\mu`, as defiend in Equation (8) @@ -204,8 +205,8 @@ def tmu(mu, data, pdf, init_pars, par_bounds, fixed_vals): >>> init_pars = model.config.suggested_init() >>> par_bounds = model.config.suggested_bounds() >>> par_bounds[model.config.poi_index] = [-10.0, 10.0] - >>> fixed_vals = [] - >>> pyhf.infer.test_statistics.tmu(test_mu, data, model, init_pars, par_bounds, []) + >>> fixed_params = model.config.suggested_fixed() + >>> pyhf.infer.test_statistics.tmu(test_mu, data, model, init_pars, par_bounds, fixed_params) array(3.9549891) Args: @@ -214,7 +215,7 @@ def tmu(mu, data, pdf, init_pars, par_bounds, fixed_vals): pdf (~pyhf.pdf.Model): The statistical model adhering to the schema model.json init_pars (`list`): Values to initialize the model parameters at for the fit par_bounds (`list` of `list`\s or `tuple`\s): The extrema of values the model parameters are allowed to reach in the fit - fixed_vals (`list`): Parameters held constant in the fit + fixed_params (`list`): Parameters held constant in the fit Returns: Float: The calculated test statistic, :math:`t_{\mu}` @@ -228,10 +229,10 @@ def tmu(mu, data, pdf, init_pars, par_bounds, fixed_vals): '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_vals) + return _tmu_like(mu, data, pdf, init_pars, par_bounds, fixed_params) -def tmu_tilde(mu, data, pdf, init_pars, par_bounds, fixed_vals): +def tmu_tilde(mu, data, pdf, init_pars, par_bounds, fixed_params): r""" The test statistic, :math:`\tilde{t}_{\mu}`, for establishing a two-sided interval on the strength parameter, :math:`\mu`, for models with @@ -264,8 +265,8 @@ def tmu_tilde(mu, data, pdf, init_pars, par_bounds, fixed_vals): >>> test_mu = 1.0 >>> init_pars = model.config.suggested_init() >>> par_bounds = model.config.suggested_bounds() - >>> fixed_vals = [] - >>> pyhf.infer.test_statistics.tmu_tilde(test_mu, data, model, init_pars, par_bounds, []) + >>> fixed_params = model.config.suggested_fixed() + >>> pyhf.infer.test_statistics.tmu_tilde(test_mu, data, model, init_pars, par_bounds, fixed_params) array(3.93824492) Args: @@ -274,7 +275,7 @@ def tmu_tilde(mu, data, pdf, init_pars, par_bounds, fixed_vals): pdf (~pyhf.pdf.Model): The statistical model adhering to the schema model.json init_pars (`list`): Values to initialize the model parameters at for the fit par_bounds (`list` of `list`\s or `tuple`\s): The extrema of values the model parameters are allowed to reach in the fit - fixed_vals (`list`): Parameters held constant in the fit + fixed_params (`list`): Parameters held constant in the fit Returns: Float: The calculated test statistic, :math:`\tilde{t}_{\mu}` @@ -288,4 +289,4 @@ def tmu_tilde(mu, data, pdf, init_pars, par_bounds, fixed_vals): '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_vals) + return _tmu_like(mu, data, pdf, init_pars, par_bounds, fixed_params) From 75fff6082e3ae7de8de71474648b49f3c5ace4cf Mon Sep 17 00:00:00 2001 From: Giordon Stark Date: Fri, 4 Sep 2020 11:15:04 -0400 Subject: [PATCH 22/24] update tests --- tests/test_backend_consistency.py | 2 +- tests/test_teststats.py | 34 +++++++++++++++++++++++-------- 2 files changed, 27 insertions(+), 9 deletions(-) diff --git a/tests/test_backend_consistency.py b/tests/test_backend_consistency.py index b42f7afaed..9831560d86 100644 --- a/tests/test_backend_consistency.py +++ b/tests/test_backend_consistency.py @@ -120,7 +120,7 @@ def test_hypotest_qmu_tilde( pdf, pdf.config.suggested_init(), pdf.config.suggested_bounds(), - [], + pdf.config.suggested_fixed(), ) test_statistic.append(qmu_tilde) diff --git a/tests/test_teststats.py b/tests/test_teststats.py index aca8f9f361..70090380c5 100644 --- a/tests/test_teststats.py +++ b/tests/test_teststats.py @@ -10,9 +10,12 @@ def test_qmu(caplog): data = [9] + model.config.auxdata init_pars = model.config.suggested_init() par_bounds = model.config.suggested_bounds() + fixed_params = model.config.suggested_fixed() with caplog.at_level(logging.WARNING, "pyhf.infer.test_statistics"): - pyhf.infer.test_statistics.qmu(mu, data, model, init_pars, par_bounds, []) + pyhf.infer.test_statistics.qmu( + mu, data, model, init_pars, par_bounds, fixed_params + ) assert "WARNING qmu test statistic used for fit" in caplog.text caplog.clear() @@ -23,10 +26,13 @@ def test_qmu_tilde(caplog): data = [9] + model.config.auxdata init_pars = model.config.suggested_init() par_bounds = model.config.suggested_bounds() + fixed_params = model.config.suggested_fixed() par_bounds[model.config.poi_index] = [-10, 10] with caplog.at_level(logging.WARNING, "pyhf.infer.test_statistics"): - pyhf.infer.test_statistics.qmu_tilde(mu, data, model, init_pars, par_bounds, []) + pyhf.infer.test_statistics.qmu_tilde( + mu, data, model, init_pars, par_bounds, fixed_params + ) assert "WARNING qmu_tilde test statistic used for fit" in caplog.text caplog.clear() @@ -37,8 +43,12 @@ def test_tmu(caplog): data = [9] + model.config.auxdata init_pars = model.config.suggested_init() par_bounds = model.config.suggested_bounds() + fixed_params = model.config.suggested_fixed() + with caplog.at_level(logging.WARNING, "pyhf.infer.test_statistics"): - pyhf.infer.test_statistics.tmu(mu, data, model, init_pars, par_bounds, []) + pyhf.infer.test_statistics.tmu( + mu, data, model, init_pars, par_bounds, fixed_params + ) assert "WARNING tmu test statistic used for fit" in caplog.text caplog.clear() @@ -49,10 +59,13 @@ def test_tmu_tilde(caplog): data = [9] + model.config.auxdata init_pars = model.config.suggested_init() par_bounds = model.config.suggested_bounds() + fixed_params = model.config.suggested_fixed() par_bounds[model.config.poi_index] = [-10, 10] with caplog.at_level(logging.WARNING, "pyhf.infer.test_statistics"): - pyhf.infer.test_statistics.tmu_tilde(mu, data, model, init_pars, par_bounds, []) + pyhf.infer.test_statistics.tmu_tilde( + mu, data, model, init_pars, par_bounds, fixed_params + ) assert "WARNING tmu_tilde test statistic used for fit" in caplog.text caplog.clear() @@ -84,9 +97,12 @@ def test_no_poi_test_stats(): data = [12] + model.config.auxdata init_pars = model.config.suggested_init() par_bounds = model.config.suggested_bounds() + fixed_params = model.config.suggested_fixed() with pytest.raises(pyhf.exceptions.UnspecifiedPOI) as excinfo: - pyhf.infer.test_statistics.qmu(test_poi, data, model, init_pars, par_bounds, []) + pyhf.infer.test_statistics.qmu( + test_poi, data, model, init_pars, par_bounds, fixed_params + ) assert ( "No POI is defined. A POI is required for profile likelihood based test statistics." in str(excinfo.value) @@ -94,7 +110,7 @@ def test_no_poi_test_stats(): with pytest.raises(pyhf.exceptions.UnspecifiedPOI) as excinfo: pyhf.infer.test_statistics.qmu_tilde( - test_poi, data, model, init_pars, par_bounds, [] + test_poi, data, model, init_pars, par_bounds, fixed_params ) assert ( "No POI is defined. A POI is required for profile likelihood based test statistics." @@ -102,7 +118,9 @@ def test_no_poi_test_stats(): ) with pytest.raises(pyhf.exceptions.UnspecifiedPOI) as excinfo: - pyhf.infer.test_statistics.tmu(test_poi, data, model, init_pars, par_bounds, []) + pyhf.infer.test_statistics.tmu( + test_poi, data, model, init_pars, par_bounds, fixed_params + ) assert ( "No POI is defined. A POI is required for profile likelihood based test statistics." in str(excinfo.value) @@ -110,7 +128,7 @@ def test_no_poi_test_stats(): with pytest.raises(pyhf.exceptions.UnspecifiedPOI) as excinfo: pyhf.infer.test_statistics.tmu_tilde( - test_poi, data, model, init_pars, par_bounds, [] + test_poi, data, model, init_pars, par_bounds, fixed_params ) assert ( "No POI is defined. A POI is required for profile likelihood based test statistics." From a3cbe9f67aa1c34784508dda0affb62ef5ae47f8 Mon Sep 17 00:00:00 2001 From: Giordon Stark Date: Fri, 4 Sep 2020 11:54:42 -0400 Subject: [PATCH 23/24] fix up notebooks --- docs/examples/notebooks/ImpactPlot.ipynb | 11 ++++++++++- docs/examples/notebooks/pullplot.ipynb | 11 ++++++++++- 2 files changed, 20 insertions(+), 2 deletions(-) diff --git a/docs/examples/notebooks/ImpactPlot.ipynb b/docs/examples/notebooks/ImpactPlot.ipynb index 9d8e6aa56e..e6f08bf467 100644 --- a/docs/examples/notebooks/ImpactPlot.ipynb +++ b/docs/examples/notebooks/ImpactPlot.ipynb @@ -82,8 +82,17 @@ " _, model, data = make_model([\"CRtt_meff\"])\n", "\n", " pyhf.set_backend(\"numpy\", pyhf.optimize.minuit_optimizer(verbose=True))\n", + " \n", + " constraints = constraints or []\n", + " init_pars = model.config.suggested_init()\n", + " fixed_params = model.config.suggested_fixed()\n", + " \n", + " for idx,fixed_val in constraints:\n", + " init_pars[idx] = fixed_val\n", + " fixed_params[idx] = True\n", + " \n", " result = pyhf.infer.mle.fit(\n", - " data, model, fixed_vals=constraints, return_uncertainties=True\n", + " data, model, init_pars=init_pars, fixed_params=fixed_params, return_uncertainties=True\n", " )\n", " bestfit = result[:, 0]\n", " errors = result[:, 1]\n", diff --git a/docs/examples/notebooks/pullplot.ipynb b/docs/examples/notebooks/pullplot.ipynb index 3b376dbb8a..67293991cc 100644 --- a/docs/examples/notebooks/pullplot.ipynb +++ b/docs/examples/notebooks/pullplot.ipynb @@ -72,8 +72,17 @@ " _, model, data = make_model([\"CRtt_meff\"])\n", "\n", " pyhf.set_backend(\"numpy\", pyhf.optimize.minuit_optimizer(verbose=True))\n", + " \n", + " constraints = constraints or []\n", + " init_pars = model.config.suggested_init()\n", + " fixed_params = model.config.suggested_fixed()\n", + " \n", + " for idx,fixed_val in constraints:\n", + " init_pars[idx] = fixed_val\n", + " fixed_params[idx] = True\n", + " \n", " result = pyhf.infer.mle.fit(\n", - " data, model, fixed_vals=constraints, return_uncertainties=True\n", + " data, model, init_pars=init_pars, fixed_params=fixed_params, return_uncertainties=True\n", " )\n", " bestfit = result[:, 0]\n", " errors = result[:, 1]\n", From 3364a79722287131c40dec65a77afa104cce64dd Mon Sep 17 00:00:00 2001 From: Giordon Stark Date: Fri, 4 Sep 2020 11:07:17 -0400 Subject: [PATCH 24/24] update fixed_poi_fit() to call fit() essentially --- src/pyhf/infer/mle.py | 32 ++++++++------------------------ 1 file changed, 8 insertions(+), 24 deletions(-) diff --git a/src/pyhf/infer/mle.py b/src/pyhf/infer/mle.py index 8a8cff69aa..a5c544ae26 100644 --- a/src/pyhf/infer/mle.py +++ b/src/pyhf/infer/mle.py @@ -49,7 +49,7 @@ def twice_nll(pars, data, pdf): def fit(data, pdf, init_pars=None, par_bounds=None, fixed_params=None, **kwargs): r""" - Run a unconstrained maximum likelihood fit. + Run a maximum likelihood fit. This is done by minimizing the objective function :func:`~pyhf.infer.mle.twice_nll` of the model parameters given the observed data. This is used to produce the maximal likelihood :math:`L\left(\hat{\mu}, \hat{\boldsymbol{\theta}}\right)` @@ -87,6 +87,7 @@ def fit(data, pdf, init_pars=None, par_bounds=None, fixed_params=None, **kwargs) pdf (~pyhf.pdf.Model): The statistical model adhering to the schema model.json init_pars (`list`): Values to initialize the model parameters at for the fit par_bounds (`list` of `list`\s or `tuple`\s): The extrema of values the model parameters are allowed to reach in the fit + fixed_params (`list`): Parameters to be held constant in the fit. kwargs: Keyword arguments passed through to the optimizer API Returns: @@ -155,6 +156,7 @@ def fixed_poi_fit( pdf (~pyhf.pdf.Model): The statistical model adhering to the schema model.json init_pars (`list`): Values to initialize the model parameters at for the fit par_bounds (`list` of `list`\s or `tuple`\s): The extrema of values the model parameters are allowed to reach in the fit + fixed_params (`list`): Parameters to be held constant in the fit. kwargs: Keyword arguments passed through to the optimizer API Returns: @@ -166,28 +168,10 @@ def fixed_poi_fit( 'No POI is defined. A POI is required to fit with a fixed POI.' ) - _, opt = get_backend() - init_pars = init_pars or pdf.config.suggested_init() - par_bounds = par_bounds or pdf.config.suggested_bounds() - fixed_params = fixed_params or pdf.config.suggested_fixed() + init_pars = [*(init_pars or pdf.config.suggested_init())] + fixed_params = [*(fixed_params or pdf.config.suggested_fixed())] - # get fixed vals from the model - fixed_vals = [ - (index, init) - for index, (init, is_fixed) in enumerate(zip(init_pars, fixed_params)) - if is_fixed - ] - # add the fixed POI - fixed_vals = fixed_vals + [(pdf.config.poi_index, poi_val)] - # de-dupe and use last-appended result for each index - fixed_vals = list(dict(fixed_vals).items()) + init_pars[pdf.config.poi_index] = poi_val + fixed_params[pdf.config.poi_index] = True - return opt.minimize( - twice_nll, - data, - pdf, - init_pars, - par_bounds, - fixed_vals, - **kwargs, - ) + return fit(data, pdf, init_pars, par_bounds, fixed_params, **kwargs)