From 5b27c8b15cae3e6e92c6bfa76218fc8af769b7bb Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Fri, 20 Dec 2024 10:27:39 -0500
Subject: [PATCH 01/32] update use

---
 docs/bib/cited.bib | 34 +++++++++++++++++++++++++---------
 1 file changed, 25 insertions(+), 9 deletions(-)

diff --git a/docs/bib/cited.bib b/docs/bib/cited.bib
index af484d5..a1fbc30 100644
--- a/docs/bib/cited.bib
+++ b/docs/bib/cited.bib
@@ -1,3 +1,15 @@
+@article{Ghosh:2024puc,
+    author        = {Ghosh, Kirtiman and Huitu, Katri and Sahu, Rameswar},
+    title         = {{Revisiting Universal Extra-Dimension Model with Gravity Mediated Decays}},
+    eprint        = {2412.09344},
+    archiveprefix = {arXiv},
+    primaryclass  = {hep-ph},
+    month         = {12},
+    year          = {2024},
+    journal       = {}
+}
+
+
 @article{Goodsell:2024aig,
     author        = {Goodsell, Mark D.},
     title         = {{HackAnalysis 2: A powerful and hackable recasting tool}},
@@ -15,25 +27,29 @@ @article{Agin:2024yfs
     eprint        = {2404.12423},
     archiveprefix = {arXiv},
     primaryclass  = {hep-ph},
-    month         = {4},
-    year          = {2024},
-    journal       = {}
+    doi           = {10.1140/epjc/s10052-024-13594-9},
+    journal       = {Eur. Phys. J. C},
+    volume        = {84},
+    number        = {11},
+    pages         = {1218},
+    year          = {2024}
 }
 
 
 @article{Feike:2024zfz,
-    archiveprefix = {arXiv},
     author        = {Feike, Alexander and Fiaschi, Juri and Fuks, Benjamin and Klasen, Michael and Puck Neuwirth, Alexander},
+    title         = {{Combination and reinterpretation of LHC SUSY searches}},
     eprint        = {2403.11715},
-    month         = {3},
+    archiveprefix = {arXiv},
     primaryclass  = {hep-ph},
     reportnumber  = {MS-TP-23-49},
-    title         = {{Combination and Reinterpretation of LHC SUSY Searches}},
-    year          = {2024},
-    journal       = {}
+    doi           = {10.1007/JHEP07(2024)122},
+    journal       = {JHEP},
+    volume        = {07},
+    pages         = {122},
+    year          = {2024}
 }
 
-
 @article{Elmer:2023wtr,
     author        = {Elmer, Nina and Madigan, Maeve and Plehn, Tilman and Schmal, Nikita},
     title         = {{Staying on Top of SMEFT-Likelihood Analyses}},

From 97a851a3d5cfac2814d20014100ef01932e67cfc Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Wed, 29 Jan 2025 10:14:45 -0500
Subject: [PATCH 02/32] allow for negative poi in UL calculation

---
 src/spey/base/hypotest_base.py              | 80 +++++++++++----------
 src/spey/hypothesis_testing/upper_limits.py |  3 +-
 2 files changed, 43 insertions(+), 40 deletions(-)

diff --git a/src/spey/base/hypotest_base.py b/src/spey/base/hypotest_base.py
index 163f450..5fbac08 100644
--- a/src/spey/base/hypotest_base.py
+++ b/src/spey/base/hypotest_base.py
@@ -840,42 +840,38 @@ def poi_upper_limit(
         low_init: Optional[float] = 1.0,
         hig_init: Optional[float] = 1.0,
         expected_pvalue: Literal["nominal", "1sigma", "2sigma"] = "nominal",
+        allow_negative_signal: bool = False,
         maxiter: int = 10000,
         optimiser_arguments: Optional[Dict[str, Any]] = None,
     ) -> Union[float, List[float]]:
         r"""
-        Compute the upper limit for the parameter of interest i.e. :math:`\mu`.
-
-        .. note::
-
-            This function uses ``"qtilde"`` test statistic which means signal values are always
-            assumed to be positive i.e. :math:`\hat\mu>0`.
+        Compute the upper limit for the parameter of interest (POI), denoted as :math:`\mu`.
 
         Args:
-            expected (~spey.ExpectationType): Sets which values the fitting algorithm should focus and
-              p-values to be computed.
-
-              * :obj:`~spey.ExpectationType.observed`: Computes the p-values with via post-fit
-                prescriotion which means that the experimental data will be assumed to be the truth
-                (default).
-              * :obj:`~spey.ExpectationType.aposteriori`: Computes the expected p-values with via
-                post-fit prescriotion which means that the experimental data will be assumed to be
-                the truth.
-              * :obj:`~spey.ExpectationType.apriori`: Computes the expected p-values with via pre-fit
-                prescription which means that the SM will be assumed to be the truth.
-
-            confidence_level (``float``, default ``0.95``): Determines the confidence level of the upper
-              limit i.e. the value of :math:`1-CL_s`. It needs to be between ``[0,1]``.
-            low_init (``Optional[float]``, default ``1.0``): Lower limit for the search algorithm to start
-              If ``None`` it the lower limit will be determined by :math:`\hat\mu + 1.5\sigma_{\hat\mu}`.
+            expected (:obj:`~spey.ExpectationType`, default :obj:`~spey.ExpectationType.observed`):
+              Specifies the type of expectation for the fitting algorithm and p-value computation.
+
+              * :obj:`~spey.ExpectationType.observed`: Computes p-values using post-fit prescription,
+                assuming experimental data as the truth (default).
+              * :obj:`~spey.ExpectationType.aposteriori`: Computes expected p-values using post-fit
+                prescription, assuming experimental data as the truth.
+              * :obj:`~spey.ExpectationType.apriori`: Computes expected p-values using pre-fit
+                prescription, assuming the Standard Model (SM) as the truth.
+
+            confidence_level (``float``, default ``0.95``): Confidence level for the upper limit,
+                representing :math:`1 - CL_s`. Must be between 0 and 1. Default is 0.95.
+            low_init (``Optional[float]``, default ``1.0``): Initial lower limit for the search
+              algorithm. If `None`, it is determined by :math:`\hat\mu + 1.5\sigma_{\hat\mu}`.
+              Default is 1.0.
 
               .. note::
 
                 :math:`\sigma_{\hat\mu}` is determined via
                 :func:`~spey.base.hypotest_base.HypothesisTestingBase.sigma_mu` function.
 
-            hig_init (``Optional[float]``, default ``1.0``): Upper limit for the search algorithm to start
-              If ``None`` it the upper limit will be determined by :math:`\hat\mu + 2.5\sigma_{\hat\mu}`.
+            hig_init (``Optional[float]``, default ``1.0``): Initial upper limit for the search
+              algorithm. If `None`, it is determined by :math:`\hat\mu + 2.5\sigma_{\hat\mu}`.
+              Default is 1.0.
 
               .. note::
 
@@ -884,30 +880,36 @@ def poi_upper_limit(
 
             expected_pvalue (``Literal["nominal", "1sigma", "2sigma"]``, default ``"nominal"``):
               In case of :obj:`~spey.ExpectationType.aposteriori` and :obj:`~spey.ExpectationType.apriori`
-              expectation, gives the choice to find excluded upper limit for statistical deviations as well.
+              expectation, specifies the type of expected p-value for upper limit calculation.
 
-              * ``"nominal"``: only find the upper limit for the central p-value. Returns a single value.
-              * ``"1sigma"``: find the upper limit for central p-value and :math:`1\sigma` fluctuation from
-                background. Returns 3 values.
-              * ``"2sigma"``: find the upper limit for central p-value and :math:`1\sigma` and
-                :math:`2\sigma` fluctuation from background. Returns 5 values.
+              * ``"nominal"``: Computes the upper limit for the central p-value. Returns a single value.
+              * ``"1sigma"``: Computes the upper limit for the central p-value and :math:`1\sigma`
+                fluctuation from background. Returns 3 values.
+              * ``"2sigma"``: Computes the upper limit for the central p-value and :math:`1\sigma`
+                and :math:`2\sigma` fluctuation from background. Returns 5 values.
 
               .. note::
 
                 For ``expected=spey.ExpectationType.observed``, ``expected_pvalue`` argument will
                 be overwritten to ``"nominal"``.
 
-            maxiter (``int``, default ``10000``): Maximum iteration limit for the optimiser.
-            optimiser_arguments (``Dict``, default ``None``): Arguments for optimiser that is used
-              to compute likelihood and its maximum.
+            allow_negative_signal (``bool``, default ``True``): Allows for negative signal values,
+                changing the computation of the test statistic. Default is False.
+            maxiter (``int``, default ``10000``): Maximum number of iterations for the optimiser.
+                Default is 10000.
+            optimiser_arguments (``Dict``, default ``None``): Additional arguments for the optimiser
+                used to compute the likelihood and its maximum. Default is None.
 
         Returns:
             ``Union[float, List[float]]``:
-            In case of nominal values it returns a single value for the upper limit. In case of
-            ``expected_pvalue="1sigma"`` or ``expected_pvalue="2sigma"`` it will return a list of
-            multiple upper limit values for fluctuations as well as the central value. The
-            output order is :math:`-2\sigma` value, :math:`-1\sigma` value, central value,
-            :math:`1\sigma` and :math:`2\sigma` value.
+
+            - A single value representing the upper limit for the nominal case.
+            - A list of values representing the upper limits for the central value and statistical
+              deviations (for "1sigma" and "2sigma" cases). The order is: :math:`-2\sigma`,
+              :math:`-1\sigma`, central value, :math:`1\sigma`, :math:`2\sigma`.
+
+        Raises:
+            AssertionError: If the confidence level is not between 0 and 1.
         """
         assert (
             0.0 <= confidence_level <= 1.0
@@ -953,7 +955,7 @@ def poi_upper_limit(
             asimov_logpdf=logpdf_asimov,
             expected=expected,
             confidence_level=confidence_level,
-            allow_negative_signal=False,
+            allow_negative_signal=allow_negative_signal,
             low_init=low_init,
             hig_init=hig_init,
             expected_pvalue=expected_pvalue,
diff --git a/src/spey/hypothesis_testing/upper_limits.py b/src/spey/hypothesis_testing/upper_limits.py
index db8073d..dd3d751 100644
--- a/src/spey/hypothesis_testing/upper_limits.py
+++ b/src/spey/hypothesis_testing/upper_limits.py
@@ -142,7 +142,8 @@ def find_poi_upper_limit(
 
         confidence_level (``float``, default ``0.95``): Determines the confidence level of the upper
               limit i.e. the value of :math:`1-CL_s`. It needs to be between ``[0,1]``.
-        allow_negative_signal (``bool``, default ``True``): _description_
+        allow_negative_signal (``bool``, default ``True``): allow for negative signal values. This will
+            change the computation of the test statistic.
         low_init (``float``, default ``None``): Lower limit for the search algorithm to start
         hig_init (``float``, default ``None``): Upper limit for the search algorithm to start
         expected_pvalue (``Text``, default ``"nominal"``): In case of :obj:`~spey.ExpectationType.aposteriori`

From c849d21e0b9db1bb468f3b2055535299c75d8036 Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Wed, 29 Jan 2025 10:44:49 -0500
Subject: [PATCH 03/32] allow upto python 3.13

---
 setup.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/setup.py b/setup.py
index 929b3ef..a09c1d0 100644
--- a/setup.py
+++ b/setup.py
@@ -48,7 +48,7 @@
     packages=find_packages(where="src"),
     entry_points={"spey.backend.plugins": backend_plugins},
     install_requires=requirements,
-    python_requires=">=3.8, <=3.12.6",
+    python_requires=">=3.8, <3.13",
     classifiers=[
         "Intended Audience :: Science/Research",
         "License :: OSI Approved :: MIT License",

From f73998c5ce10086325717f53b27e095c6a9cb758 Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Wed, 29 Jan 2025 10:46:09 -0500
Subject: [PATCH 04/32] bump the version

---
 .zenodo.json         | 6 +++---
 src/spey/_version.py | 2 +-
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/.zenodo.json b/.zenodo.json
index 4d4fab6..7340888 100644
--- a/.zenodo.json
+++ b/.zenodo.json
@@ -1,8 +1,8 @@
 {
     "description": "smooth inference for reinterpretation studies",
     "license": "MIT",
-    "title": "SpeysideHEP/spey: v0.1.11",
-    "version": "v0.1.11",
+    "title": "SpeysideHEP/spey: v0.1.12",
+    "version": "v0.1.12",
     "upload_type": "software",
     "creators": [
         {
@@ -31,7 +31,7 @@
         },
         {
             "scheme": "url",
-            "identifier": "https://github.com/SpeysideHEP/spey/tree/v0.1.11",
+            "identifier": "https://github.com/SpeysideHEP/spey/tree/v0.1.12",
             "relation": "isSupplementTo"
         },
         {
diff --git a/src/spey/_version.py b/src/spey/_version.py
index fa9697c..805b7d6 100644
--- a/src/spey/_version.py
+++ b/src/spey/_version.py
@@ -1,3 +1,3 @@
 """Version number (major.minor.patch[-label])"""
 
-__version__ = "0.1.11"
+__version__ = "0.1.12-beta"

From e2da2d05400136111f08d902013cc784f53cdafb Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Wed, 29 Jan 2025 10:57:26 -0500
Subject: [PATCH 05/32] move default_pdf to default

---
 setup.py                                    | 14 +++++++-------
 src/spey/__init__.py                        |  8 +++++++-
 src/spey/backends/default_pdf/__init__.py   | 14 +++++++-------
 src/spey/backends/default_pdf/simple_pdf.py |  6 +++---
 4 files changed, 24 insertions(+), 18 deletions(-)

diff --git a/setup.py b/setup.py
index a09c1d0..410ea74 100644
--- a/setup.py
+++ b/setup.py
@@ -17,13 +17,13 @@
 ]
 
 backend_plugins = [
-    "default_pdf.uncorrelated_background = spey.backends.default_pdf:UncorrelatedBackground",
-    "default_pdf.correlated_background = spey.backends.default_pdf:CorrelatedBackground",
-    "default_pdf.third_moment_expansion = spey.backends.default_pdf:ThirdMomentExpansion",
-    "default_pdf.effective_sigma = spey.backends.default_pdf:EffectiveSigma",
-    "default_pdf.poisson = spey.backends.default_pdf.simple_pdf:Poisson",
-    "default_pdf.normal = spey.backends.default_pdf.simple_pdf:Gaussian",
-    "default_pdf.multivariate_normal = spey.backends.default_pdf.simple_pdf:MultivariateNormal",
+    "default.uncorrelated_background = spey.backends.default_pdf:UncorrelatedBackground",
+    "default.correlated_background = spey.backends.default_pdf:CorrelatedBackground",
+    "default.third_moment_expansion = spey.backends.default_pdf:ThirdMomentExpansion",
+    "default.effective_sigma = spey.backends.default_pdf:EffectiveSigma",
+    "default.poisson = spey.backends.default_pdf.simple_pdf:Poisson",
+    "default.normal = spey.backends.default_pdf.simple_pdf:Gaussian",
+    "default.multivariate_normal = spey.backends.default_pdf.simple_pdf:MultivariateNormal",
 ]
 
 setup(
diff --git a/src/spey/__init__.py b/src/spey/__init__.py
index ffafa8e..c4e0b0a 100644
--- a/src/spey/__init__.py
+++ b/src/spey/__init__.py
@@ -150,7 +150,7 @@ def get_backend(name: str) -> Callable[[Any], StatisticalModel]:
         :linenos:
 
         >>> import spey; import numpy as np
-        >>> stat_wrapper = spey.get_backend("default_pdf.uncorrelated_background")
+        >>> stat_wrapper = spey.get_backend("default.uncorrelated_background")
 
         >>> data = np.array([1])
         >>> signal = np.array([0.5])
@@ -175,6 +175,12 @@ def get_backend(name: str) -> Callable[[Any], StatisticalModel]:
         the backend it self which is in this particular example
         :obj:`~spey.backends.simplifiedlikelihood_backend.interface.SimplifiedLikelihoodInterface`.
     """
+    if "default_pdf" in name:
+        log.warning(
+            f"`{name}` has been deprecated, please use `{name.replace('default_pdf', 'default')}` instead."
+            " This will be removed in the future."
+        )
+        name = name.replace("default_pdf", "default")
     backend = _backend_entries.get(name, False)
 
     if backend:
diff --git a/src/spey/backends/default_pdf/__init__.py b/src/spey/backends/default_pdf/__init__.py
index da8fd70..bf7a9b3 100644
--- a/src/spey/backends/default_pdf/__init__.py
+++ b/src/spey/backends/default_pdf/__init__.py
@@ -55,7 +55,7 @@ class DefaultPDFBase(BackendBase):
         :func:`autograd.numpy.array` function.
     """
 
-    name: str = "default_pdf.base"
+    name: str = "default.base"
     """Name of the backend"""
     version: str = __version__
     """Version of the backend"""
@@ -444,7 +444,7 @@ class UncorrelatedBackground(DefaultPDFBase):
     .. code:: python3
 
         >>> import spey
-        >>> stat_wrapper = spey.get_backend('default_pdf.uncorrelated_background')
+        >>> stat_wrapper = spey.get_backend('default.uncorrelated_background')
 
         >>> data = [1, 3]
         >>> signal = [0.5, 2.0]
@@ -457,7 +457,7 @@ class UncorrelatedBackground(DefaultPDFBase):
         >>> print("1-CLs : %.3f" % tuple(stat_model.exclusion_confidence_level()))
     """
 
-    name: str = "default_pdf.uncorrelated_background"
+    name: str = "default.uncorrelated_background"
     """Name of the backend"""
     version: str = DefaultPDFBase.version
     """Version of the backend"""
@@ -575,7 +575,7 @@ class CorrelatedBackground(DefaultPDFBase):
 
         >>> import spey
 
-        >>> stat_wrapper = spey.get_backend('default_pdf.correlated_background')
+        >>> stat_wrapper = spey.get_backend('default.correlated_background')
 
         >>> signal_yields = [12.0, 11.0]
         >>> background_yields = [50.0, 52.0]
@@ -586,7 +586,7 @@ class CorrelatedBackground(DefaultPDFBase):
         >>> print(statistical_model.exclusion_confidence_level())
     """
 
-    name: str = "default_pdf.correlated_background"
+    name: str = "default.correlated_background"
     """Name of the backend"""
     version: str = DefaultPDFBase.version
     """Version of the backend"""
@@ -663,7 +663,7 @@ class ThirdMomentExpansion(DefaultPDFBase):
         ``third_moment`` should also have three inputs.
     """
 
-    name: str = "default_pdf.third_moment_expansion"
+    name: str = "default.third_moment_expansion"
     """Name of the backend"""
     version: str = DefaultPDFBase.version
     """Version of the backend"""
@@ -803,7 +803,7 @@ class EffectiveSigma(DefaultPDFBase):
           * third_moments (``List[float]``): diagonal elemetns of the third moment
     """
 
-    name: str = "default_pdf.effective_sigma"
+    name: str = "default.effective_sigma"
     """Name of the backend"""
     version: str = DefaultPDFBase.version
     """Version of the backend"""
diff --git a/src/spey/backends/default_pdf/simple_pdf.py b/src/spey/backends/default_pdf/simple_pdf.py
index a8a9c1f..ff8b6c9 100644
--- a/src/spey/backends/default_pdf/simple_pdf.py
+++ b/src/spey/backends/default_pdf/simple_pdf.py
@@ -304,7 +304,7 @@ class Poisson(SimplePDFBase):
           :math:`\mu n_s^i + n_b^i + \theta^i\sigma^i`.
     """
 
-    name: str = "default_pdf.poisson"
+    name: str = "default.poisson"
     """Name of the backend"""
     version: str = SimplePDFBase.version
     """Version of the backend"""
@@ -381,7 +381,7 @@ class Gaussian(SimplePDFBase):
         absolute_uncertainties (``List[float]``): absolute uncertainties on the background
     """
 
-    name: str = "default_pdf.normal"
+    name: str = "default.normal"
     """Name of the backend"""
     version: str = SimplePDFBase.version
     """Version of the backend"""
@@ -431,7 +431,7 @@ class MultivariateNormal(SimplePDFBase):
 
     """
 
-    name: str = "default_pdf.multivariate_normal"
+    name: str = "default.multivariate_normal"
     """Name of the backend"""
     version: str = SimplePDFBase.version
     """Version of the backend"""

From f2225b9ddb17d54db1dd5fadd084e65fb269e8ce Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Wed, 29 Jan 2025 11:11:39 -0500
Subject: [PATCH 06/32] move default_pdf to default

---
 README.md                                     | 28 ++++-----
 docs/comb.rst                                 |  2 +-
 docs/exclusion.md                             |  4 +-
 docs/plugins.rst                              | 54 ++++++++--------
 docs/quick_start.rst                          | 14 ++---
 docs/tutorials/functional_tuto.ipynb          | 12 ++--
 docs/tutorials/gradients.md                   |  4 +-
 docs/tutorials/histogram.ipynb                |  6 +-
 docs/tutorials/intro_spey.ipynb               | 62 ++++++++++---------
 docs/tutorials/signal_unc.md                  |  4 +-
 src/spey/__init__.py                          |  4 +-
 .../uncorrelated_statistics_combiner.py       | 10 +--
 tests/test_combiner.py                        | 10 +--
 tests/test_default_pdf.py                     | 14 ++---
 tests/test_math.py                            |  2 +-
 tests/test_simplifiedllhd_limit.py            |  8 +--
 16 files changed, 122 insertions(+), 116 deletions(-)

diff --git a/README.md b/README.md
index 0b2c053..e16fc92 100644
--- a/README.md
+++ b/README.md
@@ -55,13 +55,13 @@ Finally, the name "Spey" originally comes from the Spey River, a river in the mi
 
 | Accessor                                | Description                                                                  |
 | --------------------------------------- | ---------------------------------------------------------------------------- |
-| ``"default_pdf.uncorrelated_background"`` | Constructs uncorrelated multi-bin statistical model.                         |
-| ``"default_pdf.correlated_background"``   | Constructs correlated multi-bin statistical model with Gaussian nuisances.   |
-| ``"default_pdf.third_moment_expansion"``  | Implements the skewness of the likelihood by using third moments.            |
-| ``"default_pdf.effective_sigma"``         | Implements the skewness of the likelihood by using asymmetric uncertainties. |
-| ``"default_pdf.poisson"``                 | Implements simple Poisson-based likelihood without uncertainties. |
-| ``"default_pdf.normal"``                  | Implements Normal distribution. |
-| ``"default_pdf.multivariate_normal"``     | Implements Multivariate normal distribution. |
+| ``"default.uncorrelated_background"`` | Constructs uncorrelated multi-bin statistical model.                         |
+| ``"default.correlated_background"``   | Constructs correlated multi-bin statistical model with Gaussian nuisances.   |
+| ``"default.third_moment_expansion"``  | Implements the skewness of the likelihood by using third moments.            |
+| ``"default.effective_sigma"``         | Implements the skewness of the likelihood by using asymmetric uncertainties. |
+| ``"default.poisson"``                 | Implements simple Poisson-based likelihood without uncertainties. |
+| ``"default.normal"``                  | Implements Normal distribution. |
+| ``"default.multivariate_normal"``     | Implements Multivariate normal distribution. |
 | ``"pyhf.uncorrelated_background"``        | Uses uncorrelated background functionality of pyhf (see [``spey-phyf`` plugin](https://github.com/SpeysideHEP/spey-pyhf)).                          |
 | ``"pyhf"``                                | Uses generic likelihood structure of pyhf (see [``spey-phyf`` plugin](https://github.com/SpeysideHEP/spey-pyhf))                                   |
 
@@ -75,17 +75,17 @@ likelihood prescriptions which can be checked via `AvailableBackends` function
 ```python
 import spey
 print(spey.AvailableBackends())
-# ['default_pdf.correlated_background',
-#  'default_pdf.effective_sigma',
-#  'default_pdf.third_moment_expansion',
-#  'default_pdf.uncorrelated_background']
+# ['default.correlated_background',
+#  'default.effective_sigma',
+#  'default.third_moment_expansion',
+#  'default.uncorrelated_background']
 ```
 
-Using ``'default_pdf.uncorrelated_background'`` one can simply create single or multi-bin
+Using ``'default.uncorrelated_background'`` one can simply create single or multi-bin
 statistical models:
 
 ```python
-pdf_wrapper = spey.get_backend('default_pdf.uncorrelated_background')
+pdf_wrapper = spey.get_backend('default.uncorrelated_background')
 
 data = [1]
 signal_yields = [0.5]
@@ -147,7 +147,7 @@ collaborations to estimate the difference from the Standard Model rather than th
 One can play the same game using the same backend for multi-bin histograms as follows;
 
 ```python
-pdf_wrapper = spey.get_backend('default_pdf.uncorrelated_background')
+pdf_wrapper = spey.get_backend('default.uncorrelated_background')
 
 data = [36, 33]
 signal = [12.0, 15.0]
diff --git a/docs/comb.rst b/docs/comb.rst
index 7f11b9d..a800672 100644
--- a/docs/comb.rst
+++ b/docs/comb.rst
@@ -39,7 +39,7 @@ Let us first import all the necessary packages and construct the data (please ad
     >>> models = {}
     >>> # loop overall data
     >>> for data in example_data["data"]:
-    >>>     pdf_wrapper = spey.get_backend("default_pdf.uncorrelated_background")
+    >>>     pdf_wrapper = spey.get_backend("default.uncorrelated_background")
 
     >>>     stat_model = pdf_wrapper(
     ...         signal_yields=data["signal_yields"],
diff --git a/docs/exclusion.md b/docs/exclusion.md
index 910dc45..514bb62 100644
--- a/docs/exclusion.md
+++ b/docs/exclusion.md
@@ -73,7 +73,7 @@ The `allow_negative_signal` keyword controls which test statistic is used and re
 
 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$.
+Below we compare the exclusion limits computed with each approach. This comparisson uses normal distribution for the likelihood ([`default.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
 ---
@@ -90,7 +90,7 @@ exclusion limit calculator comparisson for observed p-values
 ```{code-cell} ipython3
 import spey
 
-stat_wrapper = spey.get_backend("default_pdf.normal")
+stat_wrapper = spey.get_backend("default.normal")
 statistical_model = stat_wrapper(
     signal_yields=[3.0],
     background_yields=[2.0],
diff --git a/docs/plugins.rst b/docs/plugins.rst
index 7433d79..2366bb0 100644
--- a/docs/plugins.rst
+++ b/docs/plugins.rst
@@ -8,19 +8,19 @@ Plug-ins
 
    * - Keyword
      - Summary
-   * - ``'default_pdf.uncorrelated_background'``
+   * - ``'default.uncorrelated_background'``
      - :ref:`Combination of Poisson and Gaussian PDF, assuming uncorrelated bins. <uncorrelated_background>`
-   * - ``'default_pdf.correlated_background'``
+   * - ``'default.correlated_background'``
      - :ref:`Combination of Poisson and Gaussian PDF, with correlated bins. <correlated_background>`
-   * - ``'default_pdf.third_moment_expansion'``
+   * - ``'default.third_moment_expansion'``
      - :ref:`Simplified likelihood, extended with third moments of the background. <third_moment_expansion>`
-   * - ``'default_pdf.effective_sigma'``
+   * - ``'default.effective_sigma'``
      - :ref:`Simplified likelihood, extended with asymmetric uncertainties.  <effective_sigma>`
-   * - ``'default_pdf.poisson'``
+   * - ``'default.poisson'``
      - :ref:`Poisson distribution, without uncertainties.  <poisson>`
-   * - ``'default_pdf.normal'``
+   * - ``'default.normal'``
      - :ref:`Gaussian distribution.  <normal>`
-   * - ``'default_pdf.multivariate_normal'``
+   * - ``'default.multivariate_normal'``
      - :ref:`Multivariate Normal distribution.  <multinormal>`
    * - ``'pyhf'``
      - `External plug-in <https://spey-pyhf.readthedocs.io>`_ uses ``pyhf`` to construct the likelihoods.
@@ -46,10 +46,10 @@ to the available plugins can be seen using the following command:
 .. code-block:: python3
 
     >>> spey.AvailableBackends()
-    >>> # ['default_pdf.correlated_background',
-    >>> #  'default_pdf.effective_sigma',
-    >>> #  'default_pdf.third_moment_expansion',
-    >>> #  'default_pdf.uncorrelated_background']
+    >>> # ['default.correlated_background',
+    >>> #  'default.effective_sigma',
+    >>> #  'default.third_moment_expansion',
+    >>> #  'default.uncorrelated_background']
 
 where once installed without any plug-ins :func:`~spey.AvailableBackends` function
 only shows the default PDFs. In the following section, we will summarise their usability.
@@ -58,7 +58,7 @@ function e.g.
 
 .. code-block:: python3
 
-    >>> pdf_wrapper = spey.get_backend('default_pdf.uncorrelated_background')
+    >>> pdf_wrapper = spey.get_backend('default.uncorrelated_background')
 
 this will automatically create a wrapper around the likelihood prescription and allow `spey`
 to use it properly. We will demonstrate the usage for each of the default plugins below.
@@ -89,7 +89,7 @@ The first term represents the primary model, and the second represents the const
 
 .. _uncorrelated_background:
 
-``'default_pdf.uncorrelated_background'``
+``'default.uncorrelated_background'``
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 This is a basic PDF where the background is assumed to be not correlated, where the PDF has been
@@ -108,7 +108,7 @@ used as follows
 .. code-block:: python3
     :linenos:
 
-    >>> pdf_wrapper = spey.get_backend("default_pdf.uncorrelated_background")
+    >>> pdf_wrapper = spey.get_backend("default.uncorrelated_background")
     >>> statistical_model = pdf_wrapper(
     ...     signal_yields=[12.0, 15.0],
     ...     background_yields=[50.0,48.0],
@@ -140,7 +140,7 @@ API description.
 
 .. _correlated_background:
 
-``'default_pdf.correlated_background'``
+``'default.correlated_background'``
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 This plugin embeds the correlations between each bin using a covariance matrix provided by the user
@@ -158,7 +158,7 @@ defined as
 .. code-block:: python3
     :linenos:
 
-    >>> pdf_wrapper = spey.get_backend("default_pdf.correlated_background")
+    >>> pdf_wrapper = spey.get_backend("default.correlated_background")
     >>> statistical_model = pdf_wrapper(
     ...     signal_yields=[12.0, 15.0],
     ...     background_yields=[50.0,48.0],
@@ -191,7 +191,7 @@ as expected.
 
 .. _third_moment_expansion:
 
-``'default_pdf.third_moment_expansion'``
+``'default.third_moment_expansion'``
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 This plug-in implements the third-moment expansion presented in :cite:`Buckley:2018vdr`, which expands the
@@ -227,7 +227,7 @@ iterating over the same example, this PDF can be accessed as follows
 .. code-block:: python3
     :linenos:
 
-    >>> pdf_wrapper = spey.get_backend("default_pdf.third_moment_expansion")
+    >>> pdf_wrapper = spey.get_backend("default.third_moment_expansion")
     >>> statistical_model = pdf_wrapper(
     ...     signal_yields=[12.0, 15.0],
     ...     background_yields=[50.0,48.0],
@@ -257,14 +257,14 @@ reduced the exclusion limit.
    For absolute uncertainty :math:`\sigma_b`; :math:`\sigma_b = \sqrt{{\rm diag}(\Sigma)}`. The covariance matrix should be a square matrix
    and both dimensions should match the number of ``background_yields`` given as input.
  * ``third_moment``: Diagonal elements of the third moments. These can be computed using
-   :func:`~spey.backends.default_pdf.third_moment.compute_third_moments` function; however this function computes third moments using
+   :func:`~spey.backends.default.third_moment.compute_third_moments` function; however this function computes third moments using
    Bifurcated Gaussian, which may not be suitable for every case.
  * ``analysis`` (optional): Unique identifier for the analysis.
  * ``xsection`` (optional): Cross-section value for the signal hypothesis. Units determined by the user.
 
 .. _effective_sigma:
 
-``'default_pdf.effective_sigma'``
+``'default.effective_sigma'``
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 The skewness of the PDF distribution can also be captured by building an effective variance
@@ -288,7 +288,7 @@ iterating over the same example, this PDF can be utilised as follows;
 .. code-block:: python3
     :linenos:
 
-    >>> pdf_wrapper = spey.get_backend("default_pdf.effective_sigma")
+    >>> pdf_wrapper = spey.get_backend("default.effective_sigma")
     >>> statistical_model = pdf_wrapper(
     ...     signal_yields=[12.0, 15.0],
     ...     background_yields=[50.0,48.0],
@@ -324,7 +324,7 @@ Once again, the exclusion limit can be computed as
 
 .. _poisson:
 
-``'default_pdf.poisson'``
+``'default.poisson'``
 ~~~~~~~~~~~~~~~~~~~~~~~~~
 
 Simple Poisson implementation without uncertainties which can be described as follows;
@@ -338,7 +338,7 @@ It can take any number of yields.
 .. code-block:: python3
     :linenos:
 
-    >>> pdf_wrapper = spey.get_backend("default_pdf.poisson")
+    >>> pdf_wrapper = spey.get_backend("default.poisson")
     >>> statistical_model = pdf_wrapper(
     ...     signal_yields=[12.0, 15.0],
     ...     background_yields=[50.0,48.0],
@@ -357,7 +357,7 @@ It can take any number of yields.
 
 .. _normal:
 
-``'default_pdf.normal'``
+``'default.normal'``
 ~~~~~~~~~~~~~~~~~~~~~~~~~
 
 Simple Normal distribution implementation;
@@ -371,7 +371,7 @@ It can take any number of yields.
 .. code-block:: python3
     :linenos:
 
-    >>> pdf_wrapper = spey.get_backend("default_pdf.normal")
+    >>> pdf_wrapper = spey.get_backend("default.normal")
     >>> statistical_model = pdf_wrapper(
     ...     signal_yields=[12.0, 15.0],
     ...     background_yields=[50.0,48.0],
@@ -393,7 +393,7 @@ It can take any number of yields.
 
 .. _multinormal:
 
-``'default_pdf.multivariate_normal'``
+``'default.multivariate_normal'``
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 Simple Normal distribution implementation;
@@ -407,7 +407,7 @@ It can take any number of yields.
 .. code-block:: python3
     :linenos:
 
-    >>> pdf_wrapper = spey.get_backend("default_pdf.multivariate_normal")
+    >>> pdf_wrapper = spey.get_backend("default.multivariate_normal")
     >>> statistical_model = pdf_wrapper(
     ...     signal_yields=[12.0, 15.0],
     ...     background_yields=[50.0,48.0],
diff --git a/docs/quick_start.rst b/docs/quick_start.rst
index 1a451b5..e142e5d 100644
--- a/docs/quick_start.rst
+++ b/docs/quick_start.rst
@@ -50,19 +50,19 @@ function
 
     >>> import spey
     >>> print(spey.AvailableBackends())
-    >>> # ['default_pdf.correlated_background',
-    >>> #  'default_pdf.effective_sigma',
-    >>> #  'default_pdf.third_moment_expansion',
-    >>> #  'default_pdf.uncorrelated_background']
+    >>> # ['default.correlated_background',
+    >>> #  'default.effective_sigma',
+    >>> #  'default.third_moment_expansion',
+    >>> #  'default.uncorrelated_background']
 
 For details on all the backends, see `Plug-ins section <plugins.html>`_.
 
-Using ``'default_pdf.uncorrelated_background'`` one can simply create single or multi-bin
+Using ``'default.uncorrelated_background'`` one can simply create single or multi-bin
 statistical models:
 
 .. code:: python
 
-    >>> pdf_wrapper = spey.get_backend('default_pdf.uncorrelated_background')
+    >>> pdf_wrapper = spey.get_backend('default.uncorrelated_background')
 
     >>> data = [1]
     >>> signal_yields = [0.5]
@@ -128,7 +128,7 @@ One can play the same game using the same backend for multi-bin histograms as fo
 
 .. code:: python
 
-    >>> pdf_wrapper = spey.get_backend('default_pdf.uncorrelated_background')
+    >>> pdf_wrapper = spey.get_backend('default.uncorrelated_background')
 
     >>> data = [36, 33]
     >>> signal_yields = [12.0, 15.0]
diff --git a/docs/tutorials/functional_tuto.ipynb b/docs/tutorials/functional_tuto.ipynb
index 0847aa3..9602690 100644
--- a/docs/tutorials/functional_tuto.ipynb
+++ b/docs/tutorials/functional_tuto.ipynb
@@ -127,7 +127,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "pdf_wrapper = spey.get_backend(\"default_pdf.poisson\")"
+    "pdf_wrapper = spey.get_backend(\"default.poisson\")"
    ]
   },
   {
@@ -271,7 +271,7 @@
    },
    "outputs": [],
    "source": [
-    "pdf_wrapper = spey.get_backend(\"default_pdf.uncorrelated_background\")"
+    "pdf_wrapper = spey.get_backend(\"default.uncorrelated_background\")"
    ]
   },
   {
@@ -387,14 +387,14 @@
     "chi2 = []\n",
     "for c1 in np.linspace(-1, 1, 10):\n",
     "    for c2 in np.linspace(-1, 1, 10):\n",
-    "        stat_model_with_unc = spey.get_backend(\"default_pdf.uncorrelated_background\")(\n",
+    "        stat_model_with_unc = spey.get_backend(\"default.uncorrelated_background\")(\n",
     "            signal_yields=signal(c1, c2),\n",
     "            background_yields=background,\n",
     "            data=data,\n",
     "            absolute_uncertainties=bkg_unc,\n",
     "            analysis=\"with uncertainty\",\n",
     "        )\n",
-    "        stat_model = spey.get_backend(\"default_pdf.poisson\")(\n",
+    "        stat_model = spey.get_backend(\"default.poisson\")(\n",
     "            signal_yields=signal(c1, c2),\n",
     "            background_yields=background,\n",
     "            data=data,\n",
@@ -464,7 +464,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "base",
    "language": "python",
    "name": "python3"
   },
@@ -478,7 +478,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.16"
+   "version": "3.9.18"
   }
  },
  "nbformat": 4,
diff --git a/docs/tutorials/gradients.md b/docs/tutorials/gradients.md
index c1e4098..2597799 100644
--- a/docs/tutorials/gradients.md
+++ b/docs/tutorials/gradients.md
@@ -46,10 +46,10 @@ np.random.seed(14)
 
 {py:func}`spey.math.value_and_grad` returns a function that computes negative log-likelihood and its gradient for a given statistical model and {py:func}`spey.math.hessian` returns a function that computes Hessian of negative log-likelihood.
 
-Let us examine this on ``"default_pdf.uncorrelated_background"``:
+Let us examine this on ``"default.uncorrelated_background"``:
 
 ```{code-cell} ipython3
-pdf_wrapper = spey.get_backend("default_pdf.uncorrelated_background")
+pdf_wrapper = spey.get_backend("default.uncorrelated_background")
 
 data = [36, 33]
 signal_yields = [12.0, 15.0]
diff --git a/docs/tutorials/histogram.ipynb b/docs/tutorials/histogram.ipynb
index ddbaa4c..b5f85ee 100644
--- a/docs/tutorials/histogram.ipynb
+++ b/docs/tutorials/histogram.ipynb
@@ -109,7 +109,7 @@
    "id": "741b7071-56ca-42b8-9d99-abdcb88d318a",
    "metadata": {},
    "source": [
-    "Lets construct a statistical model. Since the only source of uncertainty is statistical, we can safely assume that non of the bins are correlated in systematic level. Hence we will use [`\"default_pdf.uncorrelated_background\"`](https://spey.readthedocs.io/en/main/plugins.html#default-pdf-uncorrelated-background) model."
+    "Lets construct a statistical model. Since the only source of uncertainty is statistical, we can safely assume that non of the bins are correlated in systematic level. Hence we will use [`\"default.uncorrelated_background\"`](https://spey.readthedocs.io/en/main/plugins.html#default-pdf-uncorrelated-background) model."
    ]
   },
   {
@@ -134,7 +134,7 @@
    },
    "outputs": [],
    "source": [
-    "pdf_wrapper = spey.get_backend(\"default_pdf.uncorrelated_background\")"
+    "pdf_wrapper = spey.get_backend(\"default.uncorrelated_background\")"
    ]
   },
   {
@@ -433,7 +433,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "base",
    "language": "python",
    "name": "python3"
   },
diff --git a/docs/tutorials/intro_spey.ipynb b/docs/tutorials/intro_spey.ipynb
index 21016fa..9af8446 100644
--- a/docs/tutorials/intro_spey.ipynb
+++ b/docs/tutorials/intro_spey.ipynb
@@ -10,7 +10,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "id": "8968127e-3e70-4971-8ced-81abacbaf363",
    "metadata": {
     "tags": [
@@ -29,7 +29,7 @@
    "id": "e1eea83b-de0b-4e08-b82a-b9b5eaf89c8d",
    "metadata": {},
    "source": [
-    "In this tutorial we will go over basic functionalities of Spey package. Spey is designed to be used with different plug-ins where each plug-in represents a likelihood prescription. Spey is shipped with set of default likelihood prescriptions which are identified as ```\"default_pdf\"```. The list of available likelihood prescriptions can be found using `spey.AvailableBackends()` function."
+    "In this tutorial we will go over basic functionalities of Spey package. Spey is designed to be used with different plug-ins where each plug-in represents a likelihood prescription. Spey is shipped with set of default likelihood prescriptions which are identified as ```\"default\"```. The list of available likelihood prescriptions can be found using `spey.AvailableBackends()` function."
    ]
   },
   {
@@ -41,11 +41,16 @@
     {
      "data": {
       "text/plain": [
-       "['default_pdf.correlated_background',\n",
-       " 'default_pdf.effective_sigma',\n",
-       " 'default_pdf.poisson',\n",
-       " 'default_pdf.third_moment_expansion',\n",
-       " 'default_pdf.uncorrelated_background']"
+       "['pyhf',\n",
+       " 'pyhf.simplify',\n",
+       " 'pyhf.uncorrelated_background',\n",
+       " 'default.correlated_background',\n",
+       " 'default.effective_sigma',\n",
+       " 'default.multivariate_normal',\n",
+       " 'default.normal',\n",
+       " 'default.poisson',\n",
+       " 'default.third_moment_expansion',\n",
+       " 'default.uncorrelated_background']"
       ]
      },
      "execution_count": 2,
@@ -72,7 +77,7 @@
    "source": [
     "Designing custom likelihood prescriptions are also possible, details can be found in the appendix of [arXiv:2307.06996](https://arxiv.org/abs/2307.06996) or through [this link](https://speysidehep.github.io/spey/new_plugin.html).\n",
     "\n",
-    "Each likelihood prescription is required to have certain metadata structure which will provide other users the necessary information to cite them. These metadata can be accessed via `spey.get_backend_metadata(\"<plugin>\")` function. For instance lets use it for ```'default_pdf.third_moment_expansion'```:"
+    "Each likelihood prescription is required to have certain metadata structure which will provide other users the necessary information to cite them. These metadata can be accessed via `spey.get_backend_metadata(\"<plugin>\")` function. For instance lets use it for ```'default.third_moment_expansion'```:"
    ]
   },
   {
@@ -84,12 +89,13 @@
     {
      "data": {
       "text/plain": [
-       "{'name': 'default_pdf.third_moment_expansion',\n",
+       "{'name': 'default.third_moment_expansion',\n",
        " 'author': 'SpeysideHEP',\n",
-       " 'version': '0.1.4',\n",
-       " 'spey_requires': '0.1.4',\n",
+       " 'version': '0.1.12-beta',\n",
+       " 'spey_requires': '0.1.12-beta',\n",
        " 'doi': ['10.1007/JHEP04(2019)064'],\n",
-       " 'arXiv': ['1809.05548']}"
+       " 'arXiv': ['1809.05548'],\n",
+       " 'zenodo': []}"
       ]
      },
      "execution_count": 3,
@@ -98,7 +104,7 @@
     }
    ],
    "source": [
-    "spey.get_backend_metadata(\"default_pdf.third_moment_expansion\")"
+    "spey.get_backend_metadata(\"default.third_moment_expansion\")"
    ]
   },
   {
@@ -119,7 +125,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 27,
+   "execution_count": 6,
    "id": "e5688f9e-d999-4f76-a2fb-c479722ae6ad",
    "metadata": {
     "tags": [
@@ -171,7 +177,7 @@
     "\n",
     "$$ \\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) \\cdot \\prod_{j\\in{\\rm nui}}\\mathcal{N}(\\theta^j|0, \\rho) $$\n",
     "\n",
-    "where $n,\\ n_s$ and $n_b$ stands for observed, signal and background yields, $\\sigma_b$ are background uncertainties and $\\rho$ is the correlation matrix. $\\mu$ and $\\theta$ on the other hand stands for the POI and nuisance parameters. For more details about this approximation see [CMS-NOTE-2017-001](https://cds.cern.ch/record/2242860). We will refer to this particular approximation as ```\"default_pdf.correlated_background\"``` model.\n",
+    "where $n,\\ n_s$ and $n_b$ stands for observed, signal and background yields, $\\sigma_b$ are background uncertainties and $\\rho$ is the correlation matrix. $\\mu$ and $\\theta$ on the other hand stands for the POI and nuisance parameters. For more details about this approximation see [CMS-NOTE-2017-001](https://cds.cern.ch/record/2242860). We will refer to this particular approximation as ```\"default.correlated_background\"``` model.\n",
     "\n",
     "The bottom panel above shows the exclusion limits computed by the experimental collaboration (black), correlated background model (red), third moment expansion (blue) and effective sigma (green). This plot has been retreived from [arXiv:2307.06996](https://arxiv.org/abs/2307.06996), Fig 1.\n",
     "\n",
@@ -193,7 +199,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 29,
+   "execution_count": 4,
    "id": "75a18ded-53db-4697-bfdf-be9ef2500026",
    "metadata": {
     "collapsed": true,
@@ -204,12 +210,12 @@
    },
    "outputs": [],
    "source": [
-    "pdf_wrapper = spey.get_backend(\"default_pdf.correlated_background\")"
+    "pdf_wrapper = spey.get_backend(\"default.correlated_background\")"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 30,
+   "execution_count": 7,
    "id": "16664ab1-6a3e-4b0f-ab34-554e281f02d9",
    "metadata": {},
    "outputs": [
@@ -217,7 +223,7 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "StatisticalModel(analysis='cms-sus-20-004-SL', backend=default_pdf.correlated_background, calculators=['toy', 'asymptotic', 'chi_square'])\n"
+      "StatisticalModel(analysis='cms-sus-20-004-SL', backend=default.correlated_background, calculators=['toy', 'asymptotic', 'chi_square'])\n"
      ]
     }
    ],
@@ -546,7 +552,7 @@
     "\n",
     "***Experimental collaboration provided asymmetric uncertainties but not third moments, can spey compute third moments?***\n",
     "\n",
-    "**Answer:** Yes! one can use [`compute_third_moments`](https://spey.readthedocs.io/en/main/_generated/spey.backends.default_pdf.third_moment.compute_third_moments.html#spey.backends.default_pdf.third_moment.compute_third_moments) function which computes third moments using Bifurcated Gaussian\n",
+    "**Answer:** Yes! one can use [`compute_third_moments`](https://spey.readthedocs.io/en/main/_generated/spey.backends.default.third_moment.compute_third_moments.html#spey.backends.default.third_moment.compute_third_moments) function which computes third moments using Bifurcated Gaussian\n",
     "\n",
     "$$ \n",
     "m^{(3)} = \\frac{2}{\\sigma^+ + \\sigma^-} \\left[ \\sigma^-\\int_{-\\infty}^0 x^3\\mathcal{N}(x|0,\\sigma^-)dx + \\sigma^+ \\int_0^{\\infty} x^3\\mathcal{N}(x|0,\\sigma^+)dx \\right]\n",
@@ -558,7 +564,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 18,
+   "execution_count": 8,
    "id": "e8a1aa93-5fef-4c83-95b1-0559818cb8fe",
    "metadata": {},
    "outputs": [
@@ -566,12 +572,12 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "StatisticalModel(analysis='cms-sus-20-004-TM', backend=default_pdf.third_moment_expansion, calculators=['toy', 'asymptotic', 'chi_square'])\n"
+      "StatisticalModel(analysis='cms-sus-20-004-TM', backend=default.third_moment_expansion, calculators=['toy', 'asymptotic', 'chi_square'])\n"
      ]
     }
    ],
    "source": [
-    "pdf_wrapper = spey.get_backend(\"default_pdf.third_moment_expansion\")\n",
+    "pdf_wrapper = spey.get_backend(\"default.third_moment_expansion\")\n",
     "third_mom_expansion_model = pdf_wrapper(\n",
     "    signal_yields=signal_yields,\n",
     "    background_yields=background_yields,\n",
@@ -626,7 +632,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 20,
+   "execution_count": 9,
    "id": "0e4a7660-71f3-4316-be88-6c2431c48a56",
    "metadata": {
     "tags": [
@@ -640,7 +646,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 21,
+   "execution_count": 10,
    "id": "9091d3c1-c009-41d3-9438-ae0fe5f8a6fa",
    "metadata": {},
    "outputs": [
@@ -648,12 +654,12 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "StatisticalModel(analysis='cms-sus-20-004-ES', backend=default_pdf.effective_sigma, calculators=['toy', 'asymptotic', 'chi_square'])\n"
+      "StatisticalModel(analysis='cms-sus-20-004-ES', backend=default.effective_sigma, calculators=['toy', 'asymptotic', 'chi_square'])\n"
      ]
     }
    ],
    "source": [
-    "pdf_wrapper = spey.get_backend(\"default_pdf.effective_sigma\")\n",
+    "pdf_wrapper = spey.get_backend(\"default.effective_sigma\")\n",
     "effective_sigma_model = pdf_wrapper(\n",
     "    signal_yields=signal_yields,\n",
     "    background_yields=background_yields,\n",
@@ -751,7 +757,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "base",
    "language": "python",
    "name": "python3"
   },
diff --git a/docs/tutorials/signal_unc.md b/docs/tutorials/signal_unc.md
index e5a61e7..4923c95 100644
--- a/docs/tutorials/signal_unc.md
+++ b/docs/tutorials/signal_unc.md
@@ -29,7 +29,7 @@ Note that theoretical uncertainties have different interpretations, we can inter
     import spey
     from math import sqrt
 
-    pdf_wrapper = spey.get_backend("default_pdf.uncorrelated_background")
+    pdf_wrapper = spey.get_backend("default.uncorrelated_background")
     statistical_model_sigunc = pdf_wrapper(
         signal_yields=[12.0],
         background_yields=[50.0],
@@ -57,7 +57,7 @@ $$
 where $(s)$ superscript indicates signal and $(b)$ indicates background.
 
 ```{code-cell} ipython3
-pdf_wrapper = spey.get_backend("default_pdf.uncorrelated_background")
+pdf_wrapper = spey.get_backend("default.uncorrelated_background")
 statistical_model_sigunc = pdf_wrapper(
     signal_yields=[12.0, 15.0],
     background_yields=[50.0, 48.0],
diff --git a/src/spey/__init__.py b/src/spey/__init__.py
index c4e0b0a..e52da26 100644
--- a/src/spey/__init__.py
+++ b/src/spey/__init__.py
@@ -237,13 +237,13 @@ def get_backend_metadata(name: str) -> Dict[str, Any]:
 
     .. code-block:: python3
 
-        >>> spey.get_backend_metadata("default_pdf.third_moment_expansion")
+        >>> spey.get_backend_metadata("default.third_moment_expansion")
 
     will return the following
 
     .. code-block:: python3
 
-        >>> {'name': 'default_pdf.third_moment_expansion',
+        >>> {'name': 'default.third_moment_expansion',
         ... 'author': 'SpeysideHEP',
         ... 'version': '0.0.1',
         ... 'spey_requires': '0.0.1',
diff --git a/src/spey/combiner/uncorrelated_statistics_combiner.py b/src/spey/combiner/uncorrelated_statistics_combiner.py
index ae15c4b..0dcf665 100644
--- a/src/spey/combiner/uncorrelated_statistics_combiner.py
+++ b/src/spey/combiner/uncorrelated_statistics_combiner.py
@@ -234,7 +234,7 @@ def likelihood(
 
               .. code-block:: python3
 
-                >>> statistical_model_options = {"default_pdf.uncorrelated_background" : {"init_pars" : [1., 3., 4.]}}
+                >>> statistical_model_options = {"default.uncorrelated_background" : {"init_pars" : [1., 3., 4.]}}
 
             kwargs: keyword arguments for the optimiser.
 
@@ -325,7 +325,7 @@ def generate_asimov_data(
 
               .. code-block:: python3
 
-                >>> statistical_model_options = {"default_pdf.uncorrelated_background" : {"init_pars" : [1., 3., 4.]}}
+                >>> statistical_model_options = {"default.uncorrelated_background" : {"init_pars" : [1., 3., 4.]}}
 
             kwargs: keyword arguments for the optimiser.
 
@@ -409,7 +409,7 @@ def asimov_likelihood(
 
               .. code-block:: python3
 
-                >>> statistical_model_options = {"default_pdf.uncorrelated_background" : {"init_pars" : [1., 3., 4.]}}
+                >>> statistical_model_options = {"default.uncorrelated_background" : {"init_pars" : [1., 3., 4.]}}
 
             kwargs: keyword arguments for the optimiser.
 
@@ -483,7 +483,7 @@ def maximize_likelihood(
 
               .. code-block:: python3
 
-                >>> statistical_model_options = {"default_pdf.uncorrelated_background" : {"init_pars" : [1., 3., 4.]}}
+                >>> statistical_model_options = {"default.uncorrelated_background" : {"init_pars" : [1., 3., 4.]}}
 
             kwargs: keyword arguments for the optimiser.
 
@@ -613,7 +613,7 @@ def maximize_asimov_likelihood(
 
               .. code-block:: python3
 
-                >>> statistical_model_options = {"default_pdf.uncorrelated_background" : {"init_pars" : [1., 3., 4.]}}
+                >>> statistical_model_options = {"default.uncorrelated_background" : {"init_pars" : [1., 3., 4.]}}
 
             kwargs: keyword arguments for the optimiser.
 
diff --git a/tests/test_combiner.py b/tests/test_combiner.py
index 19a661c..33c1302 100644
--- a/tests/test_combiner.py
+++ b/tests/test_combiner.py
@@ -9,14 +9,14 @@
 def test_combiner():
     """Test uncorrelated statistics combiner"""
 
-    normal = spey.get_backend("default_pdf.normal")(
+    normal = spey.get_backend("default.normal")(
         signal_yields=[1, 1, 1],
         background_yields=[2, 1, 3],
         data=[2, 2, 2],
         absolute_uncertainties=[1, 1, 1],
     )
     normal_cls = normal.exclusion_confidence_level()[0]
-    normal = spey.get_backend("default_pdf.multivariate_normal")(
+    normal = spey.get_backend("default.multivariate_normal")(
         signal_yields=[1, 1, 1],
         background_yields=[2, 1, 3],
         data=[2, 2, 2],
@@ -28,21 +28,21 @@ def test_combiner():
         normal_cls, multivar_norm_cls
     ), "Normal CLs is not the same as Multivariant normal CLs"
 
-    normal1 = spey.get_backend("default_pdf.normal")(
+    normal1 = spey.get_backend("default.normal")(
         signal_yields=[1],
         background_yields=[3],
         data=[2],
         absolute_uncertainties=[1],
         analysis="norm1",
     )
-    normal2 = spey.get_backend("default_pdf.normal")(
+    normal2 = spey.get_backend("default.normal")(
         signal_yields=[1],
         background_yields=[1],
         data=[2],
         absolute_uncertainties=[1],
         analysis="norm2",
     )
-    normal3 = spey.get_backend("default_pdf.normal")(
+    normal3 = spey.get_backend("default.normal")(
         signal_yields=[1],
         background_yields=[2],
         data=[2],
diff --git a/tests/test_default_pdf.py b/tests/test_default_pdf.py
index 4bffa28..2167d02 100644
--- a/tests/test_default_pdf.py
+++ b/tests/test_default_pdf.py
@@ -7,7 +7,7 @@
 def test_uncorrelated_background():
     """tester for uncorrelated background model"""
 
-    pdf_wrapper = spey.get_backend("default_pdf.uncorrelated_background")
+    pdf_wrapper = spey.get_backend("default.uncorrelated_background")
 
     data = [36, 33]
     signal_yields = [12.0, 15.0]
@@ -32,7 +32,7 @@ def test_uncorrelated_background():
 def test_correlated_background():
     """tester for correlated background"""
 
-    pdf_wrapper = spey.get_backend("default_pdf.correlated_background")
+    pdf_wrapper = spey.get_backend("default.correlated_background")
     statistical_model = pdf_wrapper(
         signal_yields=[12.0, 15.0],
         background_yields=[50.0, 48.0],
@@ -54,7 +54,7 @@ def test_correlated_background():
 def test_third_moment():
     """tester for the third moment"""
 
-    pdf_wrapper = spey.get_backend("default_pdf.third_moment_expansion")
+    pdf_wrapper = spey.get_backend("default.third_moment_expansion")
     statistical_model = pdf_wrapper(
         signal_yields=[12.0, 15.0],
         background_yields=[50.0, 48.0],
@@ -79,7 +79,7 @@ def test_third_moment():
 def test_effective_sigma():
     """tester for the effective sigma"""
 
-    pdf_wrapper = spey.get_backend("default_pdf.effective_sigma")
+    pdf_wrapper = spey.get_backend("default.effective_sigma")
     statistical_model = pdf_wrapper(
         signal_yields=[12.0, 15.0],
         background_yields=[50.0, 48.0],
@@ -104,7 +104,7 @@ def test_effective_sigma():
 def test_poisson():
     """tester for uncorrelated background model"""
 
-    pdf_wrapper = spey.get_backend("default_pdf.poisson")
+    pdf_wrapper = spey.get_backend("default.poisson")
 
     data = [36, 33]
     signal_yields = [12.0, 15.0]
@@ -128,7 +128,7 @@ def test_poisson():
 def test_normal():
     """tester for gaussian model"""
 
-    statistical_model = spey.get_backend("default_pdf.normal")(
+    statistical_model = spey.get_backend("default.normal")(
         signal_yields=[12.0],
         background_yields=[50.0],
         data=[36],
@@ -153,7 +153,7 @@ def test_multivariate_gauss():
     data = np.array([36, 33])
     cov = np.array([[144.0, 13.0], [25.0, 256.0]])
 
-    statistical_model = spey.get_backend("default_pdf.multivariate_normal")(
+    statistical_model = spey.get_backend("default.multivariate_normal")(
         signal_yields=signal,
         background_yields=bkg,
         data=data,
diff --git a/tests/test_math.py b/tests/test_math.py
index 9715c17..9fcbde7 100644
--- a/tests/test_math.py
+++ b/tests/test_math.py
@@ -8,7 +8,7 @@
 def test_uncorrelated_background():
     """tester for uncorrelated background model"""
 
-    pdf_wrapper = spey.get_backend("default_pdf.uncorrelated_background")
+    pdf_wrapper = spey.get_backend("default.uncorrelated_background")
 
     data = [36, 33]
     signal_yields = [12.0, 15.0]
diff --git a/tests/test_simplifiedllhd_limit.py b/tests/test_simplifiedllhd_limit.py
index e911cbc..26311f2 100644
--- a/tests/test_simplifiedllhd_limit.py
+++ b/tests/test_simplifiedllhd_limit.py
@@ -9,7 +9,7 @@
 def test_simplified_llhds_at_limit():
     """Test models at the limit"""
 
-    base_model = spey.get_backend("default_pdf.uncorrelated_background")(
+    base_model = spey.get_backend("default.uncorrelated_background")(
         signal_yields=[1, 1, 1],
         background_yields=[2, 1, 3],
         data=[2, 2, 2],
@@ -17,7 +17,7 @@ def test_simplified_llhds_at_limit():
     )
     base_model_cls = base_model.exclusion_confidence_level()[0]
 
-    correlated_model = spey.get_backend("default_pdf.correlated_background")(
+    correlated_model = spey.get_backend("default.correlated_background")(
         signal_yields=[1, 1, 1],
         background_yields=[2, 1, 3],
         data=[2, 2, 2],
@@ -29,7 +29,7 @@ def test_simplified_llhds_at_limit():
         correlated_model_cls, base_model_cls
     ), "Correlated model is not same as base model"
 
-    third_moment_model = spey.get_backend("default_pdf.third_moment_expansion")(
+    third_moment_model = spey.get_backend("default.third_moment_expansion")(
         signal_yields=[1, 1, 1],
         background_yields=[2, 1, 3],
         data=[2, 2, 2],
@@ -42,7 +42,7 @@ def test_simplified_llhds_at_limit():
         third_moment_model_cls, base_model_cls
     ), "third moment model is not same as base model"
 
-    eff_sigma_model = spey.get_backend("default_pdf.effective_sigma")(
+    eff_sigma_model = spey.get_backend("default.effective_sigma")(
         signal_yields=[1, 1, 1],
         background_yields=[2, 1, 3],
         data=[2, 2, 2],

From ecc65ed40b052b7fcd727fa269be6893075cd2f8 Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Wed, 29 Jan 2025 11:14:36 -0500
Subject: [PATCH 07/32] include latest version of python 3.12

---
 .github/workflows/tests.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index 057cfab..fdfec03 100644
--- a/.github/workflows/tests.yml
+++ b/.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.6"]
+        python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"]
         os: [ubuntu-latest, macos-latest]
 
     steps:

From 1230a427ee4ff8d493cb5972fbf5c0250bff787c Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Wed, 29 Jan 2025 16:55:35 -0500
Subject: [PATCH 08/32] change usage to normal dist instead of clipped normal

---
 src/spey/hypothesis_testing/asymptotic_calculator.py | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/src/spey/hypothesis_testing/asymptotic_calculator.py b/src/spey/hypothesis_testing/asymptotic_calculator.py
index 7e2ec90..18661e1 100644
--- a/src/spey/hypothesis_testing/asymptotic_calculator.py
+++ b/src/spey/hypothesis_testing/asymptotic_calculator.py
@@ -1,6 +1,6 @@
 """Tools for computing confidence level and pvalues at the asymptotic limit"""
 
-from typing import List, Text, Tuple
+from typing import List, Tuple
 
 from .distributions import AsymptoticTestStatisticsDistribution
 from .utils import expected_pvalues, pvalues
@@ -72,11 +72,12 @@ def compute_asymptotic_confidence_level(
         :func:`~spey.hypothesis_testing.utils.pvalues`,
         :func:`~spey.hypothesis_testing.utils.expected_pvalues`
     """
-    cutoff = -sqrt_qmuA  # use clipped normal -> normal will mean -np.inf
-    # gives more stable result for cases that \hat\mu > \mu : see eq 14, 16 :xref:`1007.1727`
+    # use normal distribution for the cutoff -> default
+    # to use clipped normal change cutoff to -sqrt_qmuA this may give more stable results
+    # for cases that \hat\mu > \mu : see eq 14, 16 :xref:`1007.1727`
 
-    sig_plus_bkg_distribution = AsymptoticTestStatisticsDistribution(-sqrt_qmuA, cutoff)
-    bkg_only_distribution = AsymptoticTestStatisticsDistribution(0.0, cutoff)
+    sig_plus_bkg_distribution = AsymptoticTestStatisticsDistribution(-sqrt_qmuA)
+    bkg_only_distribution = AsymptoticTestStatisticsDistribution(0.0)
 
     CLsb_obs, CLb_obs, CLs_obs = pvalues(
         delta_test_statistic, sig_plus_bkg_distribution, bkg_only_distribution

From ed4076d6488928a101b79dc5215a6009564ae905 Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Fri, 31 Jan 2025 13:46:32 -0500
Subject: [PATCH 09/32] update tests

---
 tests/test_cls_computer.py |  4 ++--
 tests/test_default_pdf.py  | 25 +++++++++++++++++++++----
 2 files changed, 23 insertions(+), 6 deletions(-)

diff --git a/tests/test_cls_computer.py b/tests/test_cls_computer.py
index e13f4f4..a39d069 100644
--- a/tests/test_cls_computer.py
+++ b/tests/test_cls_computer.py
@@ -18,7 +18,7 @@ def test_compute_asymptotic_confidence_level():
         0.14339349869880672,
         0.31731050786291404,
         0.5942867086725301,
-        0.5942867086725301,
+        0.8609310408460745,
     ], "expected CLs value is wrong."
 
     a, b = compute_asymptotic_confidence_level(1.0, 2.0, "q0")
@@ -29,7 +29,7 @@ def test_compute_asymptotic_confidence_level():
         0.02275013194817923,
         0.15865525393145702,
         0.5,
-        0.5,
+        0.841344746068543,
     ], "expected CLs value is wrong."
 
 
diff --git a/tests/test_default_pdf.py b/tests/test_default_pdf.py
index 2167d02..a60ef0d 100644
--- a/tests/test_default_pdf.py
+++ b/tests/test_default_pdf.py
@@ -1,6 +1,9 @@
 """Test default pdf plugin"""
 
 import numpy as np
+from scipy.optimize import minimize_scalar
+from scipy.stats import norm, multivariate_normal
+
 import spey
 
 
@@ -128,12 +131,18 @@ def test_poisson():
 def test_normal():
     """tester for gaussian model"""
 
+    sig, bkg, obs, unc = 12.0, 50.0, 36.0, 20.0
+
+    dist = norm(loc=obs, scale=unc)
+    opt = minimize_scalar(lambda x: -dist.logpdf(x * sig + bkg), bounds=(-2.0, 0.0))
+
     statistical_model = spey.get_backend("default.normal")(
-        signal_yields=[12.0],
-        background_yields=[50.0],
-        data=[36],
-        absolute_uncertainties=[20.0],
+        signal_yields=[sig],
+        background_yields=[bkg],
+        data=[obs],
+        absolute_uncertainties=[unc],
     )
+    muhat, maxnll = statistical_model.maximize_likelihood()
 
     assert np.isclose(
         statistical_model.chi2(poi_test_denominator=0),
@@ -143,6 +152,8 @@ def test_normal():
             - (0.5 * ((50.0 - 36.0) ** 2 / 20.0**2))
         ),
     ), "Gaussian chi2 is wrong"
+    assert np.isclose(muhat, opt.x), "Normal:: Muhat is wrong"
+    assert np.isclose(maxnll, opt.fun), "Normal:: MLE is wrong"
 
 
 def test_multivariate_gauss():
@@ -160,6 +171,10 @@ def test_multivariate_gauss():
         covariance_matrix=cov,
     )
 
+    dist = multivariate_normal(mean=data, cov=cov)
+    opt = minimize_scalar(lambda x: -dist.logpdf(x * signal + bkg), bounds=(-2, 0))
+    muhat, maxnll = statistical_model.maximize_likelihood()
+
     assert np.isclose(
         statistical_model.chi2(poi_test_denominator=0),
         2.0
@@ -168,3 +183,5 @@ def test_multivariate_gauss():
             - (0.5 * (bkg - data) @ np.linalg.inv(cov) @ (bkg - data))
         ),
     ), "Multivariate gauss wrong"
+    assert np.isclose(muhat, opt.x, rtol=1e-3), "MultivariateNormal:: Muhat is wrong"
+    assert np.isclose(maxnll, opt.fun, rtol=1e-3), "MultivariateNormal:: MLE is wrong"

From d633c7d2226d97329ba8dda9ba5d43304fe334ce Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Fri, 31 Jan 2025 13:48:25 -0500
Subject: [PATCH 10/32] improve wrapper

---
 src/spey/interface/statistical_model.py | 40 +++++++++++--------------
 1 file changed, 18 insertions(+), 22 deletions(-)

diff --git a/src/spey/interface/statistical_model.py b/src/spey/interface/statistical_model.py
index 02d0ff5..fb48c39 100644
--- a/src/spey/interface/statistical_model.py
+++ b/src/spey/interface/statistical_model.py
@@ -1,6 +1,7 @@
 """Statistical Model wrapper class"""
 import logging
-from typing import Any, Callable, Dict, List, Optional, Text, Tuple, Union
+from functools import wraps
+from typing import Any, Callable, Dict, List, Optional, Tuple, Union
 
 import numpy as np
 
@@ -788,6 +789,7 @@ def statistical_model_wrapper(
         * **other keyword arguments**: Backend specific keyword inputs.
     """
 
+    @wraps(func)
     def wrapper(
         *args,
         analysis: str = "__unknown_analysis__",
@@ -795,11 +797,21 @@ def wrapper(
         ntoys: int = 1000,
         **kwargs,
     ) -> StatisticalModel:
-        """
-        Statistical Model Backend wrapper.
+        return StatisticalModel(
+            backend=func(*args, **kwargs),
+            analysis=analysis,
+            xsection=xsection,
+            ntoys=ntoys,
+        )
+
+    docstring = (
+        "\n\n"
+        + "<>" * 30
+        + "\n\n"
+        + """
+        Optional arguments:
 
         Args:
-            args: Backend specific arguments.
             analysis (``Text``, default ``"__unknown_analysis__"``): Unique identifier of the
               statistical model. This attribue will be used for book keeping purposes.
             xsection (``float``, default ``nan``): cross section, unit is determined by the
@@ -807,7 +819,6 @@ def wrapper(
               cross-section value.
             ntoys (``int``, default ``1000``): Number of toy samples for hypothesis testing.
               (Only used for toy-based hypothesis testing)
-            kwargs: Backend specific keyword inputs.
 
         Raises:
             :obj:`AssertionError`: If the model input does not inherit :class:`~spey.DataBase`.
@@ -816,23 +827,8 @@ def wrapper(
             ~spey.StatisticalModel:
             Backend wraped with statistical model interface.
         """
-        return StatisticalModel(
-            backend=func(*args, **kwargs),
-            analysis=analysis,
-            xsection=xsection,
-            ntoys=ntoys,
-        )
-
-    docstring = getattr(func, "__doc__", "no docstring available")
-    if docstring is None:
-        docstring = "Documentation is not available..."
-
-    wrapper.__doc__ += (
-        "\n\t"
-        + "<>" * 30
-        + "\n\n\t Current statistical model backend properties:\n"
-        + docstring.replace("\n", "\n\t")
-        + "\n"
     )
 
+    wrapper.__doc__ += docstring
+
     return wrapper

From c379572695bf4d3208a09d15c1b150da9ed54709 Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Fri, 31 Jan 2025 15:20:00 -0500
Subject: [PATCH 11/32] test poisson

---
 tests/test_default_pdf.py | 17 +++++++++++++----
 1 file changed, 13 insertions(+), 4 deletions(-)

diff --git a/tests/test_default_pdf.py b/tests/test_default_pdf.py
index a60ef0d..ad51c54 100644
--- a/tests/test_default_pdf.py
+++ b/tests/test_default_pdf.py
@@ -2,7 +2,7 @@
 
 import numpy as np
 from scipy.optimize import minimize_scalar
-from scipy.stats import norm, multivariate_normal
+from scipy.stats import norm, multivariate_normal, poisson
 
 import spey
 
@@ -109,9 +109,9 @@ def test_poisson():
 
     pdf_wrapper = spey.get_backend("default.poisson")
 
-    data = [36, 33]
-    signal_yields = [12.0, 15.0]
-    background_yields = [50.0, 48.0]
+    data = np.array([36, 33])
+    signal_yields = np.array([12.0, 15.0])
+    background_yields = np.array([50.0, 48.0])
 
     stat_model = pdf_wrapper(
         signal_yields=signal_yields,
@@ -127,6 +127,15 @@ def test_poisson():
     ), "CLs is wrong"
     assert np.isclose(stat_model.sigma_mu(1.0), 0.5573350296644078), "Sigma mu is wrong"
 
+    opt = minimize_scalar(
+        lambda x: -sum(poisson.logpmf(data, x * signal_yields + background_yields)),
+        bounds=(-2, 0),
+    )
+    muhat, maxnll = stat_model.maximize_likelihood()
+
+    assert np.isclose(muhat, opt.x, rtol=1e-3), "Poisson:: Muhat is wrong"
+    assert np.isclose(maxnll, opt.fun), "Poisson:: MLE is wrong"
+
 
 def test_normal():
     """tester for gaussian model"""

From dfa9096d2ef769e34c6eba453232db0d0143f26a Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Fri, 31 Jan 2025 18:00:59 -0500
Subject: [PATCH 12/32] extend poiss tests

---
 tests/test_default_pdf.py | 51 ++++++++++++++++++++++++++++++---------
 1 file changed, 39 insertions(+), 12 deletions(-)

diff --git a/tests/test_default_pdf.py b/tests/test_default_pdf.py
index ad51c54..593c353 100644
--- a/tests/test_default_pdf.py
+++ b/tests/test_default_pdf.py
@@ -2,7 +2,7 @@
 
 import numpy as np
 from scipy.optimize import minimize_scalar
-from scipy.stats import norm, multivariate_normal, poisson
+from scipy.stats import multivariate_normal, norm, poisson
 
 import spey
 
@@ -109,9 +109,9 @@ def test_poisson():
 
     pdf_wrapper = spey.get_backend("default.poisson")
 
-    data = np.array([36, 33])
-    signal_yields = np.array([12.0, 15.0])
-    background_yields = np.array([50.0, 48.0])
+    data = np.array([55])
+    signal_yields = np.array([12.0])
+    background_yields = np.array([50.0])
 
     stat_model = pdf_wrapper(
         signal_yields=signal_yields,
@@ -121,21 +121,48 @@ def test_poisson():
         xsection=0.123,
     )
 
-    assert np.isclose(stat_model.poi_upper_limit(), 0.3140867496931846), "POI is wrong."
-    assert np.isclose(
-        stat_model.exclusion_confidence_level()[0], 0.9999807105228611
-    ), "CLs is wrong"
-    assert np.isclose(stat_model.sigma_mu(1.0), 0.5573350296644078), "Sigma mu is wrong"
-
     opt = minimize_scalar(
         lambda x: -sum(poisson.logpmf(data, x * signal_yields + background_yields)),
-        bounds=(-2, 0),
+        bounds=(-0.5, 1),
     )
     muhat, maxnll = stat_model.maximize_likelihood()
 
-    assert np.isclose(muhat, opt.x, rtol=1e-3), "Poisson:: Muhat is wrong"
+    assert np.isclose(
+        muhat, opt.x, rtol=1e-3
+    ), f"Poisson:: Muhat is wrong {muhat} != {opt.x}"
     assert np.isclose(maxnll, opt.fun), "Poisson:: MLE is wrong"
 
+    ## Test Exclusion limit
+    def logprob(mu, data):
+        return poisson.logpmf(data, mu * signal_yields + background_yields)
+
+    opt = minimize_scalar(
+        lambda x: -sum(poisson.logpmf(data, x * signal_yields + background_yields)),
+        bounds=(0, 1),
+    )
+    tmu = 2 * (-logprob(1, data) - opt.fun)
+
+    opt = minimize_scalar(
+        lambda x: -sum(
+            poisson.logpmf(background_yields, x * signal_yields + background_yields)
+        ),
+        bounds=(-0.5, 0.5),
+    )
+    tmuA = 2 * (-logprob(1, background_yields) - opt.fun)
+
+    sqrt_qmuA = np.sqrt(tmuA)
+    sqrt_qmu = np.sqrt(tmu)
+    delta_teststat = sqrt_qmu - sqrt_qmuA
+
+    CLsb = norm.cdf(-sqrt_qmuA - delta_teststat)
+    CLb = norm.cdf(-delta_teststat)
+    print(CLsb, CLb)
+    CLs = (CLsb / CLb)[0]
+    st_cls = stat_model.exclusion_confidence_level()[0]
+    assert np.isclose(
+        1 - CLs, st_cls
+    ), f"Poisson:: Exclusion limit is wrong {1 - CLs} != {st_cls}"
+
 
 def test_normal():
     """tester for gaussian model"""

From ba63d6b1db6f7299f555c3b60fa65d2038eb0fdc Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Fri, 31 Jan 2025 18:03:10 -0500
Subject: [PATCH 13/32] remove redundant call

---
 tests/test_default_pdf.py | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/tests/test_default_pdf.py b/tests/test_default_pdf.py
index 593c353..41b8e78 100644
--- a/tests/test_default_pdf.py
+++ b/tests/test_default_pdf.py
@@ -136,10 +136,6 @@ def test_poisson():
     def logprob(mu, data):
         return poisson.logpmf(data, mu * signal_yields + background_yields)
 
-    opt = minimize_scalar(
-        lambda x: -sum(poisson.logpmf(data, x * signal_yields + background_yields)),
-        bounds=(0, 1),
-    )
     tmu = 2 * (-logprob(1, data) - opt.fun)
 
     opt = minimize_scalar(

From 3effc7cc55a832e76ed6e8588a60246c7b81c697 Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Fri, 31 Jan 2025 20:38:46 -0500
Subject: [PATCH 14/32] simplify

---
 tests/test_default_pdf.py | 35 +++++++++++------------------------
 1 file changed, 11 insertions(+), 24 deletions(-)

diff --git a/tests/test_default_pdf.py b/tests/test_default_pdf.py
index 41b8e78..1a79567 100644
--- a/tests/test_default_pdf.py
+++ b/tests/test_default_pdf.py
@@ -109,22 +109,18 @@ def test_poisson():
 
     pdf_wrapper = spey.get_backend("default.poisson")
 
-    data = np.array([55])
-    signal_yields = np.array([12.0])
-    background_yields = np.array([50.0])
+    data = 55
+    signal_yields = 12.0
+    background_yields = 50.0
 
     stat_model = pdf_wrapper(
-        signal_yields=signal_yields,
-        background_yields=background_yields,
-        data=data,
-        analysis="multi_bin",
-        xsection=0.123,
+        signal_yields=[signal_yields], background_yields=[background_yields], data=[data]
     )
 
-    opt = minimize_scalar(
-        lambda x: -sum(poisson.logpmf(data, x * signal_yields + background_yields)),
-        bounds=(-0.5, 1),
-    )
+    def logprob(mu, data):
+        return poisson.logpmf(data, mu * signal_yields + background_yields)
+
+    opt = minimize_scalar(lambda x: -logprob(x, data), bounds=(-0.5, 1))
     muhat, maxnll = stat_model.maximize_likelihood()
 
     assert np.isclose(
@@ -133,18 +129,10 @@ def test_poisson():
     assert np.isclose(maxnll, opt.fun), "Poisson:: MLE is wrong"
 
     ## Test Exclusion limit
-    def logprob(mu, data):
-        return poisson.logpmf(data, mu * signal_yields + background_yields)
-
     tmu = 2 * (-logprob(1, data) - opt.fun)
 
-    opt = minimize_scalar(
-        lambda x: -sum(
-            poisson.logpmf(background_yields, x * signal_yields + background_yields)
-        ),
-        bounds=(-0.5, 0.5),
-    )
-    tmuA = 2 * (-logprob(1, background_yields) - opt.fun)
+    optA = minimize_scalar(lambda x: -logprob(x, background_yields), bounds=(-0.5, 0.5))
+    tmuA = 2 * (-logprob(1, background_yields) - optA.fun)
 
     sqrt_qmuA = np.sqrt(tmuA)
     sqrt_qmu = np.sqrt(tmu)
@@ -152,8 +140,7 @@ def logprob(mu, data):
 
     CLsb = norm.cdf(-sqrt_qmuA - delta_teststat)
     CLb = norm.cdf(-delta_teststat)
-    print(CLsb, CLb)
-    CLs = (CLsb / CLb)[0]
+    CLs = CLsb / CLb
     st_cls = stat_model.exclusion_confidence_level()[0]
     assert np.isclose(
         1 - CLs, st_cls

From 8fa72181c38e036565d881feb12f1152936ec74d Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Fri, 31 Jan 2025 21:48:47 -0500
Subject: [PATCH 15/32] update for numeric stability

---
 src/spey/hypothesis_testing/utils.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/src/spey/hypothesis_testing/utils.py b/src/spey/hypothesis_testing/utils.py
index 923d31c..a6213f2 100644
--- a/src/spey/hypothesis_testing/utils.py
+++ b/src/spey/hypothesis_testing/utils.py
@@ -60,6 +60,12 @@ def pvalues(
         :func:`~spey.hypothesis_testing.asymptotic_calculator.compute_asymptotic_confidence_level`,
         :func:`~spey.hypothesis_testing.utils.expected_pvalues`
     """
+    if isinstance(sig_plus_bkg_distribution, AsymptoticTestStatisticsDistribution):
+        logp_sb = sig_plus_bkg_distribution.logcdf(delta_test_statistic)
+        logp_b = bkg_only_distribution.logcdf(delta_test_statistic)
+        p_s = np.exp(logp_sb - logp_b)
+        return np.exp(logp_sb), np.exp(logp_b), p_s
+
     CLsb = sig_plus_bkg_distribution.pvalue(delta_test_statistic)
     CLb = bkg_only_distribution.pvalue(delta_test_statistic)
     with warnings.catch_warnings(record=True):

From 79675c1a07d9703a0e2e4ac4855d316dc7613b7d Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Fri, 31 Jan 2025 21:48:58 -0500
Subject: [PATCH 16/32] add logcdf

---
 src/spey/hypothesis_testing/distributions.py | 22 +++++++++++++++-----
 1 file changed, 17 insertions(+), 5 deletions(-)

diff --git a/src/spey/hypothesis_testing/distributions.py b/src/spey/hypothesis_testing/distributions.py
index 1594f56..f85f27c 100644
--- a/src/spey/hypothesis_testing/distributions.py
+++ b/src/spey/hypothesis_testing/distributions.py
@@ -5,7 +5,7 @@
 """
 
 import numpy as np
-from scipy.stats import multivariate_normal
+from scipy.stats import norm
 
 __all__ = ["AsymptoticTestStatisticsDistribution", "EmpricTestStatisticsDistribution"]
 
@@ -39,7 +39,21 @@ def pvalue(self, value: float) -> float:
             ``float``:
             p-value
         """
-        return_value = multivariate_normal.cdf(self.shift - value)
+        return_value = norm.cdf(float(self.shift - value))
+        return return_value if return_value >= self.cutoff else np.nan
+
+    def logcdf(self, value: float) -> float:
+        r"""
+        Log of cumulative distribution function
+
+        Args:
+            value (``float``): The test statistic value.
+
+        Returns:
+            ``float``:
+            log-cdf
+        """
+        return_value = norm.logcdf(float(self.shift - value))
         return return_value if return_value >= self.cutoff else np.nan
 
     def expected_value(self, nsigma: float) -> float:
@@ -95,6 +109,4 @@ def expected_value(self, nsigma: float) -> float:
             ``float``:
             expected value of test statistic.
         """
-        return np.percentile(
-            self.samples, multivariate_normal.cdf(nsigma) * 100.0, method="linear"
-        )
+        return np.percentile(self.samples, norm.cdf(nsigma) * 100.0, method="linear")

From fc905c956f2a684e3496fbd792a4c469f4126539 Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Fri, 31 Jan 2025 21:49:24 -0500
Subject: [PATCH 17/32] bugfix

---
 .../hypothesis_testing/test_statistics.py     |  4 +-
 tests/test_cls_computer.py                    | 40 +++++++++++--------
 2 files changed, 26 insertions(+), 18 deletions(-)

diff --git a/src/spey/hypothesis_testing/test_statistics.py b/src/spey/hypothesis_testing/test_statistics.py
index 2992419..90a4451 100644
--- a/src/spey/hypothesis_testing/test_statistics.py
+++ b/src/spey/hypothesis_testing/test_statistics.py
@@ -2,11 +2,11 @@
 
 import logging
 import warnings
-from typing import Callable, Text, Tuple
+from typing import Callable, Tuple
 
 import numpy as np
 
-from spey.system.exceptions import UnknownTestStatistics, AsimovTestStatZero
+from spey.system.exceptions import AsimovTestStatZero, UnknownTestStatistics
 
 __all__ = ["qmu", "qmu_tilde", "q0", "get_test_statistic", "compute_teststatistics"]
 
diff --git a/tests/test_cls_computer.py b/tests/test_cls_computer.py
index a39d069..95da95c 100644
--- a/tests/test_cls_computer.py
+++ b/tests/test_cls_computer.py
@@ -1,5 +1,7 @@
 """Test asymptotic calculator"""
 
+import numpy as np
+
 from spey.hypothesis_testing.asymptotic_calculator import (
     compute_asymptotic_confidence_level,
 )
@@ -11,26 +13,32 @@ def test_compute_asymptotic_confidence_level():
 
     a, b = compute_asymptotic_confidence_level(1.0, 2.0, test_stat="qtilde")
 
-    assert a == [0.05933583307142766], "Observed CLs value is wrong."
+    assert np.isclose(a, 0.05933583307142766), "Observed CLs value is wrong."
 
-    assert b == [
-        0.05933583307142766,
-        0.14339349869880672,
-        0.31731050786291404,
-        0.5942867086725301,
-        0.8609310408460745,
-    ], "expected CLs value is wrong."
+    assert np.isclose(
+        b,
+        [
+            0.05933583307142766,
+            0.14339349869880672,
+            0.31731050786291404,
+            0.5942867086725301,
+            0.8609310408460745,
+        ],
+    ).all(), "expected CLs value is wrong."
 
     a, b = compute_asymptotic_confidence_level(1.0, 2.0, "q0")
 
-    assert a == [0.001349898031630116], "Observed CLs value is wrong."
-    assert b == [
-        0.001349898031630116,
-        0.02275013194817923,
-        0.15865525393145702,
-        0.5,
-        0.841344746068543,
-    ], "expected CLs value is wrong."
+    assert np.isclose(a, 0.001349898031630116), "Observed CLs value is wrong."
+    assert np.isclose(
+        b,
+        [
+            0.001349898031630116,
+            0.02275013194817923,
+            0.15865525393145702,
+            0.5,
+            0.841344746068543,
+        ],
+    ).all(), "expected CLs value is wrong."
 
 
 def test_compute_toy_confidence_level():

From cbfc5155e68a6fec6e7596cf96a43c93812f5718 Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Mon, 3 Feb 2025 16:12:04 -0500
Subject: [PATCH 18/32] add analytic test

---
 tests/test_default_pdf.py | 19 +++++++++++++++++++
 1 file changed, 19 insertions(+)

diff --git a/tests/test_default_pdf.py b/tests/test_default_pdf.py
index 1a79567..ccec5c9 100644
--- a/tests/test_default_pdf.py
+++ b/tests/test_default_pdf.py
@@ -146,6 +146,25 @@ def logprob(mu, data):
         1 - CLs, st_cls
     ), f"Poisson:: Exclusion limit is wrong {1 - CLs} != {st_cls}"
 
+    # compare with analytic results
+    s_b = signal_yields + background_yields
+    sqrt_qmu = np.sqrt(-2 * (-s_b + data * np.log(s_b) - data * np.log(data) + data))
+    sqrt_qmuA = np.sqrt(
+        -2
+        * (
+            -signal_yields
+            + background_yields * np.log(1 + signal_yields / background_yields)
+        )
+    )
+    delta_teststat = sqrt_qmu - sqrt_qmuA
+
+    logp_sb = norm.logcdf(-sqrt_qmu)
+    logp_b = norm.logcdf(-sqrt_qmu + sqrt_qmuA)
+    CLs_analytic = 1 - np.exp(logp_sb - logp_b)
+    assert np.isclose(
+        CLs_analytic, st_cls
+    ), f"Poisson:: Analytic exclusion limit is wrong {CLs_analytic} != {st_cls}"
+
 
 def test_normal():
     """tester for gaussian model"""

From d5d935f41fbe5ecedf7419497b3538af88b9e18b Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Mon, 3 Feb 2025 16:12:20 -0500
Subject: [PATCH 19/32] add failsafe

---
 src/spey/hypothesis_testing/upper_limits.py | 37 ++++++++++++---------
 1 file changed, 21 insertions(+), 16 deletions(-)

diff --git a/src/spey/hypothesis_testing/upper_limits.py b/src/spey/hypothesis_testing/upper_limits.py
index dd3d751..ff40710 100644
--- a/src/spey/hypothesis_testing/upper_limits.py
+++ b/src/spey/hypothesis_testing/upper_limits.py
@@ -3,12 +3,13 @@
 import logging
 import warnings
 from functools import partial
-from typing import Callable, List, Tuple, Union, Literal
+from typing import Callable, List, Literal, Tuple, Union
 
 import numpy as np
 import scipy
 
 from spey.hypothesis_testing.test_statistics import compute_teststatistics
+from spey.system.exceptions import AsimovTestStatZero
 from spey.utils import ExpectationType
 
 from .asymptotic_calculator import compute_asymptotic_confidence_level
@@ -182,22 +183,26 @@ def find_poi_upper_limit(
 
     def computer(poi_test: float, pvalue_idx: int) -> float:
         """Compute 1 - CLs(POI) = `confidence_level`"""
-        _, sqrt_qmuA, delta_teststat = compute_teststatistics(
-            poi_test,
-            maximum_likelihood,
-            logpdf,
-            maximum_asimov_likelihood,
-            asimov_logpdf,
-            test_stat,
-        )
-        pvalue = list(
-            map(
-                lambda x: 1.0 - x,
-                compute_asymptotic_confidence_level(sqrt_qmuA, delta_teststat, test_stat)[
-                    0 if expected == ExpectationType.observed else 1
-                ],
+        try:
+            _, sqrt_qmuA, delta_teststat = compute_teststatistics(
+                poi_test,
+                maximum_likelihood,
+                logpdf,
+                maximum_asimov_likelihood,
+                asimov_logpdf,
+                test_stat,
             )
-        )
+            pvalue = list(
+                map(
+                    lambda x: 1.0 - x,
+                    compute_asymptotic_confidence_level(
+                        sqrt_qmuA, delta_teststat, test_stat
+                    )[0 if expected == ExpectationType.observed else 1],
+                )
+            )
+        except AsimovTestStatZero as err:
+            log.debug(err)
+            pvalue = [0.0] if expected == ExpectationType.observed else [0.0] * 5
         return pvalue[pvalue_idx] - confidence_level
 
     result = []

From 0e23a27f39127a0a58bd9816c1ef7f46a8d2122a Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Tue, 4 Feb 2025 15:53:33 -0500
Subject: [PATCH 20/32] generalise POI UL function

---
 src/spey/base/hypotest_base.py              |  6 ++---
 src/spey/hypothesis_testing/upper_limits.py | 27 +++++++++++++++------
 2 files changed, 23 insertions(+), 10 deletions(-)

diff --git a/src/spey/base/hypotest_base.py b/src/spey/base/hypotest_base.py
index 5fbac08..9645280 100644
--- a/src/spey/base/hypotest_base.py
+++ b/src/spey/base/hypotest_base.py
@@ -944,8 +944,8 @@ def poi_upper_limit(
                 if not np.isclose(muhat, 0.0)
                 else 1.0
             )
-            low_init = low_init or muhat + 1.5 * sigma_mu
-            hig_init = hig_init or muhat + 2.5 * sigma_mu
+            low_init = np.clip(low_init or muhat + 1.5 * sigma_mu, 1e-10, None)
+            hig_init = np.clip(hig_init or muhat + 2.5 * sigma_mu, 1e-10, None)
             log.debug(f"new low_init = {low_init}, new hig_init = {hig_init}")
 
         return find_poi_upper_limit(
@@ -955,7 +955,7 @@ def poi_upper_limit(
             asimov_logpdf=logpdf_asimov,
             expected=expected,
             confidence_level=confidence_level,
-            allow_negative_signal=allow_negative_signal,
+            test_stat="q" if allow_negative_signal else "qtilde",
             low_init=low_init,
             hig_init=hig_init,
             expected_pvalue=expected_pvalue,
diff --git a/src/spey/hypothesis_testing/upper_limits.py b/src/spey/hypothesis_testing/upper_limits.py
index ff40710..ba2ca21 100644
--- a/src/spey/hypothesis_testing/upper_limits.py
+++ b/src/spey/hypothesis_testing/upper_limits.py
@@ -112,9 +112,11 @@ def find_poi_upper_limit(
     asimov_logpdf: Callable[[float], float],
     expected: ExpectationType,
     confidence_level: float = 0.95,
-    allow_negative_signal: bool = True,
+    test_stat: str = "qtilde",
     low_init: float = 1.0,
     hig_init: float = 1.0,
+    hig_bound: float = 1e5,
+    low_bound: float = 1e-10,
     expected_pvalue: Literal["nominal", "1sigma", "2sigma"] = "nominal",
     maxiter: int = 10000,
 ) -> Union[float, List[float]]:
@@ -179,7 +181,6 @@ def find_poi_upper_limit(
     ], f"Unknown pvalue range {expected_pvalue}"
     if expected is ExpectationType.observed:
         expected_pvalue = "nominal"
-    test_stat = "q" if allow_negative_signal else "qtilde"
 
     def computer(poi_test: float, pvalue_idx: int) -> float:
         """Compute 1 - CLs(POI) = `confidence_level`"""
@@ -215,11 +216,23 @@ def computer(poi_test: float, pvalue_idx: int) -> float:
         log.debug(f"Running for p-value idx: {pvalue_idx}")
         comp = partial(computer, pvalue_idx=pvalue_idx)
         # Set an upper bound for the computation
-        hig_bound = 1e5
-        low, hig = find_root_limits(
-            comp, loc=0.0, low_ini=low_init, hig_ini=hig_init, hig_bound=hig_bound
-        )
-        log.debug(f"low: {low[-1]}, hig: {hig[-1]}")
+        hi_lo = find_root_limits(
+            comp,
+            loc=0.0,
+            low_ini=low_init,
+            hig_ini=hig_init,
+            hig_bound=hig_bound,
+            low_bound=low_bound,
+        )  # low, hig
+        if all(x > 0 for x in [low_init, hig_init, hig_bound, low_bound]):
+            low, hig = hi_lo
+            log.debug(f"low: {low[-1]}, hig: {hig[-1]}")
+        else:
+            # flip the order if the search is on the negative POI side.
+            hig, low = hi_lo
+            log.debug(
+                f"Searching for the UL on the left tail, low: {low[-1]}, hig: {hig[-1]}"
+            )
 
         # Check if its possible to find roots
         if np.sign(low[-1]) * np.sign(hig[-1]) > 0.0:

From 75dba0d59d0a5f0ab7be8b6b47bc19b4b1ed9dab Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Tue, 4 Feb 2025 15:54:30 -0500
Subject: [PATCH 21/32] bump the version

---
 .zenodo.json         | 6 +++---
 src/spey/_version.py | 2 +-
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/.zenodo.json b/.zenodo.json
index 7340888..7ec2982 100644
--- a/.zenodo.json
+++ b/.zenodo.json
@@ -1,8 +1,8 @@
 {
     "description": "smooth inference for reinterpretation studies",
     "license": "MIT",
-    "title": "SpeysideHEP/spey: v0.1.12",
-    "version": "v0.1.12",
+    "title": "SpeysideHEP/spey: v0.2.0",
+    "version": "v0.2.0",
     "upload_type": "software",
     "creators": [
         {
@@ -31,7 +31,7 @@
         },
         {
             "scheme": "url",
-            "identifier": "https://github.com/SpeysideHEP/spey/tree/v0.1.12",
+            "identifier": "https://github.com/SpeysideHEP/spey/tree/v0.2.0",
             "relation": "isSupplementTo"
         },
         {
diff --git a/src/spey/_version.py b/src/spey/_version.py
index 805b7d6..9426bdf 100644
--- a/src/spey/_version.py
+++ b/src/spey/_version.py
@@ -1,3 +1,3 @@
 """Version number (major.minor.patch[-label])"""
 
-__version__ = "0.1.12-beta"
+__version__ = "0.2.0-beta"

From e67928a560ad27596b238f13c52ca97f443856b0 Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Fri, 7 Feb 2025 11:18:26 -0500
Subject: [PATCH 22/32] first iteration on two sided limits

---
 src/spey/base/hypotest_base.py                | 98 +++++++++++++++++++
 .../hypothesis_testing/test_statistics.py     | 11 ++-
 2 files changed, 106 insertions(+), 3 deletions(-)

diff --git a/src/spey/base/hypotest_base.py b/src/spey/base/hypotest_base.py
index 9645280..f85d512 100644
--- a/src/spey/base/hypotest_base.py
+++ b/src/spey/base/hypotest_base.py
@@ -961,3 +961,101 @@ def poi_upper_limit(
             expected_pvalue=expected_pvalue,
             maxiter=maxiter,
         )
+
+    def two_sided_poi_upper_limit(
+        self,
+        expected: ExpectationType = ExpectationType.observed,
+        confidence_level: float = 0.95,
+        low_init: float = -1.0,
+        hig_init: float = -1.0,
+        maxiter: int = 10000,
+        optimiser_arguments_right_tail: Dict[str, Any] = None,
+        optimiser_arguments_left_tail: Dict[str, Any] = None,
+    ) -> List[float]:
+        r"""
+        [EXPERIMENTAL] Construct two sided POI upper limit.
+
+        .. attention::
+
+            This algorithm locates the POI upper limit on both positive and negative side.
+            Depending on the likelihood, `par_bounds` may need adjusting to be able to
+            find the correct point. `par_bounds` argument can be given through
+            `optimiser_arguments_<tail direction>` keyword argument.
+
+            Note that if there are intermediate roots, `par_bounds` needs to be further adjusted
+            to be able to extract those points.
+
+        Args:
+            expected (~spey.ExpectationType): Sets which values the fitting algorithm should
+              focus and p-values to be computed.
+
+              * :obj:`~spey.ExpectationType.observed`: Computes the p-values with via post-fit
+                  prescriotion which means that the experimental data will be assumed to be the truth
+              * :obj:`~spey.ExpectationType.aposteriori`: Computes the expected p-values with via
+                  post-fit prescriotion which means that the experimental data will be assumed to be
+                  the truth.
+              * :obj:`~spey.ExpectationType.apriori`: Computes the expected p-values with via pre-fit
+                  prescription which means that the SM will be assumed to be the truth.
+
+            confidence_level (``float``, default ``0.95``): Determines the confidence level
+              of the upper limit i.e. the value of :math:`1-CL_s`. It needs to be between ``[0,1]``.
+            low_init (``float``, default ``-1.0``): Lower limit for the search algorithm to start
+            hig_init (``float``, default ``-1.0``): Upper limit for the search algorithm to start
+            maxiter (``int``, default ``10000``): Maximum iteration limit for the optimiser.
+            optimiser_arguments_right_tail (``Dict[str, Any]``, default ``None``): Optimiser arguments
+              to find the roots on the right tail of the p-value distribution i.e. :math:`\mu>\hat\mu`
+            optimiser_arguments_left_tail (``Dict[str, Any]``, default ``None``): Optimiser arguments
+              to find the roots on the left tail of the p-value distribution i.e. :math:`\mu<\hat\mu`
+
+        Returns:
+            ``List[float]``:
+            lower and upper bounds for the POI
+
+        Raises:
+            AssertionError: If the `confidence_level` input is not between [0,1].
+        """
+        assert (
+            0.0 <= confidence_level <= 1.0
+        ), "Confidence level must be between zero and one."
+
+        # If the signal yields in all regions are zero then return inf.
+        # This means we are not able to set a bound with the given information.
+        if not self.is_alive:
+            return [np.inf, np.inf]
+
+        right_poi_ul = self.poi_upper_limit(
+            expected=expected,
+            confidence_level=confidence_level,
+            allow_negative_signal=True,
+            maxiter=maxiter,
+            optimiser_arguments=optimiser_arguments_right_tail,
+        )
+
+        optimiser_arguments = optimiser_arguments_left_tail or {}
+        (
+            maximum_likelihood,
+            logpdf,
+            maximum_asimov_likelihood,
+            logpdf_asimov,
+        ) = self._prepare_for_hypotest(
+            expected=expected,
+            test_statistics="qmu_left",
+            **optimiser_arguments,
+        )
+
+        left_poi_ul = find_poi_upper_limit(
+            maximum_likelihood=maximum_likelihood,
+            logpdf=logpdf,
+            maximum_asimov_likelihood=maximum_asimov_likelihood,
+            asimov_logpdf=logpdf_asimov,
+            expected=expected,
+            confidence_level=confidence_level,
+            test_stat="qmu_left",
+            low_init=low_init,
+            hig_init=hig_init,
+            hig_bound=-1e-10,
+            low_bound=-1e10,
+            maxiter=maxiter,
+        )
+
+        return [left_poi_ul, right_poi_ul]
diff --git a/src/spey/hypothesis_testing/test_statistics.py b/src/spey/hypothesis_testing/test_statistics.py
index 90a4451..3b613f7 100644
--- a/src/spey/hypothesis_testing/test_statistics.py
+++ b/src/spey/hypothesis_testing/test_statistics.py
@@ -83,6 +83,12 @@ def qmu(
     return 0.0 if muhat > mu else np.clip(-2.0 * (logpdf(mu) - max_logpdf), 0.0, None)
 
 
+def qmu_left(
+    mu: float, muhat: float, max_logpdf: float, logpdf: Callable[[float], float]
+) -> float:
+    return 0.0 if muhat < mu else np.clip(-2.0 * (logpdf(mu) - max_logpdf), 0.0, None)
+
+
 def q0(
     mu: float, muhat: float, max_logpdf: float, logpdf: Callable[[float], float]
 ) -> float:
@@ -112,8 +118,7 @@ def q0(
         ``float``:
         the value of :math:`q_{0}`.
     """
-    mu = 0.0
-    return 0.0 if muhat < 0.0 else np.clip(-2.0 * (logpdf(mu) - max_logpdf), 0.0, None)
+    return 0.0 if muhat < 0.0 else np.clip(-2.0 * (logpdf(0.0) - max_logpdf), 0.0, None)
 
 
 def get_test_statistic(
@@ -155,7 +160,7 @@ def get_test_statistic(
         test_stat = "q"
     elif test_stat in ["qtilde", "qmutilde"]:
         test_stat = "qmutilde"
-    options = {"qmutilde": qmu_tilde, "q": qmu, "q0": q0}
+    options = {"qmutilde": qmu_tilde, "q": qmu, "q0": q0, "qmu_left": qmu_left}
 
     if options.get(test_stat, False) is False:
         raise UnknownTestStatistics(

From 3d08109ca36190b163cfe7d3911c0aa6b93cf3fa Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Fri, 7 Feb 2025 11:43:54 -0500
Subject: [PATCH 23/32] adjust CLs

---
 src/spey/base/hypotest_base.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/src/spey/base/hypotest_base.py b/src/spey/base/hypotest_base.py
index f85d512..f34e4ac 100644
--- a/src/spey/base/hypotest_base.py
+++ b/src/spey/base/hypotest_base.py
@@ -1018,10 +1018,13 @@ def two_sided_poi_upper_limit(
             0.0 <= confidence_level <= 1.0
         ), "Confidence level must be between zero and one."
 
+        # use half of the CL due to two sided limit. Inner area should be confidence_level
+        confidence_level /= 2.0
+
         # If the signal yields in all regions are zero then return inf.
         # This means we are not able to set a bound with the given information.
-        if not self.is_alive:
-            return [np.inf, np.inf]
+        # if not self.is_alive:
+        #     return [np.inf, np.inf]
 
         right_poi_ul = self.poi_upper_limit(
             expected=expected,

From 7e5843d345aa9dbd5b413713507eaa06739e0491 Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Sat, 8 Feb 2025 10:31:35 -0500
Subject: [PATCH 24/32] first iteration of the `chi2_test` function

---
 src/spey/base/hypotest_base.py | 172 ++++++++++++++++++++++++++++++++-
 1 file changed, 168 insertions(+), 4 deletions(-)

diff --git a/src/spey/base/hypotest_base.py b/src/spey/base/hypotest_base.py
index f34e4ac..92071f7 100644
--- a/src/spey/base/hypotest_base.py
+++ b/src/spey/base/hypotest_base.py
@@ -10,6 +10,8 @@
 
 import numpy as np
 import tqdm
+from scipy.optimize import toms748
+from scipy.stats import chi2
 
 from spey.hypothesis_testing.asymptotic_calculator import (
     compute_asymptotic_confidence_level,
@@ -19,7 +21,7 @@
     get_test_statistic,
 )
 from spey.hypothesis_testing.toy_calculator import compute_toy_confidence_level
-from spey.hypothesis_testing.upper_limits import find_poi_upper_limit
+from spey.hypothesis_testing.upper_limits import find_poi_upper_limit, ComputerWrapper
 from spey.system.exceptions import (
     AsimovTestStatZero,
     CalculatorNotAvailable,
@@ -977,9 +979,9 @@ def two_sided_poi_upper_limit(
 
         .. attention::
 
-            This algorithm locates the POI upper limit on both positive and negative side.
-            Depending on the likelihood, `par_bounds` may need adjusting to be able to
-            find the correct point. `par_bounds` argument can be given through
+            This algorithm locates the POI upper limit for :math:`p_s` on both positive
+            and negative side. Depending on the likelihood, `par_bounds` may need adjusting
+            to be able to find the correct point. `par_bounds` argument can be given through
             `optimiser_arguments_<tail direction>` keyword argument.
 
             Note that if there are intermediate roots, `par_bounds` needs to be further adjusted
@@ -1062,3 +1064,165 @@ def two_sided_poi_upper_limit(
         )
 
         return [left_poi_ul, right_poi_ul]
+
+    def chi2_test(
+        self,
+        expected: ExpectationType = ExpectationType.observed,
+        confidence_level: float = 0.95,
+        limit_type: Literal["right", "left", "two-sided"] = "two-sided",
+    ) -> list[float]:
+        r"""
+        [Experimental] Find the POI value(s) that constraint :math:`\chi^2` distribution at a specific
+        confidence limit.
+
+        Args:
+            expected (~spey.ExpectationType): Sets which values the fitting algorithm should
+              focus and p-values to be computed.
+
+              * :obj:`~spey.ExpectationType.observed`: Computes the p-values with via post-fit
+                  prescriotion which means that the experimental data will be assumed to be the truth
+              * :obj:`~spey.ExpectationType.aposteriori`: Computes the expected p-values with via
+                  post-fit prescriotion which means that the experimental data will be assumed to be
+                  the truth.
+              * :obj:`~spey.ExpectationType.apriori`: Computes the expected p-values with via pre-fit
+                  prescription which means that the SM will be assumed to be the truth.
+
+            confidence_level (``float``, default ``0.95``): Determines the confidence level
+              of the upper limit. It needs to be between ``[0,1]``.
+            limit_type (``'right'``, ``'left'`` or ``'two-sided``, default ``"two-sided"``): choose
+              which side of the :math:`\chi^2` distribution should be constrained. Note that for
+              two-sided limits inner area of the :math:`\chi^2` distribution is set to `confidence_level`
+              thus the threshold becomes :math:`\alpha=(1-CL)/2` where CL corresponds to
+              `confidence_level`. For the left or right side alone :math:`\alpha=1-CL`.
+
+        Returns:
+            ``list[float]``:
+            POI value(s) that constraints the :math:`\chi^2` distribution at a given threshold.
+        """
+        assert (
+            0.0 <= confidence_level <= 1.0
+        ), "Confidence level must be between zero and one."
+        assert limit_type in [
+            "left",
+            "right",
+            "two-sided",
+        ], f"Invalid limit type: {limit_type}"
+
+        # Two sided test statistic need to be halfed, total area within
+        # two sides should be equal to confidence_level
+        alpha = (1.0 - confidence_level) * (0.5 if limit_type == "two-sided" else 1.0)
+
+        # DoF = # POI
+        chi2_threshold = chi2.isf(alpha, df=1)
+
+        muhat, mllhd = self.maximize_likelihood(
+            expected=expected, allow_negative_signal=limit_type in ["two-sided", "left"]
+        )
+
+        def computer(poi_test: float) -> float:
+            """Compute chi^2 - chi^2 threshold"""
+            llhd = self.likelihood(poi_test=poi_test, expected=expected)
+            return 2.0 * (llhd - mllhd) - chi2_threshold
+
+        try:
+            sigma_muhat = self.sigma_mu(muhat, expected=expected)
+        except:  # specify an exeption type
+            sigma_muhat = 1.0
+
+        results = []
+        if limit_type in ["left", "two-sided"]:
+            is_muhat_gt_0 = np.isclose(muhat, 0.0) or muhat > 0.0
+            low = -1.0 if is_muhat_gt_0 else muhat - 1.5 * sigma_muhat
+            hig = -1.0 if is_muhat_gt_0 else muhat - 2.5 * sigma_muhat
+            hig_bound, low_bound = -1e5, -1e-5
+
+            hig_computer = ComputerWrapper(computer)
+            while hig_computer(hig) < 0.0 and hig > hig_bound:
+                hig *= 2.0
+
+            low_computer = ComputerWrapper(computer)
+            while low_computer(low) > 0.0 and low < low_bound:
+                low *= 0.5
+
+            log.debug(
+                f"Left first attempt:: low: f({low:.5e})={low_computer[-1]:.5e},"
+                f" high: f({hig:.5e})={hig_computer[-1]:.5e}"
+            )
+            if np.sign(low_computer[-1]) == np.sign(hig_computer[-1]) and muhat > 0:
+                low_computer = ComputerWrapper(computer)
+                low = muhat
+                while low_computer(low) > 0.0 and low > 1e-5:
+                    low *= 0.5
+                log.debug(
+                    f"Left second attempt:: low: f({low:.5e})={low_computer[-1]:.5e},"
+                    f" high: f({hig:.5e})={hig_computer[-1]:.5e}"
+                )
+
+            if np.sign(low_computer[-1]) != np.sign(hig_computer[-1]):
+                x0, _ = toms748(
+                    computer,
+                    hig,
+                    low,
+                    k=2,
+                    xtol=2e-12,
+                    rtol=1e-4,
+                    full_output=True,
+                    maxiter=10000,
+                )
+                results.append(x0)
+            else:
+                log.error(
+                    "Can not find the roots on the left side."
+                    " Please check your chi^2 distribution, it might be too wide."
+                )
+                results.append(-1e5 if hig >= -1e5 else np.nan)
+
+        if limit_type in ["right", "two-sided"]:
+            is_muhat_le_0 = np.isclose(muhat, 0.0) or muhat < 0.0
+            low = 1.0 if is_muhat_le_0 else muhat + 1.5 * sigma_muhat
+            hig = 1.0 if is_muhat_le_0 else muhat + 2.5 * sigma_muhat
+            hig_bound, low_bound = 1e5, 1e-5
+
+            hig_computer = ComputerWrapper(computer)
+            while hig_computer(hig) < 0.0 and hig < hig_bound:
+                hig *= 2.0
+
+            low_computer = ComputerWrapper(computer)
+            while low_computer(low) > 0.0 and low > low_bound:
+                low *= 0.5
+
+            log.debug(
+                f"Right first attempt:: low: f({low:.5e})={low_computer[-1]:.5e},"
+                f" high: f({hig:.5e})={hig_computer[-1]:.5e}"
+            )
+
+            if np.sign(low_computer[-1]) == np.sign(hig_computer[-1]) and muhat < 0:
+                low_computer = ComputerWrapper(computer)
+                low = muhat
+                while low_computer(low) > 0.0 and low < -1e-5:
+                    low *= 0.5
+                log.debug(
+                    f"Left second attempt:: low: f({low:.5e})={low_computer[-1]:.5e},"
+                    f" high: f({hig:.5e})={hig_computer[-1]:.5e}"
+                )
+
+            if np.sign(low_computer[-1]) != np.sign(hig_computer[-1]):
+                x0, _ = toms748(
+                    computer,
+                    low,
+                    hig,
+                    k=2,
+                    xtol=2e-12,
+                    rtol=1e-4,
+                    full_output=True,
+                    maxiter=10000,
+                )
+                results.append(x0)
+            else:
+                log.error(
+                    "Can not find the roots on the right side."
+                    " Please check your chi^2 distribution, it might be too wide."
+                )
+                results.append(1e5 if hig >= 1e5 else np.nan)
+
+        return results

From 49e3039b6d6f375aac7dc191c8d4e9143cee61ca Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Sat, 8 Feb 2025 18:19:06 -0500
Subject: [PATCH 25/32] improve documentation

---
 src/spey/base/hypotest_base.py | 46 ++++++++++++++++++++--------------
 1 file changed, 27 insertions(+), 19 deletions(-)

diff --git a/src/spey/base/hypotest_base.py b/src/spey/base/hypotest_base.py
index 92071f7..aab784c 100644
--- a/src/spey/base/hypotest_base.py
+++ b/src/spey/base/hypotest_base.py
@@ -1072,32 +1072,40 @@ def chi2_test(
         limit_type: Literal["right", "left", "two-sided"] = "two-sided",
     ) -> list[float]:
         r"""
-        [Experimental] Find the POI value(s) that constraint :math:`\chi^2` distribution at a specific
-        confidence limit.
+        Determine the parameter of interest (POI) value(s) that constrain the
+        :math:`\chi^2` distribution at a specified confidence level.
+
+        .. versionadded:: 0.2.0
+
+        .. attention::
+
+            The degrees of freedom are set to one, referring to the POI. Currently, spey does not
+            support multiple POIs, but this feature is planned for future releases.
 
         Args:
-            expected (~spey.ExpectationType): Sets which values the fitting algorithm should
-              focus and p-values to be computed.
+            expected (~spey.ExpectationType): Specifies the type of expectation for the fitting
+              algorithm and p-value computation.
 
-              * :obj:`~spey.ExpectationType.observed`: Computes the p-values with via post-fit
-                  prescriotion which means that the experimental data will be assumed to be the truth
-              * :obj:`~spey.ExpectationType.aposteriori`: Computes the expected p-values with via
-                  post-fit prescriotion which means that the experimental data will be assumed to be
-                  the truth.
-              * :obj:`~spey.ExpectationType.apriori`: Computes the expected p-values with via pre-fit
-                  prescription which means that the SM will be assumed to be the truth.
+              * :obj:`~spey.ExpectationType.observed`: Computes p-values using post-fit prescription,
+                assuming experimental data as the truth.
+              * :obj:`~spey.ExpectationType.apriori`: Computes expected p-values using pre-fit
+                prescription, assuming the Standard Model (SM) as the truth.
 
-            confidence_level (``float``, default ``0.95``): Determines the confidence level
-              of the upper limit. It needs to be between ``[0,1]``.
-            limit_type (``'right'``, ``'left'`` or ``'two-sided``, default ``"two-sided"``): choose
-              which side of the :math:`\chi^2` distribution should be constrained. Note that for
-              two-sided limits inner area of the :math:`\chi^2` distribution is set to `confidence_level`
-              thus the threshold becomes :math:`\alpha=(1-CL)/2` where CL corresponds to
-              `confidence_level`. For the left or right side alone :math:`\alpha=1-CL`.
+            confidence_level (``float``, default ``0.95``): The confidence level for the upper limit.
+              Must be between 0 and 1.
+
+            limit_type (``'right'``, ``'left'`` or ``'two-sided'``, default ``"two-sided"``): Specifies
+              which side of the :math:`\chi^2` distribution should be constrained. For two-sided limits,
+              the inner area of the :math:`\chi^2` distribution is set to ``confidence_level``, making the
+              threshold :math:`\alpha=(1-CL)/2`, where CL is the `confidence_level`. For left or right
+              limits alone, :math:`\alpha=1-CL`. The :math:`\chi^2`-threshold is calculated using
+              ``chi2.isf(1.0 - confidence_level, df=1)`` for `left` and `right` limits, and
+              ``chi2.isf((1.0 - confidence_level)/2, df=1)`` for `two-sided` limits. Here, ``chi2``
+              refers to `SciPy's chi2 module <https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.chi2.html>`_.
 
         Returns:
             ``list[float]``:
-            POI value(s) that constraints the :math:`\chi^2` distribution at a given threshold.
+            POI value(s) that constrain the :math:`\chi^2` distribution at the given threshold.
         """
         assert (
             0.0 <= confidence_level <= 1.0

From 3ff81abd197d145b95c4d9c3b26c1eabf9f052b2 Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Sun, 9 Feb 2025 08:19:06 -0500
Subject: [PATCH 26/32] add chi2-test test

---
 src/spey/base/hypotest_base.py |  2 +-
 tests/test_default_pdf.py      | 13 ++++++++++++-
 2 files changed, 13 insertions(+), 2 deletions(-)

diff --git a/src/spey/base/hypotest_base.py b/src/spey/base/hypotest_base.py
index aab784c..7d9f0a2 100644
--- a/src/spey/base/hypotest_base.py
+++ b/src/spey/base/hypotest_base.py
@@ -1092,7 +1092,7 @@ def chi2_test(
                 prescription, assuming the Standard Model (SM) as the truth.
 
             confidence_level (``float``, default ``0.95``): The confidence level for the upper limit.
-              Must be between 0 and 1.
+              Must be between 0 and 1. This refers to the total inner area under the bell curve.
 
             limit_type (``'right'``, ``'left'`` or ``'two-sided'``, default ``"two-sided"``): Specifies
               which side of the :math:`\chi^2` distribution should be constrained. For two-sided limits,
diff --git a/tests/test_default_pdf.py b/tests/test_default_pdf.py
index ccec5c9..76aa7f7 100644
--- a/tests/test_default_pdf.py
+++ b/tests/test_default_pdf.py
@@ -2,7 +2,7 @@
 
 import numpy as np
 from scipy.optimize import minimize_scalar
-from scipy.stats import multivariate_normal, norm, poisson
+from scipy.stats import multivariate_normal, norm, poisson, chi2
 
 import spey
 
@@ -193,6 +193,17 @@ def test_normal():
     assert np.isclose(muhat, opt.x), "Normal:: Muhat is wrong"
     assert np.isclose(maxnll, opt.fun), "Normal:: MLE is wrong"
 
+    # add chi2-test
+
+    left_lim, right_lim = statistical_model.chi2_test()
+    left_chi2 = statistical_model.chi2(poi_test=left_lim, allow_negative_signal=True)
+    right_chi2 = statistical_model.chi2(poi_test=right_lim, allow_negative_signal=True)
+    chi2_threshold = chi2.isf((1 - 0.95) / 2, df=1)
+
+    assert np.isclose(
+        [left_chi2, right_chi2], chi2_threshold
+    ).all(), "chi2 test doesnt match chi2 threshold"
+
 
 def test_multivariate_gauss():
     """tester for multivar gauss"""

From 523c2ea9d500fdaf7b27f24402f696b62ee953b1 Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Thu, 20 Feb 2025 09:05:47 -0500
Subject: [PATCH 27/32] clean the code for the first stage of v0.2

---
 src/spey/__init__.py           |  19 ++++++-
 src/spey/base/hypotest_base.py | 101 ---------------------------------
 2 files changed, 16 insertions(+), 104 deletions(-)

diff --git a/src/spey/__init__.py b/src/spey/__init__.py
index e52da26..b835087 100644
--- a/src/spey/__init__.py
+++ b/src/spey/__init__.py
@@ -3,6 +3,7 @@
 import re
 import sys
 import textwrap
+import warnings
 from typing import Any, Callable, Dict, List, Literal, Optional, Text, Tuple, Union
 
 import numpy as np
@@ -175,11 +176,23 @@ def get_backend(name: str) -> Callable[[Any], StatisticalModel]:
         the backend it self which is in this particular example
         :obj:`~spey.backends.simplifiedlikelihood_backend.interface.SimplifiedLikelihoodInterface`.
     """
-    if "default_pdf" in name:
-        log.warning(
+    deprecated = [
+        "default_pdf.correlated_background",
+        "default_pdf.effective_sigma",
+        "default_pdf.multivariate_normal",
+        "default_pdf.normal",
+        "default_pdf.poisson",
+        "default_pdf.third_moment_expansion",
+        "default_pdf.uncorrelated_background",
+    ]
+    if name in deprecated:
+        # pylint: disable=logging-fstring-interpolation
+        message = (
             f"`{name}` has been deprecated, please use `{name.replace('default_pdf', 'default')}` instead."
-            " This will be removed in the future."
+            + " This will be removed in the future."
         )
+        warnings.warn(message, DeprecationWarning, stacklevel=2)
+        log.warning(message)
         name = name.replace("default_pdf", "default")
     backend = _backend_entries.get(name, False)
 
diff --git a/src/spey/base/hypotest_base.py b/src/spey/base/hypotest_base.py
index 7d9f0a2..b808d3f 100644
--- a/src/spey/base/hypotest_base.py
+++ b/src/spey/base/hypotest_base.py
@@ -964,107 +964,6 @@ def poi_upper_limit(
             maxiter=maxiter,
         )
 
-    def two_sided_poi_upper_limit(
-        self,
-        expected: ExpectationType = ExpectationType.observed,
-        confidence_level: float = 0.95,
-        low_init: float = -1.0,
-        hig_init: float = -1.0,
-        maxiter: int = 10000,
-        optimiser_arguments_right_tail: Dict[str, Any] = None,
-        optimiser_arguments_left_tail: Dict[str, Any] = None,
-    ) -> List[float]:
-        r"""
-        [EXPERIMENTAL] Construct two sided POI upper limit.
-
-        .. attention::
-
-            This algorithm locates the POI upper limit for :math:`p_s` on both positive
-            and negative side. Depending on the likelihood, `par_bounds` may need adjusting
-            to be able to find the correct point. `par_bounds` argument can be given through
-            `optimiser_arguments_<tail direction>` keyword argument.
-
-            Note that if there are intermediate roots, `par_bounds` needs to be further adjusted
-            to be able to extract those points.
-
-        Args:
-            expected (~spey.ExpectationType): Sets which values the fitting algorithm should
-              focus and p-values to be computed.
-
-              * :obj:`~spey.ExpectationType.observed`: Computes the p-values with via post-fit
-                  prescriotion which means that the experimental data will be assumed to be the truth
-              * :obj:`~spey.ExpectationType.aposteriori`: Computes the expected p-values with via
-                  post-fit prescriotion which means that the experimental data will be assumed to be
-                  the truth.
-              * :obj:`~spey.ExpectationType.apriori`: Computes the expected p-values with via pre-fit
-                  prescription which means that the SM will be assumed to be the truth.
-
-            confidence_level (``float``, default ``0.95``): Determines the confidence level
-              of the upper limit i.e. the value of :math:`1-CL_s`. It needs to be between ``[0,1]``.
-            low_init (``float``, default ``-1.0``): Lower limit for the search algorithm to start
-            hig_init (``float``, default ``-1.0``): Upper limit for the search algorithm to start
-            maxiter (``int``, default ``10000``): Maximum iteration limit for the optimiser.
-            optimiser_arguments_right_tail (``Dict[str, Any]``, default ``None``): Optimiser arguments
-              to find the roots on the right tail of the p-value distribution i.e. :math:`\mu>\hat\mu`
-            optimiser_arguments_left_tail (``Dict[str, Any]``, default ``None``): Optimiser arguments
-              to find the roots on the left tail of the p-value distribution i.e. :math:`\mu<\hat\mu`
-
-        Returns:
-            ``List[float]``:
-            lower and upper bounds for the POI
-
-        Raises:
-            AssertionError: If the `confidence_level` input is not between [0,1].
-        """
-        assert (
-            0.0 <= confidence_level <= 1.0
-        ), "Confidence level must be between zero and one."
-
-        # use half of the CL due to two sided limit. Inner area should be confidence_level
-        confidence_level /= 2.0
-
-        # If the signal yields in all regions are zero then return inf.
-        # This means we are not able to set a bound with the given information.
-        # if not self.is_alive:
-        #     return [np.inf, np.inf]
-
-        right_poi_ul = self.poi_upper_limit(
-            expected=expected,
-            confidence_level=confidence_level,
-            allow_negative_signal=True,
-            maxiter=maxiter,
-            optimiser_arguments=optimiser_arguments_right_tail,
-        )
-
-        optimiser_arguments = optimiser_arguments_left_tail or {}
-        (
-            maximum_likelihood,
-            logpdf,
-            maximum_asimov_likelihood,
-            logpdf_asimov,
-        ) = self._prepare_for_hypotest(
-            expected=expected,
-            test_statistics="qmu_left",
-            **optimiser_arguments,
-        )
-
-        left_poi_ul = find_poi_upper_limit(
-            maximum_likelihood=maximum_likelihood,
-            logpdf=logpdf,
-            maximum_asimov_likelihood=maximum_asimov_likelihood,
-            asimov_logpdf=logpdf_asimov,
-            expected=expected,
-            confidence_level=confidence_level,
-            test_stat="qmu_left",
-            low_init=low_init,
-            hig_init=hig_init,
-            hig_bound=-1e-10,
-            low_bound=-1e10,
-            maxiter=maxiter,
-        )
-
-        return [left_poi_ul, right_poi_ul]
-
     def chi2_test(
         self,
         expected: ExpectationType = ExpectationType.observed,

From 58267083a06982e184d288c657b521a0be0fc87c Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Thu, 20 Feb 2025 11:12:13 -0500
Subject: [PATCH 28/32] prepare CI/CD for deployment

---
 .github/workflows/check+build+deploy.yaml | 77 +++++++++++++++++++++++
 1 file changed, 77 insertions(+)
 create mode 100644 .github/workflows/check+build+deploy.yaml

diff --git a/.github/workflows/check+build+deploy.yaml b/.github/workflows/check+build+deploy.yaml
new file mode 100644
index 0000000..14819d8
--- /dev/null
+++ b/.github/workflows/check+build+deploy.yaml
@@ -0,0 +1,77 @@
+name: Check/Build/Deploy
+
+on:
+  release:
+    types: [published]
+  push:
+    branches: ["releases/**"]
+  workflow_dispatch:
+
+jobs:
+    check:
+        name: OS ${{ matrix.os }}, Python ${{ matrix.python-version }}
+        runs-on: ${{ matrix.os }}
+
+        strategy:
+          fail-fast: false
+          matrix:
+            python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"]
+            os: [ubuntu-latest, macos-latest]
+
+        steps:
+        - uses: actions/checkout@v3
+        - name: Set up Python ${{ matrix.python-version }}
+          uses: actions/setup-python@v3
+          with:
+            python-version: ${{ matrix.python-version }}
+        - name: Install dependencies
+          run: |
+            python -m pip install --upgrade pip
+            python -m pip install pytest pytest-cov
+            python -m pip install -e .
+        - name: Test with pytest
+          run: |
+            python -c 'import spey;spey.about()'
+            pytest --cov=spey tests/*py #--cov-fail-under 99
+
+    build:
+        name: Build wheel
+        runs-on: ubuntu-latest
+        needs: check
+        steps:
+          - name: Checkout git repository
+            uses: actions/checkout@v4
+
+          - name: Build Spey wheel
+            run: |
+              python -m pip install --upgrade pip
+              python -m pip install wheel
+              python setup.py bdist_wheel
+
+          - name: Upload wheels as artifacts
+            uses: actions/upload-artifact@v4
+            with:
+              name: wheel
+              path: ./dist/*.whl
+
+    deploy:
+        if:  ${{ !github.event.release.prerelease && !contains(github.ref, '-rc') && !contains(github.ref, '-beta') && !contains(github.ref, '-alpha') }}
+        runs-on: ubuntu-latest
+        needs: build
+        permissions:
+          id-token: write
+        steps:
+          - uses: actions/download-artifact@v4
+            with:
+              pattern: wheel
+              merge-multiple: true
+              path: dist
+
+          - name: Debug
+            run: ls -l dist
+
+          - uses: pypa/gh-action-pypi-publish@v1.12.4
+            with:
+              user: ${{ secrets.TWINE_USERNAME }}
+              password: ${{ secrets.TWINE_PSSWRD }}
+              repository-url: https://pypi.org/project/spey/

From 45e0dbc71c7adf68097af813152c3939a04a53be Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Thu, 20 Feb 2025 11:12:32 -0500
Subject: [PATCH 29/32] remove beta

---
 src/spey/_version.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/spey/_version.py b/src/spey/_version.py
index 9426bdf..3426ca1 100644
--- a/src/spey/_version.py
+++ b/src/spey/_version.py
@@ -1,3 +1,3 @@
 """Version number (major.minor.patch[-label])"""
 
-__version__ = "0.2.0-beta"
+__version__ = "0.2.0"

From 2c4200f4332da06543f55ae524002a366c72c447 Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Thu, 20 Feb 2025 11:17:03 -0500
Subject: [PATCH 30/32] update annotation for linux

---
 src/spey/base/hypotest_base.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/spey/base/hypotest_base.py b/src/spey/base/hypotest_base.py
index b808d3f..b5655fc 100644
--- a/src/spey/base/hypotest_base.py
+++ b/src/spey/base/hypotest_base.py
@@ -969,7 +969,7 @@ def chi2_test(
         expected: ExpectationType = ExpectationType.observed,
         confidence_level: float = 0.95,
         limit_type: Literal["right", "left", "two-sided"] = "two-sided",
-    ) -> list[float]:
+    ) -> List[float]:
         r"""
         Determine the parameter of interest (POI) value(s) that constrain the
         :math:`\chi^2` distribution at a specified confidence level.

From 85d6b3470db63bdba2688861b650baab17fee2e4 Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Thu, 20 Feb 2025 11:28:07 -0500
Subject: [PATCH 31/32] refine tests

---
 tests/test_default_pdf.py | 15 ++++++++-------
 1 file changed, 8 insertions(+), 7 deletions(-)

diff --git a/tests/test_default_pdf.py b/tests/test_default_pdf.py
index 76aa7f7..69a38ec 100644
--- a/tests/test_default_pdf.py
+++ b/tests/test_default_pdf.py
@@ -68,15 +68,16 @@ def test_third_moment():
         xsection=0.123,
     )
 
+    CLs = statistical_model.exclusion_confidence_level()[0]
+    assert np.isclose(CLs, 0.961432961), f"CLs is wrong, expected 0.961432961 got {CLs}"
+    poi_ul = statistical_model.poi_upper_limit()
     assert np.isclose(
-        statistical_model.exclusion_confidence_level()[0], 0.9614329616396733
-    ), "CLs is wrong"
-    assert np.isclose(
-        statistical_model.poi_upper_limit(), 0.9221339770245336
-    ), "POI is wrong"
+        poi_ul, 0.9221339
+    ), f"POI is wrong, expected 0.922133977 got {poi_ul}"
+    sigma_mu = statistical_model.sigma_mu(1.0)
     assert np.isclose(
-        statistical_model.sigma_mu(1.0), 0.854551194250324
-    ), "Sigma mu is wrong"
+        sigma_mu, 0.85455
+    ), f"Sigma mu is wrong, expected 0.8545511 got {sigma_mu}"
 
 
 def test_effective_sigma():

From c3b5d8ec57a8bf3176db5f4e26a4dd79c66697d3 Mon Sep 17 00:00:00 2001
From: "Jack Y. Araz" <jackaraz@gmail.com>
Date: Thu, 20 Feb 2025 11:31:48 -0500
Subject: [PATCH 32/32] refine tolerance

---
 tests/test_default_pdf.py | 12 +++++++-----
 1 file changed, 7 insertions(+), 5 deletions(-)

diff --git a/tests/test_default_pdf.py b/tests/test_default_pdf.py
index 69a38ec..bcd6894 100644
--- a/tests/test_default_pdf.py
+++ b/tests/test_default_pdf.py
@@ -69,15 +69,17 @@ def test_third_moment():
     )
 
     CLs = statistical_model.exclusion_confidence_level()[0]
-    assert np.isclose(CLs, 0.961432961), f"CLs is wrong, expected 0.961432961 got {CLs}"
+    assert np.isclose(
+        CLs, 0.961432961, rtol=1e-4
+    ), f"CLs is wrong, expected 0.961432961 got {CLs}"
     poi_ul = statistical_model.poi_upper_limit()
     assert np.isclose(
-        poi_ul, 0.9221339
-    ), f"POI is wrong, expected 0.922133977 got {poi_ul}"
+        poi_ul, 0.9221339, rtol=1e-4
+    ), f"POI is wrong, expected 0.9221339 got {poi_ul}"
     sigma_mu = statistical_model.sigma_mu(1.0)
     assert np.isclose(
-        sigma_mu, 0.85455
-    ), f"Sigma mu is wrong, expected 0.8545511 got {sigma_mu}"
+        sigma_mu, 0.85455, rtol=1e-4
+    ), f"Sigma mu is wrong, expected 0.85455 got {sigma_mu}"
 
 
 def test_effective_sigma():