We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
Cash statistics : Inconsistency between doc/code or pvalue/siginificance
Gammpay 0.17, official

Hello
In the Cash count statistics description in the "Stats-Statistics" docs
https://docs.gammapy.org/0.17/stats/index.html#cash-counts-statistic

The text says that *the significance is the square root of TS ... which is equivalent to the Li&Ma significance for Li&Ma...*

The statement is true but this is not what is in ```stat.signifcance```.

As you know the Li&Ma usual formula gives a significance which is "single-sided", i.e. that corresponds to a ```1-p_value``` obtained by the integration of the normal distribution from ```-Inf``` to the significance value (that can be negative). 

I have checked on the example that from the p_value the significance should be:
```
from scipy.stats import norm, Chi2
p_value = 1 - chi2.cdf(stat.delta_ts) # same value found : 7.365232895800901
sig_LiMa = norm.ppf(1-p_value)
```
that gives ``` 2.475649754526797```

whereas the results computed in the notebook is :```2.71389625737626``` 
I checked that this corresponds to a ``` p_value``` shared on both sides of the normal function:
```norm.ppf(1-p_value/2)```

As a consequence, if I am correct:
- the text is misleading because it mentions Li&Ma, and the situation should be explicit;
- the significance attribute should be documented as the "double-sided" significance.

However I have checked that ```CashCountsStatitics``` can return negative significances which in principle is not possible with a "double-sided" integration. So there might be another inconsistency somewhere....

Thanks!

P.S.: one should never rely on siginiicances, only on p values.


--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `gammapy/stats/counts_statistic.py`
Content:
```
1 # Licensed under a 3-clause BSD style license - see LICENSE.rst
2 import abc
3 import numpy as np
4 from scipy.optimize import brentq, newton
5 from scipy.stats import chi2
6 from .fit_statistics import cash, get_wstat_mu_bkg, wstat
7 
8 __all__ = ["WStatCountsStatistic", "CashCountsStatistic"]
9 
10 
11 class CountsStatistic(abc.ABC):
12     @property
13     def ts(self):
14         """Return stat difference (TS) of measured excess versus no excess."""
15         # Remove (small) negative TS due to error in root finding
16         ts = np.clip(self.stat_null - self.stat_max, 0, None)
17         return ts
18 
19     @property
20     def sqrt_ts(self):
21         """Return statistical significance of measured excess."""
22         return np.sign(self.n_sig) * np.sqrt(self.ts)
23 
24     @property
25     def p_value(self):
26         """Return p_value of measured excess."""
27         return chi2.sf(self.ts, 1)
28 
29     def compute_errn(self, n_sigma=1.0):
30         """Compute downward excess uncertainties.
31 
32         Searches the signal value for which the test statistics is n_sigma**2 away from the maximum.
33 
34         Parameters
35         ----------
36         n_sigma : float
37             Confidence level of the uncertainty expressed in number of sigma. Default is 1.
38         """
39         errn = np.zeros_like(self.n_on, dtype="float")
40         min_range = self.n_sig - 2 * n_sigma * (self.error + 1)
41 
42         it = np.nditer(errn, flags=["multi_index"])
43         while not it.finished:
44             try:
45                 res = brentq(
46                     self._stat_fcn,
47                     min_range[it.multi_index],
48                     self.n_sig[it.multi_index],
49                     args=(self.stat_max[it.multi_index] + n_sigma ** 2, it.multi_index),
50                 )
51                 errn[it.multi_index] = res - self.n_sig[it.multi_index]
52             except ValueError:
53                 errn[it.multi_index] = -self.n_on[it.multi_index]
54             it.iternext()
55 
56         return errn
57 
58     def compute_errp(self, n_sigma=1):
59         """Compute upward excess uncertainties.
60 
61         Searches the signal value for which the test statistics is n_sigma**2 away from the maximum.
62 
63         Parameters
64         ----------
65         n_sigma : float
66             Confidence level of the uncertainty expressed in number of sigma. Default is 1.
67         """
68         errp = np.zeros_like(self.n_on, dtype="float")
69         max_range = self.n_sig + 2 * n_sigma * (self.error + 1)
70 
71         it = np.nditer(errp, flags=["multi_index"])
72         while not it.finished:
73             errp[it.multi_index] = brentq(
74                 self._stat_fcn,
75                 self.n_sig[it.multi_index],
76                 max_range[it.multi_index],
77                 args=(self.stat_max[it.multi_index] + n_sigma ** 2, it.multi_index),
78             )
79             it.iternext()
80 
81         return errp - self.n_sig
82 
83     def compute_upper_limit(self, n_sigma=3):
84         """Compute upper limit on the signal.
85 
86         Searches the signal value for which the test statistics is n_sigma**2 away from the maximum
87         or from 0 if the measured excess is negative.
88 
89         Parameters
90         ----------
91         n_sigma : float
92             Confidence level of the upper limit expressed in number of sigma. Default is 3.
93         """
94         ul = np.zeros_like(self.n_on, dtype="float")
95 
96         min_range = np.maximum(0, self.n_sig)
97         max_range = min_range + 2 * n_sigma * (self.error + 1)
98         it = np.nditer(ul, flags=["multi_index"])
99 
100         while not it.finished:
101             TS_ref = self._stat_fcn(min_range[it.multi_index], 0.0, it.multi_index)
102 
103             ul[it.multi_index] = brentq(
104                 self._stat_fcn,
105                 min_range[it.multi_index],
106                 max_range[it.multi_index],
107                 args=(TS_ref + n_sigma ** 2, it.multi_index),
108             )
109             it.iternext()
110 
111         return ul
112 
113     def n_sig_matching_significance(self, significance):
114         """Compute excess matching a given significance.
115 
116         This function is the inverse of `significance`.
117 
118         Parameters
119         ----------
120         significance : float
121             Significance
122 
123         Returns
124         -------
125         n_sig : `numpy.ndarray`
126             Excess
127         """
128         n_sig = np.zeros_like(self.n_bkg, dtype="float")
129         it = np.nditer(n_sig, flags=["multi_index"])
130 
131         while not it.finished:
132             try:
133                 n_sig[it.multi_index] = newton(
134                     self._n_sig_matching_significance_fcn,
135                     np.sqrt(self.n_bkg[it.multi_index]) * significance,
136                     args=(significance, it.multi_index),
137                 )
138             except:
139                 n_sig[it.multi_index] = np.nan
140 
141             it.iternext()
142         return n_sig
143 
144 
145 class CashCountsStatistic(CountsStatistic):
146     """Class to compute statistics (significance, asymmetric errors , ul) for Poisson distributed variable
147     with known background.
148 
149     Parameters
150     ----------
151     n_on : int
152         Measured counts
153     mu_bkg : float
154         Known level of background
155     """
156 
157     def __init__(self, n_on, mu_bkg):
158         self.n_on = np.asanyarray(n_on)
159         self.mu_bkg = np.asanyarray(mu_bkg)
160 
161     @property
162     def n_bkg(self):
163         """Expected background counts"""
164         return self.mu_bkg
165 
166     @property
167     def n_sig(self):
168         """Excess"""
169         return self.n_on - self.n_bkg
170 
171     @property
172     def error(self):
173         """Approximate error from the covariance matrix."""
174         return np.sqrt(self.n_on)
175 
176     @property
177     def stat_null(self):
178         """Stat value for null hypothesis, i.e. 0 expected signal counts"""
179         return cash(self.n_on, self.mu_bkg + 0)
180 
181     @property
182     def stat_max(self):
183         """Stat value for best fit hypothesis, i.e. expected signal mu = n_on - mu_bkg"""
184         return cash(self.n_on, self.n_on)
185 
186     def _stat_fcn(self, mu, delta=0, index=None):
187         return cash(self.n_on[index], self.mu_bkg[index] + mu) - delta
188 
189     def _n_sig_matching_significance_fcn(self, n_sig, significance, index):
190         TS0 = cash(n_sig + self.mu_bkg[index], self.mu_bkg[index])
191         TS1 = cash(n_sig + self.mu_bkg[index], self.mu_bkg[index] + n_sig)
192         return np.sign(n_sig) * np.sqrt(np.clip(TS0 - TS1, 0, None)) - significance
193 
194 
195 class WStatCountsStatistic(CountsStatistic):
196     """Class to compute statistics (significance, asymmetric errors , ul) for Poisson distributed variable
197     with unknown background.
198 
199     Parameters
200     ----------
201     n_on : int
202         Measured counts in on region
203     n_off : int
204         Measured counts in off region
205     alpha : float
206         Acceptance ratio of on and off measurements
207     mu_sig : float
208         Expected signal counts in on region
209     """
210 
211     def __init__(self, n_on, n_off, alpha, mu_sig=None):
212         self.n_on = np.asanyarray(n_on)
213         self.n_off = np.asanyarray(n_off)
214         self.alpha = np.asanyarray(alpha)
215         if mu_sig is None:
216             self.mu_sig = np.zeros_like(self.n_on)
217         else:
218             self.mu_sig = np.asanyarray(mu_sig)
219 
220     @property
221     def n_bkg(self):
222         """Known background computed alpha * n_off"""
223         return self.alpha * self.n_off
224 
225     @property
226     def n_sig(self):
227         """Excess """
228         return self.n_on - self.n_bkg - self.mu_sig
229 
230     @property
231     def error(self):
232         """Approximate error from the covariance matrix."""
233         return np.sqrt(self.n_on + self.alpha ** 2 * self.n_off)
234 
235     @property
236     def stat_null(self):
237         """Stat value for null hypothesis, i.e. mu_sig expected signal counts"""
238         return wstat(self.n_on, self.n_off, self.alpha, self.mu_sig)
239 
240     @property
241     def stat_max(self):
242         """Stat value for best fit hypothesis, i.e. expected signal mu = n_on - alpha * n_off - mu_sig"""
243         return wstat(self.n_on, self.n_off, self.alpha, self.n_sig + self.mu_sig)
244 
245     def _stat_fcn(self, mu, delta=0, index=None):
246         return (
247             wstat(
248                 self.n_on[index],
249                 self.n_off[index],
250                 self.alpha[index],
251                 (mu + self.mu_sig[index]),
252             )
253             - delta
254         )
255 
256     def _n_sig_matching_significance_fcn(self, n_sig, significance, index):
257         stat0 = wstat(
258             n_sig + self.n_bkg[index], self.n_off[index], self.alpha[index], 0
259         )
260         stat1 = wstat(
261             n_sig + self.n_bkg[index], self.n_off[index], self.alpha[index], n_sig,
262         )
263         return np.sign(n_sig) * np.sqrt(np.clip(stat0 - stat1, 0, None)) - significance
264
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/gammapy/stats/counts_statistic.py b/gammapy/stats/counts_statistic.py
--- a/gammapy/stats/counts_statistic.py
+++ b/gammapy/stats/counts_statistic.py
@@ -18,13 +18,17 @@
 
     @property
     def sqrt_ts(self):
-        """Return statistical significance of measured excess."""
+        """Return statistical significance of measured excess.
+           The sign of the excess is applied to distinguish positive and negative fluctuations.
+        """
         return np.sign(self.n_sig) * np.sqrt(self.ts)
 
     @property
     def p_value(self):
-        """Return p_value of measured excess."""
-        return chi2.sf(self.ts, 1)
+        """Return p_value of measured excess.
+           Here the value accounts only for the positive excess significance (i.e. one-sided).
+        """
+        return 0.5*chi2.sf(self.ts, 1)
 
     def compute_errn(self, n_sigma=1.0):
         """Compute downward excess uncertainties.

{"golden_diff": "diff --git a/gammapy/stats/counts_statistic.py b/gammapy/stats/counts_statistic.py\n--- a/gammapy/stats/counts_statistic.py\n+++ b/gammapy/stats/counts_statistic.py\n@@ -18,13 +18,17 @@\n \n     @property\n     def sqrt_ts(self):\n-        \"\"\"Return statistical significance of measured excess.\"\"\"\n+        \"\"\"Return statistical significance of measured excess.\n+           The sign of the excess is applied to distinguish positive and negative fluctuations.\n+        \"\"\"\n         return np.sign(self.n_sig) * np.sqrt(self.ts)\n \n     @property\n     def p_value(self):\n-        \"\"\"Return p_value of measured excess.\"\"\"\n-        return chi2.sf(self.ts, 1)\n+        \"\"\"Return p_value of measured excess.\n+           Here the value accounts only for the positive excess significance (i.e. one-sided).\n+        \"\"\"\n+        return 0.5*chi2.sf(self.ts, 1)\n \n     def compute_errn(self, n_sigma=1.0):\n         \"\"\"Compute downward excess uncertainties.\n", "issue": "Cash statistics : Inconsistency between doc/code or pvalue/siginificance\nGammpay 0.17, official\r\n\r\nHello\r\nIn the Cash count statistics description in the \"Stats-Statistics\" docs\r\nhttps://docs.gammapy.org/0.17/stats/index.html#cash-counts-statistic\r\n\r\nThe text says that *the significance is the square root of TS ... which is equivalent to the Li&Ma significance for Li&Ma...*\r\n\r\nThe statement is true but this is not what is in ```stat.signifcance```.\r\n\r\nAs you know the Li&Ma usual formula gives a significance which is \"single-sided\", i.e. that corresponds to a ```1-p_value``` obtained by the integration of the normal distribution from ```-Inf``` to the significance value (that can be negative). \r\n\r\nI have checked on the example that from the p_value the significance should be:\r\n```\r\nfrom scipy.stats import norm, Chi2\r\np_value = 1 - chi2.cdf(stat.delta_ts) # same value found : 7.365232895800901\r\nsig_LiMa = norm.ppf(1-p_value)\r\n```\r\nthat gives ``` 2.475649754526797```\r\n\r\nwhereas the results computed in the notebook is :```2.71389625737626``` \r\nI checked that this corresponds to a ``` p_value``` shared on both sides of the normal function:\r\n```norm.ppf(1-p_value/2)```\r\n\r\nAs a consequence, if I am correct:\r\n- the text is misleading because it mentions Li&Ma, and the situation should be explicit;\r\n- the significance attribute should be documented as the \"double-sided\" significance.\r\n\r\nHowever I have checked that ```CashCountsStatitics``` can return negative significances which in principle is not possible with a \"double-sided\" integration. So there might be another inconsistency somewhere....\r\n\r\nThanks!\r\n\r\nP.S.: one should never rely on siginiicances, only on p values.\r\n\n", "before_files": [{"content": "# Licensed under a 3-clause BSD style license - see LICENSE.rst\nimport abc\nimport numpy as np\nfrom scipy.optimize import brentq, newton\nfrom scipy.stats import chi2\nfrom .fit_statistics import cash, get_wstat_mu_bkg, wstat\n\n__all__ = [\"WStatCountsStatistic\", \"CashCountsStatistic\"]\n\n\nclass CountsStatistic(abc.ABC):\n    @property\n    def ts(self):\n        \"\"\"Return stat difference (TS) of measured excess versus no excess.\"\"\"\n        # Remove (small) negative TS due to error in root finding\n        ts = np.clip(self.stat_null - self.stat_max, 0, None)\n        return ts\n\n    @property\n    def sqrt_ts(self):\n        \"\"\"Return statistical significance of measured excess.\"\"\"\n        return np.sign(self.n_sig) * np.sqrt(self.ts)\n\n    @property\n    def p_value(self):\n        \"\"\"Return p_value of measured excess.\"\"\"\n        return chi2.sf(self.ts, 1)\n\n    def compute_errn(self, n_sigma=1.0):\n        \"\"\"Compute downward excess uncertainties.\n\n        Searches the signal value for which the test statistics is n_sigma**2 away from the maximum.\n\n        Parameters\n        ----------\n        n_sigma : float\n            Confidence level of the uncertainty expressed in number of sigma. Default is 1.\n        \"\"\"\n        errn = np.zeros_like(self.n_on, dtype=\"float\")\n        min_range = self.n_sig - 2 * n_sigma * (self.error + 1)\n\n        it = np.nditer(errn, flags=[\"multi_index\"])\n        while not it.finished:\n            try:\n                res = brentq(\n                    self._stat_fcn,\n                    min_range[it.multi_index],\n                    self.n_sig[it.multi_index],\n                    args=(self.stat_max[it.multi_index] + n_sigma ** 2, it.multi_index),\n                )\n                errn[it.multi_index] = res - self.n_sig[it.multi_index]\n            except ValueError:\n                errn[it.multi_index] = -self.n_on[it.multi_index]\n            it.iternext()\n\n        return errn\n\n    def compute_errp(self, n_sigma=1):\n        \"\"\"Compute upward excess uncertainties.\n\n        Searches the signal value for which the test statistics is n_sigma**2 away from the maximum.\n\n        Parameters\n        ----------\n        n_sigma : float\n            Confidence level of the uncertainty expressed in number of sigma. Default is 1.\n        \"\"\"\n        errp = np.zeros_like(self.n_on, dtype=\"float\")\n        max_range = self.n_sig + 2 * n_sigma * (self.error + 1)\n\n        it = np.nditer(errp, flags=[\"multi_index\"])\n        while not it.finished:\n            errp[it.multi_index] = brentq(\n                self._stat_fcn,\n                self.n_sig[it.multi_index],\n                max_range[it.multi_index],\n                args=(self.stat_max[it.multi_index] + n_sigma ** 2, it.multi_index),\n            )\n            it.iternext()\n\n        return errp - self.n_sig\n\n    def compute_upper_limit(self, n_sigma=3):\n        \"\"\"Compute upper limit on the signal.\n\n        Searches the signal value for which the test statistics is n_sigma**2 away from the maximum\n        or from 0 if the measured excess is negative.\n\n        Parameters\n        ----------\n        n_sigma : float\n            Confidence level of the upper limit expressed in number of sigma. Default is 3.\n        \"\"\"\n        ul = np.zeros_like(self.n_on, dtype=\"float\")\n\n        min_range = np.maximum(0, self.n_sig)\n        max_range = min_range + 2 * n_sigma * (self.error + 1)\n        it = np.nditer(ul, flags=[\"multi_index\"])\n\n        while not it.finished:\n            TS_ref = self._stat_fcn(min_range[it.multi_index], 0.0, it.multi_index)\n\n            ul[it.multi_index] = brentq(\n                self._stat_fcn,\n                min_range[it.multi_index],\n                max_range[it.multi_index],\n                args=(TS_ref + n_sigma ** 2, it.multi_index),\n            )\n            it.iternext()\n\n        return ul\n\n    def n_sig_matching_significance(self, significance):\n        \"\"\"Compute excess matching a given significance.\n\n        This function is the inverse of `significance`.\n\n        Parameters\n        ----------\n        significance : float\n            Significance\n\n        Returns\n        -------\n        n_sig : `numpy.ndarray`\n            Excess\n        \"\"\"\n        n_sig = np.zeros_like(self.n_bkg, dtype=\"float\")\n        it = np.nditer(n_sig, flags=[\"multi_index\"])\n\n        while not it.finished:\n            try:\n                n_sig[it.multi_index] = newton(\n                    self._n_sig_matching_significance_fcn,\n                    np.sqrt(self.n_bkg[it.multi_index]) * significance,\n                    args=(significance, it.multi_index),\n                )\n            except:\n                n_sig[it.multi_index] = np.nan\n\n            it.iternext()\n        return n_sig\n\n\nclass CashCountsStatistic(CountsStatistic):\n    \"\"\"Class to compute statistics (significance, asymmetric errors , ul) for Poisson distributed variable\n    with known background.\n\n    Parameters\n    ----------\n    n_on : int\n        Measured counts\n    mu_bkg : float\n        Known level of background\n    \"\"\"\n\n    def __init__(self, n_on, mu_bkg):\n        self.n_on = np.asanyarray(n_on)\n        self.mu_bkg = np.asanyarray(mu_bkg)\n\n    @property\n    def n_bkg(self):\n        \"\"\"Expected background counts\"\"\"\n        return self.mu_bkg\n\n    @property\n    def n_sig(self):\n        \"\"\"Excess\"\"\"\n        return self.n_on - self.n_bkg\n\n    @property\n    def error(self):\n        \"\"\"Approximate error from the covariance matrix.\"\"\"\n        return np.sqrt(self.n_on)\n\n    @property\n    def stat_null(self):\n        \"\"\"Stat value for null hypothesis, i.e. 0 expected signal counts\"\"\"\n        return cash(self.n_on, self.mu_bkg + 0)\n\n    @property\n    def stat_max(self):\n        \"\"\"Stat value for best fit hypothesis, i.e. expected signal mu = n_on - mu_bkg\"\"\"\n        return cash(self.n_on, self.n_on)\n\n    def _stat_fcn(self, mu, delta=0, index=None):\n        return cash(self.n_on[index], self.mu_bkg[index] + mu) - delta\n\n    def _n_sig_matching_significance_fcn(self, n_sig, significance, index):\n        TS0 = cash(n_sig + self.mu_bkg[index], self.mu_bkg[index])\n        TS1 = cash(n_sig + self.mu_bkg[index], self.mu_bkg[index] + n_sig)\n        return np.sign(n_sig) * np.sqrt(np.clip(TS0 - TS1, 0, None)) - significance\n\n\nclass WStatCountsStatistic(CountsStatistic):\n    \"\"\"Class to compute statistics (significance, asymmetric errors , ul) for Poisson distributed variable\n    with unknown background.\n\n    Parameters\n    ----------\n    n_on : int\n        Measured counts in on region\n    n_off : int\n        Measured counts in off region\n    alpha : float\n        Acceptance ratio of on and off measurements\n    mu_sig : float\n        Expected signal counts in on region\n    \"\"\"\n\n    def __init__(self, n_on, n_off, alpha, mu_sig=None):\n        self.n_on = np.asanyarray(n_on)\n        self.n_off = np.asanyarray(n_off)\n        self.alpha = np.asanyarray(alpha)\n        if mu_sig is None:\n            self.mu_sig = np.zeros_like(self.n_on)\n        else:\n            self.mu_sig = np.asanyarray(mu_sig)\n\n    @property\n    def n_bkg(self):\n        \"\"\"Known background computed alpha * n_off\"\"\"\n        return self.alpha * self.n_off\n\n    @property\n    def n_sig(self):\n        \"\"\"Excess \"\"\"\n        return self.n_on - self.n_bkg - self.mu_sig\n\n    @property\n    def error(self):\n        \"\"\"Approximate error from the covariance matrix.\"\"\"\n        return np.sqrt(self.n_on + self.alpha ** 2 * self.n_off)\n\n    @property\n    def stat_null(self):\n        \"\"\"Stat value for null hypothesis, i.e. mu_sig expected signal counts\"\"\"\n        return wstat(self.n_on, self.n_off, self.alpha, self.mu_sig)\n\n    @property\n    def stat_max(self):\n        \"\"\"Stat value for best fit hypothesis, i.e. expected signal mu = n_on - alpha * n_off - mu_sig\"\"\"\n        return wstat(self.n_on, self.n_off, self.alpha, self.n_sig + self.mu_sig)\n\n    def _stat_fcn(self, mu, delta=0, index=None):\n        return (\n            wstat(\n                self.n_on[index],\n                self.n_off[index],\n                self.alpha[index],\n                (mu + self.mu_sig[index]),\n            )\n            - delta\n        )\n\n    def _n_sig_matching_significance_fcn(self, n_sig, significance, index):\n        stat0 = wstat(\n            n_sig + self.n_bkg[index], self.n_off[index], self.alpha[index], 0\n        )\n        stat1 = wstat(\n            n_sig + self.n_bkg[index], self.n_off[index], self.alpha[index], n_sig,\n        )\n        return np.sign(n_sig) * np.sqrt(np.clip(stat0 - stat1, 0, None)) - significance\n", "path": "gammapy/stats/counts_statistic.py"}], "after_files": [{"content": "# Licensed under a 3-clause BSD style license - see LICENSE.rst\nimport abc\nimport numpy as np\nfrom scipy.optimize import brentq, newton\nfrom scipy.stats import chi2\nfrom .fit_statistics import cash, get_wstat_mu_bkg, wstat\n\n__all__ = [\"WStatCountsStatistic\", \"CashCountsStatistic\"]\n\n\nclass CountsStatistic(abc.ABC):\n    @property\n    def ts(self):\n        \"\"\"Return stat difference (TS) of measured excess versus no excess.\"\"\"\n        # Remove (small) negative TS due to error in root finding\n        ts = np.clip(self.stat_null - self.stat_max, 0, None)\n        return ts\n\n    @property\n    def sqrt_ts(self):\n        \"\"\"Return statistical significance of measured excess.\n           The sign of the excess is applied to distinguish positive and negative fluctuations.\n        \"\"\"\n        return np.sign(self.n_sig) * np.sqrt(self.ts)\n\n    @property\n    def p_value(self):\n        \"\"\"Return p_value of measured excess.\n           Here the value accounts only for the positive excess significance (i.e. one-sided).\n        \"\"\"\n        return 0.5*chi2.sf(self.ts, 1)\n\n    def compute_errn(self, n_sigma=1.0):\n        \"\"\"Compute downward excess uncertainties.\n\n        Searches the signal value for which the test statistics is n_sigma**2 away from the maximum.\n\n        Parameters\n        ----------\n        n_sigma : float\n            Confidence level of the uncertainty expressed in number of sigma. Default is 1.\n        \"\"\"\n        errn = np.zeros_like(self.n_on, dtype=\"float\")\n        min_range = self.n_sig - 2 * n_sigma * (self.error + 1)\n\n        it = np.nditer(errn, flags=[\"multi_index\"])\n        while not it.finished:\n            try:\n                res = brentq(\n                    self._stat_fcn,\n                    min_range[it.multi_index],\n                    self.n_sig[it.multi_index],\n                    args=(self.stat_max[it.multi_index] + n_sigma ** 2, it.multi_index),\n                )\n                errn[it.multi_index] = res - self.n_sig[it.multi_index]\n            except ValueError:\n                errn[it.multi_index] = -self.n_on[it.multi_index]\n            it.iternext()\n\n        return errn\n\n    def compute_errp(self, n_sigma=1):\n        \"\"\"Compute upward excess uncertainties.\n\n        Searches the signal value for which the test statistics is n_sigma**2 away from the maximum.\n\n        Parameters\n        ----------\n        n_sigma : float\n            Confidence level of the uncertainty expressed in number of sigma. Default is 1.\n        \"\"\"\n        errp = np.zeros_like(self.n_on, dtype=\"float\")\n        max_range = self.n_sig + 2 * n_sigma * (self.error + 1)\n\n        it = np.nditer(errp, flags=[\"multi_index\"])\n        while not it.finished:\n            errp[it.multi_index] = brentq(\n                self._stat_fcn,\n                self.n_sig[it.multi_index],\n                max_range[it.multi_index],\n                args=(self.stat_max[it.multi_index] + n_sigma ** 2, it.multi_index),\n            )\n            it.iternext()\n\n        return errp - self.n_sig\n\n    def compute_upper_limit(self, n_sigma=3):\n        \"\"\"Compute upper limit on the signal.\n\n        Searches the signal value for which the test statistics is n_sigma**2 away from the maximum\n        or from 0 if the measured excess is negative.\n\n        Parameters\n        ----------\n        n_sigma : float\n            Confidence level of the upper limit expressed in number of sigma. Default is 3.\n        \"\"\"\n        ul = np.zeros_like(self.n_on, dtype=\"float\")\n\n        min_range = np.maximum(0, self.n_sig)\n        max_range = min_range + 2 * n_sigma * (self.error + 1)\n        it = np.nditer(ul, flags=[\"multi_index\"])\n\n        while not it.finished:\n            TS_ref = self._stat_fcn(min_range[it.multi_index], 0.0, it.multi_index)\n\n            ul[it.multi_index] = brentq(\n                self._stat_fcn,\n                min_range[it.multi_index],\n                max_range[it.multi_index],\n                args=(TS_ref + n_sigma ** 2, it.multi_index),\n            )\n            it.iternext()\n\n        return ul\n\n    def n_sig_matching_significance(self, significance):\n        \"\"\"Compute excess matching a given significance.\n\n        This function is the inverse of `significance`.\n\n        Parameters\n        ----------\n        significance : float\n            Significance\n\n        Returns\n        -------\n        n_sig : `numpy.ndarray`\n            Excess\n        \"\"\"\n        n_sig = np.zeros_like(self.n_bkg, dtype=\"float\")\n        it = np.nditer(n_sig, flags=[\"multi_index\"])\n\n        while not it.finished:\n            try:\n                n_sig[it.multi_index] = newton(\n                    self._n_sig_matching_significance_fcn,\n                    np.sqrt(self.n_bkg[it.multi_index]) * significance,\n                    args=(significance, it.multi_index),\n                )\n            except:\n                n_sig[it.multi_index] = np.nan\n\n            it.iternext()\n        return n_sig\n\n\nclass CashCountsStatistic(CountsStatistic):\n    \"\"\"Class to compute statistics (significance, asymmetric errors , ul) for Poisson distributed variable\n    with known background.\n\n    Parameters\n    ----------\n    n_on : int\n        Measured counts\n    mu_bkg : float\n        Known level of background\n    \"\"\"\n\n    def __init__(self, n_on, mu_bkg):\n        self.n_on = np.asanyarray(n_on)\n        self.mu_bkg = np.asanyarray(mu_bkg)\n\n    @property\n    def n_bkg(self):\n        \"\"\"Expected background counts\"\"\"\n        return self.mu_bkg\n\n    @property\n    def n_sig(self):\n        \"\"\"Excess\"\"\"\n        return self.n_on - self.n_bkg\n\n    @property\n    def error(self):\n        \"\"\"Approximate error from the covariance matrix.\"\"\"\n        return np.sqrt(self.n_on)\n\n    @property\n    def stat_null(self):\n        \"\"\"Stat value for null hypothesis, i.e. 0 expected signal counts\"\"\"\n        return cash(self.n_on, self.mu_bkg + 0)\n\n    @property\n    def stat_max(self):\n        \"\"\"Stat value for best fit hypothesis, i.e. expected signal mu = n_on - mu_bkg\"\"\"\n        return cash(self.n_on, self.n_on)\n\n    def _stat_fcn(self, mu, delta=0, index=None):\n        return cash(self.n_on[index], self.mu_bkg[index] + mu) - delta\n\n    def _n_sig_matching_significance_fcn(self, n_sig, significance, index):\n        TS0 = cash(n_sig + self.mu_bkg[index], self.mu_bkg[index])\n        TS1 = cash(n_sig + self.mu_bkg[index], self.mu_bkg[index] + n_sig)\n        return np.sign(n_sig) * np.sqrt(np.clip(TS0 - TS1, 0, None)) - significance\n\n\nclass WStatCountsStatistic(CountsStatistic):\n    \"\"\"Class to compute statistics (significance, asymmetric errors , ul) for Poisson distributed variable\n    with unknown background.\n\n    Parameters\n    ----------\n    n_on : int\n        Measured counts in on region\n    n_off : int\n        Measured counts in off region\n    alpha : float\n        Acceptance ratio of on and off measurements\n    mu_sig : float\n        Expected signal counts in on region\n    \"\"\"\n\n    def __init__(self, n_on, n_off, alpha, mu_sig=None):\n        self.n_on = np.asanyarray(n_on)\n        self.n_off = np.asanyarray(n_off)\n        self.alpha = np.asanyarray(alpha)\n        if mu_sig is None:\n            self.mu_sig = np.zeros_like(self.n_on)\n        else:\n            self.mu_sig = np.asanyarray(mu_sig)\n\n    @property\n    def n_bkg(self):\n        \"\"\"Known background computed alpha * n_off\"\"\"\n        return self.alpha * self.n_off\n\n    @property\n    def n_sig(self):\n        \"\"\"Excess \"\"\"\n        return self.n_on - self.n_bkg - self.mu_sig\n\n    @property\n    def error(self):\n        \"\"\"Approximate error from the covariance matrix.\"\"\"\n        return np.sqrt(self.n_on + self.alpha ** 2 * self.n_off)\n\n    @property\n    def stat_null(self):\n        \"\"\"Stat value for null hypothesis, i.e. mu_sig expected signal counts\"\"\"\n        return wstat(self.n_on, self.n_off, self.alpha, self.mu_sig)\n\n    @property\n    def stat_max(self):\n        \"\"\"Stat value for best fit hypothesis, i.e. expected signal mu = n_on - alpha * n_off - mu_sig\"\"\"\n        return wstat(self.n_on, self.n_off, self.alpha, self.n_sig + self.mu_sig)\n\n    def _stat_fcn(self, mu, delta=0, index=None):\n        return (\n            wstat(\n                self.n_on[index],\n                self.n_off[index],\n                self.alpha[index],\n                (mu + self.mu_sig[index]),\n            )\n            - delta\n        )\n\n    def _n_sig_matching_significance_fcn(self, n_sig, significance, index):\n        stat0 = wstat(\n            n_sig + self.n_bkg[index], self.n_off[index], self.alpha[index], 0\n        )\n        stat1 = wstat(\n            n_sig + self.n_bkg[index], self.n_off[index], self.alpha[index], n_sig,\n        )\n        return np.sign(n_sig) * np.sqrt(np.clip(stat0 - stat1, 0, None)) - significance\n", "path": "gammapy/stats/counts_statistic.py"}]}