Merge pull request #953 from ericymiao/bugfix-rebin_data_uncertainties

Fixed incorrect error propagation in rebin_data

Author: matteobachetti

docs/changes/953.bugfix.rst

@@ -0,0 +1 @@
+Fixed uncertainty propogation error in ``Lightcurve.rebin()`` 

stingray/tests/test_utils.py

@@ -88,6 +88,19 @@ def test_rebin_variable_input_mean(self):
         assert len(xbin) == 20
         assert np.allclose(ybin, counts)
 
+    def test_rebin_error_propagation(self):
+        x = np.arange(0, 4, 1.0)  # 4 points
+        y = np.ones(4) * 10
+        yerr = np.ones(4) * 2  # Each point has σ=2
+
+        xbin, ybin, ybinerr, _ = utils.rebin_data(x, y, 4.0, yerr=yerr, method="sum")
+        # Correct: sqrt(4*2**2) = sqrt(16) = 4
+        assert np.allclose(ybinerr, 4.0)
+
+        xbin, ybin, ybinerr, _ = utils.rebin_data(x, y, 4.0, yerr=yerr, method="mean")
+        # Correct: sqrt(4*2**2) / 4 = 4/4 = 1
+        assert np.allclose(ybinerr, 1.0)
+
 
 class TestRebinDataLog(object):
     @classmethod
@@ -97,7 +110,7 @@ def setup_class(cls):
         cls.xmin = 1
         cls.x = np.arange(cls.xmin, cls.xmax, cls.dx)
         cls.y = np.arange(cls.xmin, cls.xmax, cls.dx)
-        cls.y_err = np.ones_like(cls.y)
+        cls.y_err = np.ones_like(cls.y) * 2
 
         cls.true_bins = np.array([1.0, 1.0, 1.5, 1.5, 1.5, 1.5, 1.5, 2.5, 2.5, 2.5, 2.0])
 
@@ -139,26 +152,19 @@ def test_all_outputs_have_the_same_dimension_except_binx(self):
         assert binyerr.shape[0] == nsamples.shape[0]
 
     def test_binning_works_correctly(self):
-        binx, _, _, _ = utils.rebin_data_log(self.x, self.y, self.f, y_err=self.y_err, dx=self.dx)
+        binx, biny, binyerr, nsamples = utils.rebin_data_log(
+            self.x,
+            self.y,
+            self.f,
+            y_err=self.y_err,
+            dx=self.dx,
+        )
 
         assert np.allclose(np.diff(binx), self.true_bins)
-
-    def test_bin_edges_are_correct(self):
-        binx, _, _, _ = utils.rebin_data_log(self.x, self.y, self.f, y_err=self.y_err, dx=self.dx)
-
         assert np.allclose(binx, self.true_bin_edges)
-
-    def test_bin_values_are_correct(self):
-        _, biny, _, _ = utils.rebin_data_log(self.x, self.y, self.f, y_err=self.y_err, dx=self.dx)
-
         assert np.allclose(biny, self.true_values)
-
-    def test_nsamples_are_correctly_calculated(self):
-        _, _, _, nsamples = utils.rebin_data_log(
-            self.x, self.y, self.f, y_err=self.y_err, dx=self.dx
-        )
-
         assert np.allclose(nsamples, self.true_nsamples)
+        assert np.allclose(binyerr, 2 / np.sqrt(self.true_nsamples))
 
     def test_method_works_on_complex_numbers(self):
         re = np.arange(self.xmin, self.xmax, self.dx)

stingray/utils.py

@@ -534,15 +534,15 @@ def rebin_data(x, y, dx_new, yerr=None, method="sum", dx=None):
     --------
     >>> x = np.arange(0, 100, 0.01)
     >>> y = np.ones(x.size)
-    >>> yerr = np.ones(x.size)
+    >>> yerr = np.ones(x.size) * 2  # Each point has σ=2
     >>> xbin, ybin, ybinerr, step_size = rebin_data(
     ...     x, y, 4, yerr=yerr, method='sum', dx=0.01)
     >>> assert np.allclose(ybin, 400)
-    >>> assert np.allclose(ybinerr, 20)
+    >>> assert np.allclose(ybinerr, 40)  # sqrt(400 * 2**2) = 40
     >>> xbin, ybin, ybinerr, step_size = rebin_data(
     ...     x, y, 4, yerr=yerr, method='mean')
     >>> assert np.allclose(ybin, 1)
-    >>> assert np.allclose(ybinerr, 0.05)
+    >>> assert np.allclose(ybinerr, 0.1)  # sqrt(400 * 2**2) / 400 = 0.1
     """
 
     y = np.asanyarray(y)
@@ -572,7 +572,7 @@ def rebin_data(x, y, dx_new, yerr=None, method="sum", dx=None):
     xbin = np.arange(xedges[0], xedges[-1] + dx_new, dx_new)
 
     output = np.zeros(xbin.shape[0] - 1, dtype=type(y[0]))
-    outputerr = np.zeros(xbin.shape[0] - 1, dtype=type(yerr[0]))
+    outputvar = np.zeros(xbin.shape[0] - 1, dtype=type(yerr[0]))
     step_size = np.zeros(xbin.shape[0] - 1)
 
     all_x = np.searchsorted(xedges, xbin)
@@ -584,29 +584,29 @@ def rebin_data(x, y, dx_new, yerr=None, method="sum", dx=None):
         filtered_y = y[min_ind : max_ind - 1]
         filtered_yerr = yerr[min_ind : max_ind - 1]
         output[i] = np.sum(filtered_y)
-        outputerr[i] = np.sum(filtered_yerr)
+        outputvar[i] = np.sum(filtered_yerr**2)
         step_size[i] = max_ind - 1 - min_ind
 
         prev_dx = xedges[min_ind] - xedges[min_ind - 1]
         prev_frac = (xedges[min_ind] - xmin) / prev_dx
         output[i] += y[min_ind - 1] * prev_frac
-        outputerr[i] += yerr[min_ind - 1] * prev_frac
+        outputvar[i] += (yerr[min_ind - 1] * prev_frac) ** 2
         step_size[i] += prev_frac
 
         if not max_ind == xedges.size:
             dx_post = xedges[max_ind] - xedges[max_ind - 1]
             post_frac = (xmax - xedges[max_ind - 1]) / dx_post
             output[i] += y[max_ind - 1] * post_frac
-            outputerr[i] += yerr[max_ind - 1] * post_frac
+            outputvar[i] += (yerr[max_ind - 1] * post_frac) ** 2
             step_size[i] += post_frac
 
     if method in ["mean", "avg", "average", "arithmetic mean"]:
         ybin = output / step_size
-        ybinerr = np.sqrt(outputerr) / step_size
+        ybinerr = np.sqrt(outputvar) / step_size
 
     elif method == "sum":
         ybin = output
-        ybinerr = np.sqrt(outputerr)
+        ybinerr = np.sqrt(outputvar)
 
     else:
         raise ValueError(