Repository: cta-observatory/pyirf

Issue:
li_ma_significance needs float64 for large numbers
```
noff = 1e7
nexcess = 1e4

print(li_ma_significance(np.float64(noff+nexcess), np.float64(noff / 0.2), 0.2))
print(li_ma_significance(np.float32(noff+nexcess), np.float32(noff / 0.2), 0.2))

2.8861902631323018
3.7806044
```
May be this is unavoidable, I don't know... or perhaps it can be fixed just by using float64 internally in the calculation?


Discussion:

Patch:
diff --git a/pyirf/statistics.py b/pyirf/statistics.py
index 9abd7c1..b4f5b10 100644
--- a/pyirf/statistics.py
+++ b/pyirf/statistics.py
@@ -2,7 +2,7 @@ import numpy as np
 
 from .utils import is_scalar
 
-__all__  = ['li_ma_significance']
+__all__ = ["li_ma_significance"]
 
 
 def li_ma_significance(n_on, n_off, alpha=0.2):
@@ -19,7 +19,7 @@ def li_ma_significance(n_on, n_off, alpha=0.2):
     ----------
     n_on: integer or array like
         Number of events for the on observations
-    n_off: integer of array like
+    n_off: integer or array like
         Number of events for the off observations
     alpha: float
         Ratio between the on region and the off region size or obstime.
@@ -32,8 +32,11 @@ def li_ma_significance(n_on, n_off, alpha=0.2):
 
     scalar = is_scalar(n_on)
 
-    n_on = np.array(n_on, copy=False, ndmin=1)
-    n_off = np.array(n_off, copy=False, ndmin=1)
+    # Cast everything into float64 to avoid numeric instabilties
+    # when multiplying very small and very big numbers to get t1 and t2
+    n_on = np.array(n_on, copy=False, ndmin=1, dtype=np.float64)
+    n_off = np.array(n_off, copy=False, ndmin=1, dtype=np.float64)
+    alpha = np.float64(alpha)
 
     with np.errstate(divide="ignore", invalid="ignore"):
         p_on = n_on / (n_on + n_off)


Test patch:
diff --git a/pyirf/tests/test_statistics.py b/pyirf/tests/test_statistics.py
index 54ca8da..c56aac9 100644
--- a/pyirf/tests/test_statistics.py
+++ b/pyirf/tests/test_statistics.py
@@ -11,7 +11,7 @@ def test_lima():
 
 
 def test_lima_gammapy():
-    pytest.importorskip('gammapy')
+    pytest.importorskip("gammapy")
     from gammapy.stats import WStatCountsStatistic
     from pyirf.statistics import li_ma_significance
 
@@ -21,3 +21,19 @@ def test_lima_gammapy():
     for n_on, n_off, alpha in zip(n_ons, n_offs, alphas):
         sig_gammapy = WStatCountsStatistic(n_on, n_off, alpha).sqrt_ts
         assert np.isclose(li_ma_significance(n_on, n_off, alpha), sig_gammapy)
+
+
+def test_lima_accuracy():
+    from pyirf.statistics import li_ma_significance
+
+    noff = 1e7
+    nexcess = 1e4
+
+    res_f64 = li_ma_significance(
+        np.float64(noff + nexcess), np.float64(noff / 0.2), 0.2
+    )
+    res_f32 = li_ma_significance(
+        np.float32(noff + nexcess), np.float32(noff / 0.2), 0.2
+    )
+
+    assert np.isclose(res_f64, res_f32)