diff --git a/pyerrors/obs.py b/pyerrors/obs.py
index 481b71a8..03486f14 100644
--- a/pyerrors/obs.py
+++ b/pyerrors/obs.py
@@ -279,11 +279,17 @@ class Obs:
                 tmp = self.e_rho[e_name][i + 1:w_max] + np.concatenate([self.e_rho[e_name][i - 1::-1], self.e_rho[e_name][1:w_max - 2 * i]]) - 2 * self.e_rho[e_name][i] * self.e_rho[e_name][1:w_max - i]
                 self.e_drho[e_name][i] = np.sqrt(np.sum(tmp ** 2) / e_N)
 
-            # detect regular step size in irregular MC chain
-            gapsize = 1
-            for r_name in e_content[e_name][:1]:
-                if not isinstance(self.idl[r_name], range):
-                    gapsize = np.min(np.diff(self.idl[r_name]))
+            gaps = []
+            for r_name in e_content[e_name]:
+                if isinstance(self.idl[r_name], range):
+                    gaps.append(1)
+                else:
+                    gaps.append(np.min(np.diff(self.idl[r_name])))
+
+            if not np.all([gi == gaps[0] for gi in gaps]):
+                raise Exception(f"Replica for ensemble {e_name} are not equally spaced.", gaps)
+            else:
+                gapsize = gaps[0]
 
             _compute_drho(gapsize)
             if self.tau_exp[e_name] > 0:
@@ -295,7 +301,7 @@ class Obs:
                     _compute_drho(n + gapsize)
                     if (self.e_rho[e_name][n] - self.N_sigma[e_name] * self.e_drho[e_name][n]) < 0 or n >= w_max // 2 - 2:
                         # Bias correction hep-lat/0306017 eq. (49) included
-                        self.e_tauint[e_name] = self.e_n_tauint[e_name][n] * (1 + (2 * n + 1) / e_N) / (1 + 1 / e_N) + texp * np.abs(self.e_rho[e_name][n + 1])  # The absolute makes sure, that the tail contribution is always positive
+                        self.e_tauint[e_name] = self.e_n_tauint[e_name][n] * (1 + (2 * n / gapsize + 1) / e_N) / (1 + 1 / e_N) + texp * np.abs(self.e_rho[e_name][n + 1])  # The absolute makes sure, that the tail contribution is always positive
                         self.e_dtauint[e_name] = np.sqrt(self.e_n_dtauint[e_name][n] ** 2 + texp ** 2 * self.e_drho[e_name][n + 1] ** 2)
                         # Error of tau_exp neglected so far, missing term: self.e_rho[e_name][n + 1] ** 2 * d_tau_exp ** 2
                         self.e_dvalue[e_name] = np.sqrt(2 * self.e_tauint[e_name] * e_gamma[e_name][0] * (1 + 1 / e_N) / e_N)
@@ -318,7 +324,7 @@ class Obs:
                             _compute_drho(gapsize * n + gapsize)
                         if g_w[n - 1] < 0 or n >= w_max - 1:
                             n *= gapsize
-                            self.e_tauint[e_name] = self.e_n_tauint[e_name][n] * (1 + (2 * n + 1) / e_N) / (1 + 1 / e_N)  # Bias correction hep-lat/0306017 eq. (49)
+                            self.e_tauint[e_name] = self.e_n_tauint[e_name][n] * (1 + (2 * n / gapsize + 1) / e_N) / (1 + 1 / e_N)  # Bias correction hep-lat/0306017 eq. (49)
                             self.e_dtauint[e_name] = self.e_n_dtauint[e_name][n]
                             self.e_dvalue[e_name] = np.sqrt(2 * self.e_tauint[e_name] * e_gamma[e_name][0] * (1 + 1 / e_N) / e_N)
                             self.e_ddvalue[e_name] = self.e_dvalue[e_name] * np.sqrt((n + 0.5) / e_N)
diff --git a/tests/obs_test.py b/tests/obs_test.py
index d996e065..29ee67ba 100644
--- a/tests/obs_test.py
+++ b/tests/obs_test.py
@@ -649,6 +649,25 @@ def test_gamma_method_irregular():
     ol = [a, b]
     o = (ol[0] - ol[1]) / (ol[1])
 
+    N = 1000
+    dat = gen_autocorrelated_array(np.random.normal(1, .2, size=N), .8)
+
+    idl_a = list(range(0, 1001, 1))
+    idl_a.remove(101)
+
+    oa = pe.Obs([dat], ["ens1"], idl=[idl_a])
+    oa.gamma_method()
+    tau_a = oa.e_tauint["ens1"]
+
+    idl_b = list(range(0, 10010, 10))
+    idl_b.remove(1010)
+
+    ob = pe.Obs([dat], ["ens1"], idl=[idl_b])
+    ob.gamma_method()
+    tau_b = ob.e_tauint["ens1"]
+
+    assert np.isclose(tau_a, tau_b)
+
 
 def test_covariance_is_variance():
     value = np.random.normal(5, 10)