Pull Request for #19: Confidence interval for SK distance (#20)

* Methods for #19

* #19 Add get_sk_ci()

* #19 Adding unit test for get_sk_ci()

* Ceiling for #20

* #19 Apply suggestions from jms review

Co-authored-by: jmshapir <[email protected]>

* #19 Capping ub at 1

---------

Co-authored-by: jmshapir <[email protected]>

Author: rcalvo12

BootstrapReport.lyx

@@ -294,6 +294,56 @@ Here
 \end_layout
 
 \end_deeper
+\begin_layout Section*
+Confidence interval for SK distance
+\end_layout
+
+\begin_layout Itemize
+For a value 
+\begin_inset Formula $x\ge0$
+\end_inset
+
+ let 
+\begin_inset Formula 
+\begin{align*}
+\overline{x}_{\alpha} & =\min\left\{ x+\sqrt{\frac{\ln\left(2/\alpha\right)}{2J}},1\right\} \\
+\underline{x}_{\alpha} & =\max\left\{ x-\sqrt{\frac{\ln\left(2/\alpha\right)}{2J}},0\right\} 
+\end{align*}
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Itemize
+We calculate the 
+\begin_inset CommandInset href
+LatexCommand href
+name "DKW"
+target "https://en.wikipedia.org/wiki/Dvoretzky%E2%80%93Kiefer%E2%80%93Wolfowitz_inequality"
+literal "false"
+
+\end_inset
+
+ confidence interval whose bounds are
+\begin_inset Formula 
+\[
+\min\left\{ \overline{SK}_{\alpha}^{+}\left(N\left(\mu,\sigma^{2}\right),\hat{\eta}\right)+\overline{SK}_{\alpha}^{-}\left(N\left(\mu,\sigma^{2}\right),\hat{\eta}\right),1\right\} 
+\]
+
+\end_inset
+
+and
+\begin_inset Formula 
+\[
+\underline{SK}_{\alpha}^{+}\left(N\left(\mu,\sigma^{2}\right),\hat{\eta}\right)+\underline{SK}_{\alpha}^{-}\left(N\left(\mu,\sigma^{2}\right),\hat{\eta}\right)
+\]
+
+\end_inset
+
+
+\end_layout
+
 \begin_layout Section*
 Number of changes in direction in the difference in CDFs
 \end_layout

src/BootstrapReport/helpers.py

@@ -138,6 +138,22 @@ def get_sk_dist(rep, normal, sep = False):
     else:
         return -neg_dist + pos_dist
 
+def get_sk_ci(num_replicates, neg_dist, pos_dist, alpha = 0.05):
+    '''calculates confidence interval for sk distance
+    :param num_replicates: number of replicates
+    :param neg_dist: maximum negative distance 
+    :param post_dist: maximum positive distance
+    :alpha: 1 - alpha = confidence level for confidence bands
+    '''
+    sk_ci_lb = max(pos_dist - np.sqrt(np.log(2/alpha)/(2 * num_replicates)), 0) + \
+               max(neg_dist - np.sqrt(np.log(2/alpha)/(2 * num_replicates)), 0)
+                   
+    sk_ci_ub = min(pos_dist + np.sqrt(np.log(2/alpha)/(2 * num_replicates)), 1) + \
+               min(neg_dist + np.sqrt(np.log(2/alpha)/(2 * num_replicates)), 1)
+    sk_ci_ub = min(sk_ci_ub, 1)
+    
+    return sk_ci_lb, sk_ci_ub
+
 def select_bandwidth(rot_se, max_grid, min_grid, num_gridpoints, 
                      estimate, se, num_sets, num_replicates, lbound, rbound, first_seed,
                      num_expansions, max_subdivisions):

src/BootstrapReport/main.py

@@ -25,7 +25,9 @@ def __init__(self, estimate, se, replicates):
         self.best_bandwidth_value = np.NaN
         self.bias_value = None
         self.crossings = None
-        self.sk_dist = helpers.get_sk_dist(self.replicates, norm(loc = estimate, scale = se), sep = False)
+        self.sk_neg, self.sk_pos = helpers.get_sk_dist(self.replicates, norm(loc = self.estimate, scale = self.se), sep = True)
+        self.sk_dist = self.sk_neg + self.sk_pos
+        self.sk_ci_lb, self.sk_ci_ub = helpers.get_sk_ci(len(self.replicates), self.sk_neg, self.sk_pos)
 
     def get_bias_corrected_tvd(self, num_gridpoints=10, num_sets=100, min_grid=-1, max_grid=1,
                                bounds_of_integration=np.inf, first_seed=42, second_seed=11, 

tests/test_BootstrapReport.py

@@ -17,6 +17,15 @@ def test_sk_distance():
     assert test.neg_dist <= 1 and test.neg_dist >= 0
     assert test.pos_dist <= 1 and test.pos_dist >= 0
 
+def test_sk_ci():
+    test_replicates = pd.read_csv('examples/gamma_replicates.csv')['replicate_value'].values
+    estimate, standard_error = 0, 1
+    test = ObjectOfInterest(estimate = estimate, se = standard_error, replicates = test_replicates)
+        
+    assert test.sk_ci_lb <= test.sk_dist and test.sk_ci_ub >= test.sk_dist
+    assert test.sk_ci_ub <= 1
+    assert test.sk_ci_lb >= 0
+
 def test_density_plot():
     test_replicates = pd.read_csv('examples/gamma_replicates.csv')['replicate_value'].values
     estimate, standard_error = 0, 1