Bayesian Ordered Probit

Model

The cumulative ordered probit model places a latent continuous quality score $z_{ig}$ behind each discrete rating. A rater from group $g$ (movie) gives rating $k$ whenever the latent score falls in the $k$ -th interval defined by shared thresholds $\boldsymbol{\tau}$ .

z_{ig} \sim \mathcal{N}(\mu_g,\, \sigma_g^2), \qquad y_{ig} = k \iff \tau_{k-1} < z_{ig} \le \tau_k

$\mu_g$ — latent quality (location) for movie $g$
$\sigma_g$ — latent polarisation (scale) for movie $g$
$\boldsymbol{\tau}$ — shared ordered thresholds across all movies

The proportional-odds model ( $\texttt{MASS::polr}$ in R) estimates only a location shift per group, conflating quality and polarisation. The heteroscedastic extension adds a per-group scale $\sigma_g$ : two movies with the same mean rating can differ fundamentally — one beloved by all, one split between fans and detractors.

Identification requires pinning two outer thresholds. For $K=5$ stars: $\tau_1 = 1.5$ and $\tau_4 = 4.5$ are fixed; $\tau_2$ and $\tau_3$ are estimated.

Priors

\mu_g \sim \mathcal{N}(m_0,\, s_0^2), \qquad \sigma_g^{-2} \sim \text{Gamma}(a_0,\, b_0), \qquad \tau_k \sim U(\tau_{k-1},\, \tau_{k+1})

4-Block Gibbs — Albert & Chib (1993)

The 4-block Albert & Chib (1993) Gibbs sampler cycles through exact full conditionals. All blocks are conjugate or have tractable closed-form distributions.

Block	Draw	Distribution
(1)	$z_{ig} \mid \mu_g, \sigma_g, y_{ig}, \boldsymbol{\tau}$	Truncated normal $\text{TN}(\mu_g, \sigma_g^2;\; \tau_{y-1}, \tau_y)$
(2)	$\mu_g \mid \mathbf{z}_g, \sigma_g$	Normal conjugate
(3)	$\sigma_g^{-2} \mid \mathbf{z}_g, \mu_g$	Gamma conjugate
(4)	$\tau_k \mid \mathbf{z}, \mathbf{y}$	Uniform slice (exact full conditional)

Notebooks

Applied to 36 Amazon Prime movies from Liddell & Kruschke (2018), each rated on a 1–5 star scale by a varying number of reviewers. The heteroscedastic model reveals that several movies with similar average ratings differ sharply in polarisation — a distinction invisible to the metric mean or the proportional-odds model.

R Ordered probit / logit — MASS::polr (R) Open → Python Ordered probit — A&C Gibbs, heteroscedastic (Python) Open → Python Ordered probit — heteroscedastic (PyMC) Open →

Downloads

ord_probit_gibbs.py A&C Gibbs sampler module MoviesData.csv 36 Amazon Prime movies — 1–5 star rating counts

Gibbs Sampler — Source Code

"""
Cumulative ordered probit (Albert & Chib 1993) — Gibbs sampler.

Model — heteroscedastic specification (Liddell & Kruschke 2018):
  z_ig ~ N(mu_g, sigma_g^2),   y_ig = k  iff  tau[k-1] < z_ig <= tau[k]

  mu_g     : latent quality location of group g
  sigma_g  : latent quality scale (audience polarisation)

Priors
------
  mu_g        ~ N(m0, s0^2)
  1/sigma_g^2 ~ Gamma(a0, b0)   [shape / rate parameterisation]
  tau_k (free): uniform slice update per A&C §3.3

Identification
--------------
  tau[0]   = LO = 1.5      (pinned)
  tau[K-2] = HI = K - 0.5  (pinned, equals K-0.5 for K categories)
  For K=5: tau in R^4 = (1.5, tau1_free, tau2_free, 4.5)

4-block Gibbs
-------------
  (1)  z_ig | mu_g, sigma_g, y_ig, tau  ->  TN(mu_g, sigma_g^2; lo, hi)
  (2)  mu_g | z_g, sigma_g              ->  Normal conjugate
  (3)  1/sigma_g^2 | z_g, mu_g         ->  Gamma conjugate
  (4)  tau_k (free) | z, y             ->  U(lo_k, hi_k)  slice update

Usage
-----
    from ord_probit_gibbs import fit_ordered_probit, posterior_summary
    out = fit_ordered_probit(counts, R=5000, burn=1000)
    # out['mu']    : (R_kept, G)   latent quality draws
    # out['sigma'] : (R_kept, G)   latent scale draws
    # out['tau']   : (R_kept, K-1) threshold draws (outer two pinned)
"""

import time
import numpy as np
from scipy.special import ndtr, ndtri     # fast standard-normal CDF / PPF


# ── private helpers ───────────────────────────────────────────────────────────

def _rtn(lo, hi, mu, sigma, rng):
    """
    Vectorised draw from N(mu, sigma^2) truncated to (lo_i, hi_i).
    Falls back to interval midpoint when probability mass is negligible.
    """
    lo_s = (lo - mu) / sigma
    hi_s = (hi - mu) / sigma
    p_lo = ndtr(lo_s)
    p_hi = ndtr(hi_s)
    gap  = p_hi - p_lo
    u    = rng.uniform(size=len(lo))
    p    = np.clip(u * gap + p_lo, 1e-15, 1.0 - 1e-15)
    z    = ndtri(p)
    return np.where(gap > 1e-12, z * sigma + mu, (lo + hi) * 0.5)


# ── public helpers ────────────────────────────────────────────────────────────

def posterior_summary(draws, param_names=None, true_vals=None):
    """
    Posterior mean, SD, 95% CrI for each column of a (R, p) draws array.

    Parameters
    ----------
    draws       : (R, p) ndarray
    param_names : length-p list of strings (optional)
    true_vals   : length-p array of true values; adds 'true' and 'covered' cols

    Returns
    -------
    pandas.DataFrame indexed by param_names
    """
    import pandas as pd
    p = draws.shape[1]
    if param_names is None:
        param_names = [f'p{j}' for j in range(p)]
    lo = np.quantile(draws, 0.025, axis=0)
    hi = np.quantile(draws, 0.975, axis=0)
    df = pd.DataFrame({
        'mean':   draws.mean(0),
        'sd':     draws.std(0),
        'q025':   lo,
        'q975':   hi,
    }, index=param_names)
    if true_vals is not None:
        tv = np.asarray(true_vals)
        df['true']    = tv
        df['covered'] = (tv >= lo) & (tv <= hi)
    return df.round(4)


def simulate_ordinal(G, K, mu_true, sigma_true, tau_true, N_per, seed=0):
    """
    Simulate ordinal rating counts from the heteroscedastic model.

    Parameters
    ----------
    G, K         : number of groups, number of categories
    mu_true      : (G,) latent quality locations
    sigma_true   : (G,) latent quality scales
    tau_true     : (K-1,) shared thresholds (including pinned outer two)
    N_per        : ratings per group (int or list of ints)

    Returns
    -------
    counts : (G, K) int array
    """
    rng = np.random.default_rng(seed)
    mu_true    = np.asarray(mu_true,    float)
    sigma_true = np.asarray(sigma_true, float)
    tau_true   = np.asarray(tau_true,   float)
    Ns = [N_per] * G if np.isscalar(N_per) else list(N_per)

    counts = np.zeros((G, K), dtype=int)
    for g in range(G):
        z = rng.normal(mu_true[g], sigma_true[g], Ns[g])
        y = np.searchsorted(tau_true, z) + 1          # 1-based categories
        y = np.clip(y, 1, K)
        counts[g] = np.bincount(y, minlength=K + 1)[1:]
    return counts


# ── main sampler ─────────────────────────────────────────────────────────────

def fit_ordered_probit(
    counts,
    R=5000,
    burn=1000,
    m0=3.0,
    s0=1.5,
    a0=1.0,
    b0=0.44,
    seed=42,
    nprint=500,
):
    """
    4-block Gibbs sampler for the A&C / L&K heteroscedastic ordered probit.

    Parameters
    ----------
    counts  : (G, K) int array — rating counts, one row per group/movie
    R, burn : total iterations; burn-in (R_kept = R - burn draws stored)
    m0, s0  : Normal(m0, s0^2) prior for mu_g
    a0, b0  : Gamma(a0, b0) prior for 1/sigma_g^2  [shape, rate]
              b0=0.44 gives mode sigma_g ~ 1.5
    seed    : int — reproducibility
    nprint  : print progress every nprint iterations (0 = silent)

    Returns
    -------
    dict with:
      'mu'    : (R_kept, G)    posterior draws — latent quality
      'sigma' : (R_kept, G)    posterior draws — latent scale
      'tau'   : (R_kept, K-1)  posterior draws — thresholds
      'G', 'K', 'n_g', 'R_kept'
    """
    rng    = np.random.default_rng(seed)
    counts = np.asarray(counts, dtype=int)
    G, K   = counts.shape
    assert K >= 3, "Need at least 3 ordered categories"

    LO = 1.5
    HI = float(K) - 0.5

    # ── expand counts to per-group observation vectors ────────────────────────
    cats   = np.arange(1, K + 1)
    y_list = [np.repeat(cats, counts[g]) for g in range(G)]
    n_g    = counts.sum(axis=1)

    # ── initialise ────────────────────────────────────────────────────────────
    tau     = np.linspace(LO, HI, K - 1)           # evenly spaced start
    mu      = np.full(G, m0)
    sigma   = np.full(G, 1.5)

    tau_ext = np.r_[-np.inf, tau, np.inf]
    z_list  = [(tau_ext[y - 1] + tau_ext[y]) / 2.0 for y in y_list]

    # ── storage ───────────────────────────────────────────────────────────────
    R_kept  = R - burn
    mu_d    = np.empty((R_kept, G))
    sigma_d = np.empty((R_kept, G))
    tau_d   = np.empty((R_kept, K - 1))
    kept    = 0
    t0      = time.time()

    for r in range(1, R + 1):
        tau_ext = np.r_[-np.inf, tau, np.inf]

        # ── Block 1: z_ig | mu_g, sigma_g, y, tau  (truncated normal) ────────
        for g in range(G):
            y         = y_list[g]
            z_list[g] = _rtn(tau_ext[y - 1], tau_ext[y], mu[g], sigma[g], rng)

        # ── Block 2: mu_g | z_g, sigma_g  (Normal-Normal conjugate) ──────────
        for g in range(G):
            prec  = 1.0 / s0 ** 2 + n_g[g] / sigma[g] ** 2
            mn    = (m0 / s0 ** 2 + z_list[g].sum() / sigma[g] ** 2) / prec
            mu[g] = rng.normal(mn, 1.0 / np.sqrt(prec))

        # ── Block 3: 1/sigma_g^2 | z_g, mu_g  (Normal-Gamma conjugate) ───────
        for g in range(G):
            ss       = np.sum((z_list[g] - mu[g]) ** 2)
            tau_g    = rng.gamma(a0 + n_g[g] / 2.0, 1.0 / (b0 + ss / 2.0))
            sigma[g] = 1.0 / np.sqrt(tau_g)

        # ── Block 4: free thresholds | z, y  (A&C uniform slice) ─────────────
        if K >= 4:
            z_all = np.concatenate(z_list)
            y_all = np.concatenate(y_list)
            # tau is 0-indexed; pinned: tau[0]=LO, tau[K-2]=HI
            # free: indices 1 .. K-3  (for K=5: indices 1 and 2)
            for idx in range(1, K - 2):
                cat_lo = idx + 1          # 1-based category whose upper bound is tau[idx]
                cat_hi = idx + 2          # 1-based category whose lower bound is tau[idx]
                lo_k   = tau[idx - 1]
                hi_k   = tau[idx + 1]
                m_lo   = y_all == cat_lo
                m_hi   = y_all == cat_hi
                if m_lo.any():
                    lo_k = max(lo_k, z_all[m_lo].max())
                if m_hi.any():
                    hi_k = min(hi_k, z_all[m_hi].min())
                if lo_k < hi_k:
                    tau[idx] = rng.uniform(lo_k, hi_k)

        if r > burn:
            mu_d[kept]    = mu
            sigma_d[kept] = sigma
            tau_d[kept]   = tau
            kept += 1

        if nprint > 0 and r % nprint == 0:
            elapsed = time.time() - t0
            print(f'  Iter {r:5d}/{R}  ({elapsed:.0f}s elapsed)')

    elapsed = time.time() - t0
    print(f'Done in {elapsed / 60:.1f} min  |  kept draws: {R_kept}')

    return {
        'mu': mu_d, 'sigma': sigma_d, 'tau': tau_d,
        'G': G, 'K': K, 'n_g': n_g, 'R_kept': R_kept,
    }

References

Albert, J. H. & Chib, S. (1993). Bayesian analysis of binary and polychotomous response data. Journal of the American Statistical Association 88(422), 669–679. — the data-augmentation Gibbs sampler for the ordered probit, including the threshold full conditional
Liddell, T. M. & Kruschke, J. K. (2018). Analyzing ordinal data with metric models: what could possibly go wrong? Journal of Experimental Social Psychology 79, 328–348. — the movie-ratings example and the case for a per-group location and scale
McCullagh, P. (1980). Regression models for ordinal data. Journal of the Royal Statistical Society: Series B 42(2), 109–142. — the cumulative-link / proportional-odds model that the R polr fit represents
Venables, W. N. & Ripley, B. D. (2002). Modern Applied Statistics with S (4th ed.). Springer. — MASS::polr, the homoscedastic proportional-odds MLE used as the R comparison
Johnson, V. E. & Albert, J. H. (1999). Ordinal Data Modeling. Springer. — book-length treatment of Bayesian ordinal models and latent-scale identification