Merge pull request #38 from SpeysideHEP/sigunc_bugfix

Bugfix in signal uncertainties

.zenodo.json

@@ -1,8 +1,8 @@
 {
     "description": "smooth inference for reinterpretation studies",
     "license": "MIT",
-    "title": "SpeysideHEP/spey: v0.1.7",
-    "version": "v0.1.7",
+    "title": "SpeysideHEP/spey: v0.1.8",
+    "version": "v0.1.8",
     "upload_type": "software",
     "creators": [
         {
@@ -29,7 +29,7 @@
         },
         {
             "scheme": "url",
-            "identifier": "https://github.com/SpeysideHEP/spey/tree/v0.1.7",
+            "identifier": "https://github.com/SpeysideHEP/spey/tree/v0.1.8",
             "relation": "isSupplementTo"
         },
         {

docs/releases/changelog-v0.1.md

@@ -66,6 +66,9 @@ Specific upgrades for the latest release can be found [here](https://github.com/
 * Bug fix in signal uncertainty synthesizer
   ([#34](https://github.com/SpeysideHEP/spey/pull/34))
 
+* Signal uncertainties were causing a narrower $\chi^2$ distribution due to the weight of the constraint term. 
+  ([#38](https://github.com/SpeysideHEP/spey/pull/38))
+
 ## Contributors
 
 This release contains contributions from (in alphabetical order):

docs/tutorials/intro_spey.ipynb

@@ -257,11 +257,10 @@
    "id": "99dcc3b1",
    "metadata": {},
    "source": [
-    "````{margin}\n",
     "```{attention}\n",
+    ":class: dropdown\n",
     "Notice that we started POI scan from -0.3, this is because this function is not defined for $\\mu<-0.3$ otherwise the Poisson will get negative values. This value can be checked via ```corr_background_model.backend.config().minimum_poi``` attribute.\n",
-    "```\n",
-    "````"
+    "```"
    ]
   },
   {
@@ -542,8 +541,8 @@
    "id": "e0f22242",
    "metadata": {},
    "source": [
-    "````{margin}\n",
     "```{admonition} **Question:** \n",
+    ":class: dropdown\n",
     "\n",
     "***Experimental collaboration provided asymmetric uncertainties but not third moments, can spey compute third moments?***\n",
     "\n",
@@ -554,8 +553,7 @@
     "$$\n",
     "\n",
     "But this is an assumption, if collaboration provides exact third moments, please always use those. \n",
-    "```\n",
-    "````"
+    "```"
    ]
   },
   {
@@ -767,7 +765,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.16"
+   "version": "3.9.18"
   }
  },
  "nbformat": 4,

docs/tutorials/signal_unc.md

@@ -30,35 +30,71 @@ Note that theoretical uncertainties have different interpretations, we can inter
 ```{code-cell} ipython3
 :tags: [hide-cell]
 import spey
+import numpy as np
+import matplotlib.pyplot as plt
 ```
 
 We can add uncorrelated signal uncertainties just like in uncorrelated background case
 
 $$
-    \mathcal{L}(\mu, \theta) = \prod_{i\in{\rm bins}}{\rm Poiss}(n^i|\mu n_s^i + n_b^i + \theta^i\sigma_b^i + \theta_i^{(s)}\sigma_s^i) \cdot \prod_{j\in{\rm nui}}\mathcal{N}(\theta^j|0, 1) \cdot \prod_{j\in{\rm nui}}\mathcal{N}(\theta_{(s)}^j|0, 1)\ ,
+    \mathcal{L}(\mu, \theta) = \prod_{i\in{\rm bins}}{\rm Poiss}(n_i|\mu n^{(s)}_i + \theta_i^{(s)}\sigma^{(s)}_i + n^{(b)}_i + \theta_i^{(b)}\sigma^{(b)}_i) \cdot \prod_{j\in{\rm nui}}\mathcal{N}(\theta_j^{(b)}|0, 1) \cdot \prod_{j\in{\rm nui}}\mathcal{N}(\theta^{(s)}_j|0, 1)\ ,
 $$
 
+where $(s)$ superscript indicates signal and $(b)$ indicates background.
+
 ```{code-cell} ipython3
 pdf_wrapper = spey.get_backend("default_pdf.uncorrelated_background")
-statistical_model = pdf_wrapper(
+statistical_model_sigunc = pdf_wrapper(
     signal_yields=[12.0, 15.0],
-    background_yields=[50.0,48.0],
+    background_yields=[50.0, 48.0],
     data=[36, 33],
-    absolute_uncertainties=[12.0,16.0],
+    absolute_uncertainties=[12.0, 16.0],
     signal_uncertainty_configuration={"absolute_uncertainties": [3.0, 4.0]},
-    analysis="example",
-    xsection=0.123,
 )
 ```
 
 Similarly, we can construct signal uncertainties using ``"absolute_uncertainty_envelops"`` keyword which accepts upper and lower uncertainties as ``[(upper, lower)]``. We can also add a correlation matrix with ``"correlation_matrix"`` keyword and third moments with ``"third_moments"`` keyword. Notice that these are completely independent of background. Now we can simply compute the limits as follows
 
 ```{code-cell} ipython3
-print(f"1 - CLs: {statistical_model.exclusion_confidence_level()[0]:.5f}")
-print(f"POI upper limit: {statistical_model.poi_upper_limit():.5f}")
+print(f"1 - CLs: {statistical_model_sigunc.exclusion_confidence_level()[0]:.5f}")
+print(f"POI upper limit: {statistical_model_sigunc.poi_upper_limit():.5f}")
 ```
 
 ```python
-1 - CLs: 0.99563
-POI upper limit: 0.51504
-```
+1 - CLs: 0.96607
+POI upper limit: 0.88808
+```
+
+Let us also check the $\chi^2$ distribution with respect to POI which we expect the distribution should get wider with signal uncertainties. For this comparison we first need to define the model without signal uncertainties:
+
+```{code-cell} ipython3
+statistical_model = pdf_wrapper(
+    signal_yields=[12.0, 15.0],
+    background_yields=[50.0, 48.0],
+    data=[36, 33],
+    absolute_uncertainties=[12.0, 16.0],
+)
+```
+
+Using ``statistical_model`` and ``statistical_model_sigunc`` we can compute the $\chi^2$ distribution
+
+```{code-cell} ipython3
+:tags: [hide-cell]
+poi = np.linspace(-3,2,20)
+plt.plot(poi, [statistical_model.chi2(p, allow_negative_signal=True) for p in poi], color="b", label="no signal uncertainties")
+plt.plot(poi, [statistical_model_sigunc.chi2(p, allow_negative_signal=True) for p in poi], color="r", label="with signal uncertainties")
+plt.legend()
+plt.xlabel("$\mu$")
+plt.ylabel("$\chi^2(\mu)$")
+plt.show()
+```
+
+```{figure} ../figs/sig_unc_chi2.png
+---
+width: 60%
+figclass: caption
+alt: chi-square distribution
+name: fig1
+---
+$\chi^2(\mu)$ distribution comparisson for statistical model with and without signal uncertainties.
+```

src/spey/_version.py

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

src/spey/backends/default_pdf/uncertainty_synthesizer.py

@@ -12,21 +12,32 @@
 def constraint_from_corr(
     correlation_matrix: List[List[float]], size: int, domain: slice
 ) -> List[Dict[Text, Any]]:
+    """
+    Derive constraints from inputs
+
+    Args:
+        correlation_matrix (``List[List[float]]``): correlation matrix
+        size (``int``): size of the signal vector
+        domain (``slice``): domain of the nuisances
+
+    Returns:
+        ``List[Dict[Text, Any]]``:
+    """
     if correlation_matrix is not None:
         corr = np.array(correlation_matrix)
         constraint_term = [
             {
                 "distribution_type": "multivariatenormal",
                 "args": [np.zeros(size), corr],
-                "kwargs": {"domain": domain, "weight": lambda pars: pars[0]},
+                "kwargs": {"domain": domain},
             }
         ]
     else:
         constraint_term = [
             {
                 "distribution_type": "normal",
                 "args": [np.zeros(size), np.ones(size)],
-                "kwargs": {"domain": domain, "weight": lambda pars: pars[0]},
+                "kwargs": {"domain": domain},
             }
         ]
 
@@ -67,7 +78,7 @@ def signal_uncertainty_synthesizer(
         absolute_uncertainties = np.array(absolute_uncertainties)
 
         def lam_signal(pars: np.ndarray) -> np.ndarray:
-            return pars[0] * absolute_uncertainties * pars[domain]
+            return absolute_uncertainties * pars[domain]
 
         constraint_term = constraint_from_corr(
             correlation_matrix, len(absolute_uncertainties), domain
@@ -86,12 +97,12 @@ def effective_sigma(pars: np.ndarray) -> np.ndarray:
             """Compute effective sigma"""
             return np.sqrt(
                 sigma_plus * sigma_minus
-                + (sigma_plus - sigma_minus) * (pars - signal_yields)
+                + (sigma_plus - sigma_minus) * (pars[domain] - signal_yields)
             )
 
         def lam_signal(pars: np.ndarray) -> np.ndarray:
             """Compute lambda for Main model"""
-            return pars[0] * effective_sigma(pars[domain]) * pars[domain]
+            return effective_sigma(pars) * pars[domain]
 
         constraint_term = constraint_from_corr(
             correlation_matrix, len(sigma_plus), domain
@@ -120,14 +131,13 @@ def lam_signal(pars: np.ndarray) -> np.ndarray:
                 expectation value of the poisson distribution with respect to
                 nuisance parameters.
             """
-            nI = A + B * pars[domain] + C * np.square(pars[domain])
-            return pars[0] * nI
+            return A + B * pars[domain] + C * np.square(pars[domain])
 
         constraint_term = [
             {
                 "distribution_type": "multivariatenormal",
                 "args": [np.zeros(len(signal_yields)), corr],
-                "kwargs": {"domain": domain, "weight": lambda pars: pars[0]},
+                "kwargs": {"domain": domain},
             }
         ]
 

src/spey/backends/distributions.py

@@ -1,5 +1,6 @@
 """Autograd based differentiable distribution classes"""
 
+import logging
 import warnings
 from typing import Any, Callable, Dict, List, Text, Union
 
@@ -8,7 +9,9 @@
 from autograd.scipy.stats.poisson import logpmf
 from scipy.stats import multivariate_normal, norm, poisson
 
-# pylint: disable=E1101
+# pylint: disable=E1101, W1203
+
+log = logging.getLogger("Spey")
 
 __all__ = ["Poisson", "Normal", "MultivariateNormal", "MainModel", "ConstraintModel"]
 
@@ -239,11 +242,13 @@ def __init__(self, pdf_descriptions: List[Dict[Text, Any]]):
         self._pdfs = []
         distributions = {"normal": Normal, "multivariatenormal": MultivariateNormal}
 
+        log.debug("Adding constraint terms:")
         for desc in pdf_descriptions:
             assert desc["distribution_type"].lower() in [
                 "normal",
                 "multivariatenormal",
             ], f"Unknown distribution type: {desc['distribution_type']}"
+            log.debug(f"{desc}")
             self._pdfs.append(
                 distributions[desc["distribution_type"]](
                     *desc.get("args", []), **desc.get("kwargs", {})

src/spey/system/webutils.py

@@ -2,9 +2,10 @@
 from typing import Literal, Text
 
 import requests
-from pkg_resources import get_distribution
 from semantic_version import Version
 
+from spey._version import __version__
+
 log = logging.getLogger("Spey")
 
 __all__ = ["get_bibtex", "check_updates", "ConnectionError"]
@@ -101,7 +102,7 @@ def check_updates() -> None:
         response.encoding = "utf-8"
         pypi_info = response.json()
         pypi_version = pypi_info.get("info", {}).get("version", False)
-        version = get_distribution("spey").version
+        version = __version__
         if pypi_version:
             log.debug(f"Curernt version {version}, latest version {pypi_version}.")
             if Version(version) < Version(pypi_version):