<jupyter_start><jupyter_text>The Wilson score interval and weighted histogramsThe Wilson Score interval (WSI) was won the LHCb internal challenge for the most suitable interval for efficiencies. The question what to do when events are weighted was left open at the time. It is addressed here, the WSI can be used also with weighted histograms if ordinary counts in unweighted histograms are replaced with effective counts in weighted histograms.<jupyter_code>import numpy as np
import boost_histogram as bh
from matplotlib import pyplot as plt<jupyter_output><empty_output><jupyter_text>https://en.wikipedia.org/wiki/Binomial_proportion_confidence_intervalWilson score interval:$$p_\text{min} = \frac{n_S + \frac 12 z^2}{n_S + n_F + z^2} - \frac{z}{n_S + n_F + z^2} \sqrt{\frac{n_S n_F}{n_S + n_F} + \frac{z^2}{4}}$$$$p_\text{max} = \frac{n_S + \frac 12 z^2}{n_S + n_F + z^2} + \frac{z}{n_S + n_F + z^2} \sqrt{\frac{n_S n_F}{n_S + n_F} + \frac{z^2}{4}}$$with $n_S$ counts of selected events, $n_F$ counts of rejected events, $z$ number of sigmas ($z = 1$ for standard intervals).In case of weighted events, the $n_S,n_F$ are the effective counts, computed as$$n_\text{effective} = \frac{\big(\sum_i w_i\big)^2}{\sum_i w^2_i}.$$In boost-histogram, this is computed as `n_eff = h.values() ** 2 / h.variances()` for a histogram `h` with `WeightedSum` storage.<jupyter_code>def wilson_score_interval(n_s, n_f, z=1):
    n = n_s + n_f
    z2 = z * z
    center = (n_s + 0.5 * z2) / (n + z2)
    delta = z / (n + z2) * np.sqrt(n_s * n_f / n + 0.25 * z2)
    return center - delta, center + delta
p_truth = 0.2

rng = np.random.default_rng(1)

w_s = bh.accumulators.WeightedSum()
w_f = bh.accumulators.WeightedSum()

w = rng.exponential(10, size=100)
m = rng.uniform(size=len(w)) < p_truth

w_s.fill(w[m])
w_f.fill(w[~m])
n_s = w_s.value ** 2 / w_s.variance if w_s.variance > 0 else 0
n_f = w_f.value ** 2 / w_f.variance if w_f.variance > 0 else 0

p = n_s / (n_s + n_f)
p_min, p_max = wilson_score_interval(n_s, n_f)

plt.errorbar([1], [p], [[p-p_min], [p_max-p]], fmt="ok")

print(f"efficiency = {p:.2f} - {p-p_min:.2f} + {p_max-p:.2f}")<jupyter_output>efficiency = 0.16 - 0.05 + 0.06<jupyter_text>Note that the error bar is asymmetric in general. This is especially noticable when p is close to 0 or 1.Let's check with a parametric bootstrap simulation whether this uncertainty is correct. It is parametric because I put in we know that the weights are exponentially distributed.<jupyter_code>rng = np.random.default_rng(1)

w_mean = np.mean(w)

p_boot = []
for iboot in range(1000):
    wi = rng.exponential(w_mean, size=len(w))
    mi = rng.uniform(size=len(w)) < p
    t_s = bh.accumulators.WeightedSum()
    t_f = bh.accumulators.WeightedSum()
    t_s.fill(wi[mi])
    t_f.fill(wi[~mi])
    p_boot.append(t_s.value / (t_s.value + t_f.value))
    
center = np.mean(p_boot)
delta = np.std(p_boot)
p_min_boot = center - delta
p_max_boot = center + delta

plt.errorbar([1], [p], [[p-p_min], [p_max-p]], fmt="ok")
plt.errorbar([2], [p], [[p-p_min_boot], [p_max_boot-p]], fmt="sr")
plt.axhline(p_truth, ls="--", color="0.5")

print(f"efficiency = {p:.2f} - {p-p_min:.2f} + {p_max-p:.2f}")
print(f"efficiency(boot) = {p:.2f} - {p-p_min_boot:.2f} + {p_max_boot-p:.2f}")<jupyter_output>efficiency = 0.16 - 0.05 + 0.06
efficiency(boot) = 0.16 - 0.05 + 0.05
{'fraction_non_alphanumeric': 0.1592946802151823, 'fraction_numerical': 0.027196652719665274, 'mean_word_length': 4.254317111459969, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 13, 'https://': 1, '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': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '508609', 'license': 'MIT', 'n_tokens_mistral': 1462, 'n_tokens_neox': 1354, 'n_words': 347, 'path': 'Wilson Score Interval with Weighted Histograms.ipynb', 'repo_name': 'HDembinski/essays', 'size': '3346'}