|
| 1 | +# Copyright (c) Meta Platforms, Inc. and affiliates. |
| 2 | + |
| 3 | +import numpy as np |
| 4 | +import pandas as pd |
| 5 | +import plotly.graph_objects as go |
| 6 | +from plotly import express as px |
| 7 | +from scipy.stats import beta, gamma |
| 8 | + |
| 9 | + |
| 10 | +##### Initialize constants |
| 11 | + |
| 12 | +# fix the random seed for reproducibility |
| 13 | +np.random.seed(0) |
| 14 | + |
| 15 | +# How many rows are in the population table |
| 16 | +population_size = 1000 |
| 17 | + |
| 18 | +# What % of population rows will get sampled |
| 19 | +sample_rate = 0.01 |
| 20 | + |
| 21 | +# The number of population rows that get sampled |
| 22 | +sample_size = int(population_size * sample_rate) |
| 23 | +sample_size = 100 |
| 24 | + |
| 25 | +# The beta parameters that are used for sampling ground truth labels |
| 26 | +# for each row. |
| 27 | +beta_prior_negative = 50 |
| 28 | +beta_prior_positive = 10 |
| 29 | + |
| 30 | + |
| 31 | +###### Create the population |
| 32 | + |
| 33 | +# Sample the p(positive) score for each population row |
| 34 | +pop_positive_probs = np.random.beta( |
| 35 | + beta_prior_positive, beta_prior_negative, size=(population_size) |
| 36 | +) |
| 37 | + |
| 38 | +# Sample a ground truth label for each population row |
| 39 | +pop_labels = np.random.rand((pop_positive_probs.size)) < pop_positive_probs |
| 40 | +pop_mean = pop_labels.mean() |
| 41 | +print("Population mean label:", pop_mean) |
| 42 | + |
| 43 | +num_trials = 100 |
| 44 | + |
| 45 | +# Add some noise to the probabilities to simulate importance sampling with an imperfect classifier |
| 46 | +noise_mean = 0 |
| 47 | +noise_std_dev = 0.1 |
| 48 | +classifier_preds = (pop_positive_probs) + np.random.normal( |
| 49 | + noise_mean, noise_std_dev, pop_positive_probs.shape |
| 50 | +) |
| 51 | +classifier_preds = np.clip(classifier_preds, 0, 1) + 1e-9 |
| 52 | + |
| 53 | + |
| 54 | +# take the sqrt because that's what we do today |
| 55 | +# TODO verify necessity |
| 56 | +sample_weights = np.sqrt(classifier_preds) |
| 57 | + |
| 58 | +pop_sampling_probs = sample_weights / sample_weights.sum() |
| 59 | + |
| 60 | + |
| 61 | +sampled_indices = np.random.choice( |
| 62 | + np.arange(population_size), |
| 63 | + size=(num_trials, sample_size), |
| 64 | + replace=True, |
| 65 | + p=pop_sampling_probs, |
| 66 | +) |
| 67 | + |
| 68 | +# the probability each sample is picked |
| 69 | +sample_probs = pop_sampling_probs[sampled_indices] |
| 70 | +sample_labels = pop_labels[sampled_indices] |
| 71 | +sample_rep_weights = (sample_size / sample_probs) / (sample_size / sample_probs).sum( |
| 72 | + -1, keepdims=True |
| 73 | +) |
| 74 | + |
| 75 | + |
| 76 | +# This function implements the conjugate beta estimator algorithm. |
| 77 | +# Its forumala is |
| 78 | +# alpha = mu*n + alpha_prior |
| 79 | +# beta = (1-mu)*n + beta_prior |
| 80 | +# ci = [ppf(0.05, alpha, beta), ppf(0.95, alpha, beta] |
| 81 | +# where mu is the mean of the (weighted) sample labels and n is the sample size in bits |
| 82 | +# (see below). |
| 83 | +def cbe(sample_mean, sample_size): |
| 84 | + alpha_param = sample_mean * sample_size + 1e-9 # optional: add an alpha prior |
| 85 | + beta_param = (1 - sample_mean) * sample_size + 1e-9 # optional: add a beta prior |
| 86 | + lower_bound = beta.ppf(0.025, alpha_param, beta_param) |
| 87 | + median = beta.ppf(0.5, alpha_param, beta_param) |
| 88 | + upper_bound = beta.ppf(0.975, alpha_param, beta_param) |
| 89 | + return lower_bound.mean(), median.mean(), upper_bound.mean() |
| 90 | + |
| 91 | + |
| 92 | +### Perfect labels section |
| 93 | +sample_mean = (sample_labels * sample_rep_weights).sum(-1) |
| 94 | +res = cbe(sample_mean, sample_size) |
| 95 | +print("CI with perfect labels:", res) |
| 96 | + |
| 97 | + |
| 98 | +### Noisy labels section |
| 99 | + |
| 100 | +# Generate random noise |
| 101 | +rns = np.random.uniform(0, 1, size=sample_labels.shape) |
| 102 | +noisy_labels = sample_labels.copy().astype(int) |
| 103 | + |
| 104 | +sensitivity = 0.85 |
| 105 | +specificity = 0.93 |
| 106 | + |
| 107 | +# Add random noise to negative samples according to sensitivity |
| 108 | +noisy_labels[sample_labels == 0] += rns[sample_labels == 0] > sensitivity |
| 109 | + |
| 110 | +# Add random noise to positive samples according to specificity |
| 111 | +noisy_labels[sample_labels == 1] = ( |
| 112 | + noisy_labels[sample_labels == 1] |
| 113 | + + (rns[sample_labels == 1] > specificity).astype(int) |
| 114 | +) % 2 |
| 115 | + |
| 116 | + |
| 117 | +# Rogan Gladen estimator for the sample mean. See: https://en.wikipedia.org/wiki/Beth_Gladen |
| 118 | +# To derive it, we start with the following equation, where o is the observed mean, p is the true mean, |
| 119 | +# and tpr and fpr are the true positive and false positive rates. |
| 120 | +# o = tpr*p + fpr*(1-p) |
| 121 | +# o = tpr*p + fpr - fpr*p |
| 122 | +# p*(tpr - fpr) = o - fpr |
| 123 | +# p = (o - fpr)/(tpr - fpr) |
| 124 | +def rg(mean): |
| 125 | + return (mean + specificity - 1) / (sensitivity + specificity - 1) |
| 126 | + |
| 127 | + |
| 128 | +noisy_mean = (noisy_labels * sample_rep_weights).sum(-1) |
| 129 | +print("Noisy mean:", noisy_mean.mean()) |
| 130 | +res = cbe(noisy_mean, sample_size) |
| 131 | +print("CI with noisy labels:", res) |
| 132 | + |
| 133 | +rg_mean = rg((noisy_mean)) |
| 134 | +print("RG sample mean:", rg_mean.mean()) |
| 135 | + |
| 136 | +res = cbe(rg_mean, sample_size) |
| 137 | +print("CI with RG mean and standard sample size", res) |
| 138 | + |
| 139 | + |
| 140 | +def entropy(val): |
| 141 | + return -val * np.log2(val) - (1 - val) * np.log2(1 - val) |
| 142 | + |
| 143 | + |
| 144 | +# When the mean of sensitivity and specificity approachs 0.5, the denominator |
| 145 | +# in the Rogan Gladen formula, (mean + specificity - 1)/(sensitivity + specificity - 1) |
| 146 | +# approaches 0. This causes the RG mean estimate to be highly unstable to small |
| 147 | +# fluctuations in sensitivity and/or specificity. To counter this instability, we want |
| 148 | +# to widen the CI. The entropy in this case is 1 (because entropy(0.5) = 1). |
| 149 | +# In the formula below, an entropy value of 1 leads to 0 bits per sample, which |
| 150 | +# amounts to having no samples when those samples are basically random noise. |
| 151 | +# However, as the mean of sensitivity and specificity approaches 1, the entropy |
| 152 | +# approaches 0, leading to about 1 bit per sample, which is the behavior |
| 153 | +# of the original CBE algorithm. |
| 154 | +num_bits = (1 - entropy((sensitivity + specificity) / 2)) * sample_size |
| 155 | +res = cbe(rg_mean, num_bits) |
| 156 | +print("CI with RG mean and entropy weighted sample size:", res) |
0 commit comments