adjust test

src/pyhf/infer/calculators.py

@@ -0,0 +1,66 @@
+from .mle import fixed_poi_fit
+from .. import get_backend
+from .test_statistics import qmu
+
+
+def generate_asimov_data(asimov_mu, data, pdf, init_pars, par_bounds):
+    """Compute Asimov Dataset (expected yields at best-fit values) for a given POI value."""
+    bestfit_nuisance_asimov = fixed_poi_fit(asimov_mu, data, pdf, init_pars, par_bounds)
+    return pdf.expected_data(bestfit_nuisance_asimov)
+
+
+class AsymptoticTestStatDistribution(object):
+    def __init__(self, shift):
+        self.shift = shift
+
+    def pvalue(self, value):
+        tensorlib, _ = get_backend()
+        return 1 - tensorlib.normal_cdf(value - self.shift)
+
+    def expected_value(self, nsigma):
+        return nsigma
+
+
+class AsymptoticCalculator(object):
+    def __init__(self, data, pdf, init_pars=None, par_bounds=None, qtilde=False):
+        self.data = data
+        self.pdf = pdf
+        self.init_pars = init_pars or pdf.config.suggested_init()
+        self.par_bounds = par_bounds or pdf.config.suggested_bounds()
+        self.qtilde = qtilde
+
+    def distributions(self, poi_test):
+        tensorlib, _ = get_backend()
+        asimov_mu = 0.0
+        asimov_data = generate_asimov_data(
+            asimov_mu, self.data, self.pdf, self.init_pars, self.par_bounds
+        )
+        qmuA_v = qmu(poi_test, asimov_data, self.pdf, self.init_pars, self.par_bounds)
+        self.sqrtqmuA_v = tensorlib.sqrt(qmuA_v)
+
+        sb_dist = AsymptoticTestStatDistribution(-self.sqrtqmuA_v)
+        b_dist = AsymptoticTestStatDistribution(0.0)
+        return sb_dist, b_dist
+
+    def teststatistic(self, poi_test):
+        tensorlib, _ = get_backend()
+        qmu_v = qmu(poi_test, self.data, self.pdf, self.init_pars, self.par_bounds)
+        sqrtqmu_v = tensorlib.sqrt(qmu_v)
+        if not self.qtilde:  # qmu
+            teststat = sqrtqmu_v - self.sqrtqmuA_v
+        else:  # qtilde
+
+            def _true_case():
+                teststat = sqrtqmu_v - self.sqrtqmuA_v
+                return teststat
+
+            def _false_case():
+                qmu = tensorlib.power(sqrtqmu_v, 2)
+                qmu_A = tensorlib.power(self.sqrtqmuA_v, 2)
+                teststat = (qmu - qmu_A) / (2 * self.sqrtqmuA_v)
+                return teststat
+
+            teststat = tensorlib.conditional(
+                (sqrtqmu_v < self.sqrtqmuA_v), _true_case, _false_case
+            )
+        return teststat