Add summed_flag to power confidence limits

docs/changes/952.feature.rst

@@ -0,0 +1 @@
+Added ``summed_flag`` to ``power_confidence_limits``, to make it work with averaged powers too.

stingray/stats.py

@@ -906,7 +906,7 @@ def _pavnosigfun(power, nspec):
     return sum
 
 
-def power_confidence_limits(preal, n=1, c=0.95):
+def power_confidence_limits(preal, n=1, c=0.95, summed_flag=True):
     """Confidence limits on power, given a (theoretical) signal power.
 
     This is to be used when we *expect* a given power (e.g. from the pulsed
@@ -930,6 +930,9 @@ def power_confidence_limits(preal, n=1, c=0.95):
         power in a QPO, or the n in Z^2_n
     c: float
         The confidence level (e.g. 0.95=95%)
+    summed_flag: bool
+        Whether the power is (if True) summed or (if False) averaged. Only
+        relevant if n > 1
 
     Returns
     -------
@@ -939,11 +942,22 @@ def power_confidence_limits(preal, n=1, c=0.95):
 
     Examples
     --------
-    >>> cl = power_confidence_limits(150, c=0.68)
-    >>> assert np.allclose(cl, [128, 176], atol=1)
+    >>> cl = power_confidence_limits(150, c=0.68, n=1, summed_flag=True)
+    >>> assert np.allclose(cl, [127, 176], atol=1)
+    >>> cl = power_confidence_limits(150, c=0.68, n=8, summed_flag=True)
+    >>> assert np.allclose(cl, [141, 190], atol=1)
+    >>> cl = power_confidence_limits(150, c=0.68, n=1, summed_flag=False)
+    >>> assert np.allclose(cl, [127, 176], atol=1)
+    >>> cl = power_confidence_limits(150, c=0.68, n=8, summed_flag=False)
+    >>> assert np.allclose(cl, [143, 160], atol=1)
     """
-    rv = stats.ncx2(2 * n, preal)
-    return rv.ppf([(1 - c) / 2, (1 + c) / 2])
+    if summed_flag:
+        rv = stats.ncx2(2 * n, preal)
+        ints = rv.ppf([(1 - c) / 2, (1 + c) / 2])
+    else:
+        rv = stats.ncx2(2 * n, preal * n)
+        ints = rv.ppf([(1 - c) / 2, (1 + c) / 2]) / n
+    return ints
 
 
 def power_upper_limit(pmeas, n=1, c=0.95, summed_flag=True):

stingray/tests/test_stats.py

@@ -270,6 +270,7 @@ def test_power_upper_limit_averaging(self):
 
     def test_power_confidence_limits(self):
         from scipy.stats import ncx2
+
         preal = 25.0
         n = 2
         c = 0.95