<jupyter_start><jupyter_text>Conditional Median Inference Simulations This notebook replicates the simulations seen in Section 4 of [1]. Specifically, it tests the performance of the general Conditional Quantile Inference algorithm (Algorithm 2) on 3 distributions with 4 different conformity scores. There are 4 sections:1.   **Distributions** - Describes the distribution class and covers the 3 relevant distributions2.   **Conformity Scores** - Describes the conformity score class and covers the 4 relevant scores3.   **Conditional Quantile Algorithm** - Contains an implementation of Algorithm 24.   **Testing Algorithm Performance** - Tests algorithm performance given each combination of distribution and conformity score[1] <NAME> and <NAME>, "Distribution-Free Conditional Median Inference." 2021.<jupyter_code>!pip install scikit-garden
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from datetime import datetime
from sklearn.utils import shuffle
import pickle
import tqdm
import os
import pandas as pd

# Used to save images to Google Drive
from google.colab import drive
drive.mount('/content/gdrive')

from skgarden import RandomForestQuantileRegressor, RandomForestRegressor


metadata_dir = '/content/gdrive/My Drive/Conditional Median Inference/Metadata/'
images_dir = '/content/gdrive/My Drive/Conditional Median Inference/Images/'<jupyter_output><empty_output><jupyter_text>Distributions This section contains different distributions that can be passed into the conditional quantile inference algorithm to test performance. Each distribution contains a sampling method and a conditional quantile calculation method.  Distribution 1<jupyter_code>class NormalCorrelatedDistribution:

  def __init__(self, d = 3, rho = 0, sigma = 1):
    assert d >= 3, f"Need d at least 3 but got {d}"
    self.mu = np.zeros(d)
    self.Sigma = (1. - rho) * np.identity(d) + rho * np.ones((d,d))
    self.resid_sigma = sigma

  def get_stdev(self, X, n):
    return 0.1 + 0.25 * np.square(np.linalg.norm(X, axis = 1)).reshape((n,1))

  def sample(self, n):
    X = np.random.multivariate_normal(mean = self.mu, cov = self.Sigma, size = n)
    Y_mean = np.square(X[:,[0]] + X[:,[1]]) - X[:,[2]]
    Y_resid = np.multiply(self.get_stdev(X, n), np.random.normal(loc = 0., scale = self.resid_sigma, size = (n,1)))
    Y = Y_mean + Y_resid
    return X, Y

  def get_quantile(self, X, q):
    Y_mean = np.square(X[:,[0]] + X[:,[1]]) - X[:,[2]]
    Y_quantile = self.get_stdev(X, X.shape[0]) * norm.ppf(q, loc = 0., scale = self.resid_sigma)
    return Y_mean + Y_quantile<jupyter_output><empty_output><jupyter_text>Distribution 2<jupyter_code>class IncreasingVarianceDistribution:

  def __init__(self, scale = np.pi, power = 1.):
    assert scale > 0, f"Need positive scale but got {scale}"
    assert power > 0, f"Need positive power but got {power}"
    self.scale = scale
    self.power = power

  def f(self, X):
    return 1 + np.abs(X) * np.square(np.sin(X))

  def sample(self, n):
    X = np.random.uniform(low = -1. * self.scale, high = self.scale, size = (n,1))
    Y = self.f(X) * (np.random.uniform(low = 0., high = 1., size = (n,1)) ** self.power)
    return X, Y

  def get_quantile(self, X, q):
    return self.f(X) * (np.full(shape = (X.shape[0], 1), fill_value = q) ** self.power)<jupyter_output><empty_output><jupyter_text>Distribution 3<jupyter_code>class ConfusingDistribution:

  def __init__(self, delta = 0.01, M = 10):
    assert delta > 0 and delta <= 0.25, f"Need delta in (0, 0.25] but got {delta}"
    assert isinstance(M, int) and M > 0, f"Need M to be a positive integer but got {M}"
    self.delta = delta
    self.M = M
    self.gamma = 1. / M

  def f(self, X):
    integral_part, fractional_part = np.divmod(X * self.M, 1.)
    return self.gamma * fractional_part - self.gamma / 2. - ((-1) ** (integral_part)) * (1. - self.gamma / 2.)

  def sample(self, n):
    X = np.random.uniform(low = -1., high = 1., size = (n,1))
    Y = self.f(X) * np.random.binomial(n = 1, p = 0.5 + 2 * self.delta, size = (n,1))
    return X, Y

  def get_quantile(self, X, q):
    nonzero_values = self.f(X)
    return (np.sign(nonzero_values) * (q - 0.5) >= -2. * self.delta) * nonzero_values<jupyter_output><empty_output><jupyter_text>Distribution Plots We plot the distributions, as well as the conditional medians, of Distributions 2 and 3. For Distribution 1, we cannot plot the full distribution due to dimensionality; we instead plot the distribution of the datapoints as well as the distribution of the conditional median.<jupyter_code>n = 10000
q = 0.5

bins = 150

save_figs = False
dpi = 2000
distribution1 = NormalCorrelatedDistribution(d = 10, rho = 0.25, sigma = 1)
X, Y = distribution1.sample(n)
Y_median = distribution1.get_quantile(X, q)

plt.hist(Y, density = True, bins = bins, alpha = 1, color = 'c', label = 'Datapoints')
plt.hist(Y_median, density = True, bins = bins, alpha = 0.5, color = 'r', label = 'Conditional Median')
plt.legend()
plt.xlabel('Y')
plt.ylabel('Density')
plt.title(r'$P_1$ Y Distribution and Conditional Median Distribution')

if save_figs:
  filename = "Distribution 1.pdf"
  plt.savefig(images_dir + filename)
plt.show()
datapoint_size = 5
datapoint_color = 'c'
datapoint_scale = 2.

median_linewidth = 2
median_color = 'k'

title_fontsize = 16
label_fontsize = 14
scale = 4 * np.pi
power = 0.25

distribution2 = IncreasingVarianceDistribution(scale = scale, power = power)
X, Y = distribution2.sample(n)

X_range = np.linspace(start = -1. * scale, stop = scale, num = n)
Y_median = distribution2.get_quantile(X_range.reshape((n,1)), 0.5)

plt.figure(figsize = (5,4))
plt.tight_layout()
plt.scatter(X,Y, s = datapoint_size, color = datapoint_color, label = 'Datapoints')
plt.plot(X_range, Y_median, color = median_color, linewidth = median_linewidth, label = 'Conditional Median')
plt.legend(markerscale= datapoint_scale)

plt.xlabel(r"$X$", fontsize = label_fontsize)
plt.ylabel(r"$Y$", fontsize = label_fontsize)
plt.title(r"$P_2$", fontsize = title_fontsize)

save_figs = True
if save_figs:
  filename = "Distribution 2.pdf"
  plt.savefig(images_dir + filename,bbox_inches='tight')

plt.show()
M = 25
delta = 0.0001
gamma = 1./M

distribution3 = ConfusingDistribution(delta = delta, M = M)
X, Y = distribution3.sample(n)

plt.figure(figsize = (5,4))
plt.tight_layout()
plt.scatter(X,Y, s = 8, color = datapoint_color, label = 'Datapoints')
plt.plot([0, 0], [-1, -1], color = median_color, linewidth = median_linewidth, label = 'Conditional Median')

for i in range(-M, M):
  if (i % 2) == 0:
    plt.plot([gamma * i, gamma * i + gamma], [-1, -1 + gamma], color = median_color, linewidth = median_linewidth)
  if (i % 2) == 1:
    plt.plot([gamma * i, gamma * i + gamma], [1 - gamma, 1], color = median_color, linewidth = median_linewidth)


plt.ylim((-1.1, 1.6))
plt.legend(markerscale= 1.5, loc = 'upper center')

plt.xlabel(r"$X$", fontsize = label_fontsize)
plt.ylabel(r"$Y$", fontsize = label_fontsize)
plt.title(r"$P_3$", fontsize = title_fontsize)

save_figs = True
if save_figs:
  filename = "Distribution 3.pdf"
  plt.savefig(images_dir + filename,bbox_inches='tight')

plt.show()<jupyter_output><empty_output><jupyter_text>Conformity Scores This section contains different conformity scores that can be passed into the conditional median inference algorithm. Each score contains a training method, a scoring method, and an inverse method. Note that we use a random forest predictor for all 4 conformity scores; for the first two, we train a random forest to predict the conditional mean, while the last two train a random forest to predict the conditional quantile.  Conformity Score 1<jupyter_code>class LinearDistanceConformityScore:

  def __init__(self, **kwargs):
    self.mu = RandomForestRegressor(**kwargs)

  def train(self, X, Y):
    self.mu.fit(X, Y)

  def score(self, X, Y):
    n = X.shape[0]
    prediction = self.mu.predict(X).reshape((n,1))
    return Y - prediction

  def inverse(self, X, score):
    n = X.shape[0]
    prediction = self.mu.predict(X).reshape((n,1))
    return prediction + score<jupyter_output><empty_output><jupyter_text>Conformity Score 2<jupyter_code>class NormalizedDistanceConformityScore: 

  def __init__(self, **kwargs):
    self.mu = RandomForestRegressor(**kwargs)

  def train(self, X, Y):
    self.mu.fit(X, Y)

  def score(self, X, Y):
    n = X.shape[0]
    prediction, std = self.mu.predict(X, return_std = True)
    prediction = prediction.reshape((n,1))
    std = std.reshape((n,1))
    return (Y - prediction) / std

  def inverse(self, X, score):
    n = X.shape[0]
    prediction, std = self.mu.predict(X, return_std = True)
    prediction = prediction.reshape((n,1))
    std = std.reshape((n,1))
    return prediction + score * std<jupyter_output><empty_output><jupyter_text>Conformity Score 3<jupyter_code>class QuantileDistanceConformityScore:

  def __init__(self, quantile = None, **kwargs):
    self.mu = RandomForestQuantileRegressor(**kwargs)
    self.quantile = 100 * quantile

  def train(self, X, Y):
    self.mu.fit(X, Y)

  def score(self, X, Y):
    n = X.shape[0]
    prediction = self.mu.predict(X, quantile = self.quantile).reshape((n,1))
    return Y - prediction

  def inverse(self, X, score):
    n = X.shape[0]
    prediction = self.mu.predict(X, quantile = self.quantile).reshape((n,1))
    return prediction + score<jupyter_output><empty_output><jupyter_text>Conformity Score 4<jupyter_code>class InverseCDFConformityScore:

  quantile_list = [0,1,2,3,4,5,10,20,30,40,50,60,70,80,90,95,96,97,98,99,100]

  def __init__(self, quantiles = None, **kwargs):
    self.mu = RandomForestQuantileRegressor(**kwargs)
    if quantiles:
      self.quantiles = quantiles
    else:
      self.quantiles = self.quantile_list

  def train(self, X, Y):
    self.mu.fit(X, Y)

  def score_individual(self, X_i, Y_i):
    Y_value = Y_i[0,0]
    quantile_predictions = {}

    low_index = 0
    high_index = len(self.quantiles) - 1
    mid_index = 0

    while low_index + 1 < high_index:
      mid_index = int((low_index + high_index) / 2)
      mid_value = self.quantiles[mid_index]
      quantile_predictions[mid_value] = self.mu.predict(X_i, quantile = mid_value)[0]

      if quantile_predictions[mid_value] < Y_value:
        low_index = mid_index
      elif quantile_predictions[mid_value] > Y_value:
        high_index = mid_index
      else:
        return mid_value

    low_value = self.quantiles[low_index]
    high_value = self.quantiles[high_index]
    if low_value not in quantile_predictions:
      quantile_predictions[low_value] = self.mu.predict(X_i, quantile = low_value)[0]
    if high_value not in quantile_predictions:
      quantile_predictions[high_value] = self.mu.predict(X_i, quantile = high_value)[0]

    if Y_value < quantile_predictions[low_value]:
      return self.quantiles[0] + Y_value - quantile_predictions[low_value]
    if Y_value > quantile_predictions[high_value]:
      return self.quantiles[-1] + Y_value - quantile_predictions[high_value]
    return np.interp(Y_value, [quantile_predictions[low_value], quantile_predictions[high_value]], [low_value, high_value])

  def score(self, X, Y):
    n = X.shape[0]
    scores = np.zeros((n, 1))

    for i in range(n):
      X_i = X[[i],:]
      Y_i = Y[[i],:]
      scores[i, 0] = self.score_individual(X_i, Y_i)
    return scores

  def inverse_individual(self, X_i, score_i):
    score_value = score_i[0,0]
    quantile_predictions = {}

    if score_value < self.quantiles[0]:
      quantile_predictions[self.quantiles[0]] = self.mu.predict(X_i, quantile = self.quantiles[0])[0]
      return quantile_predictions[self.quantiles[0]] + score_value - self.quantiles[0]
    if score_value > self.quantiles[-1]:
      quantile_predictions[self.quantiles[-1]] = self.mu.predict(X_i, quantile = self.quantiles[-1])[0]
      return quantile_predictions[self.quantiles[-1]] + score_value - self.quantiles[-1]

    low_index = np.searchsorted(self.quantiles, score_value, side = 'right') - 1
    high_index = np.searchsorted(self.quantiles, score_value, side = 'left')
    low_value = self.quantiles[low_index]
    high_value = self.quantiles[high_index]
    if low_value not in quantile_predictions:
      quantile_predictions[low_value] = self.mu.predict(X_i, quantile = low_value)[0]
    if high_value not in quantile_predictions:
      quantile_predictions[high_value] = self.mu.predict(X_i, quantile = high_value)[0]
    return np.interp(score_value, [low_value,high_value], [quantile_predictions[low_value], quantile_predictions[high_value]])
  
  def inverse(self, X, score):
    n = X.shape[0]
    inverses = np.zeros((n, 1))

    for i in range(n):
      X_i = X[[i],:]
      score_i = score[[i],:]
      inverses[i, 0] = self.inverse_individual(X_i, score_i)
    return inverses<jupyter_output><empty_output><jupyter_text>Conditional Quantile Algorithm This section contains the conditional quantile algorithm (Algorithm 2). It contains a train and test method.<jupyter_code>class ConditionalMedianAlgorithm:

  def __init__(self, conformity_score_low, conformity_score_high, 
               alpha = 0.1, quantile = 0.5, 
               data_split_proportion = 0.5, probability_split_proportion = 0.5):
    
    assert alpha > 0 and alpha < 1, f"Need alpha to be in (0,1) but got {alpha}"
    assert quantile > 0 and quantile < 1, f"Need quantile to be in (0,1) but got {quantile}"
    assert data_split_proportion > 0 and quantile < 1, f"Need data_split_proportion to be in (0,1) but got {data_split_proportion}"
    assert probability_split_proportion > 0 and probability_split_proportion < 1, f"Need probability_split_proportion to be in (0,1) but got {probability_split_proportion}"


    self.f_lo = conformity_score_low
    self.f_hi = conformity_score_high
    self.alpha = alpha
    self.q = quantile
    self.data_split_proportion = data_split_proportion
    self.probability_split_proportion = probability_split_proportion

    self.r = alpha * probability_split_proportion
    self.s = alpha * (1. - probability_split_proportion)

  def train(self, X, Y):
    n = X.shape[0]
    assert Y.shape[0] == n, f"Need X and Y to be 2-dimensional arrays with equal number of rows but got X.shape={X.shape} and Y.shape={Y.shape}"

    n_1 = int(n * self.data_split_proportion)
    n_2 = n - n_1
    X_shuffled, Y_shuffled = shuffle(X, Y)
    X_I_1 = X_shuffled[:n_1, :]
    X_I_2 = X_shuffled[n_1:, :]
    Y_I_1 = Y_shuffled[:n_1, :]
    Y_I_2 = Y_shuffled[n_1:, :]

    self.f_lo.train(X_I_1, Y_I_1)
    self.f_hi.train(X_I_1, Y_I_1)

    f_lo_scores = self.f_lo.score(X_I_2, Y_I_2).flatten()
    f_hi_scores = self.f_hi.score(X_I_2, Y_I_2).flatten()

    lower_quantile = self.r * self.q * (1. + 1. / n_2) - 1. / n_2
    upper_quantile = (1 - self.s * (1 - self.q)) * (1. + 1. / n_2)

    self.lower_score_bound = np.quantile(f_lo_scores, lower_quantile)
    self.upper_score_bound = np.quantile(f_hi_scores, upper_quantile)

  def test(self, X):
    n = X.shape[0]
    CI_lower = self.f_lo.inverse(X, np.full((n, 1), fill_value = self.lower_score_bound))
    CI_upper = self.f_hi.inverse(X, np.full((n, 1), fill_value = self.upper_score_bound))
    return np.concatenate((CI_lower, CI_upper), axis = 1)<jupyter_output><empty_output><jupyter_text>Quantile Forest Regression Algorithm This section contains a nonconformalized quantile regression algorithm. It is intended to be used to compare performance against.<jupyter_code>class RawMedianRandomForest:

  def __init__(self, alpha = 0.1, quantile = 0.5, **kwargs):
    assert alpha > 0 and alpha < 1, f"Need alpha to be in (0,1) but got {alpha}"
    
    self.alpha = alpha
    self.lower_quantile = 100 * (self.alpha / 2)
    self.upper_quantile = 100 * (1 - self.alpha / 2)
    self.mu = RandomForestQuantileRegressor(**kwargs)
    
  def train(self, X, Y):
    self.mu.fit(X, Y)

  def test(self, X):
    n = X.shape[0]
    CI_lower = self.mu.predict(X, quantile = self.lower_quantile).reshape((n,1))
    CI_upper = self.mu.predict(X, quantile = self.upper_quantile).reshape((n,1))
    return np.concatenate((CI_lower, CI_upper), axis = 1)<jupyter_output><empty_output><jupyter_text>Testing Algorithm Performance We now go over the performance of the Conditional Median Algorithm by testing its performance on different distributions and conformity scores.  Simulation Parameters<jupyter_code>alpha = 0.1

num_trials = 500

n_train = 5000
n_test = 5000<jupyter_output><empty_output><jupyter_text>Helper Functions<jupyter_code>def get_indices_accurate(values, CI):
  return np.logical_and(np.less_equal(CI[:,[0]], values), np.less_equal(values, CI[:,[1]]))

def get_proportion_accurate(values, CI):
  return np.average(get_indices_accurate(values, CI))

def get_widths(CI):
  return CI[:,[1]] - CI[:,[0]]

def get_avg_width(CI):
  return np.average(get_widths(CI))<jupyter_output><empty_output><jupyter_text>Distribution Parameter Specifications<jupyter_code>def get_dist(distribution):
  if distribution == 1:
    return NormalCorrelatedDistribution(d = 10, rho = 0.25, sigma = 1)
  if distribution == 2:
    return IncreasingVarianceDistribution(scale = 4 * np.pi, power = 0.25)
  if distribution == 3:
    return ConfusingDistribution(delta = 0.0001, M = 25)
  assert f"Need distribution to be 1, 2, or 3, but got {distribution}"<jupyter_output><empty_output><jupyter_text>Distribution Test Points<jupyter_code>dist_1_test_data = get_dist(1).sample(n_test)[0]
dist_2_test_data = get_dist(2).sample(n_test)[0]
dist_3_test_data = get_dist(3).sample(n_test)[0]

def get_test_datapoints(distribution, n):
  if distribution == 1:
    return dist_1_test_data[:n, :]
  if distribution == 2:
    return dist_2_test_data[:n, :]
  if distribution == 3:
    return dist_3_test_data[:n, :]
  assert f"Need distribution to be 1, 2, or 3, but got {distribution}"

for dist_number in [1,2,3]:
  test_datapoints = get_test_datapoints(dist_number, n_test)
  output = open(metadata_dir + f"Test Points/Distribution {dist_number}.pkl", 'wb')
  pickle.dump(test_datapoints, output)
  output.close()

def get_test_dist(dist_number):
  with open(metadata_dir + f"Test Points/Distribution {dist_number}.pkl", 'rb') as handle:
    return pickle.load(handle)<jupyter_output><empty_output><jupyter_text>Random Forest Hyperparameters<jupyter_code>kwargs = {"n_estimators": 20,
          "min_samples_split": 80,
          "min_samples_leaf": 40}

def get_scores(score, alpha):
  if score == 1:
    return LinearDistanceConformityScore(**kwargs), LinearDistanceConformityScore(**kwargs)
  if score == 2:
    return NormalizedDistanceConformityScore(**kwargs), NormalizedDistanceConformityScore(**kwargs)
  if score == 3:
    return QuantileDistanceConformityScore(quantile = alpha/2., criterion = "mae", **kwargs), QuantileDistanceConformityScore(quantile = 1. - alpha/2., criterion = "mae", **kwargs)
  if score == 4:
    return InverseCDFConformityScore(quantiles = None, criterion = "mae", **kwargs), InverseCDFConformityScore(quantiles = None, criterion = "mae", **kwargs)
  assert f"Need score to be 1, 2, 3, or 4, but got {score}"<jupyter_output><empty_output><jupyter_text>Singular Trial<jupyter_code>def run_singular_trial(distribution_number, score_number, make_fig = True, save_fig = False, trial_num = 0, print_stats = False, save_results = False):
  setup = f"Distribution {distribution_number}; Conformity Score {score_number}; Trial {trial_num}"
  print(setup)
  dist = get_dist(distribution_number)

  ## Raw quantile regression algorithm
  if score_number == 5:
    algorithm = RawMedianRandomForest(alpha = alpha, n_estimators = 20, min_samples_split = 160, min_samples_leaf = 80)
  else:
    f_lo, f_hi = get_scores(score_number, alpha = alpha)
    algorithm = ConditionalMedianAlgorithm(conformity_score_low = f_lo, conformity_score_high = f_hi, alpha = alpha)

  X_train, Y_train = dist.sample(n_train)
  algorithm.train(X_train, Y_train)

  X_test = get_test_dist(distribution_number)
  Y_median_true = dist.get_quantile(X_test, 0.5)
  Y_median_CI = algorithm.test(X_test)
  prop_accurate = get_proportion_accurate(Y_median_true, Y_median_CI)
  avg_width = get_avg_width(Y_median_CI)

  if print_stats:
    print(f"Proportion Accurate: {prop_accurate}")
    print(f"Average Width: {avg_width}")

  if make_fig and distribution_number in [2,3]:
    plt.figure(figsize = (5,4))
    plt.tight_layout()
    plt.plot(X_test, Y_median_true, color = median_color, linewidth = median_linewidth, label = 'Conditional Median')
    plt.fill_between(X_test.ravel(), Y_median_CI[:,[0]].ravel(), Y_median_CI[:,[1]].ravel(), label = 'Confidence Interval', color = 'r', alpha = 0.2)
    plt.scatter(X_train, Y_train, s = 0.5, color = datapoint_color, alpha = 0.25, label = 'Datapoints')
    plt.legend(markerscale = 5)

    plt.xlabel(r"$X$", fontsize = label_fontsize)
    plt.ylabel(r"$Y$", fontsize = label_fontsize)
    plt.title(f"Conformity Score {score_number}", fontsize = title_fontsize)
    if save_fig:
      filename = f"{setup}.{trial_num}.pdf"
      plt.savefig(images_dir + filename,bbox_inches='tight')
    plt.show()

  if save_results:
    output = open(metadata_dir + f"Performance/Distribution{distribution_number}.Score{score_number}.Trial{trial_num}.pkl", 'wb')
    output_data = {'Coverage': get_indices_accurate(Y_median_true, Y_median_CI), 'Width': get_widths(Y_median_CI)}
    pickle.dump(output_data, output)
    output.close()

  return get_indices_accurate(Y_median_true, Y_median_CI), get_widths(Y_median_CI)<jupyter_output><empty_output><jupyter_text>Get Images of Singular Trial on Distribution 2 for all conformity scores<jupyter_code>distribution_number = 2

for score_number in tqdm.tqdm([1,2,3,4]):
  indices_acc, widths = run_singular_trial(distribution_number, score_number, make_fig = True, save_fig = True, trial_num = 0)

  prop_accurate = np.average(indices_acc)
  avg_width = np.average(widths)
  performance = {"Proportion Accurate": prop_accurate,
                "Average Width": avg_width}
  output = open(images_dir + f"Performance Distribution {distribution_number} Score {score_number}.pkl", 'wb')
  pickle.dump(performance, output)
  output.close()<jupyter_output><empty_output><jupyter_text>Run 500 trials for each combination of distribution and conformity score to test performance<jupyter_code>for score_number in [1,2,3,4,5]:
  for dist_number in [1,2,3]:
    for trial in tqdm.tqdm(num_trials):
      run_singular_trial(dist_number, score_number, make_fig = False, save_fig = False, trial_num = trial + 1, save_results = True)<jupyter_output><empty_output><jupyter_text>Calculating Simulation Metrics We calculate the average rate of coverage, the standard deviation of rate of coverage, the minimal conditional coverage, the average confidence interval width, and the standard deviation of average confidence interval width.<jupyter_code>for score_number in [1,2,3,4,5]:
  for dist_number in [1,2,3]:

    statistics = {}

    coverages = []
    widths = []
    for trial_num in tqdm.tqdm(range(1, num_trials)):
      with open(metadata_dir + f"Performance/Distribution{dist_number}.Score{score_number}.Trial{trial_num}.pkl", 'rb') as handle:
        output =  pickle.load(handle)
      coverages.append(output['Coverage'])
      widths.append(output['Width'])
    coverages = np.array(coverages)
    widths = np.array(widths)

    statistics['AC'] = np.mean(np.mean(coverages, axis = 1))
    statistics['SDAC'] = np.std(np.mean(coverages, axis = 1), ddof = 1)
    statistics['MCC'] = np.min(np.mean(coverages, axis = 0))
    statistics['AW'] = np.mean(np.mean(widths, axis = 1))
    statistics['SDAW'] = np.std(np.mean(widths, axis = 1), ddof = 1)

    output = open(metadata_dir + f"Distribution{dist_number}.Score{score_number}.pkl", 'wb')
    pickle.dump(statistics, output)
    output.close()<jupyter_output><empty_output><jupyter_text>View Performance for any Distribution / Score<jupyter_code>dist_number = 1
score_number = 1

with open(metadata_dir + f"Distribution{dist_number}.Score{score_number}.pkl", 'rb') as handle:
  statistics =  pickle.load(handle)

statistics<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.123980626733452, 'fraction_numerical': 0.020242579790536903, 'mean_word_length': 4.514266149280986, 'pattern_counts': {'":': 5, '<': 97, '<?xml version=': 0, '>': 101, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '635687', 'license': 'MIT', 'n_tokens_mistral': 8994, 'n_tokens_neox': 8354, 'n_words': 2022, 'path': 'Conditional Median Inference Simulations.ipynb', 'repo_name': 'dhruvm2/conditional-median-inference', 'size': '24177'}