#!/usr/bin/env python3
"""
.. module:: basicStats
:synopsis: a module thats collects various statistical codes, that \
are shared among other submodules.
.. moduleauthor:: Wolfgang Waltenberger <wolfgang.waltenberger@gmail.com>
"""
from scipy import stats
from smodels.base.smodelsLogging import logger
import numpy as np
from smodels.statistics.exceptions import SModelSStatisticsError as SModelSError
from typing import Text
__all__ = [ "CLsfromNLL", "determineBrentBracket", "chi2FromLmax" ]
[docs]def CLsfromNLL(
nllA: float, nll0A: float, nll: float, nll0: float,
return_type: Text = "CLs-alpha" ) -> float:
"""
compute the CLs - alpha from the NLLs
TODO: following needs explanation
:param nllA:
:param nll0A:
:param nll:
:param nll0:
:param return_type: (Text) can be "CLs-alpha", "1-CLs", "CLs" \
CLs-alpha: returns CLs - 0.05 \
1-CLs: returns 1-CLs value \
CLs: returns CLs value
:return: Cls-type value, see above
"""
assert return_type in ["CLs-alpha", "1-CLs", "CLs"], f"Unknown return type: {return_type}."
qmu = 0.0 if 2 * (nll - nll0) < 0.0 else 2 * (nll - nll0)
sqmu = np.sqrt(qmu)
qA = 2 * (nllA - nll0A)
if qA < 0.0:
qA = 0.0
sqA = np.sqrt(qA)
if qA >= qmu:
CLsb = 1.0 - stats.multivariate_normal.cdf(sqmu)
CLb = stats.multivariate_normal.cdf(sqA - sqmu)
else:
CLsb = 1.0 if qA == 0.0 else 1.0 - stats.multivariate_normal.cdf((qmu + qA) / (2 * sqA))
CLb = 1.0 if qA == 0.0 else 1.0 - stats.multivariate_normal.cdf((qmu - qA) / (2 * sqA))
CLs = CLsb / CLb if CLb > 0 else 0.0
if return_type == "1-CLs":
return 1.0 - CLs
elif return_type == "CLs":
return CLs
return CLs - 0.05
[docs]def determineBrentBracket(mu_hat, sigma_mu, rootfinder,
allowNegative = True ):
"""find a, b for brent bracketing
:param mu_hat: mu that maximizes likelihood
:param sigm_mu: error on mu_hat (not too reliable)
:param rootfinder: function that finds the root (usually root_func)
:param allowNegative: if False, then do not allow a or b to become negative
:returns: the interval a,b
"""
sigma_mu = max(sigma_mu, 0.5) # there is a minimum on sigma_mu
sigma_mu = min(sigma_mu, 100.) # there is a maximum on sigma_mu
# the root should be roughly at mu_hat + 2*sigma_mu
a = mu_hat + 1.5 * sigma_mu
ra = rootfinder(a)
ntrials = 20
i = 0
foundExtra = False
while ra < 0.0:
# if this is negative, we move it to the left
i += 1
a -= (i**2.0) * sigma_mu
ra = rootfinder(a)
if i > ntrials or a < -10000.0 or ra is None or ( a < 0 and not allowNegative ):
avalues = [0.0, 1.0, -1.0, 3.0, -3.0, 10.0, -10.0, 0.1, -0.1, 0.01, -0.01, .001, -.001, 100., -100., .0001, -.0001 ]
if not allowNegative:
avalues = [0.0, 1.0, 3.0, 10.0, 0.1, 0.01, .001, 100., .0001 ]
for a in avalues:
ra = rootfinder(a)
if ra is None: # if cls computation failed, try with next a value
continue
if ra > 0.0:
foundExtra = True
break
if not foundExtra:
logger.error(
f"cannot find an a that is left of the root. last attempt, a={a:.2f}, root = {ra:.2f}."
)
logger.error(f"mu_hat={mu_hat:.2f}, sigma_mu={sigma_mu:.2f}")
raise SModelSError(
f"cannot find an a that is left of the root. last attempt, a={a:.2f}, root = {ra:.2f}."
)
i = 0
foundExtra = False
b = mu_hat + 2.5 * sigma_mu
rb = rootfinder(b)
while rb > 0.0:
# if this is positive, we move it to the right
i += 1
b += (i**2.0) * sigma_mu
rb = rootfinder(b)
if rb is None: # if cls computation failed, try with a bigger b value
continue
closestr, closest = float("inf"), None
if i > ntrials or ( b < 0 and not allowNegative ):
bvalues = [1.0, 0.0, 3.0, -1.0, 10.0, -3.0, 0.1, -0.1, -10.0, 100.0, -100.0, 1000.0, .01, -.01, .001, -.001, 10000.0, 100000.0, 1000000.0 ]
if not allowNegative:
bvalues = [1.0, 0.0, 3.0, 10.0, 0.1, 100.0, 1000.0, .01, .001, 10000.0, 100000.0, 1000000.0 ]
for b in bvalues:
rb = rootfinder(b)
if rb is None: # if cls computation failed, try with next b value
continue
if rb < 0.0:
foundExtra = True
break
if rb < closestr:
closestr = rb
closest = b
if not foundExtra:
logger.error(f"cannot find a b that is right of the root (i.e. rootfinder(b) < 0).")
logger.error(f"closest to zero rootfinder({closest})={closestr}")
logger.error(f"mu_hat was at {mu_hat:.2f} sigma_mu at {sigma_mu:.2f}")
raise SModelSError()
return a, b
def deltaChi2FromLlhd(likelihood):
"""compute the delta chi2 value from a likelihood (convenience function)"""
if likelihood == 0.0:
return 1e10 # a very high number is good
elif likelihood is None:
return None
return -2.0 * np.log(likelihood)
[docs]def chi2FromLmax(llhd, lmax):
"""compute the chi2 from likelihood and lmax"""
if llhd is None or lmax is None:
return None
chi2 = 0.0
if llhd > 1e-200:
from math import log
chi2 = 2 * log(lmax / llhd)
if chi2 < 0.0 and llhd < 1e-100:
# numerical inaccuracy
chi2 = 0.0
return chi2