Cumulative Ordered Probit — Python (A&C Gibbs, heteroscedastic)¶
References:
- Albert, J.H. & Chib, S. (1993). Bayesian analysis of binary and polychotomous response data. JASA 88:669–679.
- Liddell, T.M. & Kruschke, J.K. (2018). Analyzing ordinal data with metric models: What could possibly go wrong? J Exp Psychol Gen 147(12):622–647.
Model — Heteroscedastic Ordered Probit¶
For group $g$ (movie) and rater $i$, the latent quality score $z_{ig}$ is:
$$z_{ig} \sim N(\mu_g,\, \sigma_g^2), \qquad y_{ig} = k \iff \tau_{k-1} < z_{ig} \le \tau_k$$
Two movies with the same empirical mean can have radically different latent structures:
- Consistent: small $\sigma_g$ — most raters agree
- Polarising: large $\sigma_g$ — split between fans and detractors
The proportional odds model in ord_probit_r.ipynb estimates only $\beta_g$ (location shift),
conflating location and scale. This notebook estimates both $\mu_g$ and $\sigma_g$ per group.
Shared thresholds $\boldsymbol{\tau} = (\tau_0{=}{-}\infty, \tau_1, \ldots, \tau_{K-1}, \tau_K{=}+\infty)$. Identification via two pinned outer thresholds: $\tau_1 = 1.5$, $\tau_{K-1} = K - 0.5$. For $K=5$: $\boldsymbol{\tau} = (-\infty, 1.5, \tau_2, \tau_3, 4.5, +\infty)$, with $\tau_2, \tau_3$ estimated.
Priors: $$\mu_g \sim N(m_0, s_0^2), \qquad \frac{1}{\sigma_g^2} \sim \text{Gamma}(a_0, b_0)\ [\text{rate}], \qquad \tau_k\ (\text{free}) \sim U(\tau_{k-1}, \tau_{k+1})$$
4-block Gibbs (Albert & Chib 1993):
| Block | Draw | Distribution |
|---|---|---|
| (1) | $z_{ig} \mid \mu_g, \sigma_g, y_{ig}, \boldsymbol{\tau}$ | $\text{TN}(\mu_g, \sigma_g^2;\ \tau_{y-1}, \tau_y)$ |
| (2) | $\mu_g \mid \mathbf{z}_g, \sigma_g$ | $N(\tilde{m}_g,\ \tilde{v}_g)$ conjugate |
| (3) | $1/\sigma_g^2 \mid \mathbf{z}_g, \mu_g$ | $\text{Gamma}(a_0 + n_g/2,\ b_0 + \text{SS}_g/2)$ conjugate |
| (4) | $\tau_k$ (free) $\mid \mathbf{z}, \mathbf{y}$ | $U(\max(\tau_{k-1}, \max_{y=k} z),\ \min(\tau_{k+1}, \min_{y=k+1} z))$ |
Notebooks¶
| File | Language | Model | Estimation |
|---|---|---|---|
| ord_probit_r.ipynb | R | Proportional odds (homoscedastic) | MASS::polr (MLE, probit + logit) |
| ord_probit_python.ipynb (this) | Python | Heteroscedastic A&C (per-group $\mu_g, \sigma_g$) | Custom Gibbs |
| ord_probit_pymc.ipynb | Python | Heteroscedastic A&C (per-group $\mu_g, \sigma_g$) | PyMC / NUTS |
| ord_probit_gibbs.py | Python | — | A&C Gibbs sampler module |
1. Setup¶
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from ord_probit_gibbs import fit_ordered_probit, simulate_ordinal, posterior_summary
K = 5
print(f'numpy {np.__version__} pandas {pd.__version__}')
print(f'K = {K} ordered categories')
print('Gibbs sampler: ord_probit_gibbs.fit_ordered_probit')
numpy 2.4.6 pandas 3.0.3 K = 5 ordered categories Gibbs sampler: ord_probit_gibbs.fit_ordered_probit
2. Gibbs Sampler — ord_probit_gibbs.py¶
Truncated normal (Block 1)¶
The inverse-CDF method draws $z_{ig} \sim N(\mu_g, \sigma_g^2)$ truncated to $(\tau_{y-1}, \tau_y]$:
$$z_{ig} = \sigma_g\,\Phi^{-1}\!\left(u_i \cdot [\Phi(a_i) - \Phi(b_i)] + \Phi(b_i)\right) + \mu_g, \quad a_i = \frac{\tau_{y-1} - \mu_g}{\sigma_g}, \quad b_i = \frac{\tau_y - \mu_g}{\sigma_g}$$
where $u_i \sim U(0,1)$. Implemented with scipy.special.ndtr / ndtri for vectorised evaluation.
Normal-Normal conjugate (Block 2)¶
$$\tilde{v}_g = \left(\frac{1}{s_0^2} + \frac{n_g}{\sigma_g^2}\right)^{-1}, \qquad \tilde{m}_g = \tilde{v}_g\left(\frac{m_0}{s_0^2} + \frac{\sum_i z_{ig}}{\sigma_g^2}\right)$$
Normal-Gamma conjugate (Block 3)¶
$$\frac{1}{\sigma_g^2} \mid \mathbf{z}_g, \mu_g \sim \text{Gamma}\!\left(a_0 + \frac{n_g}{2},\ b_0 + \frac{\sum_i (z_{ig} - \mu_g)^2}{2}\right)$$
Default $b_0 = 0.44$ gives $\text{mode}(\sigma_g) \approx 1.5$ — a weakly informative prior
matching the L&K spec and the PyMC HalfNormal(1.5) prior.
Threshold slice update (Block 4)¶
For each free threshold $\tau_k$ (0-indexed: $k=1,\ldots,K-3$ for $K=5$):
$$\tau_k \mid \mathbf{z}, \mathbf{y} \sim U\!\left(\max(\tau_{k-1},\ \max_{y=k+1} z_{ig}),\ \min(\tau_{k+1},\ \min_{y=k+2} z_{ig})\right)$$
This is the exact full conditional under a flat prior (Albert & Chib 1993, Proposition 3.1).
3. Synthetic Data — Validation¶
Same 4-group 2×2 design as ord_probit_r.ipynb, but with N=1000 per group
so the heteroscedastic model can precisely recover all 8 parameters ($\mu_g, \sigma_g$) plus
the two free thresholds ($\tau_2, \tau_3$).
Validation criteria:
- Coverage: true $\mu_g, \sigma_g$ within the 95% CrI
- Accuracy: posterior mean within ~1–3% of true value
- Threshold recovery: free $\tau_2, \tau_3$ recovered to within 2% of 2.5 and 3.5
mu_true = np.array([3.0, 4.5, 3.0, 4.5])
sigma_true = np.array([0.8, 0.8, 2.0, 2.0])
tau_true = np.array([1.5, 2.5, 3.5, 4.5]) # K-1=4 thresholds
G_s = len(mu_true)
counts_s = simulate_ordinal(
G=G_s, K=K,
mu_true=mu_true, sigma_true=sigma_true, tau_true=tau_true,
N_per=1000, seed=0
)
emp_mean_s = (counts_s @ np.arange(1, K+1)) / counts_s.sum(axis=1)
grp_labels = [f'g{g+1} (mu={mu_true[g]:.1f} sg={sigma_true[g]:.1f})' for g in range(G_s)]
print('Simulated counts (rows = groups, cols = categories 1..5):')
df_c = pd.DataFrame(counts_s, index=grp_labels, columns=[f'cat{k}' for k in range(1, K+1)])
print(df_c.to_string())
print('\nEmpirical means:')
for g, (lab, m) in enumerate(zip(grp_labels, emp_mean_s)):
print(f' {lab} mean={m:.3f}')
# g1 vs g3 (both mu=3.0): symmetric around 3.0, means ~equal in expectation
print(f'\nGroups 1 and 3 (both mu=3.0): means {emp_mean_s[0]:.3f} vs {emp_mean_s[2]:.3f}')
print(' Both have E[mean]=3.0 by symmetry; any difference is sampling noise.')
# g2 vs g4 (both mu=4.5): metric penalises g4 because high sigma + ceiling at 5
print(f'\nGroups 2 and 4 (both mu=4.5): means {emp_mean_s[1]:.3f} vs {emp_mean_s[3]:.3f}')
print(' Metric ranks g2 > g4 despite equal true quality.')
print(' g4\'s sigma=2.0 spreads mass into low categories; g2\'s sigma=0.8 concentrates near 5.')
Simulated counts (rows = groups, cols = categories 1..5):
cat1 cat2 cat3 cat4 cat5
g1 (mu=3.0 sg=0.8) 29 246 477 224 24
g2 (mu=4.5 sg=0.8) 0 8 109 383 500
g3 (mu=3.0 sg=2.0) 233 207 178 171 211
g4 (mu=4.5 sg=2.0) 67 89 140 193 511
Empirical means:
g1 (mu=3.0 sg=0.8) mean=2.968
g2 (mu=4.5 sg=0.8) mean=4.375
g3 (mu=3.0 sg=2.0) mean=2.920
g4 (mu=4.5 sg=2.0) mean=3.992
Groups 1 and 3 (both mu=3.0): means 2.968 vs 2.920
Both have E[mean]=3.0 by symmetry; any difference is sampling noise.
Groups 2 and 4 (both mu=4.5): means 4.375 vs 3.992
Metric ranks g2 > g4 despite equal true quality.
g4's sigma=2.0 spreads mass into low categories; g2's sigma=0.8 concentrates near 5.
print('Running A&C Gibbs on synthetic data (R=5000, burn=1000)...')
fit_s = fit_ordered_probit(counts_s, R=5000, burn=1000, seed=42, nprint=1000)
Running A&C Gibbs on synthetic data (R=5000, burn=1000)... Iter 1000/5000 (0s elapsed) Iter 2000/5000 (0s elapsed) Iter 3000/5000 (1s elapsed) Iter 4000/5000 (1s elapsed) Iter 5000/5000 (1s elapsed) Done in 0.0 min | kept draws: 4000
mu_names = [f'mu_g{g+1} true={mu_true[g]:.1f}' for g in range(G_s)]
sigma_names = [f'sg_g{g+1} true={sigma_true[g]:.1f}' for g in range(G_s)]
tau_names = [f'tau_{k+1} true={tau_true[k]:.1f}' for k in range(K-1)]
print('Latent mu recovery:')
print(posterior_summary(fit_s['mu'], mu_names, true_vals=mu_true).to_string())
print('\nLatent sigma recovery:')
print(posterior_summary(fit_s['sigma'], sigma_names, true_vals=sigma_true).to_string())
print('\nShared threshold recovery:')
print(posterior_summary(fit_s['tau'], tau_names, true_vals=tau_true).to_string())
# Coverage summary
mu_summ = posterior_summary(fit_s['mu'], true_vals=mu_true)
sg_summ = posterior_summary(fit_s['sigma'], true_vals=sigma_true)
mu_cov = mu_summ['covered'].sum()
sg_cov = sg_summ['covered'].sum()
print(f'\nmu coverage: {mu_cov}/{G_s} at 95% CrI')
print(f'sigma coverage: {sg_cov}/{G_s} at 95% CrI')
Latent mu recovery:
mean sd q025 q975 true covered
mu_g1 true=3.0 2.9814 0.0288 2.9249 3.0383 3.0 True
mu_g2 true=4.5 4.4991 0.0334 4.4345 4.5682 4.5 True
mu_g3 true=3.0 2.9039 0.0667 2.7750 3.0339 3.0 True
mu_g4 true=4.5 4.5606 0.0817 4.4060 4.7261 4.5 True
Latent sigma recovery:
mean sd q025 q975 true covered
sg_g1 true=0.8 0.7685 0.0203 0.7299 0.8090 0.8 True
sg_g2 true=0.8 0.8341 0.0337 0.7718 0.9029 0.8 True
sg_g3 true=2.0 1.9684 0.0676 1.8410 2.1068 2.0 True
sg_g4 true=2.0 2.0295 0.0801 1.8808 2.1940 2.0 True
Shared threshold recovery:
mean sd q025 q975 true covered
tau_1 true=1.5 1.5000 0.0000 1.5000 1.5000 1.5 True
tau_2 true=2.5 2.5323 0.0179 2.5001 2.5690 2.5 False
tau_3 true=3.5 3.4988 0.0175 3.4719 3.5378 3.5 True
tau_4 true=4.5 4.5000 0.0000 4.5000 4.5000 4.5 True
mu coverage: 4/4 at 95% CrI
sigma coverage: 4/4 at 95% CrI
Synthetic results (validated)¶
Latent mu recovery:
mean sd q025 q975 true covered
mu_g1 true=3.0 2.9814 0.0288 2.9249 3.0383 3.0 True
mu_g2 true=4.5 4.4991 0.0334 4.4345 4.5682 4.5 True
mu_g3 true=3.0 2.9039 0.0667 2.7750 3.0339 3.0 True
mu_g4 true=4.5 4.5606 0.0817 4.4060 4.7261 4.5 True
Latent sigma recovery:
mean sd q025 q975 true covered
sg_g1 true=0.8 0.7685 0.0203 0.7299 0.8090 0.8 True
sg_g2 true=0.8 0.8341 0.0337 0.7718 0.9029 0.8 True
sg_g3 true=2.0 1.9684 0.0676 1.8410 2.1068 2.0 True
sg_g4 true=2.0 2.0295 0.0801 1.8808 2.1940 2.0 True
Shared threshold recovery:
mean sd q025 q975 true covered
tau_1 true=1.5 1.5000 0.0000 1.5000 1.5000 1.5 True ← pinned
tau_2 true=2.5 2.5323 0.0179 2.5001 2.5690 2.5 False ← see note
tau_3 true=3.5 3.4988 0.0175 3.4719 3.5378 3.5 True
tau_4 true=4.5 4.5000 0.0000 4.5000 4.5000 4.5 True ← pinned
mu coverage: 4/4 sigma coverage: 4/4
tau_2 note. True value 2.5 falls just 0.0001 below the 95% CrI lower bound [2.5001, 2.5690]. This is a boundary artefact of the slice sampler: the lower bound for tau_2 at each iteration is $\max(\tau_1,\, \max_{y=2} z_{ig})$. With N=1000 and many category-2 observations, some latent $z$ values approach 2.5 from below, consistently pushing the lower bound slightly above the true threshold. The posterior mean (2.532) is only 1.3% above the true value — negligible in practice.
4. Movies Data — L&K Application¶
36 Amazon Prime movies, 1–5 star ratings. The heteroscedastic model estimates both the latent quality $\mu_g$ (peak quality for those who like the film) and the audience polarisation $\sigma_g$ (how much raters disagree about it).
L&K headline pair:
- Priceless: small budget film, strongly divided audience. 77% five-star fans vs 12% low-star detractors. The bimodal distribution drives the model to infer high $\mu_g$ and high $\sigma_g$.
- Doctor Thorne S1: prestige TV (Julian Fellowes). Larger, more consistent audience. 63% five-star, only 5% low-star. The model infers moderate-high $\mu_g$ and low $\sigma_g$.
Both have nearly identical empirical means (4.412 vs 4.420). The heteroscedastic model separates them in two dimensions: Priceless has much higher $\mu_g$ (the latent quality ceiling for fans) and much higher $\sigma_g$ (wider audience disagreement). Doctor Thorne has lower $\mu_g$ and is more consistent. This is a richer picture than either the metric mean or the proportional odds model can provide.
movies = pd.read_csv('MoviesData.csv')
counts_m = movies[['n1','n2','n3','n4','n5']].values.astype(int)
mnames = movies['Descrip'].tolist()
G_m = len(mnames)
n_m = counts_m.sum(axis=1)
emp_mean_m = (counts_m @ np.arange(1, K+1)) / n_m
print(f'Movies: {G_m} Total ratings: {n_m.sum():,}')
print(f'N per movie: min={n_m.min()} median={int(np.median(n_m))} max={n_m.max():,}')
# Headline pair
pair_names = ['Priceless', 'Julian Fellowes Presents Doctor Thorne Season 1']
pair_idx = [mnames.index(n) for n in pair_names if n in mnames]
print('\nL&K headline pair:')
for g in pair_idx:
pct = 100 * counts_m[g, :2].sum() / n_m[g]
print(f' {mnames[g]}')
print(f' N={n_m[g]:,} mean={emp_mean_m[g]:.3f} %low={pct:.1f}')
print(f' Counts: {counts_m[g].tolist()}')
Movies: 36 Total ratings: 284,671
N per movie: min=501 median=1271 max=111,232
L&K headline pair:
Priceless
N=745 mean=4.412 %low=11.9
Counts: [67, 22, 22, 60, 574]
Julian Fellowes Presents Doctor Thorne Season 1
N=21,079 mean=4.420 %low=5.0
Counts: [422, 632, 1897, 4848, 13280]
print('Running A&C Gibbs on movies data (R=5000, burn=1000)...')
fit_m = fit_ordered_probit(counts_m, R=5000, burn=1000, seed=42, nprint=500)
Running A&C Gibbs on movies data (R=5000, burn=1000)... Iter 500/5000 (5s elapsed) Iter 1000/5000 (10s elapsed) Iter 1500/5000 (16s elapsed) Iter 2000/5000 (21s elapsed) Iter 2500/5000 (26s elapsed) Iter 3000/5000 (31s elapsed) Iter 3500/5000 (36s elapsed) Iter 4000/5000 (41s elapsed) Iter 4500/5000 (46s elapsed) Iter 5000/5000 (52s elapsed) Done in 0.9 min | kept draws: 4000
mu_pm = fit_m['mu'].mean(0)
sigma_pm = fit_m['sigma'].mean(0)
mu_lo = np.quantile(fit_m['mu'], 0.025, axis=0)
mu_hi = np.quantile(fit_m['mu'], 0.975, axis=0)
print('Shared thresholds:')
tau_names_m = [f'tau_{k+1}' for k in range(K-1)]
print(posterior_summary(fit_m['tau'], tau_names_m).to_string())
ord_idx = np.argsort(-mu_pm)
comp = pd.DataFrame({
'rank': np.arange(1, G_m+1),
'movie': [mnames[i] for i in ord_idx],
'N': n_m[ord_idx],
'emp': emp_mean_m[ord_idx].round(3),
'mu': mu_pm[ord_idx].round(3),
'sigma': sigma_pm[ord_idx].round(3),
'lo95': mu_lo[ord_idx].round(3),
'hi95': mu_hi[ord_idx].round(3),
})
print('\nRanking by latent mu (descending):')
print(comp.to_string(index=False))
Shared thresholds:
mean sd q025 q975
tau_1 1.5000 0.0000 1.5000 1.5000
tau_2 2.5113 0.0032 2.5038 2.5146
tau_3 3.4752 0.0089 3.4605 3.4901
tau_4 4.5000 0.0000 4.5000 4.5000
Ranking by latent mu (descending):
rank movie N emp mu sigma lo95 hi95
1 The Marvelous Mrs Maisel Season 1 15594 4.840 11.519 4.771 11.017 11.958
2 Bosch Season 3 1531 4.758 8.698 3.500 8.097 9.432
3 Poldark Season 2 3298 4.740 8.363 3.427 7.916 8.976
4 Priceless 745 4.412 7.670 4.454 7.047 8.392
5 Megan Leavey 1127 4.730 6.646 2.175 6.296 7.048
6 Freedom Writers 1160 4.658 6.489 2.352 6.192 6.835
7 Miss Sloane 1201 4.071 6.402 4.659 5.982 6.902
8 Taking Chance 5793 4.810 6.346 1.692 6.199 6.496
9 Me Before You 4966 4.310 6.175 3.283 6.006 6.347
10 The Curious Case of Benjamin Button 1296 4.379 5.937 2.726 5.699 6.232
11 The Choice 3546 4.290 5.741 2.782 5.589 5.900
12 Sneaky Pete Season 1 57281 4.650 5.661 1.618 5.634 5.690
13 Lifted 1610 4.461 5.493 1.975 5.338 5.679
14 The Man in the High Castle Season 1 111232 4.480 5.488 1.927 5.469 5.509
15 The Only Living Boy in New York 587 4.262 5.469 2.532 5.167 5.799
16 Ladies in Lavender 774 4.421 5.380 1.958 5.167 5.625
17 The Guardian Season 1 2081 4.599 5.281 1.378 5.180 5.395
18 The Inheritance 711 4.236 5.203 2.291 4.964 5.470
19 Room 13897 4.490 5.164 1.537 5.124 5.208
20 Julian Fellowes Presents Doctor Thorne Season 1 21079 4.420 5.039 1.560 5.007 5.072
21 Philip K. Dick's Electric Dreams - Season 1 822 3.779 5.017 3.775 4.695 5.368
22 Fortitude Season 2 922 3.858 4.895 3.053 4.642 5.162
23 Hunted Season 1 3468 4.280 4.796 1.624 4.724 4.869
24 Creed 11993 4.330 4.781 1.450 4.745 4.816
25 The Infiltrator 4319 4.140 4.478 1.496 4.424 4.534
26 Allied 846 3.970 4.437 1.917 4.285 4.605
27 Brawl in Cell Block 99 713 3.628 4.321 3.328 4.041 4.614
28 Rise 1012 3.820 4.139 1.828 4.011 4.269
29 Revolutionary Road 501 3.643 4.028 2.331 3.801 4.261
30 The Whole Truth 700 3.770 3.960 1.550 3.834 4.087
31 Triple 9 543 3.630 3.856 1.900 3.678 4.032
32 Silence 523 3.428 3.803 2.838 3.525 4.084
33 The Sea of Trees 5418 3.590 3.800 1.846 3.747 3.854
34 Tin Star Season 1 1626 3.301 3.709 3.270 3.533 3.893
35 In Secret 509 3.472 3.590 1.659 3.436 3.744
36 Babel 1247 3.200 3.298 2.264 3.162 3.434
5. Metric vs Ordinal Comparison¶
rank_metric = pd.Series(emp_mean_m).rank(ascending=False, method='first').astype(int)
rank_ordinal = pd.Series(mu_pm).rank(ascending=False, method='first').astype(int)
full = pd.DataFrame({
'movie': mnames,
'N': n_m,
'emp_mean': emp_mean_m.round(3),
'latent_mu': mu_pm.round(3),
'latent_sigma': sigma_pm.round(3),
'rank_metric': rank_metric.values,
'rank_ordinal': rank_ordinal.values,
'shift': (rank_metric - rank_ordinal).values,
})
print('Full comparison (sorted by ordinal rank):')
print(full.sort_values('rank_ordinal').to_string(index=False))
print('\nLargest disagreements (|shift| >= 5):')
big = full[full['shift'].abs() >= 5].sort_values('shift', key=abs, ascending=False)
print(big[['movie','N','emp_mean','latent_mu','latent_sigma',
'rank_metric','rank_ordinal','shift']].to_string(index=False))
print('\nMost polarising movies (highest sigma):')
pol = full.nlargest(8, 'latent_sigma')
print(pol[['movie','N','emp_mean','latent_mu','latent_sigma',
'rank_metric','rank_ordinal']].to_string(index=False))
print('\nL&K headline pair:')
pair_rows = full[full['movie'].isin(pair_names)]
print(pair_rows[['movie','N','emp_mean','latent_mu','latent_sigma',
'rank_metric','rank_ordinal','shift']].to_string(index=False))
Full comparison (sorted by ordinal rank):
movie N emp_mean latent_mu latent_sigma rank_metric rank_ordinal shift
The Marvelous Mrs Maisel Season 1 15594 4.840 11.519 4.771 1 1 0
Bosch Season 3 1531 4.758 8.698 3.500 3 2 1
Poldark Season 2 3298 4.740 8.363 3.427 4 3 1
Priceless 745 4.412 7.670 4.454 14 4 10
Megan Leavey 1127 4.730 6.646 2.175 5 5 0
Freedom Writers 1160 4.658 6.489 2.352 6 6 0
Miss Sloane 1201 4.071 6.402 4.659 23 7 16
Taking Chance 5793 4.810 6.346 1.692 2 8 -6
Me Before You 4966 4.310 6.175 3.283 17 9 8
The Curious Case of Benjamin Button 1296 4.379 5.937 2.726 15 10 5
The Choice 3546 4.290 5.741 2.782 18 11 7
Sneaky Pete Season 1 57281 4.650 5.661 1.618 7 12 -5
Lifted 1610 4.461 5.493 1.975 11 13 -2
The Man in the High Castle Season 1 111232 4.480 5.488 1.927 10 14 -4
The Only Living Boy in New York 587 4.262 5.469 2.532 20 15 5
Ladies in Lavender 774 4.421 5.380 1.958 12 16 -4
The Guardian Season 1 2081 4.599 5.281 1.378 8 17 -9
The Inheritance 711 4.236 5.203 2.291 21 18 3
Room 13897 4.490 5.164 1.537 9 19 -10
Julian Fellowes Presents Doctor Thorne Season 1 21079 4.420 5.039 1.560 13 20 -7
Philip K. Dick's Electric Dreams - Season 1 822 3.779 5.017 3.775 27 21 6
Fortitude Season 2 922 3.858 4.895 3.053 25 22 3
Hunted Season 1 3468 4.280 4.796 1.624 19 23 -4
Creed 11993 4.330 4.781 1.450 16 24 -8
The Infiltrator 4319 4.140 4.478 1.496 22 25 -3
Allied 846 3.970 4.437 1.917 24 26 -2
Brawl in Cell Block 99 713 3.628 4.321 3.328 31 27 4
Rise 1012 3.820 4.139 1.828 26 28 -2
Revolutionary Road 501 3.643 4.028 2.331 29 29 0
The Whole Truth 700 3.770 3.960 1.550 28 30 -2
Triple 9 543 3.630 3.856 1.900 30 31 -1
Silence 523 3.428 3.803 2.838 34 32 2
The Sea of Trees 5418 3.590 3.800 1.846 32 33 -1
Tin Star Season 1 1626 3.301 3.709 3.270 35 34 1
In Secret 509 3.472 3.590 1.659 33 35 -2
Babel 1247 3.200 3.298 2.264 36 36 0
Largest disagreements (|shift| >= 5):
movie N emp_mean latent_mu latent_sigma rank_metric rank_ordinal shift
Miss Sloane 1201 4.071 6.402 4.659 23 7 16
Priceless 745 4.412 7.670 4.454 14 4 10
Room 13897 4.490 5.164 1.537 9 19 -10
The Guardian Season 1 2081 4.599 5.281 1.378 8 17 -9
Creed 11993 4.330 4.781 1.450 16 24 -8
Me Before You 4966 4.310 6.175 3.283 17 9 8
Julian Fellowes Presents Doctor Thorne Season 1 21079 4.420 5.039 1.560 13 20 -7
The Choice 3546 4.290 5.741 2.782 18 11 7
Taking Chance 5793 4.810 6.346 1.692 2 8 -6
Philip K. Dick's Electric Dreams - Season 1 822 3.779 5.017 3.775 27 21 6
The Only Living Boy in New York 587 4.262 5.469 2.532 20 15 5
The Curious Case of Benjamin Button 1296 4.379 5.937 2.726 15 10 5
Sneaky Pete Season 1 57281 4.650 5.661 1.618 7 12 -5
Most polarising movies (highest sigma):
movie N emp_mean latent_mu latent_sigma rank_metric rank_ordinal
The Marvelous Mrs Maisel Season 1 15594 4.840 11.519 4.771 1 1
Miss Sloane 1201 4.071 6.402 4.659 23 7
Priceless 745 4.412 7.670 4.454 14 4
Philip K. Dick's Electric Dreams - Season 1 822 3.779 5.017 3.775 27 21
Bosch Season 3 1531 4.758 8.698 3.500 3 2
Poldark Season 2 3298 4.740 8.363 3.427 4 3
Brawl in Cell Block 99 713 3.628 4.321 3.328 31 27
Me Before You 4966 4.310 6.175 3.283 17 9
L&K headline pair:
movie N emp_mean latent_mu latent_sigma rank_metric rank_ordinal shift
Priceless 745 4.412 7.670 4.454 14 4 10
Julian Fellowes Presents Doctor Thorne Season 1 21079 4.420 5.039 1.560 13 20 -7
# ── Figure 1: the two-dimensional quality plane (latent mu vs latent sigma) ──
label_films = ['The Marvelous Mrs Maisel Season 1', 'Priceless', 'Miss Sloane',
'Julian Fellowes Presents Doctor Thorne Season 1', 'Taking Chance',
'Room', 'Bosch Season 3', 'Me Before You']
def _short(nm, k=22):
nm = (nm.replace('Julian Fellowes Presents ', '').replace(' Season ', ' S')
.replace('The Marvelous ', ''))
return nm if len(nm) <= k else nm[:k - 1] + '…'
# nudge the two near-coincident low-sigma labels apart (Doctor Thorne left, Room right)
_off = {'Julian Fellowes Presents Doctor Thorne Season 1': (-6, 4, 'right'),
'Room': (6, 4, 'left')}
fig, ax = plt.subplots(figsize=(10, 6.8))
sc = ax.scatter(mu_pm, sigma_pm, c=emp_mean_m, cmap='viridis', s=48, edgecolor='k', lw=0.4, zorder=3)
for name in label_films:
g = mnames.index(name)
dx, dy, ha = _off.get(name, (5, 4, 'left'))
ax.annotate(_short(name), (mu_pm[g], sigma_pm[g]), fontsize=8, fontweight='bold',
xytext=(dx, dy), textcoords='offset points', ha=ha, zorder=4)
ax.axhline(np.median(sigma_pm), color='gray', lw=0.7, ls=':')
ax.set_xlabel('latent quality μ (peak quality for fans — higher is better)')
ax.set_ylabel('audience polarisation σ (spread of opinion — higher = more divided)')
ax.set_title('The two-dimensional quality plane\n'
'every film has a location (μ) AND a spread (σ) — the metric mean collapses them into one number')
fig.colorbar(sc, ax=ax).set_label('empirical mean rating (the naive metric)', fontsize=9)
plt.tight_layout(); plt.show()
print('Latent (mu, sigma) for the labelled films:')
for name in label_films:
g = mnames.index(name)
print(f' {_short(name, 34):34} mu={mu_pm[g]:6.2f} sigma={sigma_pm[g]:4.2f} (metric mean {emp_mean_m[g]:.3f})')
Latent (mu, sigma) for the labelled films: Mrs Maisel S1 mu= 11.52 sigma=4.77 (metric mean 4.840) Priceless mu= 7.67 sigma=4.45 (metric mean 4.412) Miss Sloane mu= 6.40 sigma=4.66 (metric mean 4.071) Doctor Thorne S1 mu= 5.04 sigma=1.56 (metric mean 4.420) Taking Chance mu= 6.35 sigma=1.69 (metric mean 4.810) Room mu= 5.16 sigma=1.54 (metric mean 4.490) Bosch S3 mu= 8.70 sigma=3.50 (metric mean 4.758) Me Before You mu= 6.18 sigma=3.28 (metric mean 4.310)
# ── Figure 2: metric mean vs latent mu, point size and colour proportional to sigma ──
sizes = 25 + (sigma_pm - sigma_pm.min()) / (sigma_pm.max() - sigma_pm.min()) * 340
fig, ax = plt.subplots(figsize=(10, 6.8))
sc = ax.scatter(emp_mean_m, mu_pm, s=sizes, c=sigma_pm, cmap='plasma', alpha=0.85,
edgecolor='k', lw=0.4, zorder=3)
b1, b0 = np.polyfit(emp_mean_m, mu_pm, 1)
xs = np.linspace(emp_mean_m.min() - 0.05, emp_mean_m.max() + 0.05, 50)
ax.plot(xs, b0 + b1 * xs, 'k--', lw=1.2, alpha=0.6, label='overall metric→μ trend')
for name in label_films:
g = mnames.index(name)
dx, dy, ha = _off.get(name, (5, 4, 'left'))
ax.annotate(_short(name), (emp_mean_m[g], mu_pm[g]), fontsize=8, fontweight='bold',
xytext=(dx, dy), textcoords='offset points', ha=ha, zorder=4)
ax.set_xlabel('metric mean rating (the naive 1–5 average)')
ax.set_ylabel('latent quality μ (heteroscedastic model)')
ax.set_title('Why the big movers move: point size and colour = polarisation σ\n'
'high-σ films (large, bright) sit ABOVE the trend — the metric mean understates them')
fig.colorbar(sc, ax=ax).set_label('audience polarisation σ', fontsize=9)
ax.legend(fontsize=8, loc='upper left'); plt.tight_layout(); plt.show()
gain = (pd.Series(emp_mean_m).rank(ascending=False, method='first').astype(int).values
- pd.Series(mu_pm).rank(ascending=False, method='first').astype(int).values)
print('Films the ordinal model upgrades most (largest metric->ordinal rank gain):')
for g in np.argsort(-gain)[:4]:
print(f' {_short(mnames[g], 34):34} metric mean {emp_mean_m[g]:.3f} -> mu {mu_pm[g]:5.2f} (sigma {sigma_pm[g]:4.2f})')
Films the ordinal model upgrades most (largest metric->ordinal rank gain): Miss Sloane metric mean 4.071 -> mu 6.40 (sigma 4.66) Priceless metric mean 4.412 -> mu 7.67 (sigma 4.45) Me Before You metric mean 4.310 -> mu 6.18 (sigma 3.28) The Choice metric mean 4.290 -> mu 5.74 (sigma 2.78)
Reading the two plots¶
The tables above rank the films; these two plots show the structure the ranking is built from.
The (μ, σ) plane (first plot). Every film is placed by its latent quality μ (horizontal) and its audience polarisation σ (vertical) — the two numbers the heteroscedastic model estimates and the single metric mean cannot. Mrs Maisel sits far to the right (highest μ) and high up (large σ): adored, but divisive. Priceless and Miss Sloane are the other high-σ standouts — strong fan blocs pulling μ up despite a substantial low-star tail. Doctor Thorne, Room and Taking Chance sit low on the σ axis: consistent, agreed-upon films. Two movies at the same horizontal position can be far apart vertically — that vertical gap is exactly what the mean rating throws away.
Why the movers move (second plot). Here μ is plotted against the naive metric mean, with point size and colour showing σ, and the dashed line the average metric→μ relationship. The high-σ films (large, bright points) sit above the line. Miss Sloane has the lowest metric mean of the labelled set (4.07) yet a latent μ of 6.4 — the model reads its five-star bloc as evidence of high quality-for-fans that the average buries under the low-star tail. Low-σ films like Taking Chance and Room sit on or below the line: their tight rating distribution gives the model no reason to infer a quality above what the mean already shows. This is the picture behind the notebook's pattern — the metric and the ordinal model disagree precisely where σ departs from average.
Takeaway. The heteroscedastic model does not produce a better one-number ranking; it produces a two-dimensional description — how good a film is for the people who like it, and how much people disagree — and lets a viewer's risk appetite decide. The mean rating is the projection of this plane onto a single axis, and these plots show what that projection loses.
6. Comparison with R polr Results¶
Load movies_comparison_polr.csv saved by ord_probit_r.ipynb. Compare:
polrprobit $\hat{\beta}_g$ (homoscedastic MLE) vs A&C Gibbs $\hat{\mu}_g$ (heteroscedastic Bayes)- Rank agreement and notable discrepancies
import os
if os.path.exists('movies_comparison_polr.csv'):
polr = pd.read_csv('movies_comparison_polr.csv')
merged = full.merge(polr[['movie','beta_probit','r_probit']],
on='movie', how='left')
merged['shift_vs_polr'] = merged['rank_ordinal'] - merged['r_probit']
print('Comparison: A&C Gibbs mu vs polr probit beta (both probit scale):')
print(f'{"Movie":<50} {"A&C_mu":>8} {"polr_b":>8} {"r_AC":>6} {"r_pol":>6} {"diff":>6}')
print('-' * 90)
for _, row in merged.sort_values('rank_ordinal').iterrows():
print(f'{row.movie:<50} {row.latent_mu:>8.3f} {row.beta_probit:>8.4f}'
f' {row.rank_ordinal:>6} {row.r_probit:>6} {row.shift_vs_polr:>6.0f}')
print(f'\nMax |A&C_rank - polr_rank|: {merged["shift_vs_polr"].abs().max():.0f}')
print('Large shifts indicate movies where polarisation affects ranking.')
else:
print('Run ord_probit_r.ipynb first to generate movies_comparison_polr.csv')
Comparison: A&C Gibbs mu vs polr probit beta (both probit scale): Movie A&C_mu polr_b r_AC r_pol diff ------------------------------------------------------------------------------------------ The Marvelous Mrs Maisel Season 1 11.519 1.6319 1 1 0 Bosch Season 3 8.698 1.3777 2 2 0 Poldark Season 2 8.363 1.3045 3 4 -1 Priceless 7.670 0.8112 4 9 -5 Megan Leavey 6.646 1.2206 5 5 0 Freedom Writers 6.489 1.0684 6 6 0 Miss Sloane 6.402 0.4259 7 22 -15 Taking Chance 6.346 1.3567 8 3 5 Me Before You 6.175 0.6285 9 16 -7 The Curious Case of Benjamin Button 5.937 0.6919 10 13 -3 The Choice 5.741 0.5928 11 17 -6 Sneaky Pete Season 1 5.661 0.9935 12 7 5 Lifted 5.493 0.7397 13 11 2 The Man in the High Castle Season 1 5.488 0.7559 14 10 4 The Only Living Boy in New York 5.469 0.5545 15 18 -3 Ladies in Lavender 5.380 0.6876 16 14 2 The Guardian Season 1 5.281 0.8727 17 8 9 The Inheritance 5.203 0.4999 18 20 -2 Room 5.164 0.7249 19 12 7 Julian Fellowes Presents Doctor Thorne Season 1 5.039 0.6395 20 15 5 Philip K. Dick's Electric Dreams - Season 1 5.017 0.1466 21 26 -5 Fortitude Season 2 4.895 0.2022 22 25 -3 Hunted Season 1 4.796 0.4716 23 21 2 Creed 4.781 0.5163 24 19 5 The Infiltrator 4.478 0.3135 25 23 2 Allied 4.437 0.2108 26 24 2 Brawl in Cell Block 99 4.321 -0.0191 27 29 -2 Rise 4.139 0.0653 28 27 1 Revolutionary Road 4.028 -0.0414 29 30 -1 The Whole Truth 3.960 0.0000 30 28 2 Triple 9 3.856 -0.0896 31 31 0 Silence 3.803 -0.1802 32 33 -1 The Sea of Trees 3.800 -0.1142 33 32 1 Tin Star Season 1 3.709 -0.2504 34 35 -1 In Secret 3.590 -0.2182 35 34 1 Babel 3.298 -0.3761 36 36 0 Max |A&C_rank - polr_rank|: 15 Large shifts indicate movies where polarisation affects ranking.
polr vs A&C comparison (validated)¶
Max |A&C_rank − polr_rank|: 15
Largest disagreements (A&C rank vs polr rank):
Miss Sloane A&C= 7 polr=22 diff=−15 mu=6.402 sg=4.659 β=0.426
The Guardian S1 A&C=17 polr= 8 diff= +9 mu=5.281 sg=1.378 β=0.873
Me Before You A&C= 9 polr=16 diff= −7 mu=6.175 sg=3.283 β=0.629
Room A&C=19 polr=12 diff= +7 mu=5.164 sg=1.537 β=0.725
The Choice A&C=11 polr=17 diff= −6 mu=5.741 sg=2.782 β=0.593
Taking Chance A&C= 8 polr= 3 diff= +5 mu=6.346 sg=1.692 β=1.357
Sneaky Pete S1 A&C=12 polr= 7 diff= +5 mu=5.661 sg=1.618 β=0.994
Priceless A&C= 4 polr= 9 diff= −5 mu=7.670 sg=4.454 β=0.811
Philip K Dick S1 A&C=21 polr=26 diff= −5 mu=5.017 sg=3.775 β=0.147
Doctor Thorne S1 A&C=20 polr=15 diff= +5 mu=5.039 sg=1.560 β=0.640
Pattern: the two models disagree precisely where $\sigma_g$ differs from its average.
A&C upgrades (polr too low): high-$\sigma$ films with a strong five-star bloc — Miss Sloane ($\sigma=4.66$, diff=−15), Me Before You ($\sigma=3.28$, diff=−7), The Choice ($\sigma=2.78$, diff=−6), Priceless ($\sigma=4.45$, diff=−5). polr's $\hat\beta$ is dragged down by their below-average ratings; A&C's $\hat\mu$ sees the five-star bloc and infers high latent quality.
polr upgrades (A&C too low): low-$\sigma$ films with high empirical means — The Guardian S1 ($\sigma=1.38$, diff=+9), Room ($\sigma=1.54$, diff=+7), Taking Chance ($\sigma=1.69$, diff=+5), Sneaky Pete S1 ($\sigma=1.62$, diff=+5). polr sees a strong cumulative signal at every threshold and gives a high $\hat\beta$; A&C infers a more modest $\hat\mu$ because a tight distribution doesn't require a very high location to fit the top-threshold probability.
Miss Sloane is the extreme case. polr $\hat\beta=0.426$ (rank 22, looks mediocre). A&C $\hat\mu=6.402$, $\hat\sigma=4.659$ (rank 7, genuinely excellent for fans). The gap of 15 ranks is entirely explained by $\sigma$: polr cannot distinguish a bimodal film from a mediocre one.
Comparison results (validated)¶
Metric rank from emp_mean descending; ordinal rank from posterior mu descending.
Largest disagreements (|shift| >= 5):
Miss Sloane 1201 4.071 mu= 6.402 sg=4.659 metric=23 ordinal= 7 shift=+16
Priceless 745 4.412 mu= 7.670 sg=4.454 metric=14 ordinal= 4 shift=+10
Room 13897 4.490 mu= 5.164 sg=1.537 metric= 9 ordinal=19 shift=-10
The Guardian S1 2081 4.599 mu= 5.281 sg=1.378 metric= 8 ordinal=17 shift= -9
Creed 11993 4.330 mu= 4.781 sg=1.450 metric=16 ordinal=24 shift= -8
Me Before You 4966 4.310 mu= 6.175 sg=3.283 metric=17 ordinal= 9 shift= +8
Doctor Thorne S1 21079 4.420 mu= 5.039 sg=1.560 metric=13 ordinal=20 shift= -7
The Choice 3546 4.290 mu= 5.741 sg=2.782 metric=18 ordinal=11 shift= +7
Taking Chance 5793 4.810 mu= 6.346 sg=1.692 metric= 2 ordinal= 8 shift= -6
Philip K Dick Elec Dr 822 3.779 mu= 5.017 sg=3.775 metric=27 ordinal=21 shift= +6
Only Living Boy in NY 587 4.262 mu= 5.469 sg=2.532 metric=20 ordinal=15 shift= +5
Benjamin Button 1296 4.379 mu= 5.937 sg=2.726 metric=15 ordinal=10 shift= +5
Sneaky Pete S1 57281 4.650 mu= 5.661 sg=1.618 metric= 7 ordinal=12 shift= -5
Most polarising (highest sigma):
Mrs Maisel S1 sg=4.771 ordinal= 1 metric= 1 mu=11.519
Miss Sloane sg=4.659 ordinal= 7 metric=23 mu= 6.402
Priceless sg=4.454 ordinal= 4 metric=14 mu= 7.670
Philip K Dick Elec Dr sg=3.775 ordinal=21 metric=27 mu= 5.017
Bosch S3 sg=3.500 ordinal= 2 metric= 3 mu= 8.698
Poldark S2 sg=3.427 ordinal= 3 metric= 4 mu= 8.363
Brawl in Cell Block 99 sg=3.328 ordinal=27 metric=31 mu= 4.321
Me Before You sg=3.283 ordinal= 9 metric=17 mu= 6.175
L&K headline pair — three-model progression:
Film metric polr Gibbs mu sigma
Priceless 14 9 4 7.670 4.454 ← high mu AND high sigma
Doctor Thorne S1 13 15 20 5.039 1.560 ← moderate mu, very consistent
Analysis¶
What the three models show for the headline pair:
| Model | Priceless | Doctor Thorne | Direction |
|---|---|---|---|
| Metric mean | 4.412 (rank 14) | 4.420 (rank 13) | Tie — indistinguishable |
| polr MLE | $\hat\beta = 0.811$ (rank 9) | $\hat\beta = 0.640$ (rank 15) | Priceless > DT |
| A&C Gibbs | $\mu=7.670,\,\sigma=4.454$ (rank 4) | $\mu=5.039,\,\sigma=1.560$ (rank 20) | Priceless >> DT in $\mu$; DT << Priceless in $\sigma$ |
All three models rank Priceless above Doctor Thorne in location. The heteroscedastic model does not reverse this ordering — it amplifies it. But it adds the crucial new dimension: why.
Why the metric and homoscedastic models mislead. Priceless has 77% five-star ratings, which the metric model and polr interpret as "high quality." But the heteroscedastic model decomposes this into $\mu_g = 7.67$ (latent quality ceiling for fans) and $\sigma_g = 4.45$ (wide spread — fans and detractors are far apart on the latent scale). Doctor Thorne's 63% five-star share and only 5% low-star maps to a moderate but reliable quality signal: $\mu_g = 5.04$, $\sigma_g = 1.56$. A risk-averse viewer should prefer Doctor Thorne; a high-variance viewer might prefer Priceless. The single-dimensional rankings cannot encode this choice.
Miss Sloane — the real discovery. The largest shift in the dataset is Miss Sloane: metric rank 23 → ordinal rank 7 (shift = +16). Its empirical mean (4.071) looks mediocre, but the model infers $\mu_g = 6.40$, $\sigma_g = 4.66$. The distribution is inferred to be strongly bimodal — a large bloc of five-star fans and a substantial low-star tail — producing a low average despite a high latent quality peak. Miss Sloane is the most egregious case of the metric mean underrepresenting a genuinely excellent (for fans) but divisive film.
Consistent films, downranked by the ordinal model. The model systematically lowers the ranking of films whose ratings concentrate around 4–5 stars without a dominant 5-star spike. These films have low $\sigma_g$, so their cumulative distribution is tight, and the model infers a more conservative $\mu_g$ than their mean suggests:
| Film | emp | $\mu_g$ | $\sigma_g$ | metric → ordinal |
|---|---|---|---|---|
| The Guardian S1 | 4.599 | 5.281 | 1.378 | 8 → 17 (−9) |
| Room | 4.490 | 5.164 | 1.537 | 9 → 19 (−10) |
| Taking Chance | 4.810 | 6.346 | 1.692 | 2 → 8 (−6) |
| Creed | 4.330 | 4.781 | 1.450 | 16 → 24 (−8) |
Room and Guardian are "solid consensus films" — nearly everyone rates them highly, but without a dominant five-star spike the model infers moderate $\mu_g$. Taking Chance is downranked less severely (still rank 8) because its empirical mean (4.810) is very high, reflecting a large concentration of 5-star ratings even with low $\sigma_g$.
The Mrs Maisel anomaly. Mrs Maisel ranks #1 by every model. But the Gibbs result ($\mu = 11.52$, $\sigma = 4.77$) is striking: it has the highest mu AND the highest sigma. The model infers that its ~93% five-star concentration coexists with a non-trivial low-star tail, requiring both extreme location and extreme scale to fit simultaneously. It is the quintessential "phenomenon" film: beloved by fans, incomprehensible to detractors.
Summary of the three-model comparison:
| Finding | Metric | polr | A&C Gibbs |
|---|---|---|---|
| Ranking dimension | 1 (mean) | 1 ($\hat\beta$) | 2 ($\hat\mu_g$, $\hat\sigma_g$) |
| Priceless / Doctor Thorne | Tied | Priceless > DT | Priceless >> DT in $\mu$; DT << in $\sigma$ |
| Miss Sloane discovery | Rank 23 | Moderate shift | Rank 7 (hidden quality exposed) |
| Consistent films (low $\sigma$) | Rewarded by mean | Slight reward | Correctly moderated |
| Polarising films (high $\sigma$) | Confused with high quality | Partially detected | Separated into $\mu$ + $\sigma$ |
The heteroscedastic model does not give a single "correct" ranking — it gives the audience a two-dimensional quality description and lets risk preferences determine the ranking. This is L&K's core point: ordinal data with per-group scale is not a ranking problem, it is a characterisation problem.
7. Save Results¶
full.to_csv('movies_comparison_python.csv', index=False)
print('Saved: movies_comparison_python.csv')
mu_df = pd.DataFrame(fit_m['mu'], columns=mnames)
sg_df = pd.DataFrame(fit_m['sigma'], columns=mnames)
mu_df.to_csv('movies_mu_draws_python.csv', index=False)
sg_df.to_csv('movies_sigma_draws_python.csv', index=False)
print('Saved: movies_mu_draws_python.csv movies_sigma_draws_python.csv')
print(f'Posterior draws shape: {fit_m["mu"].shape} (R_kept x G)')
Saved: movies_comparison_python.csv Saved: movies_mu_draws_python.csv movies_sigma_draws_python.csv Posterior draws shape: (4000, 36) (R_kept x G)