"""Interval estimation"""
import numpy as np
from scipy.optimize import toms748
from pyhf import get_backend
from pyhf.infer import hypotest
__all__ = ["upper_limit", "linear_grid_scan", "toms748_scan"]
def __dir__():
return __all__
def _interp(x, xp, fp):
tb, _ = get_backend()
return tb.astensor(np.interp(x, xp.tolist(), fp.tolist()))
[docs]
def toms748_scan(
data,
model,
bounds_low,
bounds_up,
level=0.05,
atol=2e-12,
rtol=1e-4,
from_upper_limit_fn=False,
**hypotest_kwargs,
):
"""
Calculate an upper limit interval ``(0, poi_up)`` for a single
Parameter of Interest (POI) using an automatic scan through
POI-space, using the :func:`~scipy.optimize.toms748` algorithm.
Example:
>>> import numpy as np
>>> import pyhf
>>> pyhf.set_backend("numpy")
>>> model = pyhf.simplemodels.uncorrelated_background(
... signal=[12.0, 11.0], bkg=[50.0, 52.0], bkg_uncertainty=[3.0, 7.0]
... )
>>> observations = [51, 48]
>>> data = pyhf.tensorlib.astensor(observations + model.config.auxdata)
>>> obs_limit, exp_limits = pyhf.infer.intervals.upper_limits.toms748_scan(
... data, model, 0., 5., rtol=0.01
... )
>>> obs_limit
array(1.01156939)
>>> exp_limits
[array(0.5600747), array(0.75702605), array(1.06234693), array(1.50116923), array(2.05078912)]
Args:
data (:obj:`tensor`): The observed data.
model (~pyhf.pdf.Model): The statistical model adhering to the schema ``model.json``.
bounds_low (:obj:`float`): Lower boundary of search interval.
bounds_up (:obj:`float`): Upper boundary of search interval.
level (:obj:`float`): The threshold value to evaluate the interpolated results at.
Defaults to ``0.05``.
atol (:obj:`float`): Absolute tolerance.
The iteration will end when the result is within absolute
*or* relative tolerance of the true limit.
rtol (:obj:`float`): Relative tolerance.
For optimal performance this argument should be set
to the highest acceptable relative tolerance.
hypotest_kwargs (:obj:`string`): Kwargs for the calls to
:class:`~pyhf.infer.hypotest` to configure the fits.
Returns:
Tuple of Tensors:
- Tensor: The observed upper limit on the POI.
- Tensor: The expected upper limits on the POI.
.. versionadded:: 0.7.0
"""
cache = {}
def f_cached(poi):
if poi not in cache:
cache[poi] = hypotest(
poi,
data,
model,
return_expected_set=True,
**hypotest_kwargs,
)
return cache[poi]
def f(poi, level, limit=0):
# Use integers for limit so we don't need a string comparison
# limit == 0: Observed
# else: expected
return (
f_cached(poi)[0] - level
if limit == 0
else f_cached(poi)[1][limit - 1] - level
)
def best_bracket(limit):
# return best bracket
ks = np.asarray(list(cache))
vals = np.asarray(
[
value[0] - level if limit == 0 else value[1][limit - 1] - level
for value in cache.values()
]
)
pos = vals >= 0
neg = vals < 0
lower = ks[pos][np.argmin(vals[pos])]
upper = ks[neg][np.argmax(vals[neg])]
return (lower, upper)
# extend bounds_low and bounds_up if they don't bracket CLs level
lower_results = f_cached(bounds_low)
# {lower,upper}_results[0] is an array and {lower,upper}_results[1] is a
# list of arrays so need to turn {lower,upper}_results[0] into list to
# concatenate them
while np.any(np.asarray([lower_results[0]] + lower_results[1]) < level):
bounds_low /= 2
lower_results = f_cached(bounds_low)
upper_results = f_cached(bounds_up)
while np.any(np.asarray([upper_results[0]] + upper_results[1]) > level):
bounds_up *= 2
upper_results = f_cached(bounds_up)
tb, _ = get_backend()
obs = tb.astensor(
toms748(f, bounds_low, bounds_up, args=(level, 0), k=2, xtol=atol, rtol=rtol)
)
exp = [
tb.astensor(
toms748(f, *best_bracket(idx), args=(level, idx), k=2, xtol=atol, rtol=rtol)
)
for idx in range(1, 6)
]
if from_upper_limit_fn:
return obs, exp, (list(cache), list(cache.values()))
return obs, exp
[docs]
def linear_grid_scan(
data, model, scan, level=0.05, return_results=False, **hypotest_kwargs
):
"""
Calculate an upper limit interval ``(0, poi_up)`` for a single
Parameter of Interest (POI) using a linear scan through POI-space.
Example:
>>> import numpy as np
>>> import pyhf
>>> pyhf.set_backend("numpy")
>>> model = pyhf.simplemodels.uncorrelated_background(
... signal=[12.0, 11.0], bkg=[50.0, 52.0], bkg_uncertainty=[3.0, 7.0]
... )
>>> observations = [51, 48]
>>> data = pyhf.tensorlib.astensor(observations + model.config.auxdata)
>>> scan = np.linspace(0, 5, 21)
>>> obs_limit, exp_limits, (scan, results) = pyhf.infer.intervals.upper_limits.upper_limit(
... data, model, scan, return_results=True
... )
>>> obs_limit
array(1.01764089)
>>> exp_limits
[array(0.59576921), array(0.76169166), array(1.08504773), array(1.50170482), array(2.06654952)]
Args:
data (:obj:`tensor`): The observed data.
model (~pyhf.pdf.Model): The statistical model adhering to the schema ``model.json``.
scan (:obj:`iterable`): Iterable of POI values.
level (:obj:`float`): The threshold value to evaluate the interpolated results at.
return_results (:obj:`bool`): Whether to return the per-point results.
hypotest_kwargs (:obj:`string`): Kwargs for the calls to
:class:`~pyhf.infer.hypotest` to configure the fits.
Returns:
Tuple of Tensors:
- Tensor: The observed upper limit on the POI.
- Tensor: The expected upper limits on the POI.
- Tuple of Tensors: The given ``scan`` along with the
:class:`~pyhf.infer.hypotest` results at each test POI.
Only returned when ``return_results`` is ``True``.
.. versionadded:: 0.7.0
"""
tb, _ = get_backend()
results = [
hypotest(mu, data, model, return_expected_set=True, **hypotest_kwargs)
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])
result_array = tb.concatenate([obs, exp], axis=1).T
# observed limit and the (0, +-1, +-2)sigma expected limits
limits = [_interp(level, result_array[idx][::-1], scan[::-1]) for idx in range(6)]
obs_limit, exp_limits = limits[0], limits[1:]
if return_results:
return obs_limit, exp_limits, (scan, results)
return obs_limit, exp_limits
[docs]
def upper_limit(
data, model, scan=None, level=0.05, return_results=False, **hypotest_kwargs
):
"""
Calculate an upper limit interval ``(0, poi_up)`` for a single Parameter of
Interest (POI) using root-finding or a linear scan through POI-space.
Example:
>>> import numpy as np
>>> import pyhf
>>> pyhf.set_backend("numpy")
>>> model = pyhf.simplemodels.uncorrelated_background(
... signal=[12.0, 11.0], bkg=[50.0, 52.0], bkg_uncertainty=[3.0, 7.0]
... )
>>> observations = [51, 48]
>>> data = pyhf.tensorlib.astensor(observations + model.config.auxdata)
>>> scan = np.linspace(0, 5, 21)
>>> obs_limit, exp_limits, (scan, results) = pyhf.infer.intervals.upper_limits.upper_limit(
... data, model, scan, return_results=True
... )
>>> obs_limit
array(1.01764089)
>>> exp_limits
[array(0.59576921), array(0.76169166), array(1.08504773), array(1.50170482), array(2.06654952)]
Args:
data (:obj:`tensor`): The observed data.
model (~pyhf.pdf.Model): The statistical model adhering to the schema ``model.json``.
scan (:obj:`iterable` or ``None``): Iterable of POI values or ``None`` to use
:class:`~pyhf.infer.intervals.upper_limits.toms748_scan`.
level (:obj:`float`): The threshold value to evaluate the interpolated results at.
return_results (:obj:`bool`): Whether to return the per-point results.
Returns:
Tuple of Tensors:
- Tensor: The observed upper limit on the POI.
- Tensor: The expected upper limits on the POI.
- Tuple of Tensors: The given ``scan`` along with the
:class:`~pyhf.infer.hypotest` results at each test POI.
Only returned when ``return_results`` is ``True``.
.. versionadded:: 0.7.0
"""
if scan is not None:
return linear_grid_scan(
data, model, scan, level, return_results, **hypotest_kwargs
)
# else:
bounds = model.config.suggested_bounds()[
model.config.par_slice(model.config.poi_name).start
]
obs_limit, exp_limit, results = toms748_scan(
data,
model,
bounds[0],
bounds[1],
from_upper_limit_fn=True,
**hypotest_kwargs,
)
if return_results:
return obs_limit, exp_limit, results
return obs_limit, exp_limit