Fix Issue ACCLAB/dabestr#48 : bug in calculating 95CIs for paired diffs

Author: josesho

dabest/_stats_tools/confint_1group.py

@@ -2,8 +2,28 @@
 # -*-coding: utf-8 -*-
 # Author: Joses Ho
 # Email : joseshowh@gmail.com
+"""
+A range of functions to compute bootstraps for a single sample.
+"""
 
+def create_bootstrap_indexes(array, resamples=5000, random_seed=12345):
+    """Given an array-like, returns a generator of bootstrap indexes
+    to be used for resampling.
+    """
+    import numpy as np
+    
+    # Set seed.
+    np.random.seed(random_seed)
+    
+    indexes = range(0, len(array))
 
+    out = (np.random.choice(indexes, len(indexes), replace=True)
+            for i in range(0, resamples))
+    
+    # Reset seed
+    np.random.seed()
+    
+    return out
 
 def compute_1group_jackknife(x, func, *args, **kwargs):
     """
@@ -20,22 +40,6 @@ def compute_1group_jackknife(x, func, *args, **kwargs):
 def compute_1group_acceleration(jack_dist):
     from . import confint_2group_diff as ci_2g
     return ci_2g._calc_accel(jack_dist)
-    
-    
-    
-def _create_bootstrap_indexes(array, resamples=5000):
-    """Given an array-like, returns a generator of bootstrap indexes
-    to be used for resampling.
-    """
-    import numpy as np
-    
-    indexes = range(0, len(array))
-
-    out = (np.random.choice(indexes, len(indexes), replace=True)
-            for i in range(0, resamples))
-    
-    return out
-
 
 
 
@@ -49,8 +53,9 @@ def compute_1group_bootstraps(x, func, resamples=5000, random_seed=12345,
     np.random.seed(random_seed)
 
     # Create bootstrap indexes.
-    boot_indexes = _create_bootstrap_indexes(x, resamples)
-    
+    boot_indexes = create_bootstrap_indexes(x, resamples=resamples,
+                                            random_seed=random_seed)
+
     out = [func(x[b], *args, **kwargs) for b in boot_indexes]
 
     del boot_indexes
@@ -123,11 +128,13 @@ def summary_ci_1group(x, func, resamples=5000, alpha=0.05, random_seed=12345,
     from . import confint_2group_diff as ci2g
     from numpy import sort as npsort
 
-    boots = compute_1group_bootstraps(x, func, resamples, random_seed) 
-    bias  = compute_1group_bias_correction(x, boots, func)
+    boots = compute_1group_bootstraps(x, func, resamples=resamples,
+                                      random_seed=random_seed,
+                                      *args, **kwargs)
+    bias = compute_1group_bias_correction(x, boots, func)
 
-    jk = compute_1group_jackknife(x, func)
-    accel = ci2g._calc_accel(jk)
+    jk = compute_1group_jackknife(x, func, *args, **kwargs)
+    accel = compute_1group_acceleration(jk)
     del jk
 
     ci_idx = ci2g.compute_interval_limits(bias, accel, resamples, alpha)

dabest/_stats_tools/confint_2group_diff.py

@@ -2,8 +2,10 @@
 # -*-coding: utf-8 -*-
 # Author: Joses Ho
 # Email : joseshowh@gmail.com
-
-
+"""
+A range of functions to compute bootstraps for the mean difference 
+between two groups.
+"""
 
 def create_jackknife_indexes(data):
     """
@@ -103,9 +105,34 @@ def _calc_accel(jack_dist):
 
 
 
+# def compute_bootstrapped_diff(x0, x1, is_paired, effect_size,
+#                                 resamples=5000, random_seed=12345):
+#     """Bootstraps the effect_size for 2 groups."""
+#     from . import effsize as __es
+#     import numpy as np
+# 
+#     np.random.seed(random_seed)
+# 
+#     out = np.repeat(np.nan, resamples)
+#     x0_len = len(x0)
+#     x1_len = len(x1)
+# 
+#     for i in range(int(resamples)):
+#         x0_boot = np.random.choice(x0, x0_len, replace=True)
+#         x1_boot = np.random.choice(x1, x1_len, replace=True)
+#         out[i] = __es.two_group_difference(x0_boot, x1_boot,
+#                                           is_paired, effect_size)
+# 
+#     # reset seed
+#     np.random.seed()
+# 
+#     return out
+
+
 def compute_bootstrapped_diff(x0, x1, is_paired, effect_size,
                                 resamples=5000, random_seed=12345):
     """Bootstraps the effect_size for 2 groups."""
+    
     from . import effsize as __es
     import numpy as np
 
@@ -114,11 +141,20 @@ def compute_bootstrapped_diff(x0, x1, is_paired, effect_size,
     out = np.repeat(np.nan, resamples)
     x0_len = len(x0)
     x1_len = len(x1)
-
+    
     for i in range(int(resamples)):
-        x0_boot = np.random.choice(x0, x0_len, replace=True)
-        x1_boot = np.random.choice(x1, x1_len, replace=True)
-        out[i] = __es.two_group_difference(x0_boot, x1_boot,
+        
+        if is_paired:
+            if x0_len != x1_len:
+                raise ValueError("The two arrays do not have the same length.")
+            random_idx = np.random.choice(x0_len, x0_len, replace=True)
+            x0_sample = x0[random_idx]
+            x1_sample = x1[random_idx]
+        else:
+            x0_sample = np.random.choice(x0, x0_len, replace=True)
+            x1_sample = np.random.choice(x1, x1_len, replace=True)
+            
+        out[i] = __es.two_group_difference(x0_sample, x1_sample,
                                           is_paired, effect_size)
 
     # reset seed
@@ -128,6 +164,7 @@ def compute_bootstrapped_diff(x0, x1, is_paired, effect_size,
 
 
 
+
 def compute_meandiff_bias_correction(bootstraps, effsize):
     """
     Computes the bias correction required for the BCa method