fix: Use MLEs of NPs to create sampling distributions in ToyCalculator (

#1610)

* Use the best fit parameter values of pyhf.infer.mle.fixed_poi_fit as the values of the nuisance parameters when creating the sampling distributions to generate pseudo-experiments with the ToyCalculator
   - This implements the recommendations from ATLAS and CMS in "Procedure for the LHC Higgs boson search combination in Summer 2011" (Section 2.1 'Observed limits')
   - https://inspirehep.net/literature/1196797
   - Also matches the description of RooStats::FrequentistCalculator where "The nuisance parameters are fixed to their MLEs"
* Update test values in response to the difference in values from pseudo-experiments
* Add Mason Proffitt to the list of contributors

docs/contributors.rst

@@ -24,3 +24,4 @@ Contributors include:
 - Henry Schreiner
 - Saransh Chopra
 - Sviatoslav Sydorenko
+- Mason Proffitt

src/pyhf/infer/calculators.py

@@ -667,7 +667,16 @@ def __init__(
 
     def distributions(self, poi_test, track_progress=None):
         """
-        Probability Distributions of the test statistic value under the signal + background and background-only hypothesis.
+        Probability distributions of the test statistic value under the signal + background and background-only hypotheses.
+
+        These distributions are produced by generating pseudo-data ("toys")
+        with the nuisance parameters set to their conditional maximum likelihood
+        estimators at the corresponding value of the parameter of interest for
+        each hypothesis, following the joint recommendations of the ATLAS and CMS
+        experiments in |LHC Higgs search combination procedure|_.
+
+        .. _LHC Higgs search combination procedure: https://inspirehep.net/literature/1196797
+        .. |LHC Higgs search combination procedure| replace:: *Procedure for the LHC Higgs boson search combination in Summer 2011*
 
         Example:
 
@@ -686,7 +695,7 @@ def distributions(self, poi_test, track_progress=None):
             ... )
             >>> sig_plus_bkg_dist, bkg_dist = toy_calculator.distributions(mu_test)
             >>> sig_plus_bkg_dist.pvalue(mu_test), bkg_dist.pvalue(mu_test)
-            (array(0.14), array(0.76))
+            (array(0.14), array(0.79))
 
         Args:
             poi_test (:obj:`float` or :obj:`tensor`): The value for the parameter of interest.
@@ -699,14 +708,26 @@ def distributions(self, poi_test, track_progress=None):
         tensorlib, _ = get_backend()
         sample_shape = (self.ntoys,)
 
-        signal_pars = self.pdf.config.suggested_init()
-        signal_pars[self.pdf.config.poi_index] = poi_test
-        signal_pdf = self.pdf.make_pdf(tensorlib.astensor(signal_pars))
+        signal_pars = fixed_poi_fit(
+            poi_test,
+            self.data,
+            self.pdf,
+            self.init_pars,
+            self.par_bounds,
+            self.fixed_params,
+        )
+        signal_pdf = self.pdf.make_pdf(signal_pars)
         signal_sample = signal_pdf.sample(sample_shape)
 
-        bkg_pars = self.pdf.config.suggested_init()
-        bkg_pars[self.pdf.config.poi_index] = 1.0 if self.test_stat == 'q0' else 0.0
-        bkg_pdf = self.pdf.make_pdf(tensorlib.astensor(bkg_pars))
+        bkg_pars = fixed_poi_fit(
+            1.0 if self.test_stat == 'q0' else 0.0,
+            self.data,
+            self.pdf,
+            self.init_pars,
+            self.par_bounds,
+            self.fixed_params,
+        )
+        bkg_pdf = self.pdf.make_pdf(bkg_pars)
         bkg_sample = bkg_pdf.sample(sample_shape)
 
         teststat_func = utils.get_test_stat(self.test_stat)
@@ -774,7 +795,7 @@ def pvalues(self, teststat, sig_plus_bkg_distribution, bkg_only_distribution):
             >>> sig_plus_bkg_dist, bkg_dist = toy_calculator.distributions(mu_test)
             >>> CLsb, CLb, CLs = toy_calculator.pvalues(q_tilde, sig_plus_bkg_dist, bkg_dist)
             >>> CLsb, CLb, CLs
-            (array(0.01), array(0.41), array(0.02439024))
+            (array(0.03), array(0.37), array(0.08108108))
 
         Args:
             teststat (:obj:`tensor`): The test statistic.
@@ -820,7 +841,7 @@ def expected_pvalues(self, sig_plus_bkg_distribution, bkg_only_distribution):
             >>> sig_plus_bkg_dist, bkg_dist = toy_calculator.distributions(mu_test)
             >>> CLsb_exp_band, CLb_exp_band, CLs_exp_band = toy_calculator.expected_pvalues(sig_plus_bkg_dist, bkg_dist)
             >>> CLs_exp_band
-            [array(0.), array(0.), array(0.06186224), array(0.28450033), array(1.)]
+            [array(0.), array(0.), array(0.08403955), array(0.21892596), array(0.86072977)]
 
         Args:
             sig_plus_bkg_distribution (~pyhf.infer.calculators.EmpiricalDistribution):

src/pyhf/infer/utils.py

@@ -57,7 +57,7 @@ def create_calculator(calctype, *args, **kwargs):
         ... )
         >>> qmu_sig, qmu_bkg = toy_calculator.distributions(mu_test)
         >>> qmu_sig.pvalue(mu_test), qmu_bkg.pvalue(mu_test)
-        (array(0.14), array(0.76))
+        (array(0.14), array(0.79))
 
     Args:
         calctype (:obj:`str`): The calculator to create. Choose either

tests/test_infer.py

@@ -365,30 +365,30 @@ def test_toy_calculator(tmpdir, hypotest_args):
     assert qtilde_mu_sig.samples.tolist() == pytest.approx(
         [
             0.0,
-            0.13298492825293806,
-            0.0,
-            0.7718560148925349,
-            1.814884694401428,
+            0.017350013494649374,
             0.0,
+            0.2338008822475217,
+            0.020328779776718875,
+            0.8911134903562186,
+            0.04408274703718007,
             0.0,
+            0.03977591672014569,
             0.0,
-            0.0,
-            0.06586643485326249,
         ],
         1e-07,
     )
     assert qtilde_mu_bkg.samples.tolist() == pytest.approx(
         [
-            2.2664625749100082,
-            1.081660887453154,
-            2.7570218408936853,
-            1.3835691388297846,
-            0.4707467005909507,
-            0.0,
-            3.7166483705294127,
-            3.8021896732709592,
-            5.114135391143066,
-            1.3511153731000718,
+            5.642956861215396,
+            0.37581364290284114,
+            4.875367689039649,
+            3.4299006094989295,
+            1.0161021805475343,
+            0.03345317321810626,
+            0.21984803001140563,
+            1.274869119189077,
+            9.368264062021098,
+            3.0716486684082156,
         ],
         1e-07,
     )

tests/test_validation.py

@@ -125,8 +125,8 @@ def expected_result_1bin_shapesys_q0():
 @pytest.fixture(scope='module')
 def expected_result_1bin_shapesys_q0_toys():
     expected_result = {
-        "exp": [0.0, 0.0005, 0.0145, 0.1205, 0.402761],
-        "obs": 0.005,
+        "exp": [0.0, 0.0, 0.0135, 0.1365, 0.39854497],
+        "obs": 0.004,
     }
     return expected_result