fix stats

StFroese · Jan 27, 2025 · 79eaf94 · 79eaf94
1 parent 15499cf
commit 79eaf94
Show file tree

Hide file tree

Showing 2 changed files with 19 additions and 7 deletions.
diff --git a/titrate/datasets.py b/titrate/datasets.py
@@ -56,7 +56,7 @@ def fake(self):
         npred = self.npred()
         data = np.nan_to_num(npred.data, copy=True, nan=0.0, posinf=0.0, neginf=0.0)
         npred.data = data
-        self.counts = npred.data.copy()
+        self.counts.data = npred.data.copy()
 
     @classmethod
     def from_SpectralDataset(self, dataset):

diff --git a/titrate/statistics.py b/titrate/statistics.py
@@ -2,7 +2,7 @@
 
 import numpy as np
 from gammapy.modeling import Fit
-from scipy.stats import kstwo, ncx2, norm
+from scipy.stats import kstwo, norm
 
 from titrate.datasets import AsimovMapDataset, AsimovSpectralDataset
 
@@ -260,8 +260,6 @@ def asympotic_approximation_pdf(
                     / (2 * np.sqrt(2 * np.pi) * np.sqrt(ts_val))
                     * np.exp(-0.5 * (np.sqrt(ts_val)) ** 2)
                 ),
-                #     0,
-                # )
                 (
                     1
                     / (np.sqrt(2 * np.pi) * (2 * np.sqrt(mu_sigma)))
@@ -280,8 +278,6 @@ def asympotic_approximation_pdf(
                 / (2 * np.sqrt(2 * np.pi) * np.sqrt(ts_val))
                 * np.exp(-0.5 * (np.sqrt(ts_val) - np.sqrt(nc)) ** 2)
             ),
-            #     0,
-            # )
             (
                 1
                 / (np.sqrt(2 * np.pi) * (2 * np.sqrt(mu_sigma)))
@@ -293,7 +289,23 @@ def asympotic_approximation_cdf(
         self, ts_val, poi_val, same=True, poi_true_val=None
     ):
         nc = self.evaluate(poi_val)
-        return ncx2.cdf(ts_val, nc=nc, df=1)
+
+        if same:
+            sigma = np.sqrt(self.fit_result.covariance_result.matrix[0, 0])
+            mu_sigma = poi_val**2 / sigma**2
+            return np.where(
+                (ts_val > 0) & (ts_val <= mu_sigma),
+                norm.cdf(np.sqrt(ts_val)),
+                norm.cdf((ts_val + mu_sigma) / (2 * np.sqrt(mu_sigma))),
+            )
+
+        sigma = poi_val / np.sqrt(nc)
+        mu_sigma = poi_val**2 / sigma**2
+        return np.where(
+            (ts_val > 0) & (ts_val <= mu_sigma),
+            norm.cdf(np.sqrt(ts_val) - np.sqrt(nc)),
+            norm.cdf((ts_val - nc) / (2 * np.sqrt(mu_sigma))),
+        )
 
     def pvalue(self, poi_val, same=True, poi_true_val=None, ts_val=None):
         if ts_val is None: