diff --git a/.zenodo.json b/.zenodo.json index 5d2c7d4..cd0af80 100644 --- a/.zenodo.json +++ b/.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" }, { diff --git a/docs/figs/sig_unc_chi2.png b/docs/figs/sig_unc_chi2.png new file mode 100644 index 0000000..227a7e2 Binary files /dev/null and b/docs/figs/sig_unc_chi2.png differ diff --git a/docs/releases/changelog-v0.1.md b/docs/releases/changelog-v0.1.md index a46edc2..26ae0e3 100644 --- a/docs/releases/changelog-v0.1.md +++ b/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): diff --git a/docs/tutorials/intro_spey.ipynb b/docs/tutorials/intro_spey.ipynb index bf47651..21016fa 100644 --- a/docs/tutorials/intro_spey.ipynb +++ b/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, diff --git a/docs/tutorials/signal_unc.md b/docs/tutorials/signal_unc.md index 312e5ca..6cc3bea 100644 --- a/docs/tutorials/signal_unc.md +++ b/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 -``` \ No newline at end of file +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. +``` diff --git a/src/spey/_version.py b/src/spey/_version.py index 9a13fa2..e83e13c 100644 --- a/src/spey/_version.py +++ b/src/spey/_version.py @@ -1,3 +1,3 @@ """Version number (major.minor.patch[-label])""" -__version__ = "0.1.7" +__version__ = "0.1.8" diff --git a/src/spey/backends/default_pdf/uncertainty_synthesizer.py b/src/spey/backends/default_pdf/uncertainty_synthesizer.py index 8b54a9a..99520a9 100644 --- a/src/spey/backends/default_pdf/uncertainty_synthesizer.py +++ b/src/spey/backends/default_pdf/uncertainty_synthesizer.py @@ -12,13 +12,24 @@ 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: @@ -26,7 +37,7 @@ def constraint_from_corr( { "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}, } ] diff --git a/src/spey/backends/distributions.py b/src/spey/backends/distributions.py index 8588b5f..265eeee 100644 --- a/src/spey/backends/distributions.py +++ b/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", {}) diff --git a/src/spey/system/webutils.py b/src/spey/system/webutils.py index b67c5ac..51b77b2 100644 --- a/src/spey/system/webutils.py +++ b/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):