Bayesian rate estimation

Deep dive

How it works

In: s of n successes, prior · Out: posterior, interval, P(B beats A) · State: two counters per arm · Cost: O(1) update

1. Prior on the rate p: Beta(a0, b0), density proportional to p^(a0-1) (1-p)^(b0-1). Beta(1,1) is uniform, Beta(0.5,0.5) is Jeffreys; a0 + b0 is the prior strength in pseudo-trials.
2. Data: s successes in n trials (binomial likelihood). The posterior is conjugate: Beta(a0 + s, b0 + n - s). No integral is needed.
3. Summaries: mean a / (a + b); equal-tailed credible interval = the t and 1 - t quantiles of the posterior (t = 0.025 for 95 %). The quantiles come from the incomplete beta function: CDF(x) = I_x(a, b).
4. Comparing arms: P(p_B > p_A) is an integral over two Betas. For integer parameters it is a finite sum (Evan Miller): sum_{i=0}^{aB-1} B(aA+i, bA+bB) / ((bB+i) B(1+i, bB) B(aA, bA)), with B the Beta function.
5. Count rates: for k events over exposure t, a Gamma(alpha, beta) prior (shape, rate) gives Gamma(alpha + k, beta + t); mean (alpha + k) / (beta + t).
6. Empirical Bayes: when you have many comparable groups, fit (a0, b0) to their spread (method of moments, below) so the prior reflects the population instead of a guess. Small groups then shrink strongly, large ones barely.

raw     0/3      1/5     30/90   110/600
        |         |       |        |        prior Beta(2.6, 20.0), mean 0.116
shrunk  0.103   0.131    0.290    0.181     (small n pulled to prior mean, large n kept)

Code

/// ln Gamma(x), x > 0 (Lanczos, g = 7, n = 9), about 1e-15 relative error.
fn ln_gamma(x: f64) -> f64 {
    const C: [f64; 9] = [
        0.999_999_999_999_809_9, 676.520_368_121_885_1, -1_259.139_216_722_402_8,
        771.323_428_777_653_1, -176.615_029_162_140_6, 12.507_343_278_686_905,
        -0.138_571_095_265_720_12, 9.984_369_578_019_572e-6, 1.505_632_735_149_311_6e-7,
    ];
    if x < 0.5 {
        return (std::f64::consts::PI / (std::f64::consts::PI * x).sin()).ln() - ln_gamma(1.0 - x);
    }
    let x = x - 1.0;
    let t = x + 7.5;
    let s = C[0] + (1..9).map(|i| C[i] / (x + i as f64)).sum::<f64>();
    0.5 * (2.0 * std::f64::consts::PI).ln() + (x + 0.5) * t.ln() - t + s.ln()
}
fn ln_beta(a: f64, b: f64) -> f64 { ln_gamma(a) + ln_gamma(b) - ln_gamma(a + b) }

/// Regularized incomplete beta I_x(a, b) = Beta(a,b) CDF (modified Lentz continued fraction).
pub fn beta_cdf(x: f64, a: f64, b: f64) -> f64 {
    if x <= 0.0 { return 0.0; }
    if x >= 1.0 { return 1.0; }
    let front = (a * x.ln() + b * (1.0 - x).ln() - ln_beta(a, b)).exp();
    if x < (a + 1.0) / (a + b + 2.0) {
        front * cf(x, a, b) / a
    } else {
        1.0 - front * cf(1.0 - x, b, a) / b // symmetry: converges faster
    }
}
fn cf(x: f64, a: f64, b: f64) -> f64 {
    const TINY: f64 = 1e-300;
    let (mut c, mut d) = (1.0, 1.0 - (a + b) * x / (a + 1.0));
    if d.abs() < TINY { d = TINY; }
    d = 1.0 / d;
    let mut h = d;
    for m in 1..=300 {
        let m = m as f64;
        for step in 0..2 {
            let num = if step == 0 {
                m * (b - m) * x / ((a + 2.0 * m - 1.0) * (a + 2.0 * m))
            } else {
                -(a + m) * (a + b + m) * x / ((a + 2.0 * m) * (a + 2.0 * m + 1.0))
            };
            d = 1.0 + num * d; if d.abs() < TINY { d = TINY; }
            c = 1.0 + num / c; if c.abs() < TINY { c = TINY; }
            d = 1.0 / d;
            let delta = c * d;
            h *= delta;
            if step == 1 && (delta - 1.0).abs() < 1e-14 { return h; }
        }
    }
    h
}

/// Quantile of Beta(a, b) by bisection on the CDF (monotone, so it cannot fail).
pub fn beta_quantile(p: f64, a: f64, b: f64) -> f64 {
    let (mut lo, mut hi) = (0.0, 1.0);
    for _ in 0..60 {
        let mid = 0.5 * (lo + hi);
        if beta_cdf(mid, a, b) < p { lo = mid } else { hi = mid }
    }
    0.5 * (lo + hi)
}

/// Posterior of a rate after `s` successes in `n` trials with a Beta(a0, b0) prior.
#[derive(Clone, Copy, Debug)]
pub struct Post { pub a: f64, pub b: f64 }
impl Post {
    pub fn new(a0: f64, b0: f64, s: u64, n: u64) -> Post {
        Post { a: a0 + s as f64, b: b0 + (n - s) as f64 }
    }
    pub fn mean(&self) -> f64 { self.a / (self.a + self.b) }
    /// Equal-tailed credible interval.
    pub fn interval(&self, level: f64) -> (f64, f64) {
        let t = 0.5 * (1.0 - level);
        (beta_quantile(t, self.a, self.b), beta_quantile(1.0 - t, self.a, self.b))
    }
}

/// P(p_b > p_a), exact for integer a and b parameters (sum of at most b_b terms).
pub fn prob_b_beats_a(a: Post, b: Post) -> f64 {
    let (aa, ba, ab, bb) = (a.a, a.b, b.a, b.b);
    let (mut sum, base) = (0.0, ln_beta(aa, ba));
    for i in 0..(ab as u64) {
        let i = i as f64;
        sum += (ln_beta(aa + i, ba + bb) - (bb + i).ln() - ln_beta(1.0 + i, bb) - base).exp();
    }
    sum
}

/// Empirical-Bayes Beta prior by method of moments over many (successes, trials) groups.
/// Moments are trial-weighted; returns (a0, b0). Needs real between-group spread.
pub fn fit_beta_prior(groups: &[(u64, u64)]) -> Option<(f64, f64)> {
    let nt: f64 = groups.iter().map(|g| g.1 as f64).sum();
    let m: f64 = groups.iter().map(|g| g.0 as f64).sum::<f64>() / nt;
    let v: f64 = groups.iter().map(|g| g.1 as f64 * (g.0 as f64 / g.1 as f64 - m).powi(2)).sum::<f64>() / nt;
    let nbar = nt / groups.len() as f64;
    let v = v - m * (1.0 - m) / nbar; // remove sampling noise; approximate for unequal n
    if v <= 0.0 || v >= m * (1.0 - m) { return None; } // no usable spread
    let k = m * (1.0 - m) / v - 1.0; // prior strength a0 + b0
    Some((m * k, (1.0 - m) * k))
}

Checked by running: Beta CDF against closed forms (I_0.2(5,3) = 0.004672), the 0/100 interval against 1 - 0.05^(1/101) (to 1e-9), and P(B>A) against a 20000-draw Monte Carlo. A 120/1000 vs 150/1000 test with uniform priors gives means 0.121 and 0.151, 95 % intervals [0.101, 0.142] and [0.129, 0.173], P(B>A) = 0.975 (Monte Carlo 0.976). The empirical-Bayes fit on ten groups (n from 3 to 600, raw rates 0 to 0.33) gave Beta(2.64, 20.03), and shrank 0/3 to 0.103 and 30/90 to 0.290.

Complexity and constants

• Update: two integer additions, 16 bytes per arm. Memory is independent of the number of observations.
• CDF: continued fraction, order of tens of iterations (about 1e-14 tolerance) at a few flops each plus 3 ln_gamma evaluations; the quantile here is bisection with 60 CDF calls.
• Measured here (Apple M-series, release build, scalar): a 95 % interval (two quantiles) about 25-30 us; P(B>A) at 1000-trial arms about 9 us (a ~150-term sum). Order of magnitude only.
• Cost drivers at scale are not the maths but the counting: bytes read to produce (s, n), and the number of arms you evaluate per request.

Efficiency

• Sufficient statistics: keep (s, n) per arm and discard raw events. Merge shards, days or replicas by adding; no per-event storage, no re-scan to change the prior (re-derive from s, n).
• Closed form beats sampling: P(B>A) as a finite sum costs order of a thousand ln_gamma evaluations (the terms also follow a cheap recurrence, so a production version needs far fewer); a Monte Carlo estimate needs O(1/eps^2) draws for accuracy eps (0.0035 worst case with 20000 draws). MCMC is only needed for non-conjugate or hierarchical models.
• Cheaper quantiles: seed Newton or a few Halley steps with the normal approximation instead of 60 bisections, or call a tuned library (Boost.Math, SciPy betaincinv) once per refresh, not per request. Cache the interval per arm and recompute only when the counters change.
• Batching: accumulate counts in memory (sharded atomics) and flush periodically; one write per interval, not per event. Skip the interval maths entirely for arms outside the decision range.
• Early stopping of work: once P(B>A) is far from the decision threshold, stop collecting for that arm; fewer samples is fewer bytes logged, shipped and stored.
• Cheaper alternative: for a single rate with large n, a normal or Wilson interval needs no special functions; the Bayesian posterior mostly matters when n is small or you want a probability statement.

Implementation notes

• Work in log space for ln_beta; the incomplete-beta series swaps (x, a, b) to (1-x, b, a) when x >= (a+1)/(a+b+2) for fast convergence (done above). In production use a vetted implementation.
• The P(B>A) sum needs integer aB and does aB terms. When aB is large, use the mirror form: p_B > p_A iff 1-p_B < 1-p_A, i.e. P(Beta(bA,aA) > Beta(bB,aB)), which loops over bA terms instead; pick the cheaper form, or use a normal approximation or numerical integration for thousands of terms.
• Prior defaults: Beta(1,1) when you know nothing; a weak informative prior (strength of order 10-50 trials, centred at the historical rate) for rate-ranking; never a=b=0 (improper, posterior breaks at s=0 or s=n).
• Empirical-Bayes moment fit: needs real between-group variance beyond binomial noise; if the corrected variance is <= 0 the groups look alike and the fit is undefined (the code returns None). Unequal n makes the moment estimate approximate; a marginal-likelihood (beta-binomial) fit is better.
• Count the right unit: successes and trials must be the same unit (users, not clicks per user).

Pitfalls

• Reporting the posterior mean as if it were the rate: show the interval too, and the sample size.
• A strong prior from the wrong population (last year's conversion on a new product) biases small samples for a long time.
• Peeking: stopping an experiment the first time P(B>A) > 0.95 does not preserve a 5 % false-positive rate; if error rates matter use a sequential test or a decision threshold with a minimum sample and loss budget.
• P(B>A) ignores effect size: 0.99 for a 0.01 % lift is not useful; also report the interval of the difference or the expected loss.
• Empirical-Bayes treats the fitted prior as exact; with few groups its own uncertainty is large (use a hierarchical model).
• Dependent observations (same user many times) make n too large and intervals too narrow.
• Mixing Beta-Binomial with rates that vary by day of week or traffic source: stratify or model it.

Variants and alternatives

In production

Conversion reporting in A/B tests (Miller); click rates and bandits in online advertising (Chapelle and Li); ranking items with few ratings; batting averages in sports with empirical Bayes (Robinson).

Implementations

Sources

• Andrew Gelman, John B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari, Donald B. Rubin, Bayesian Data Analysis, 3rd ed., CRC Press 2013 (Beta-Binomial, Gamma-Poisson, hierarchical models).
• Evan Miller, Formulas for Bayesian A/B Testing, 2015. article
• Lawrence D. Brown, T. Tony Cai, Anirban DasGupta, Interval Estimation for a Binomial Proportion, Statistical Science 2001. paper
• Bradley Efron, Carl Morris, Data Analysis Using Stein's Estimator and Its Generalizations, JASA 1975. paper
• Herbert Robbins, An Empirical Bayes Approach to Statistics, Proc. Third Berkeley Symposium on Mathematical Statistics and Probability, vol. 1, 1956.
• NIST, Digital Library of Mathematical Functions, 8.17 (incomplete beta function, continued fraction). reference