Bayesian average rating

Deep dive

How it works

In: per item n ratings with sum S; global mean m, prior strength C · Out: one sortable score per item · State: n, S, and S2 (sum of squares, to fit C) · Cost: O(1) per rating and per score

1. Model: rating x ~ N(q, σ²) around item quality q; qualities q ~ N(m, τ²) across items.
2. The posterior of q is normal with mean (m/τ² + S/σ²) / (1/τ² + n/σ²). Multiply through by σ² and set C = σ²/τ²: mean (C·m + S) / (C + n), variance σ² / (C + n).
3. C is the prior strength in ratings; own data gets weight n/(C+n). With C = 8.8: one rating 10 %, 100 ratings 92 %.
4. Stars: a Dirichlet prior with pseudo-counts α_k per star gives the same formula with C = Σα_k, m = Σ k·α_k / C. Up/down: a Beta(a0, b0) prior gives C = a0 + b0, m = a0/C (Bayesian rate estimation).
5. Fit by method of moments: m = mean of item means; σ² = pooled within-item variance; τ² = mean of (x̄_i − m)² − σ²/n_i (spread minus noise); C = σ²/τ². If τ² ≤ 0 the items are indistinguishable and all scores collapse to m.
6. Rank by the score, or by score − z·σ/√(C+n) for a cautious list. Refit m, C on a schedule, not per rating.

Why it works: shrinking noisy estimates toward a common mean lowers total squared error for three or more items (James and Stein 1961, equal noise; Efron and Morris 1975 on batting averages).

Code

// Damped (Bayesian) average: (C*m + sum)/(C + n), with C and m fitted from the data.
struct Rng(u64);
impl Rng {
    fn next(&mut self) -> u64 { // splitmix64
        self.0 = self.0.wrapping_add(0x9E37_79B9_7F4A_7C15);
        let mut z = self.0;
        z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
        z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
        z ^ (z >> 31)
    }
    fn unit(&mut self) -> f64 { ((self.next() >> 11) as f64 + 0.5) / (1u64 << 53) as f64 }
    fn normal(&mut self) -> f64 { // Box-Muller
        (-2.0 * self.unit().ln()).sqrt() * (2.0 * std::f64::consts::PI * self.unit()).cos()
    }
}

/// Per-item state: just two numbers plus the sum of squares (for the variance fit).
#[derive(Clone, Copy, Default)]
struct Item { n: f64, sum: f64, sumsq: f64 }
impl Item {
    fn add(&mut self, x: f64) { self.n += 1.0; self.sum += x; self.sumsq += x * x; }
    /// Exponential decay with half-life `h`: call before add() when `dt` time has passed.
    fn decay(&mut self, dt: f64, h: f64) {
        let f = (-dt / h * std::f64::consts::LN_2).exp();
        self.n *= f; self.sum *= f; self.sumsq *= f;
    }
    fn raw(&self) -> f64 { if self.n > 0.0 { self.sum / self.n } else { f64::NAN } }
    /// Posterior mean with prior mean m and strength c (pseudo-ratings).
    fn damped(&self, m: f64, c: f64) -> f64 { (c * m + self.sum) / (c + self.n) }
    /// Conservative score: posterior mean minus z posterior standard deviations.
    fn lower(&self, m: f64, c: f64, sigma2: f64, z: f64) -> f64 {
        self.damped(m, c) - z * (sigma2 / (c + self.n)).sqrt()
    }
}

/// Method of moments: returns (m, C, sigma2) with C = sigma2 / tau2.
fn fit_prior(items: &[Item]) -> Option<(f64, f64, f64)> {
    let (mut ss, mut df) = (0.0, 0.0);
    for it in items.iter().filter(|i| i.n >= 2.0) {
        ss += it.sumsq - it.sum * it.sum / it.n; // within-item sum of squares
        df += it.n - 1.0;
    }
    let sigma2 = ss / df; // pooled noise variance of one rating
    let m = items.iter().map(|i| i.raw()).sum::<f64>() / items.len() as f64;
    // E[(mean_i - m)^2] = tau2 + sigma2 / n_i, so subtract the noise part.
    let tau2 = items.iter().map(|i| (i.raw() - m).powi(2) - sigma2 / i.n).sum::<f64>() / items.len() as f64;
    if tau2 <= 0.0 { return None; } // items look identical: shrink everything to m
    Some((m, sigma2 / tau2, sigma2))
}

fn top_k(score: &[f64], k: usize) -> Vec<usize> {
    let mut idx: Vec<usize> = (0..score.len()).collect();
    idx.sort_by(|&a, &b| score[b].total_cmp(&score[a]));
    idx.truncate(k);
    idx
}

fn main() {
    let (n_items, k, (mu, tau, sigma)) = (5000, 100, (3.6, 0.35, 1.0));
    let mut rng = Rng(42);
    let mut truth = vec![];
    let mut items = vec![];
    for _ in 0..n_items {
        let q = mu + tau * rng.normal(); // true quality drawn from the population
        let n = 1 + (2.0 + 1.3 * rng.normal()).exp().min(5000.0) as usize; // heavy-tailed counts
        let mut it = Item::default();
        for _ in 0..n { it.add(q + sigma * rng.normal()); }
        truth.push(q);
        items.push(it);
    }
    let (m, c, s2) = fit_prior(&items).unwrap();
    println!("fitted m = {m:.3}, C = {c:.2} (true C = {:.2}), sigma = {:.3}", sigma * sigma / (tau * tau), s2.sqrt());
    let raw: Vec<f64> = items.iter().map(|i| i.raw()).collect();
    let damp: Vec<f64> = items.iter().map(|i| i.damped(m, c)).collect();
    let low: Vec<f64> = items.iter().map(|i| i.lower(m, c, s2, 1.0)).collect();
    let mse = |s: &[f64]| s.iter().zip(&truth).map(|(a, b)| (a - b).powi(2)).sum::<f64>() / s.len() as f64;
    let best = top_k(&truth, k);
    let hit = |s: &[f64]| top_k(s, k).iter().filter(|i| best.contains(i)).count();
    println!("MSE   raw {:.4}  damped {:.4}", mse(&raw), mse(&damp));
    println!("top-{k} hits  raw {}  damped {}  damped-minus-1sd {}", hit(&raw), hit(&damp), hit(&low));
    assert!(mse(&damp) < mse(&raw) && hit(&damp) > hit(&raw));
    // decay: after one half-life the old evidence counts half
    let mut it = Item::default();
    it.add(5.0); it.decay(30.0, 30.0);
    assert!((it.n - 0.5).abs() < 1e-12);
}

Output of this code, release build, seed 42 (measured on this Mac (Apple M4 Pro)): m = 3.583, C = 8.81 (true 1.0²/0.35² = 8.16); MSE against true quality raw 0.2113, damped 0.0616; top-100 hits raw 15, damped 37, damped minus 1 sd 34. Seeds 1 to 10: damped hits 27 to 41, raw 10 to 20.

Complexity and constants

• Update: three float adds, 24 bytes per item. Score: one multiply-add, one divide; the cost drivers are the sort and the counter reads.
• Measured on this Mac (Apple M4 Pro), release build, 1 million synthetic items, separate harness around the same functions: fit_prior 2.0 ms, scoring 0.6 ms per pass, scoring plus a full sort for the top 100 about 23 ms.

Efficiency

• Cost: O(1) per rating and per score; about 0.6 ns per item to score, 2 ns per item to refit (measured on this Mac (Apple M4 Pro)). Bound: memory bandwidth over the item table, or the counter read in a service.
• Levers:
  • Keep (n, S, S2), not raw ratings — constant memory, merge by addition — when raw ratings are archived elsewhere.
  • Refit C offline on a schedule (2 ms per million items, measured) — always.
  • Keep scores in a sorted index, update only changed items, batch writes — when n is large and single ratings barely move a score.
  • select_nth_unstable for top-k: linear time instead of a full O(N log N) sort (Rust std docs) — when showing one page.
  • Decay lazily with a stored timestamp at read time — when ratings expire.
• Naive → standard → state of the art:
  • raw mean, hide items under a minimum count · O(1), cliff at the threshold
  • damped mean, fitted m, C · O(1), 24 bytes per item
  • hierarchical model with category, rater and covariate effects · heavy offline fit, O(1) lookup

Implementation notes

• Thin data: too few items to fit C? Set it from a guess of σ²/τ² (rating noise over quality spread; 8.8 in the run above); m from the catalogue or category, never 0.
• Use f64 for S2: S2 − S²/n cancels badly in f32.
• Time decay: multiply n, S, S2 by 2^(−dt/H) (half-life H) before adding a rating (Miller 2012 for up/down counts). Unrated items drift back to m.

Pitfalls

• Skewed stars. Stars are bounded and J-shaped; the normal model is a convenience; a Dirichlet suits polarising items.
• Different raters. Fans rate one item, everyone another: scores are not comparable and one global m hides it. Shrink to a category mean, correct rater leniency (item response theory style rater and item effects), risk-adjust league tables for case mix.
• Mean versus rank. If wrongly promoting costs more than wrongly burying, rank by a lower quantile: Miller's loss factor L gives the 1/(1+L) quantile; L = 5 is the 16.7 % quantile, about mean − 0.97 sd. Here the plain mean found more of the true top 100 (37 vs 34): the bound is a risk preference, not accuracy.
• Manipulation. Fake ratings move thin items most: with m = 3.6, C = 8.8, an item with 10 real ratings averaging 3.0 scores 3.28; 14 fake 5-star ratings lift it to 4.0, while 1,000 real ratings need 1,004. Damping raises the price but does not check who rates. Defences sit around the formula: verified raters, one rating per account, rate limits, burst detection (bot detection, reputation scores).
• Feedback loop. Low score, low exposure, few new ratings. Add exploration (bandits).

Variants and alternatives

In production

• Film and book lists: IMDb publishes weighted averages and does not disclose the weighting (help page).
• Sports: empirical Bayes batting averages (Efron and Morris; Robinson's book).
• League tables: the same shrinkage, called random-effects or hierarchical models.

Implementations

Sources

• Evan Miller, Bayesian Average Ratings, 2012. article
• Evan Miller, How Not To Sort By Average Rating, 2009. article
• Bradley Efron, Carl Morris, Data Analysis Using Stein's Estimator and Its Generalizations, JASA 70(350), 1975, pp. 311-319. paper
• W. James, Charles Stein, Estimation with Quadratic Loss, Proc. Fourth Berkeley Symposium on Mathematical Statistics and Probability, vol. 1, 1961; reprinted in Breakthroughs in Statistics, Springer 1992, pp. 443-460. reprint
• IMDb Help, Ratings FAQ: weighted average ratings. page
• David Robinson, Introduction to Empirical Bayes: Examples from Baseball Statistics, 2017. book