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/
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:
diff --git a/.zenodo.json b/.zenodo.json
index 4d4fab6..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.11",
-    "version": "v0.1.11",
+    "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.11",
+            "identifier": "https://github.com/SpeysideHEP/spey/tree/v0.2.0",
             "relation": "isSupplementTo"
         },
         {
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/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}},
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/setup.py b/setup.py
index 929b3ef..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(
@@ -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",
diff --git a/src/spey/__init__.py b/src/spey/__init__.py
index ffafa8e..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
@@ -150,7 +151,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 +176,24 @@ 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`.
     """
+    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."
+        )
+        warnings.warn(message, DeprecationWarning, stacklevel=2)
+        log.warning(message)
+        name = name.replace("default_pdf", "default")
     backend = _backend_entries.get(name, False)
 
     if backend:
@@ -231,13 +250,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/_version.py b/src/spey/_version.py
index fa9697c..3426ca1 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.2.0"
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"""
diff --git a/src/spey/base/hypotest_base.py b/src/spey/base/hypotest_base.py
index 163f450..b5655fc 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,
@@ -840,42 +842,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 +882,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
@@ -942,8 +946,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(
@@ -953,9 +957,179 @@ def poi_upper_limit(
             asimov_logpdf=logpdf_asimov,
             expected=expected,
             confidence_level=confidence_level,
-            allow_negative_signal=False,
+            test_stat="q" if allow_negative_signal else "qtilde",
             low_init=low_init,
             hig_init=hig_init,
             expected_pvalue=expected_pvalue,
             maxiter=maxiter,
         )
+
+    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"""
+        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): 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.
+              * :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``): The confidence level for the upper limit.
+              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,
+              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 constrain the :math:`\chi^2` distribution at the 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
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/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
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")
diff --git a/src/spey/hypothesis_testing/test_statistics.py b/src/spey/hypothesis_testing/test_statistics.py
index 2992419..3b613f7 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"]
 
@@ -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(
diff --git a/src/spey/hypothesis_testing/upper_limits.py b/src/spey/hypothesis_testing/upper_limits.py
index db8073d..ba2ca21 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
@@ -111,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]]:
@@ -142,7 +145,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`
@@ -177,26 +181,29 @@ 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`"""
-        _, 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 = []
@@ -209,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:
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):
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
diff --git a/tests/test_cls_computer.py b/tests/test_cls_computer.py
index e13f4f4..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.5942867086725301,
-    ], "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.5,
-    ], "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():
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..bcd6894 100644
--- a/tests/test_default_pdf.py
+++ b/tests/test_default_pdf.py
@@ -1,13 +1,16 @@
 """Test default pdf plugin"""
 
 import numpy as np
+from scipy.optimize import minimize_scalar
+from scipy.stats import multivariate_normal, norm, poisson, chi2
+
 import spey
 
 
 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 +35,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 +57,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],
@@ -65,21 +68,24 @@ def test_third_moment():
         xsection=0.123,
     )
 
+    CLs = statistical_model.exclusion_confidence_level()[0]
     assert np.isclose(
-        statistical_model.exclusion_confidence_level()[0], 0.9614329616396733
-    ), "CLs is wrong"
+        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(
-        statistical_model.poi_upper_limit(), 0.9221339770245336
-    ), "POI is wrong"
+        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(
-        statistical_model.sigma_mu(1.0), 0.854551194250324
-    ), "Sigma mu is wrong"
+        sigma_mu, 0.85455, rtol=1e-4
+    ), f"Sigma mu is wrong, expected 0.85455 got {sigma_mu}"
 
 
 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,36 +110,80 @@ 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]
-    background_yields = [50.0, 48.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]
     )
 
-    assert np.isclose(stat_model.poi_upper_limit(), 0.3140867496931846), "POI is wrong."
+    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(
-        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"
+        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
+    tmu = 2 * (-logprob(1, data) - 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)
+    delta_teststat = sqrt_qmu - sqrt_qmuA
+
+    CLsb = norm.cdf(-sqrt_qmuA - delta_teststat)
+    CLb = norm.cdf(-delta_teststat)
+    CLs = CLsb / CLb
+    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}"
+
+    # 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"""
 
-    statistical_model = spey.get_backend("default_pdf.normal")(
-        signal_yields=[12.0],
-        background_yields=[50.0],
-        data=[36],
-        absolute_uncertainties=[20.0],
+    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=[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 +193,19 @@ 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"
+
+    # 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():
@@ -153,13 +216,17 @@ 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,
         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 +235,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"
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],