Merge pull request #42 from SpeysideHEP/chi2-bugfix

Bugfix in exclusion limit calculator

.github/workflows/tests.yml

@@ -18,7 +18,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"]
+        python-version: ["3.8", "3.9", "3.10", "3.11", "3.12.6"]
         os: [ubuntu-latest, macos-latest]
 
     steps:

.zenodo.json

@@ -6,7 +6,7 @@
     "upload_type": "software",
     "creators": [
         {
-            "affiliation": "Thomas Jefferson National Accelerator Facility",
+            "affiliation": "Stony Brook University",
             "name": "Araz, Jack Y.",
             "orcid": "0000-0001-8721-8042"
         }

docs/exclusion.md

@@ -0,0 +1,76 @@
+---
+myst:
+  html_meta:
+    "property=og:title": "Exclusion limits"
+    "property=og:description": "How does the exclusion limits work in spey"
+    "property=og:image": "https://spey.readthedocs.io/en/main/_static/spey-logo.png"
+    "property=og:url": "https://spey.readthedocs.io/en/main/exclusion.html"
+jupytext:
+  formats: ipynb,md:myst
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.12
+    jupytext_version: 1.8.2
+kernelspec:
+  display_name: Python 3
+  language: python
+  name: python3
+---
+
+# Exclusion limits
+
+Any Spey statistical model can compute the exclusion confidence level using three different calculators. The availability of these options depends on the functions used in the likelihood construction (see {ref}`this section <sec_new_plugin>` for details).
+
+The exclusion limits can be calculated for two different test statistic measure i.e. test statistic
+
+$$
+q_\mu = \begin{cases}
+        0 & {\rm if\ } \hat{\mu}>\mu\ , \\
+        -2\log\frac{\mathcal{L}(\mu, \theta_\mu)}{\mathcal{L}(\hat{\mu}, \hat{\theta})} & {\rm otherwise}
+    \end{cases} \quad ,
+$$
+
+or alternate test statistic
+
+$$
+\tilde{q}_\mu = \begin{cases}
+        0 & {\rm if\ } \hat{\mu}>\mu\ , \\
+        -2\log\frac{\mathcal{L}(\mu, \theta_\mu)}{\mathcal{L}(\hat{\mu}, \hat{\theta})} & {\rm if\ } 0\leq\hat{\mu}\leq\mu\ ,\\
+        -2\log\frac{\mathcal{L}(\mu, \theta_\mu)}{\mathcal{L}(0, \theta_0)} & {\rm otherwise}
+    \end{cases}\quad .
+$$
+
+The main distinction being, alternate test statistic assumes that the signal can not be negative, for instance in case of negative interference effects in EFT one should use test statistic. By default {func}`~spey.StatisticalModel.exclusion_confidence_level` function assumes alternate test statistic being use and this can be changed with ``allow_negative_signal`` boolean argument.
+
+````{margin}
+```{note}
+To check which calculators are available for a given model, users can call the {obj}`~spey.StatisticalModel.available_calculators` function.
+```
+````
+
+The {func}`~spey.StatisticalModel.exclusion_confidence_level` function allows users to specify the calculator through the ``calculator`` keyword, offering three options:
+
+- **"asymptotic"**: Uses asymptotic formulae to compute p-values (see ref. {cite}`Cowan:2010js` for details). This method requires the likelihood to support sampling capabilities, as it compares the signal-like test statistic to the Asimov test statistic using the formula $t_\mu = \sqrt{\tilde{q}_\mu} - \sqrt{\tilde{q}_{\mu,A}}$. This is the most commonly used method in LHC analyses.
+- **"toy"**: Relies on the likelihood's sampling functionality. It calculates p-values by generating samples from both the signal-plus-background and background-only distributions.
+- **"chi_square"**: Compares the signal hypothesis ($\tilde{q}_{\mu=1}$) to the null hypothesis, with the test statistic $t_\mu = \sqrt{\chi^2(\mu=0)} - \sqrt{\tilde{q}_{\mu=1}}$. This approach was widely used during the Tevatron era.
+
+The `expected` keyword allows users to select between computing observed or expected exclusion limits. It also supports the calculation of prefit expected exclusion limits, which can be enabled by setting `expected=spey.ExpectationType.apriori`. This option ignores experimental data and computes the expected exclusion limit based solely on the simulated Standard Model (SM) background yields. On the other hand, {obj}`~spey.ExpectationType.observed` (for observed limits) and {obj}`~spey.ExpectationType.aposteriori` (for post-fit expected limits) compute the exclusion confidence limits after fitting the model.
+
+In both expected cases, the exclusion limits are returned with $\pm1\sigma$ and $\pm2\sigma$ variations around the background model, resulting in five values: $[-2\sigma, -1\sigma, 0, 1\sigma, 2\sigma]$. However, the observed exclusion limit returns a single value. The ``"chi_square"`` calculator is an exception—it only provides one value for both observed and expected limits.
+
+The `allow_negative_signal` keyword controls which test statistic is used and restricts the values that $\mu$ (the signal strength) can take when computing the maximum likelihood. When `allow_negative_signal=True`, the $q_\mu$ test statistic is applied; otherwise, the $\tilde{q}_\mu$ statistic is used (for further details, see {cite}`Cowan:2010js, Araz:2023bwx`).
+
+For complex statistical models, optimizing the likelihood can be challenging and depends on the choice of optimizer. Spey uses SciPy for optimization and fitting tasks. Any additional keyword arguments not explicitly covered in the {func}`~spey.StatisticalModel.exclusion_confidence_level` function description are passed directly to the optimizer, allowing users to customize its behavior through the interface.
+
+Below we compare the exclusion limits computed with each approach. This comparisson uses normal distribution for the likelihood ([`default_pdf.normal`](#normal)) background yields are set to $n_b$, uncertainties are shown with $\sigma$ and observations are given with $n$.
+
+```{figure} ./figs/comparisson_observed.png
+---
+width: 100%
+figclass: caption
+alt: exclusion limit calculator comparisson
+name: fig2
+---
+exclusion limit calculator comparisson for observed p-values
+```

docs/releases/changelog-v0.1.md

@@ -79,6 +79,9 @@ Specific upgrades for the latest release can be found [here](https://github.com/
 * Compatibility issues have been fixed in Python 3.12
   ([#39](https://github.com/SpeysideHEP/spey/pull/39))
 
+* p-value computation in $\chi^2$ and toy calculators have been fixed.
+  ([#42](https://github.com/SpeysideHEP/spey/pull/42))
+
 ## Contributors
 
 This release contains contributions from (in alphabetical order):

setup.py

@@ -42,13 +42,13 @@
     },
     download_url=f"https://github.com/SpeysideHEP/spey/archive/refs/tags/v{version}.tar.gz",
     author="Jack Y. Araz",
-    author_email=("[email protected]"),
+    author_email=("[email protected]"),
     license="MIT",
     package_dir={"": "src"},
     packages=find_packages(where="src"),
     entry_points={"spey.backend.plugins": backend_plugins},
     install_requires=requirements,
-    python_requires=">=3.8, <=3.12.4",
+    python_requires=">=3.8, <=3.12.6",
     classifiers=[
         "Intended Audience :: Science/Research",
         "License :: OSI Approved :: MIT License",

src/spey/__init__.py

@@ -53,6 +53,14 @@ def set_log_level(level: Literal[0, 1, 2, 3]) -> None:
     """
     Set log level for spey
 
+    Log level can also be set through terminal by the following command
+
+    .. code::
+
+        export SPEY_LOGLEVEL=3
+
+    value corresponds to the levels shown below.
+
     Args:
         level (``int``): log level
 
@@ -332,5 +340,8 @@ def cite() -> List[Text]:
         return ""
 
 
+if int(os.environ.get("SPEY_LOGLEVEL", -1)) >= 0:
+    set_log_level(int(os.environ.get("SPEY_LOGLEVEL")))
+
 if os.environ.get("SPEY_CHECKUPDATE", "ON").upper() != "OFF":
     check_updates()

src/spey/_version.py

@@ -1,3 +1,3 @@
 """Version number (major.minor.patch[-label])"""
 
-__version__ = "0.1.10-beta"
+__version__ = "0.1.10"