from pyhf import get_backend
from pyhf.infer.mle import fixed_poi_fit, fit
from pyhf.exceptions import UnspecifiedPOI
import logging
log = logging.getLogger(__name__)
__all__ = ["q0", "qmu", "qmu_tilde", "tmu", "tmu_tilde"]
def __dir__():
return __all__
def _qmu_like(
mu, data, pdf, init_pars, par_bounds, fixed_params, return_fitted_pars=False
):
"""
Clipped version of _tmu_like where the returned test statistic
is 0 if muhat > 0 else tmu_like_stat.
If the lower bound of the POI is 0 this automatically implements
qmu_tilde. Otherwise this is qmu (no tilde).
"""
tensorlib, optimizer = get_backend()
tmu_like_stat, (mubhathat, muhatbhat) = _tmu_like(
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
)
if return_fitted_pars:
return qmu_like_stat, (mubhathat, muhatbhat)
return qmu_like_stat
def _tmu_like(
mu, data, pdf, init_pars, par_bounds, fixed_params, return_fitted_pars=False
):
"""
Basic Profile Likelihood test statistic.
If the lower bound of the POI is 0 this automatically implements
tmu_tilde. Otherwise this is tmu (no tilde).
"""
tensorlib, optimizer = get_backend()
mubhathat, fixed_poi_fit_lhood_val = fixed_poi_fit(
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_params, return_fitted_val=True
)
log_likelihood_ratio = fixed_poi_fit_lhood_val - unconstrained_fit_lhood_val
tmu_like_stat = tensorlib.astensor(
tensorlib.clip(log_likelihood_ratio, 0.0, max_value=None)
)
if return_fitted_pars:
return tmu_like_stat, (mubhathat, muhatbhat)
return tmu_like_stat
[docs]
def qmu(mu, data, pdf, init_pars, par_bounds, fixed_params, return_fitted_pars=False):
r"""
The test statistic, :math:`q_{\mu}`, for establishing an upper
limit on the strength parameter, :math:`\mu`, as defined in
Equation (14) in :xref:`arXiv:1007.1727`
.. math::
:nowrap:
\begin{equation}
q_{\mu} = \left\{\begin{array}{ll}
-2\ln\lambda\left(\mu\right), &\hat{\mu} \leq \mu,\\
0, & \hat{\mu} > \mu
\end{array}\right.
\end{equation}
where :math:`\lambda\left(\mu\right)` is the profile likelihood ratio as defined in Equation (7)
.. math::
\lambda\left(\mu\right) = \frac{L\left(\mu, \hat{\hat{\boldsymbol{\theta}}}\right)}{L\left(\hat{\mu}, \hat{\boldsymbol{\theta}}\right)}\,.
Example:
>>> 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)
>>> test_mu = 1.0
>>> init_pars = model.config.suggested_init()
>>> par_bounds = model.config.suggested_bounds()
>>> par_bounds[model.config.poi_index] = [-10.0, 10.0]
>>> fixed_params = model.config.suggested_fixed()
>>> pyhf.infer.test_statistics.qmu(
... test_mu, data, model, init_pars, par_bounds, fixed_params
... )
array(3.9549891)
>>> pyhf.infer.test_statistics.qmu(
... test_mu, data, model, init_pars, par_bounds, fixed_params, return_fitted_pars=True
... )
(array(3.9549891), (array([1. , 0.97224597, 0.87553894]), array([-0.06679525, 1.00555369, 0.96930896])))
Args:
mu (Number or Tensor): The signal strength parameter
data (Tensor): The data to be considered
pdf (~pyhf.pdf.Model): The HistFactory statistical model used in the likelihood ratio calculation
init_pars (:obj:`list` of :obj:`float`): The starting values of the model parameters for minimization.
par_bounds (:obj:`list` of :obj:`list`/:obj:`tuple`): The extrema of values the model parameters
are allowed to reach in the fit.
The shape should be ``(n, 2)`` for ``n`` model parameters.
fixed_params (:obj:`tuple` or :obj:`list` of :obj:`bool`): The flag to set a parameter
constant to its starting value during minimization.
return_fitted_pars (:obj:`bool`): Return the best-fit parameter tensors
the fixed-POI and unconstrained fits have converged on
(i.e. :math:`\mu, \hat{\hat{\theta}}` and :math:`\hat{\mu}, \hat{\theta}`)
Returns:
Tuple of a Float and a Tuple of Tensors:
- The calculated test statistic, :math:`q_{\mu}`
- The parameter tensors corresponding to the constrained and unconstrained best fit,
:math:`\mu, \hat{\hat{\theta}}` and :math:`\hat{\mu}, \hat{\theta}`.
Only returned if ``return_fitted_pars`` is ``True``.
"""
if pdf.config.poi_index is None:
raise UnspecifiedPOI(
'No POI is defined. A POI is required for profile likelihood based test statistics.'
)
if par_bounds[pdf.config.poi_index][0] == 0:
log.warning(
'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.\n'
+ 'If you called this from pyhf.infer.mle or pyhf.infer.hypotest, set test_stat="qtilde".'
)
return _qmu_like(
mu,
data,
pdf,
init_pars,
par_bounds,
fixed_params,
return_fitted_pars=return_fitted_pars,
)
[docs]
def qmu_tilde(
mu, data, pdf, init_pars, par_bounds, fixed_params, return_fitted_pars=False
):
r"""
The "alternative" test statistic, :math:`\tilde{q}_{\mu}`, for establishing
an upper limit on the strength parameter, :math:`\mu`, for models with
bounded POI, as defined in Equation (16) in :xref:`arXiv:1007.1727`
.. math::
:nowrap:
\begin{equation}
\tilde{q}_{\mu} = \left\{\begin{array}{ll}
-2\ln\tilde{\lambda}\left(\mu\right), &\hat{\mu} \leq \mu,\\
0, & \hat{\mu} > \mu
\end{array}\right.
\end{equation}
where :math:`\tilde{\lambda}\left(\mu\right)` is the constrained profile likelihood ratio as defined in Equation (10)
.. math::
:nowrap:
\begin{equation}
\tilde{\lambda}\left(\mu\right) = \left\{\begin{array}{ll}
\frac{L\left(\mu, \hat{\hat{\boldsymbol{\theta}}}(\mu)\right)}{L\left(\hat{\mu}, \hat{\hat{\boldsymbol{\theta}}}(0)\right)}, &\hat{\mu} < 0,\\
\frac{L\left(\mu, \hat{\hat{\boldsymbol{\theta}}}(\mu)\right)}{L\left(\hat{\mu}, \hat{\boldsymbol{\theta}}\right)}, &\hat{\mu} \geq 0.
\end{array}\right.
\end{equation}
Example:
>>> 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)
>>> test_mu = 1.0
>>> init_pars = model.config.suggested_init()
>>> par_bounds = model.config.suggested_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)
>>> pyhf.infer.test_statistics.qmu_tilde(
... test_mu, data, model, init_pars, par_bounds, fixed_params, return_fitted_pars=True
... )
(array(3.93824492), (array([1. , 0.97224597, 0.87553894]), array([0. , 1.0030512 , 0.96266961])))
Args:
mu (Number or Tensor): The signal strength parameter
data (:obj:`tensor`): The data to be considered
pdf (~pyhf.pdf.Model): The statistical model adhering to the schema model.json
init_pars (:obj:`list` of :obj:`float`): The starting values of the model parameters for minimization.
par_bounds (:obj:`list` of :obj:`list`/:obj:`tuple`): The extrema of values the model parameters
are allowed to reach in the fit.
The shape should be ``(n, 2)`` for ``n`` model parameters.
fixed_params (:obj:`tuple` or :obj:`list` of :obj:`bool`): The flag to set a parameter
constant to its starting value during minimization.
return_fitted_pars (:obj:`bool`): Return the best-fit parameter tensors
the fixed-POI and unconstrained fits have converged on
(i.e. :math:`\mu, \hat{\hat{\theta}}` and :math:`\hat{\mu}, \hat{\theta}`)
Returns:
Tuple of a Float and a Tuple of Tensors:
- The calculated test statistic, :math:`\tilde{q}_{\mu}`
- The parameter tensors corresponding to the constrained best fit,
:math:`\mu, \hat{\hat{\theta}}`, and the unconstrained best fit,
:math:`\hat{\mu}, \hat{\theta}`.
Only returned if ``return_fitted_pars`` is ``True``.
"""
if pdf.config.poi_index is None:
raise UnspecifiedPOI(
'No POI is defined. A POI is required for profile likelihood based test statistics.'
)
if par_bounds[pdf.config.poi_index][0] != 0:
log.warning(
'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.\n'
+ 'If you called this from pyhf.infer.mle or pyhf.infer.hypotest, set test_stat="q".'
)
return _qmu_like(
mu,
data,
pdf,
init_pars,
par_bounds,
fixed_params,
return_fitted_pars=return_fitted_pars,
)
[docs]
def tmu(mu, data, pdf, init_pars, par_bounds, fixed_params, return_fitted_pars=False):
r"""
The test statistic, :math:`t_{\mu}`, for establishing a two-sided
interval on the strength parameter, :math:`\mu`, as defined in Equation (8)
in :xref:`arXiv:1007.1727`
.. math::
t_{\mu} = -2\ln\lambda\left(\mu\right)
where :math:`\lambda\left(\mu\right)` is the profile likelihood ratio as defined in Equation (7)
.. math::
\lambda\left(\mu\right) = \frac{L\left(\mu, \hat{\hat{\boldsymbol{\theta}}}\right)}{L\left(\hat{\mu}, \hat{\boldsymbol{\theta}}\right)}\,.
Example:
>>> 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)
>>> test_mu = 1.0
>>> init_pars = model.config.suggested_init()
>>> par_bounds = model.config.suggested_bounds()
>>> par_bounds[model.config.poi_index] = [-10.0, 10.0]
>>> fixed_params = model.config.suggested_fixed()
>>> pyhf.infer.test_statistics.tmu(
... test_mu, data, model, init_pars, par_bounds, fixed_params
... )
array(3.9549891)
>>> pyhf.infer.test_statistics.tmu(
... test_mu, data, model, init_pars, par_bounds, fixed_params, return_fitted_pars=True
... )
(array(3.9549891), (array([1. , 0.97224597, 0.87553894]), array([-0.06679525, 1.00555369, 0.96930896])))
Args:
mu (Number or Tensor): The signal strength parameter
data (Tensor): The data to be considered
pdf (~pyhf.pdf.Model): The statistical model adhering to the schema model.json
init_pars (:obj:`list` of :obj:`float`): The starting values of the model parameters for minimization.
par_bounds (:obj:`list` of :obj:`list`/:obj:`tuple`): The extrema of values the model parameters
are allowed to reach in the fit.
The shape should be ``(n, 2)`` for ``n`` model parameters.
fixed_params (:obj:`tuple` or :obj:`list` of :obj:`bool`): The flag to set a parameter
constant to its starting value during minimization.
return_fitted_pars (:obj:`bool`): Return the best-fit parameter tensors
the fixed-POI and unconstrained fits have converged on
(i.e. :math:`\mu, \hat{\hat{\theta}}` and :math:`\hat{\mu}, \hat{\theta}`)
Returns:
Tuple of a Float and a Tuple of Tensors:
- The calculated test statistic, :math:`t_{\mu}`
- The parameter tensors corresponding to the constrained best fit,
:math:`\mu, \hat{\hat{\theta}}`, and the unconstrained best fit,
:math:`\hat{\mu}, \hat{\theta}`.
Only returned if ``return_fitted_pars`` is ``True``.
"""
if pdf.config.poi_index is None:
raise UnspecifiedPOI(
'No POI is defined. A POI is required for profile likelihood based test statistics.'
)
if par_bounds[pdf.config.poi_index][0] == 0:
log.warning(
'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_params,
return_fitted_pars=return_fitted_pars,
)
[docs]
def tmu_tilde(
mu, data, pdf, init_pars, par_bounds, fixed_params, return_fitted_pars=False
):
r"""
The test statistic, :math:`\tilde{t}_{\mu}`, for establishing a two-sided
interval on the strength parameter, :math:`\mu`, for models with
bounded POI, as defined in Equation (11) in :xref:`arXiv:1007.1727`
.. math::
\tilde{t}_{\mu} = -2\ln\tilde{\lambda}\left(\mu\right)
where :math:`\tilde{\lambda}\left(\mu\right)` is the constrained profile likelihood ratio as defined in Equation (10)
.. math::
:nowrap:
\begin{equation}
\tilde{\lambda}\left(\mu\right) = \left\{\begin{array}{ll}
\frac{L\left(\mu, \hat{\hat{\boldsymbol{\theta}}}(\mu)\right)}{L\left(\hat{\mu}, \hat{\hat{\boldsymbol{\theta}}}(0)\right)}, &\hat{\mu} < 0,\\
\frac{L\left(\mu, \hat{\hat{\boldsymbol{\theta}}}(\mu)\right)}{L\left(\hat{\mu}, \hat{\boldsymbol{\theta}}\right)}, &\hat{\mu} \geq 0.
\end{array}\right.
\end{equation}
Example:
>>> 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)
>>> test_mu = 1.0
>>> init_pars = model.config.suggested_init()
>>> par_bounds = model.config.suggested_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)
>>> pyhf.infer.test_statistics.tmu_tilde(
... test_mu, data, model, init_pars, par_bounds, fixed_params, return_fitted_pars=True
... )
(array(3.93824492), (array([1. , 0.97224597, 0.87553894]), array([0. , 1.0030512 , 0.96266961])))
Args:
mu (Number or Tensor): The signal strength parameter
data (:obj:`tensor`): The data to be considered
pdf (~pyhf.pdf.Model): The statistical model adhering to the schema model.json
init_pars (:obj:`list` of :obj:`float`): The starting values of the model parameters for minimization.
par_bounds (:obj:`list` of :obj:`list`/:obj:`tuple`): The extrema of values the model parameters
are allowed to reach in the fit.
The shape should be ``(n, 2)`` for ``n`` model parameters.
fixed_params (:obj:`tuple` or :obj:`list` of :obj:`bool`): The flag to set a parameter
constant to its starting value during minimization.
return_fitted_pars (:obj:`bool`): Return the best-fit parameter tensors
the fixed-POI and unconstrained fits have converged on
(i.e. :math:`\mu, \hat{\hat{\theta}}` and :math:`\hat{\mu}, \hat{\theta}`)
Returns:
Tuple of a Float and a Tuple of Tensors:
- The calculated test statistic, :math:`\tilde{t}_{\mu}`
- The parameter tensors corresponding to the constrained best fit,
:math:`\mu, \hat{\hat{\theta}}`, and the unconstrained best fit,
:math:`\hat{\mu}, \hat{\theta}`.
Only returned if ``return_fitted_pars`` is ``True``.
"""
if pdf.config.poi_index is None:
raise UnspecifiedPOI(
'No POI is defined. A POI is required for profile likelihood based test statistics.'
)
if par_bounds[pdf.config.poi_index][0] != 0:
log.warning(
'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_params,
return_fitted_pars=return_fitted_pars,
)
[docs]
def q0(mu, data, pdf, init_pars, par_bounds, fixed_params, return_fitted_pars=False):
r"""
The test statistic, :math:`q_{0}`, for discovery of a positive signal
as defined in Equation (12) in :xref:`arXiv:1007.1727`, for :math:`\mu=0`.
.. math::
:nowrap:
\begin{equation}
q_{0} = \left\{\begin{array}{ll}
-2\ln\lambda\left(0\right), &\hat{\mu} \ge 0,\\
0, & \hat{\mu} < 0,
\end{array}\right.
\end{equation}
Example:
>>> 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 = [60, 65]
>>> data = pyhf.tensorlib.astensor(observations + model.config.auxdata)
>>> test_mu = 0.0
>>> init_pars = model.config.suggested_init()
>>> par_bounds = model.config.suggested_bounds()
>>> fixed_params = model.config.suggested_fixed()
>>> pyhf.infer.test_statistics.q0(test_mu, data, model, init_pars, par_bounds, fixed_params)
array(2.98339447)
>>> pyhf.infer.test_statistics.q0(
... test_mu, data, model, init_pars, par_bounds, fixed_params, return_fitted_pars=True
... )
(array(2.98339447), (array([0. , 1.03050845, 1.12128752]), array([0.95260667, 0.99635345, 1.02140172])))
Args:
mu (Number or Tensor): The signal strength parameter (must be set to zero)
data (Tensor): The data to be considered
pdf (~pyhf.pdf.Model): The HistFactory statistical model used in the likelihood ratio calculation
init_pars (:obj:`list` of :obj:`float`): The starting values of the model parameters for minimization.
par_bounds (:obj:`list` of :obj:`list`/:obj:`tuple`): The extrema of values the model parameters
are allowed to reach in the fit.
The shape should be ``(n, 2)`` for ``n`` model parameters.
fixed_params (:obj:`tuple` or :obj:`list` of :obj:`bool`): The flag to set a parameter
constant to its starting value during minimization.
return_fitted_pars (:obj:`bool`): Return the best-fit parameter tensors
the fixed-POI and unconstrained fits have converged on
(i.e. :math:`\mu, \hat{\hat{\theta}}` and :math:`\hat{\mu}, \hat{\theta}`)
Returns:
Tuple of a Float and a Tuple of Tensors:
- The calculated test statistic, :math:`q_{0}`
- The parameter tensors corresponding to the constrained best fit,
:math:`\mu, \hat{\hat{\theta}}`, and the unconstrained best fit,
:math:`\hat{\mu}, \hat{\theta}`.
Only returned if ``return_fitted_pars`` is ``True``.
"""
if pdf.config.poi_index is None:
raise UnspecifiedPOI(
'No POI is defined. A POI is required for profile likelihood based test statistics.'
)
if mu != 0.0:
log.warning(
'q0 test statistic only used for fit configuration with POI set to zero. Setting mu=0.'
)
mu = 0.0
tensorlib, optimizer = get_backend()
tmu_like_stat, (mubhathat, muhatbhat) = _tmu_like(
mu, data, pdf, init_pars, par_bounds, fixed_params, return_fitted_pars=True
)
q0_stat = tensorlib.where(
muhatbhat[pdf.config.poi_index] < 0, tensorlib.astensor(0.0), tmu_like_stat
)
if return_fitted_pars:
return q0_stat, (mubhathat, muhatbhat)
return q0_stat