Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

feat: Test statistic specified by string instead of using qtilde kwarg #1231

Merged
merged 20 commits into from
Dec 19, 2020
Merged
Show file tree
Hide file tree
Changes from 12 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 1 addition & 3 deletions src/pyhf/cli/infer.py
Original file line number Diff line number Diff line change
Expand Up @@ -192,8 +192,6 @@ def cls(

ws = Workspace(spec)

is_qtilde = teststat == 'qtilde'

patches = [json.loads(click.open_file(pfile, 'r').read()) for pfile in patch]
model = ws.model(
measurement_name=measurement,
Expand Down Expand Up @@ -226,7 +224,7 @@ def cls(
testpoi,
ws.data(model),
model,
qtilde=is_qtilde,
test_stat=teststat,
calctype=calctype,
return_expected_set=True,
)
Expand Down
6 changes: 6 additions & 0 deletions src/pyhf/exceptions/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -91,6 +91,12 @@ class InvalidInterpCode(Exception):
"""


class InvalidTestStatistic(Exception):
"""
InvalidTestStatistic is raised when an invalid/unimplemented test statistic is requested.
"""


class ImportBackendError(Exception):
"""
MissingLibraries is raised when something is imported by sustained an import error due to missing additional, non-default libraries.
Expand Down
2 changes: 1 addition & 1 deletion src/pyhf/infer/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ def hypotest(
>>> data = pyhf.tensorlib.astensor(observations + model.config.auxdata)
>>> mu_test = 1.0
>>> CLs_obs, CLs_exp_band = pyhf.infer.hypotest(
... mu_test, data, model, return_expected_set=True, qtilde=True
... mu_test, data, model, return_expected_set=True, test_stat="qtilde"
... )
>>> CLs_obs
array(0.05251497)
Expand Down
53 changes: 42 additions & 11 deletions src/pyhf/infer/calculators.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,9 +9,13 @@
"""
from .mle import fixed_poi_fit
from .. import get_backend
from .test_statistics import qmu, qmu_tilde
from . import test_statistics
import tqdm

import logging

log = logging.getLogger(__name__)


def generate_asimov_data(asimov_mu, data, pdf, init_pars, par_bounds, fixed_params):
"""
Expand Down Expand Up @@ -159,17 +163,22 @@ def __init__(
init_pars=None,
par_bounds=None,
fixed_params=None,
qtilde=True,
test_stat="qtilde",
**kwargs,
):
r"""
Asymptotic Calculator.

.. deprecated:: 0.6
``qtilde=True`` has been deprecated in favor of ``test_stat="qtilde"``.

Args:
data (:obj:`tensor`): The observed data.
pdf (~pyhf.pdf.Model): The statistical model adhering to the schema ``model.json``.
init_pars (:obj:`tensor`): The initial parameter values to be used for fitting.
par_bounds (:obj:`tensor`): The parameter value bounds to be used for fitting.
fixed_params (:obj:`tensor`): Whether to fix the parameter to the init_pars value during minimization
test_stat (:obj:`str`): The test statistic to use as a numerical summary of the data.
qtilde (:obj:`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`
Expand All @@ -186,7 +195,15 @@ def __init__(
self.par_bounds = par_bounds or pdf.config.suggested_bounds()
self.fixed_params = fixed_params or pdf.config.suggested_fixed()

self.qtilde = qtilde
# handle deprecation for qtilde gracefully
qtilde = kwargs.pop('qtilde', None)
if qtilde is not None:
test_stat = "qtilde" if qtilde else "q"
log.warning(
f"Setting qtilde={qtilde} is deprecated. Use test_stat='{test_stat}' instead."
)

self.test_stat = test_stat
self.sqrtqmuA_v = None

def distributions(self, poi_test):
Expand All @@ -205,7 +222,7 @@ def distributions(self, poi_test):
>>> observations = [51, 48]
>>> data = observations + model.config.auxdata
>>> mu_test = 1.0
>>> asymptotic_calculator = pyhf.infer.calculators.AsymptoticCalculator(data, model, qtilde=True)
>>> asymptotic_calculator = pyhf.infer.calculators.AsymptoticCalculator(data, model, test_stat="qtilde")
>>> _ = asymptotic_calculator.teststatistic(mu_test)
>>> qmu_sig, qmu_bkg = asymptotic_calculator.distributions(mu_test)
>>> qmu_sig.pvalue(mu_test), qmu_bkg.pvalue(mu_test)
Expand Down Expand Up @@ -238,7 +255,7 @@ def teststatistic(self, poi_test):
>>> observations = [51, 48]
>>> data = observations + model.config.auxdata
>>> mu_test = 1.0
>>> asymptotic_calculator = pyhf.infer.calculators.AsymptoticCalculator(data, model, qtilde=True)
>>> asymptotic_calculator = pyhf.infer.calculators.AsymptoticCalculator(data, model, test_stat="qtilde")
>>> asymptotic_calculator.teststatistic(mu_test)
0.14043184405388176

Expand All @@ -251,7 +268,7 @@ def teststatistic(self, poi_test):
"""
tensorlib, _ = get_backend()

teststat_func = qmu_tilde if self.qtilde else qmu
teststat_func = test_statistics.get(self.test_stat)

qmu_v = teststat_func(
poi_test,
Expand Down Expand Up @@ -282,7 +299,7 @@ def teststatistic(self, poi_test):
)
self.sqrtqmuA_v = tensorlib.sqrt(qmuA_v)

if not self.qtilde: # qmu
if self.test_stat == "q":
teststat = sqrtqmu_v - self.sqrtqmuA_v
else: # qtilde

Expand Down Expand Up @@ -449,19 +466,24 @@ def __init__(
init_pars=None,
par_bounds=None,
fixed_params=None,
qtilde=True,
test_stat="qtilde",
ntoys=2000,
track_progress=True,
**kwargs,
):
r"""
Toy-based Calculator.

.. deprecated:: 0.6
``qtilde=True`` has been deprecated in favor of ``test_stat="qtilde"``.

Args:
data (:obj:`tensor`): The observed data.
pdf (~pyhf.pdf.Model): The statistical model adhering to the schema ``model.json``.
init_pars (:obj:`tensor`): The initial parameter values to be used for fitting.
par_bounds (:obj:`tensor`): The parameter value bounds to be used for fitting.
fixed_params (:obj:`tensor`): Whether to fix the parameter to the init_pars value during minimization
test_stat (:obj:`str`): The test statistic to use as a numerical summary of the data.
qtilde (:obj:`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`
Expand All @@ -480,7 +502,16 @@ def __init__(
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()
self.qtilde = qtilde

# handle deprecation for qtilde gracefully
qtilde = kwargs.pop('qtilde', None)
if qtilde is not None:
test_stat = "qtilde" if qtilde else "q"
log.warning(
f"Setting qtilde={qtilde} is deprecated. Use test_stat='{test_stat}' instead."
)

self.test_stat = test_stat
self.track_progress = track_progress

def distributions(self, poi_test, track_progress=None):
Expand Down Expand Up @@ -527,7 +558,7 @@ def distributions(self, poi_test, track_progress=None):
bkg_pdf = self.pdf.make_pdf(tensorlib.astensor(bkg_pars))
bkg_sample = bkg_pdf.sample(sample_shape)

teststat_func = qmu_tilde if self.qtilde else qmu
teststat_func = test_statistics.get(self.test_stat)

tqdm_options = dict(
total=self.ntoys,
Expand Down Expand Up @@ -597,7 +628,7 @@ def teststatistic(self, poi_test):
Float: The value of the test statistic.

"""
teststat_func = qmu_tilde if self.qtilde else qmu
teststat_func = test_statistics.get(self.test_stat)
teststat = teststat_func(
poi_test,
self.data,
Expand Down
3 changes: 2 additions & 1 deletion src/pyhf/infer/intervals.py
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,8 @@ def upperlimit(data, model, scan, level=0.05, return_results=False):
"""
tb, _ = get_backend()
results = [
hypotest(mu, data, model, qtilde=True, return_expected_set=True) for mu in scan
hypotest(mu, data, model, test_stat="qtilde", return_expected_set=True)
for mu in scan
]
obs = tb.astensor([[r[0]] for r in results])
exp = tb.astensor([[r[1][idx] for idx in range(5)] for r in results])
Expand Down
31 changes: 30 additions & 1 deletion src/pyhf/infer/test_statistics.py
Original file line number Diff line number Diff line change
@@ -1,12 +1,41 @@
from .. import get_backend
from .mle import fixed_poi_fit, fit
from ..exceptions import UnspecifiedPOI
from ..exceptions import UnspecifiedPOI, InvalidTestStatistic

import logging

log = logging.getLogger(__name__)


def get(name):
"""
Get the test statistic function by name.

Example:

>>> import pyhf
>>> pyhf.infer.test_statistics.get("q")
<function qmu at 0x...>
>>> pyhf.infer.test_statistics.get("qtilde")
<function qmu_tilde at 0x...>

Args:
name (:obj:`str`): The name of the test statistic to retrieve


Returns:
callable: The test statistic function
"""
_mapping = {
"q": qmu,
"qtilde": qmu_tilde,
}
try:
return _mapping[name]
except KeyError:
raise InvalidTestStatistic


def _qmu_like(mu, data, pdf, init_pars, par_bounds, fixed_params):
"""
Clipped version of _tmu_like where the returned test statistic
Expand Down
2 changes: 1 addition & 1 deletion src/pyhf/infer/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ def create_calculator(calctype, *args, **kwargs):
>>> data = observations + model.config.auxdata
>>> mu_test = 1.0
>>> toy_calculator = pyhf.infer.utils.create_calculator(
... "toybased", data, model, ntoys=100, qtilde=True, track_progress=False
... "toybased", data, model, ntoys=100, test_stat="qtilde", track_progress=False
... )
>>> qmu_sig, qmu_bkg = toy_calculator.distributions(mu_test)
>>> qmu_sig.pvalue(mu_test), qmu_bkg.pvalue(mu_test)
Expand Down
12 changes: 12 additions & 0 deletions tests/test_calculator.py
Original file line number Diff line number Diff line change
@@ -1,7 +1,19 @@
import pyhf
import pyhf.infer.calculators
import pytest
import logging


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("calculator", ["asymptotics", "toybased"])
@pytest.mark.parametrize("qtilde", [True, False])
def test_deprecated_qtilde(caplog, mocker, calculator, qtilde):
with caplog.at_level(logging.WARNING):
pyhf.infer.utils.create_calculator(
calculator, ['fake data'], mocker.Mock(), qtilde=qtilde
)
assert "is deprecated. Use test_stat" in caplog.text
12 changes: 6 additions & 6 deletions tests/test_infer.py
Original file line number Diff line number Diff line change
Expand Up @@ -25,9 +25,9 @@ def test_upperlimit(tmpdir, hypotest_args):
Check that the default return structure of pyhf.infer.hypotest is as expected
"""
_, data, model = hypotest_args
observed_limit, expected_limits = pyhf.infer.intervals.upperlimit(
data, model, scan=np.linspace(0, 5, 11)
)
results = pyhf.infer.intervals.upperlimit(data, model, scan=np.linspace(0, 5, 11))
assert len(results) == 2
observed_limit, expected_limits = results
assert observed_limit == pytest.approx(1.0262704738584554)
assert expected_limits == pytest.approx(
[0.65765653, 0.87999725, 1.12453992, 1.50243428, 2.09232927]
Expand Down Expand Up @@ -230,10 +230,10 @@ def logpdf(self, pars, data):
assert np.isclose(cls, 0.7267836451638846)


@pytest.mark.parametrize("qtilde", [True, False])
def test_calculator_distributions_without_teststatistic(qtilde):
@pytest.mark.parametrize("test_stat", ["qtilde", "q"])
def test_calculator_distributions_without_teststatistic(test_stat):
calc = pyhf.infer.calculators.AsymptoticCalculator(
[0.0], {}, [1.0], [(0.0, 10.0)], [False], qtilde=qtilde
[0.0], {}, [1.0], [(0.0, 10.0)], [False], test_stat=test_stat
)
with pytest.raises(RuntimeError):
calc.distributions(1.0)
Expand Down
2 changes: 1 addition & 1 deletion tests/test_regression.py
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ def calculate_CLs(bkgonly_json, signal_patch_json):
},
)
result = pyhf.infer.hypotest(
1.0, workspace.data(model), model, qtilde=True, return_expected_set=True
1.0, workspace.data(model), model, test_stat="qtilde", return_expected_set=True
)
return result[0].tolist(), result[-1]

Expand Down
10 changes: 10 additions & 0 deletions tests/test_teststats.py
Original file line number Diff line number Diff line change
Expand Up @@ -134,3 +134,13 @@ def test_no_poi_test_stats():
"No POI is defined. A POI is required for profile likelihood based test statistics."
in str(excinfo.value)
)


@pytest.mark.parametrize("test_stat", ["qtilde", "q"])
def test_get_teststat_by_name(test_stat):
assert pyhf.infer.test_statistics.get(test_stat)


def test_get_teststat_error():
with pytest.raises(pyhf.exceptions.InvalidTestStatistic):
pyhf.infer.test_statistics.get("look at me i'm not real")
2 changes: 1 addition & 1 deletion tests/test_validation.py
Original file line number Diff line number Diff line change
Expand Up @@ -576,7 +576,7 @@ def validate_hypotest(pdf, data, mu_test, expected_result, tolerance=1e-6):
init_pars,
par_bounds,
return_expected_set=True,
qtilde=False,
test_stat="q",
)
assert abs(CLs_obs - expected_result['obs']) / expected_result['obs'] < tolerance
for result, expected in zip(CLs_exp_set, expected_result['exp']):
Expand Down