<jupyter_start><jupyter_text>Libraries<jupyter_code># default libraries
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
# plotting configuration
sns.set(style = 'darkgrid',
        font_scale = 1.4,
        rc = {"lines.linewidth": 2,
          "figure.figsize": (12, 8)})

sns.set_palette('Set2')
# load implemented packages
from tools import get_alpha_quantile_student, jackknife_cookie_bucket_mean,\
jackknife_cookie_bucket_quantile<jupyter_output><empty_output><jupyter_text>Test computation of quantiles<jupyter_code>assert np.abs(get_alpha_quantile_student(ci_level = 0.9) - 1.73) < 0.01
assert np.abs(get_alpha_quantile_student(ci_level = 0.95) - 2.09) < 0.01
assert np.abs(get_alpha_quantile_student(ci_level = 0.99) - 2.86) < 0.01<jupyter_output><empty_output><jupyter_text>Jackknife cookie bucket test for means  DescriptionLet $M$ be some metric of interest. Let subscript 0 correspond to the control group and 1 - to the treatment group. Then the Jackknife Cookie Bucket test is defined as follows:* Data is split randomly in $N$ buckets (we use $N=20$)* Let $M_{0,-j}$ is the mean value of the metric observed in the control $M_{0}$ across all buckets **except** $j$. $M_{1,-j}$ is defined similarly. The second index omited in case one refers to the total change in the metric value across all buckets* Overall percent change in all buckets is defined as: $$ PC= 100\times \left(  \frac{M_1}{M_0}-1 \right) $$* Overall percent change in all buckets, except $j$, is defined as: $$ PC_{-j}= 100\times \left(  \frac{M_{1,-j}}{M_{0,-j}}-1 \right)$$* Then define: $$ D_{-j} = N*PC-(N-1)*PC_{-j} $$ and $$\overline{D}=\frac{1}{N}\sum_{j=1}^N D_{-j}, \qquad SE_{D}^2=\frac{1}{N-1}\sum_{j=1}^{N}(\overline{D}-D_{-j})^2$$ * The final confidence interval for performing the test is defined as: $$ \left[ \overline{D}- q_{N-1,\alpha/2}\cdot SE_{D}/\sqrt{N};\ \overline{D}+ q_{N-1,\alpha}\cdot SE_{D}/\sqrt{N} \right] $$where $q_{N-1,\alpha}$ corresponds to the corresponding quantile of the Student t-distribution with $N-1$ degrees of freedom. If $N=20$, then $q_{19,0.9}\approx 1.73, q_{19,0.95}\approx 2.09, q_{19,0.99}\approx 2.86$  Test for the Normal case with the same mean and same variance (Null)<jupyter_code>mean_metric_control = 50
mean_metric_treatement = 50
num_of_clients = 5000
num_of_simulations = 1000
# 1 will correspond to rejecting the null
coverage_results = list()

for cur_sim in range(num_of_simulations):
  # simulate data
  metric_treatment = np.random.normal(loc = mean_metric_treatement, scale=5, size=[num_of_clients, 1])
  metric_control = np.random.normal(loc = mean_metric_control, scale=5, size=[num_of_clients, 1])
  client_id = np.arange(1, num_of_clients + 1).reshape(num_of_clients, 1) 

  # create required data frames
  df_treatment = pd.DataFrame(data = np.concatenate([client_id, metric_treatment],axis=1), columns=['Client ID', 'Metric Value'] )
  df_control = pd.DataFrame(data = np.concatenate([client_id, metric_control],axis=1), columns=['Client ID', 'Metric Value'] ) 

  # perform test
  ci_center, ci_width = jackknife_cookie_bucket_mean(df_treatment, df_control, 20, 0.95)

  if (ci_center - ci_width < 0) and  (ci_center + ci_width > 0):
    coverage_results += [0]
  else:
    coverage_results += [1]

print('Rejection rate: {:2.2%}'.format(np.mean(coverage_results)))
ci_center
ci_width
plt.hist(metric_treatment, bins=50, label='Treatment group', color='darkslateblue')
plt.axvline(np.mean(metric_treatment), c='darkslateblue')
plt.hist(metric_control, bins=50, alpha=0.4, label='Control group', color='b')
plt.axvline(np.mean(metric_control), c='b')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Test for the Normal case with the same mean and different variance (Null)<jupyter_code># 1 will correspond to rejecting the null
coverage_results = list()

for cur_sim in range(num_of_simulations):
  # simulate data
  metric_treatment = np.random.normal(loc = mean_metric_treatement, scale=5, size=[num_of_clients, 1])
  metric_control = np.random.normal(loc = mean_metric_control, scale=10, size=[num_of_clients, 1])
  client_id = np.arange(1, num_of_clients + 1).reshape(num_of_clients, 1) 

  # create required data frames
  df_treatment = pd.DataFrame(data = np.concatenate([client_id, metric_treatment],axis=1), columns=['Client ID', 'Metric Value'] )
  df_control = pd.DataFrame(data = np.concatenate([client_id, metric_control],axis=1), columns=['Client ID', 'Metric Value'] ) 

  # perform test
  ci_center, ci_width = jackknife_cookie_bucket_mean(df_treatment, df_control, 20, 0.95)

  if (ci_center - ci_width < 0) and  (ci_center + ci_width > 0):
    coverage_results += [0]
  else:
    coverage_results += [1]

print('Rejection rate: {:2.2%}'.format(np.mean(coverage_results)))
plt.hist(metric_treatment, bins=50, label='Treatment group', color='darkslateblue')
plt.axvline(np.mean(metric_treatment), c='darkslateblue')
plt.hist(metric_control, bins=50, alpha=0.4, label='Control group', color='b')
plt.axvline(np.mean(metric_control), c='b')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Test for the Normal case with the different mean and same variance (Alternative)<jupyter_code>mean_metric_control = 50
mean_metric_treatement = 45
# 1 will correspond to rejecting the null
coverage_results = list()

for cur_sim in range(num_of_simulations):
  # simulate data
  metric_treatment = np.random.normal(loc = mean_metric_treatement, scale=5, size=[num_of_clients, 1])
  metric_control = np.random.normal(loc = mean_metric_control, scale=5, size=[num_of_clients, 1])
  client_id = np.arange(1, num_of_clients + 1).reshape(num_of_clients, 1) 

  # create required data frames
  df_treatment = pd.DataFrame(data = np.concatenate([client_id, metric_treatment],axis=1), columns=['Client ID', 'Metric Value'] )
  df_control = pd.DataFrame(data = np.concatenate([client_id, metric_control],axis=1), columns=['Client ID', 'Metric Value'] ) 

  # perform test
  ci_center, ci_width = jackknife_cookie_bucket_mean(df_treatment, df_control, 20, 0.95)

  if (ci_center-ci_width < 0) and  (ci_center+ci_width > 0):
    coverage_results += [0]
  else:
    coverage_results += [1]

print('Rejection rate: {:2.2%}'.format(np.mean(coverage_results)))
plt.hist(metric_treatment, bins=50, label='Treatment group', color='darkslateblue')
plt.axvline(np.mean(metric_treatment), c='darkslateblue')
plt.hist(metric_control, bins=50, alpha=0.4, label='Control group', color='b')
plt.axvline(np.mean(metric_control), c='b')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Test for the Normal case with the different mean and different variance (Alternative)<jupyter_code># 1 will correspond to rejecting the null
coverage_results = list()

for cur_sim in range(num_of_simulations):
  # simulate data
  metric_treatment = np.random.normal(loc = mean_metric_treatement, scale=5, size=[num_of_clients, 1])
  metric_control = np.random.normal(loc = mean_metric_control, scale=10, size=[num_of_clients, 1])
  client_id = np.arange(1, num_of_clients + 1).reshape(num_of_clients, 1) 

  # create required data frames
  df_treatment = pd.DataFrame(data = np.concatenate([client_id, metric_treatment],axis=1), columns=['Client ID', 'Metric Value'] )
  df_control = pd.DataFrame(data = np.concatenate([client_id, metric_control],axis=1), columns=['Client ID', 'Metric Value'] ) 

  # perform test
  ci_center, ci_width = jackknife_cookie_bucket_mean(df_treatment, df_control, 20, 0.95)

  if (ci_center-ci_width < 0) and  (ci_center+ci_width > 0):
    coverage_results += [0]
  else:
    coverage_results += [1]

print('Rejection rate: {:2.2%}'.format(np.mean(coverage_results)))<jupyter_output>Rejection rate: 100.00%<jupyter_text>Jackknife Cookie bucket test for quantiles  DescriptionLet $M$ be some metric of interest. Let subscript 0 correspond to the control group and 1 - to the treatment group. Then the Jackknife Cookie Bucket test is defined as follows:* Data is split randomly in $N$ buckets (we use $N=20$)* Let $M_{0,-j}$ is the chosen quantile of the metric observed in the control $M_{0}$ across all buckets **except** $j$. $M_{1,-j}$ is defined similarly. The second index omited in case one refers to the total change in the metric value across all buckets* Overall percent change in all buckets is defined as: $$ PC= 100\times \left(  \frac{M_1}{M_0}-1 \right) $$* Overall percent change in all buckets except $j$-th is defined as: $$ PC_{-j}= 100\times \left(  \frac{M_{1,-j}}{M_{0,-j}}-1 \right)$$* Then define the pseudovalue of the percent change for bucket $j$: $$ps_{-j} =  N \cdot PC - (N-1)\cdot PC_{-j}$$* Define: $$ps = \sum_{i=1}^Nps_{-j}, \quad SE_{ps}^2 =  \frac{1}{N-1}\sum_{j=1}^{N}(ps-ps_{-j})^2$$* The final confidence interval for performing the test is defined as: $$ \left[ ps- q_{N-1,\alpha/2}\cdot SE_{ps}/\sqrt{N};\ ps+ q_{N-1,\alpha}\cdot SE_{ps}/\sqrt{N} \right] $$where $q_{N-1,\alpha}$ corresponds to the corresponding quantile of the Student t-distribution with $N-1$ degrees of freedom. If $N=20$, then $q_{19,0.9}\approx 1.73, q_{19,0.95}\approx 2.09, q_{19,0.99}\approx 2.86$  Problems with Jackknife for quantilesGiven a sample $X=(X_1, \dots, X_n)$, we estimates parameter $\theta$ using an estimator $\phi(X)=\phi_n(X_1,\dots,X_n)$. $i$-th pseudovalue of $\phi_n(X)$ is:$$ps_i(X)=n\phi_n(X_1,\dots,X_n)-(n-1)\phi_{n-1}((X_1,\dots,X_n)_{[i]}) = \phi_n(X)+(n-1)(\phi_n(X)-\phi_{n-1}(X_{[i]}))$$where $X_{[i]}=(X_1,\dots,X_n)_{[i]}$ means that $i$-th value in a sample has been deleted. Then:$$ ps(X)=\frac{1}{n}\sum_{i=1}^nps_i(X), \quad V_{ps}(X)=\frac{1}{n-1}\sum_{i=1}^n(ps_i(X)-ps(X))^2 $$with the consequent construction of the CI for $\theta$. Generalization, which is used actually here is delete-$k$ or block jackknife.Problems arise when the estimator is too discontinuous as a function of $X_i$ or if $\phi(X)$ depends on one or a few values in $X$, which is actually the case for quantiles.**Example** Assume that $\phi_n(X)$ is the sample median. Assume that $X_i$'s are distinct and $n=2m$ is even. Then:* there are two values $X_{(m)}$ and $X_{(m+1)}$ such that:$$\phi_{n-1}(X_{[i]}) = X_{(m+1)}, \quad i\leq m$$and  $$\phi_{n-1}(X_{[i]}) = X_{(m)}, \quad i\geq m+1$$Hence,$$ps(X) = \frac{1}{2}(X_{(m)}+X_{(m+1)}) = \phi_n(X)$$Thus,$$\quad V_{ps}(X)=\frac{n(n-1)}{4}(X_{(m)}-X_{(m+1)})^2 \rightarrow \sigma^2 (\chi^2_2/2)^2$$where $\sigma^2$ is the asymptotic variance of $n^{1/2}(\phi_n(X)- med(F_X))$ meaning that variance estimator is inconsistent. For delete-$k$ Jackknife k should satisfy $n^{1/2}/k\rightarrow 0$ and $n-k\rightarrow\infty$ to achieve consistency.In our case, we have 20 buckets meaning that we take $k=n/20$  Test for the Normal case with the same mean and same variance (Null)<jupyter_code>mean_metric_control = 50
mean_metric_treatement = 50
num_of_clients=5000
# 1 will correspond to rejecting the null
coverage_results = list()
tracking_samples = list()
for cur_sim in range(num_of_simulations):
  # simulate data
  metric_treatment = np.random.normal(loc = mean_metric_treatement, scale=5, size=[num_of_clients, 1])
  metric_control = np.random.normal(loc = mean_metric_control, scale=5, size=[num_of_clients, 1])
  client_id = np.arange(1, num_of_clients + 1).reshape(num_of_clients, 1) 

  # create required data frames
  df_treatment = pd.DataFrame(data = np.concatenate([client_id, metric_treatment],axis=1),
    columns=['Client ID', 'Metric Value'] )
  df_control = pd.DataFrame(data = np.concatenate([client_id, metric_control],axis=1),
    columns=['Client ID', 'Metric Value'] ) 

  # perform test
  ci_center, ci_width = jackknife_cookie_bucket_quantile(df_treatment, df_control, 20, 0.95,
    quantile_to_test=0.5)

  if (ci_center-ci_width < 0) and  (ci_center+ci_width > 0):
    coverage_results += [0]
  else:
    coverage_results += [1]
    tracking_samples+= [[[metric_treatment],[metric_control]]]

print('Rejection rate: {:2.2%}'.format(np.mean(coverage_results)))
len(tracking_samples)
tracking_samples[cur_sample][0][0]
cur_sample=5
plt.hist(tracking_samples[cur_sample][0][0], bins=50, label='Treatment group', color='darkslateblue')
plt.axvline(np.quantile(tracking_samples[cur_sample][0][0], q=0.5), c='darkslateblue')
plt.hist(tracking_samples[cur_sample][1][0], bins=50, alpha=0.4, label='Control group', color='b')
plt.axvline(np.quantile(tracking_samples[cur_sample][1][0], q=0.5), c='b')
plt.legend()
plt.show()
tracking_samples[cur_sample][0][0].mean()
np.quantile(tracking_samples[cur_sample][0][0],q=0.5)
tracking_samples[cur_sample][1][0].mean()
np.quantile(tracking_samples[cur_sample][1][0],q=0.5)<jupyter_output><empty_output><jupyter_text>Test for the Normal case with the same mean and different variance (Null)<jupyter_code># 1 will correspond to rejecting the null
coverage_results = list()

for cur_sim in range(num_of_simulations):
  # simulate data
  metric_treatment = np.random.normal(loc = mean_metric_treatement, scale=5, size=[num_of_clients, 1])
  metric_control = np.random.normal(loc = mean_metric_control, scale=10, size=[num_of_clients, 1])
  client_id = np.arange(1, num_of_clients + 1).reshape(num_of_clients, 1) 

  # create required data frames
  df_treatment = pd.DataFrame(data = np.concatenate([client_id, metric_treatment], axis=1),
    columns=['Client ID', 'Metric Value'] )
  df_control = pd.DataFrame(data = np.concatenate([client_id, metric_control], axis=1),
    columns=['Client ID', 'Metric Value'] ) 

  # perform test
  ci_center, ci_width = jackknife_cookie_bucket_quantile(df_treatment, df_control, 20, 0.95,
    quantile_to_test=0.5)

  if (ci_center-ci_width < 0) and  (ci_center+ci_width > 0):
    coverage_results += [0]
  else:
    coverage_results += [1]

print('Rejection rate: {:2.2%}'.format(np.mean(coverage_results)))<jupyter_output>Rejection rate: 17.20%<jupyter_text>Test for the Normal case with different mean and same variance (Alternative)<jupyter_code>mean_metric_control = 50
mean_metric_treatement = 45
# 1 will correspond to rejecting the null
coverage_results = list()

for cur_sim in range(num_of_simulations):
  # simulate data
  metric_treatment = np.random.normal(loc = mean_metric_treatement, scale=5, size=[num_of_clients, 1])
  metric_control = np.random.normal(loc = mean_metric_control, scale=5, size=[num_of_clients, 1])
  client_id = np.arange(1, num_of_clients + 1).reshape(num_of_clients, 1) 

  # create required data frames
  df_treatment = pd.DataFrame(data = np.concatenate([client_id, metric_treatment],axis=1),
    columns=['Client ID', 'Metric Value'] )
  df_control = pd.DataFrame(data = np.concatenate([client_id, metric_control],axis=1),
    columns=['Client ID', 'Metric Value'] ) 

  # perform test
  ci_center, ci_width = jackknife_cookie_bucket_quantile(df_treatment, df_control, 20, 0.95,
    quantile_to_test=0.5)

  if (ci_center-ci_width < 0) and  (ci_center+ci_width > 0):
    coverage_results += [0]
  else:
    coverage_results += [1]

print('Rejection rate: {:2.2%}'.format(np.mean(coverage_results)))<jupyter_output>Rejection rate: 100.00%<jupyter_text>Test for the Normal case with different mean and different variance (Alternative)<jupyter_code># 1 will correspond to rejecting the null
coverage_results = list()

for cur_sim in range(num_of_simulations):
  # simulate data
  metric_treatment = np.random.normal(loc = mean_metric_treatement, scale=5, size=[num_of_clients, 1])
  metric_control = np.random.normal(loc = mean_metric_control, scale=10, size=[num_of_clients, 1])
  client_id = np.arange(1, num_of_clients + 1).reshape(num_of_clients, 1) 

  # create required data frames
  df_treatment = pd.DataFrame(data = np.concatenate([client_id, metric_treatment],axis=1),
    columns=['Client ID', 'Metric Value'] )
  df_control = pd.DataFrame(data = np.concatenate([client_id, metric_control],axis=1),
    columns=['Client ID', 'Metric Value'] ) 

  # perform test
  ci_center, ci_width = jackknife_cookie_bucket_quantile(df_treatment, df_control, 20, 0.95,
    quantile_to_test = 0.5)

  if (ci_center-ci_width < 0) and  (ci_center+ci_width > 0):
    coverage_results += [0]
  else:
    coverage_results += [1]

print('Rejection rate: {:2.2%}'.format(np.mean(coverage_results)))<jupyter_output>Rejection rate: 100.00%<jupyter_text>Extreme quantiles in case of the same mean and variance (Null)  More observations<jupyter_code>mean_metric_control = 50
mean_metric_treatement = 50
num_of_clients = 20000
# 1 will correspond to rejecting the null
coverage_results = list()

for cur_sim in range(num_of_simulations):
  # simulate data
  metric_treatment = np.random.normal(loc = mean_metric_treatement, scale=5, size=[num_of_clients, 1])
  metric_control = np.random.normal(loc = mean_metric_control, scale=5, size=[num_of_clients, 1])
  client_id = np.arange(1, num_of_clients + 1).reshape(num_of_clients, 1) 

  # create required data frames
  df_treatment = pd.DataFrame(data = np.concatenate([client_id, metric_treatment],axis=1),
    columns=['Client ID', 'Metric Value'] )
  df_control = pd.DataFrame(data = np.concatenate([client_id, metric_control],axis=1),
    columns=['Client ID', 'Metric Value'] ) 

  # perform test
  ci_center, ci_width = jackknife_cookie_bucket_quantile(df_treatment, df_control, 20, 0.95,
    quantile_to_test=0.95)

  if (ci_center-ci_width < 0) and  (ci_center+ci_width > 0):
    coverage_results += [0]
  else:
    coverage_results += [1]

print('Rejection rate: {:2.2%}'.format(np.mean(coverage_results)))
plt.hist(metric_treatment, bins=50, label='Treatment group', color='darkslateblue')
plt.axvline(np.quantile(metric_treatment, q=0.95), c='darkslateblue')
plt.hist(metric_control, bins=50, alpha=0.4, label='Control group', color='b')
plt.axvline(np.quantile(metric_control, q=0.95), c='b')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Less observations<jupyter_code>mean_metric_control = 50
mean_metric_treatement = 50
num_of_clients = 1000
# 1 will correspond to rejecting the null
coverage_results = list()

for cur_sim in range(num_of_simulations):
  # simulate data
  metric_treatment = np.random.normal(loc = mean_metric_treatement, scale=5, size=[num_of_clients, 1])
  metric_control = np.random.normal(loc = mean_metric_control, scale=5, size=[num_of_clients, 1])
  client_id = np.arange(1, num_of_clients + 1).reshape(num_of_clients, 1) 

  # create required data frames
  df_treatment = pd.DataFrame(data = np.concatenate([client_id, metric_treatment],axis=1),
    columns=['Client ID', 'Metric Value'] )
  df_control = pd.DataFrame(data = np.concatenate([client_id, metric_control],axis=1),
    columns=['Client ID', 'Metric Value'] ) 

  # perform test
  ci_center, ci_width = jackknife_cookie_bucket_quantile(df_treatment, df_control, 20, 0.95,
    quantile_to_test=0.95)

  if (ci_center-ci_width < 0) and  (ci_center+ci_width > 0):
    coverage_results += [0]
  else:
    coverage_results += [1]

print('Rejection rate: {:2.2%}'.format(np.mean(coverage_results)))
plt.hist(metric_treatment, bins=50, label='Treatment group', color='darkslateblue')
plt.axvline(np.quantile(metric_treatment, q=0.95), c='darkslateblue')
plt.hist(metric_control, bins=50, alpha=0.4, label='Control group', color='b')
plt.axvline(np.quantile(metric_control, q=0.95), color='b')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Extreme quantile in case of the same mean, but different variance (Alternative)<jupyter_code>mean_metric_control = 50
mean_metric_treatement = 50
num_of_clients = 10000
# 1 will correspond to rejecting the null
coverage_results = list()

for cur_sim in range(num_of_simulations):
  # simulate data
  metric_treatment = np.random.normal(loc = mean_metric_treatement, scale=5, size=[num_of_clients, 1])
  metric_control = np.random.normal(loc = mean_metric_control, scale=10, size=[num_of_clients, 1])
  client_id = np.arange(1, num_of_clients + 1).reshape(num_of_clients, 1) 

  # create required data frames
  df_treatment = pd.DataFrame(data = np.concatenate([client_id, metric_treatment],axis=1),
    columns=['Client ID', 'Metric Value'] )
  df_control = pd.DataFrame(data = np.concatenate([client_id, metric_control],axis=1),
    columns=['Client ID', 'Metric Value'] ) 

  # perform test
  ci_center, ci_width = jackknife_cookie_bucket_quantile(df_treatment, df_control, 20, 0.95,
    quantile_to_test=0.95)

  if (ci_center-ci_width < 0) and  (ci_center+ci_width > 0):
    coverage_results += [0]
  else:
    coverage_results += [1]

print('Rejection rate: {:2.2%}'.format(np.mean(coverage_results)))
plt.hist(metric_treatment, bins=50, label='Treatment group', color='darkslateblue')
plt.axvline(np.quantile(metric_treatment, q=0.95), c='darkslateblue')
plt.hist(metric_control, bins=50, alpha=0.4, label='Control group', color='b')
plt.axvline(np.quantile(metric_control, q=0.95), color='b')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Testing for median and for the mean<jupyter_code>mean_metric_control = 50
mean_metric_treatement = 50
num_of_clients=20000
# 1 will correspond to rejecting the null
coverage_results_mean = list()
coverage_results_median = list()
tracking_samples = list()
n_buckets = 20
width_median = list()
width_mean = list()


for cur_sim in range(num_of_simulations):
  # simulate data
  metric_treatment = np.random.normal(loc = mean_metric_treatement, scale=5, size=[num_of_clients, 1])
  metric_control = np.random.normal(loc = mean_metric_control, scale=5, size=[num_of_clients, 1])
  client_id = np.arange(1, num_of_clients + 1).reshape(num_of_clients, 1) 

  # create required data frames
  df_treatment = pd.DataFrame(data = np.concatenate([client_id, metric_treatment],axis=1),
    columns=['Client ID', 'Metric Value'] )
  df_control = pd.DataFrame(data = np.concatenate([client_id, metric_control],axis=1),
    columns=['Client ID', 'Metric Value'] ) 

  # perform test for median
  ci_center_med, ci_width_med = jackknife_cookie_bucket_quantile(df_treatment, df_control, n_buckets, 0.95,
    quantile_to_test=0.5)

  ci_center_mean, ci_width_mean = jackknife_cookie_bucket_mean(df_treatment, df_control, n_buckets, 0.95)


  if (ci_center_mean-ci_width_mean < 0) and  (ci_center_mean+ci_width_mean > 0):
    test_mean = 0
  else:
    test_mean = 1

  if (ci_center_med-ci_width_med < 0) and  (ci_center_med+ci_width_med > 0):
    test_med = 0
  else:
    test_med = 1
  
  width_median+=[ci_width_med]
  width_mean+=[ci_width_mean]
  coverage_results_mean+=[test_mean]
  coverage_results_median+=[test_med]
  if test_mean==0 and test_med ==1:
    tracking_samples+= [[metric_treatment,metric_control]]

print('Rejection rate: {:2.2%}'.format(np.mean(coverage_results_mean)))
print('Rejection rate: {:2.2%}'.format(np.mean(coverage_results_median)))
plt.hist(width_median,label='Width of CI for the median', bins=30)
plt.hist(width_mean,label='Width of CI for the mean', bins=30)
plt.legend()
plt.show()
len(tracking_samples)
cur_sample=10
plt.hist(tracking_samples[cur_sample][0], bins=50, label='Treatment group', color='darkslateblue')
plt.axvline(np.quantile(tracking_samples[cur_sample][0], q=0.5), c='darkslateblue')
plt.hist(tracking_samples[cur_sample][1], bins=50, alpha=0.4, label='Control group', color='b')
plt.axvline(np.quantile(tracking_samples[cur_sample][1], q=0.5), c='b')
plt.legend()
plt.show()
tracking_samples[cur_sample][0].mean()
np.quantile(tracking_samples[cur_sample][0],q=0.5)
tracking_samples[cur_sample][1].mean()
np.quantile(tracking_samples[cur_sample][1],q=0.5)
mean_metric_control = 50
mean_metric_treatement = 50
num_of_clients=10000
# 1 will correspond to rejecting the null
coverage_results_mean = list()
coverage_results_median = list()
tracking_samples = list()
n_buckets = 20
width_median = list()
width_mean = list()


for cur_sim in range(num_of_simulations):
  # simulate data
  metric_treatment = np.random.normal(loc = mean_metric_treatement, scale=5, size=[num_of_clients, 1])
  metric_control = np.random.normal(loc = mean_metric_control, scale=5, size=[num_of_clients, 1])
  client_id = np.arange(1, num_of_clients + 1).reshape(num_of_clients, 1) 

  # create required data frames
  df_treatment = pd.DataFrame(data = np.concatenate([client_id, metric_treatment],axis=1),
    columns=['Client ID', 'Metric Value'] )
  df_control = pd.DataFrame(data = np.concatenate([client_id, metric_control],axis=1),
    columns=['Client ID', 'Metric Value'] ) 

  # perform test for median
  ci_center_med, ci_width_med = jackknife_cookie_bucket_quantile(df_treatment, df_control, n_buckets, 0.95,
    quantile_to_test=0.5)

  ci_center_mean, ci_width_mean = jackknife_cookie_bucket_mean(df_treatment, df_control, n_buckets, 0.95)


  if (ci_center_mean-ci_width_mean < 0) and  (ci_center_mean+ci_width_mean > 0):
    test_mean = 0
  else:
    test_mean = 1

  if (ci_center_med-ci_width_med < 0) and  (ci_center_med+ci_width_med > 0):
    test_med = 0
  else:
    test_med = 1
  
  width_median+=[ci_width_med]
  width_mean+=[ci_width_mean]
  coverage_results_mean+=[test_mean]
  coverage_results_median+=[test_med]
  if test_mean==0 and test_med ==1:
    tracking_samples+= [[metric_treatment,metric_control]]

print('Rejection rate: {:2.2%}'.format(np.mean(coverage_results_mean)))
print('Rejection rate: {:2.2%}'.format(np.mean(coverage_results_median)))
plt.hist(width_median,label='Width of CI for the median', bins=30)
plt.hist(width_mean,label='Width of CI for the mean', bins=30)
plt.legend()
plt.show()
df_treatment = pd.DataFrame(data = np.concatenate([client_id, tracking_samples[0][0]],axis=1),
    columns=['Client ID', 'Metric Value'] )
df_control = pd.DataFrame(data = np.concatenate([client_id, tracking_samples[0][0]],axis=1),
    columns=['Client ID', 'Metric Value'] ) 
df_treatment['Bucket id'] = (
        df_treatment['Client ID'] % n_buckets).astype('int')
df_control['Bucket id'] = (
    df_control['Client ID'] % n_buckets).astype('int')

# list of all bucket ids
bucket_ids = np.arange(n_buckets).tolist()


# for stacking the results for both treatment and control groups
treatment_group_mv = list()
control_group_mv = list()

for cur_bucket in bucket_ids:
    # leave one bucket out and compute the mean
    treatment_group_mv += [df_treatment.loc[df_treatment['Bucket id']
                                                != cur_bucket]['Metric Value'].mean()]
    control_group_mv += [df_control.loc[df_control['Bucket id']
                                            != cur_bucket]['Metric Value'].mean()]

# removing temporary column
#treatment_group = treatment_group.drop(axis=1, columns=['Bucket id'])
#control_group = control_group.drop(axis=1, columns=['Bucket id'])

# convert to array for further purposes
#treatment_group_mv = np.stack(treatment_group_mv)
#control_group_mv = np.stack(control_group_mv)
plt.hist(df_treatment.loc[df_treatment['Bucket id'] != 0]['Metric Value'], bins=40, alpha=0.7)
plt.hist(df_treatment.loc[df_treatment['Bucket id'] != 1]['Metric Value'], bins=40, alpha=0.7)
#plt.hist(df_treatment.loc[df_treatment['Bucket id'] != 2]['Metric Value'], bins=40, alpha=0.5)
plt.show()
for cur_bucket in range(20):
    plt.axvline(np.quantile(df_treatment.loc[df_treatment['Bucket id'] != cur_bucket]['Metric Value'], q=0.5))
plt.xlim([49.5,50.5])
for cur_bucket in range(20):
    plt.axvline(np.mean(df_treatment.loc[df_treatment['Bucket id'] != cur_bucket]['Metric Value']))
plt.xlim([49.5,50.5])
plt.axvline(np.mean(df_treatment.loc[df_treatment['Bucket id'] != 0]['Metric Value']))
plt.axvline(np.mean(df_treatment.loc[df_treatment['Bucket id'] != 4]['Metric Value']))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14347259430888423, 'fraction_numerical': 0.027622909275215407, 'mean_word_length': 5.900574569073195, 'pattern_counts': {'":': 2, '<': 71, '<?xml version=': 0, '>': 68, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '309721', 'license': 'Apache-2.0', 'n_tokens_mistral': 10610, 'n_tokens_neox': 9719, 'n_words': 2155, 'path': 'jackknife_cookie_bucket.ipynb', 'repo_name': 'sashapodkopaev/chrome-robust-hypothesis-testing', 'size': '27622'}