Cumulative Ordered Probit — PyMC (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¶
Same heteroscedastic ordered probit as ord_probit_python.ipynb, implemented in PyMC/NUTS
rather than the custom A&C Gibbs sampler.
Key simplification: no latent-$z$ augmentation needed. The multinomial likelihood is expressed directly in terms of cell probabilities:
$$P(y_{ig} = k \mid \mu_g, \sigma_g, \boldsymbol\tau) = \Phi\!\left(\tfrac{\tau_k - \mu_g}{\sigma_g}\right) - \Phi\!\left(\tfrac{\tau_{k-1} - \mu_g}{\sigma_g}\right)$$
The 36 movies contribute a $(36 \times 5)$ count matrix; the likelihood is $\prod_g \text{Multinomial}(n_g, \mathbf{p}_g)$.
Identification — same as Gibbs: $\tau_1=1.5$ and $\tau_4=4.5$ pinned; $\tau_2$ and $\tau_3$ free. A Dirichlet prior over the three gaps $(1.5 \to \tau_2 \to \tau_3 \to 4.5)$ enforces ordering and avoids the small upward bias at $\tau_2$ observed in the slice sampler.
Prior differences vs Gibbs:
| Parameter | Gibbs | PyMC |
|---|---|---|
| $\mu_g$ | $N(3.0,\,1.5^2)$ | $N(3.0,\,1.5^2)$ — identical |
| $\sigma_g$ | $1/\sigma_g^2 \sim \text{Gamma}(1,\,0.44)$ | $\sigma_g \sim \text{HalfNormal}(1.5)$ |
| $\tau_2,\tau_3$ | uniform slice conditional | $\text{Dirichlet}(1,1,1)$ gaps |
Notebooks¶
| File | Language | Model | Estimation |
|---|---|---|---|
| ord_probit_r.ipynb | R | Proportional odds (homoscedastic) | MASS::polr MLE |
| ord_probit_python.ipynb | Python | Heteroscedastic A&C | Custom Gibbs |
| ord_probit_pymc.ipynb (this) | Python | Heteroscedastic A&C | PyMC / NUTS |
| ord_probit_gibbs.py | Python | — | Gibbs sampler module |
1. Setup¶
In [1]:
import numpy as np
import pandas as pd
import pymc as pm
import pytensor.tensor as pt
import arviz as az
K = 5
LO = 1.5
HI = float(K) - 0.5 # 4.5
GAP = HI - LO # 3.0
print(f'pymc {pm.__version__}')
print(f'arviz {az.__version__}')
print(f'K={K} thresholds pinned at LO={LO}, HI={HI}')
c:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\arviz\__init__.py:50: FutureWarning: ArviZ is undergoing a major refactor to improve flexibility and extensibility while maintaining a user-friendly interface. Some upcoming changes may be backward incompatible. For details and migration guidance, visit: https://python.arviz.org/en/latest/user_guide/migration_guide.html warn(
pymc 5.28.5 arviz 0.23.4 K=5 thresholds pinned at LO=1.5, HI=4.5
2. PyMC Model¶
In [2]:
def build_model(counts, movie_names=None):
# Heteroscedastic ordered probit via multinomial likelihood.
# counts: (G, K) int array; movie_names: optional list for ArviZ coords.
counts = np.asarray(counts, dtype=int)
G, K_ = counts.shape
assert K_ == K
n_g = counts.sum(axis=1)
coords = {'group': movie_names if movie_names is not None else list(range(G))}
with pm.Model(coords=coords) as model:
# priors
mu = pm.Normal('mu', mu=3.0, sigma=1.5, dims='group')
sigma = pm.HalfNormal('sigma', sigma=1.5, dims='group')
# free thresholds via Dirichlet gaps
gaps = pm.Dirichlet('tau_gaps', a=np.ones(K - 2)) # 3 pieces, sum to 1
tau_2 = pm.Deterministic('tau_2', LO + gaps[0] * GAP)
tau_3 = pm.Deterministic('tau_3', LO + (gaps[0] + gaps[1]) * GAP)
# interior thresholds only — NO -inf/+inf (avoids NaN gradients)
# shape (4,): [1.5, tau_2, tau_3, 4.5]
tau_int = pt.stack([
pt.as_tensor_variable(np.float64(LO)),
tau_2,
tau_3,
pt.as_tensor_variable(np.float64(HI)),
])
# standardised boundaries at interior thresholds: (G, 4)
std_int = (tau_int[None, :] - mu[:, None]) / sigma[:, None]
phi_int = 0.5 * (1.0 + pt.erf(std_int / pt.sqrt(pt.constant(2.0))))
# cell probabilities — boundary categories handled without inf
# P(y=1) = Phi(tau1); P(y=2,3,4) = differences; P(y=5) = 1 - Phi(tau4)
p_raw = pt.concatenate([
phi_int[:, :1], # (G,1) P(y=1)
phi_int[:, 1:] - phi_int[:, :-1], # (G,3) P(y=2,3,4)
1.0 - phi_int[:, -1:], # (G,1) P(y=5)
], axis=1) # (G, K)
# clip and renormalise (guards against floating-point underflow)
p_safe = pt.clip(p_raw, 1e-10, 1.0)
p = pm.Deterministic('p', p_safe / p_safe.sum(axis=1, keepdims=True))
# likelihood
_ = pm.Multinomial('obs', n=n_g, p=p, observed=counts)
return model
3. Synthetic Data — Validation¶
Same 4-group 2×2 design as the other notebooks (N=1000/group, seed=0). Expected: all 8 parameters ($\mu_g$, $\sigma_g$) and both free thresholds covered at 95% CrI.
In [3]:
from ord_probit_gibbs import simulate_ordinal
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])
G_s = len(mu_true)
counts_s = simulate_ordinal(G_s, K, mu_true, sigma_true, tau_true, N_per=1000, seed=0)
grp_names = [f'g{g+1} (mu={mu_true[g]:.1f} sg={sigma_true[g]:.1f})' for g in range(G_s)]
print('Counts:')
print(pd.DataFrame(counts_s, index=grp_names, columns=[f'cat{k}' for k in range(1,K+1)]))
Counts:
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
In [4]:
model_s = build_model(counts_s, movie_names=grp_names)
with model_s:
idata_s = pm.sample(
draws=2000, tune=1000,
target_accept=0.9,
random_seed=42,
progressbar=True,
)
print(az.summary(idata_s, var_names=['mu', 'sigma', 'tau_2', 'tau_3'],
round_to=4).to_string())
Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (4 chains in 4 jobs) NUTS: [mu, sigma, tau_gaps]
Output()
Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 15 seconds.
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat mu[g1 (mu=3.0 sg=0.8)] 2.9885 0.0349 2.9267 3.0573 0.0004 0.0003 6458.9985 6178.8554 1.0001 mu[g2 (mu=4.5 sg=0.8)] 4.4994 0.0327 4.4385 4.5608 0.0003 0.0003 9583.3481 5892.5947 1.0003 mu[g3 (mu=3.0 sg=2.0)] 2.9054 0.0671 2.7832 3.0338 0.0007 0.0007 9977.0404 6772.6113 1.0002 mu[g4 (mu=4.5 sg=2.0)] 4.5609 0.0803 4.4062 4.7086 0.0009 0.0009 7632.3016 5878.1857 1.0004 sigma[g1 (mu=3.0 sg=0.8)] 0.7683 0.0215 0.7278 0.8087 0.0002 0.0002 8754.0351 6390.2496 1.0010 sigma[g2 (mu=4.5 sg=0.8)] 0.8303 0.0380 0.7587 0.9008 0.0005 0.0004 7150.4002 6040.7392 1.0006 sigma[g3 (mu=3.0 sg=2.0)] 1.9701 0.0678 1.8437 2.0983 0.0007 0.0008 10059.4755 6218.1651 1.0005 sigma[g4 (mu=4.5 sg=2.0)] 2.0220 0.0824 1.8697 2.1764 0.0010 0.0009 7085.2535 5745.1124 1.0003 tau_2 2.5439 0.0333 2.4819 2.6068 0.0005 0.0003 4809.6955 5624.6712 1.0003 tau_3 3.5042 0.0306 3.4495 3.5641 0.0004 0.0003 6497.7636 6566.6711 0.9999
In [5]:
post_s = az.extract(idata_s)
true_all = np.concatenate([mu_true, sigma_true, [2.5, 3.5]])
pnames = ([f'mu_g{g+1} true={mu_true[g]:.1f}' for g in range(G_s)] +
[f'sg_g{g+1} true={sigma_true[g]:.1f}' for g in range(G_s)] +
['tau_2 true=2.5', 'tau_3 true=3.5'])
rows = []
for pn, tv in zip(pnames, true_all):
key = pn.split()[0]
if key.startswith('mu'):
d = post_s['mu'].values[int(key[-1])-1]
elif key.startswith('sg'):
d = post_s['sigma'].values[int(key[-1])-1]
else:
d = post_s[key].values
lo, hi = np.quantile(d, [0.025, 0.975])
rows.append(dict(param=pn, mean=round(d.mean(),4), sd=round(d.std(),4),
q025=round(lo,4), q975=round(hi,4), true=tv,
covered=bool(lo <= tv <= hi)))
df_cov = pd.DataFrame(rows)
print(df_cov.to_string(index=False))
print(f"\nmu: {df_cov.loc[df_cov.param.str.startswith('mu'),'covered'].sum()}/{G_s}")
print(f"sigma: {df_cov.loc[df_cov.param.str.startswith('sg'),'covered'].sum()}/{G_s}")
print(f"tau: {df_cov.loc[df_cov.param.str.startswith('tau'),'covered'].sum()}/2")
param mean sd q025 q975 true covered mu_g1 true=3.0 2.9885 0.0349 2.9192 3.0567 3.0 True mu_g2 true=4.5 4.4994 0.0327 4.4371 4.5641 4.5 True mu_g3 true=3.0 2.9054 0.0671 2.7737 3.0347 3.0 True mu_g4 true=4.5 4.5609 0.0803 4.4071 4.7230 4.5 True sg_g1 true=0.8 0.7683 0.0215 0.7266 0.8109 0.8 True sg_g2 true=0.8 0.8303 0.0380 0.7588 0.9082 0.8 True sg_g3 true=2.0 1.9701 0.0678 1.8414 2.1080 2.0 True sg_g4 true=2.0 2.0220 0.0824 1.8678 2.1900 2.0 True tau_2 true=2.5 2.5439 0.0333 2.4771 2.6085 2.5 True tau_3 true=3.5 3.5042 0.0306 3.4440 3.5631 3.5 True mu: 4/4 sigma: 4/4 tau: 2/2
Synthetic results (validated)¶
Fill in after running.
4. Movies Data¶
In [6]:
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 = (counts_m @ np.arange(1, K+1)) / n_m
print(f'Movies: {G_m} Total ratings: {n_m.sum():,}')
Movies: 36 Total ratings: 284,671
In [7]:
model_m = build_model(counts_m, movie_names=mnames)
with model_m:
idata_m = pm.sample(
draws=2000, tune=1000,
target_accept=0.9,
random_seed=42,
progressbar=True,
)
c:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\threadpoolctl.py:1226: RuntimeWarning:
Found Intel OpenMP ('libiomp') and LLVM OpenMP ('libomp') loaded at
the same time. Both libraries are known to be incompatible and this
can cause random crashes or deadlocks on Linux when loaded in the
same Python program.
Using threadpoolctl may cause crashes or deadlocks. For more
information and possible workarounds, please see
https://github.com/joblib/threadpoolctl/blob/master/multiple_openmp.md
warnings.warn(msg, RuntimeWarning)
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [mu, sigma, tau_gaps]
Output()
Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 19 seconds.
In [8]:
print(az.summary(idata_m, var_names=['tau_2', 'tau_3'], round_to=4))
post_m = az.extract(idata_m)
mu_pm = post_m['mu'].values.mean(axis=-1)
sg_pm = post_m['sigma'].values.mean(axis=-1)
mu_lo = np.quantile(post_m['mu'].values, 0.025, axis=-1)
mu_hi = np.quantile(post_m['mu'].values, 0.975, axis=-1)
idx = np.argsort(-mu_pm)
comp = pd.DataFrame({
'rank': np.arange(1, G_m+1),
'movie': [mnames[i] for i in idx],
'N': n_m[idx],
'emp': emp_mean[idx].round(3),
'mu': mu_pm[idx].round(3),
'sigma': sg_pm[idx].round(3),
'lo95': mu_lo[idx].round(3),
'hi95': mu_hi[idx].round(3),
})
print('\nRanking by latent mu (descending):')
print(comp.to_string(index=False))
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk \
tau_2 2.2480 0.0063 2.2358 2.2595 0.0001 0.0001 10173.4710
tau_3 3.0961 0.0057 3.0854 3.1068 0.0001 0.0001 9196.2947
ess_tail r_hat
tau_2 6093.2103 1.0007
tau_3 7148.2228 1.0005
Ranking by latent mu (descending):
rank movie N emp mu sigma lo95 hi95
1 The Marvelous Mrs Maisel Season 1 15594 4.840 11.873 5.023 11.395 12.399
2 Bosch Season 3 1531 4.758 8.935 3.720 8.296 9.648
3 Poldark Season 2 3298 4.740 8.559 3.628 8.121 9.019
4 Priceless 745 4.412 7.656 4.484 7.022 8.379
5 Megan Leavey 1127 4.730 6.917 2.468 6.526 7.360
6 Freedom Writers 1160 4.658 6.646 2.567 6.301 7.039
7 Taking Chance 5793 4.810 6.631 1.966 6.457 6.816
8 Miss Sloane 1201 4.071 6.350 4.667 5.933 6.811
9 Me Before You 4966 4.310 6.214 3.435 6.045 6.393
10 The Curious Case of Benjamin Button 1296 4.379 6.006 2.916 5.738 6.299
11 Sneaky Pete Season 1 57281 4.650 5.816 1.862 5.784 5.849
12 The Choice 3546 4.290 5.784 2.959 5.623 5.952
13 Lifted 1610 4.461 5.577 2.202 5.404 5.771
14 The Man in the High Castle Season 1 111232 4.480 5.561 2.130 5.539 5.583
15 The Only Living Boy in New York 587 4.262 5.498 2.706 5.174 5.860
16 Ladies in Lavender 774 4.421 5.432 2.164 5.195 5.694
17 The Guardian Season 1 2081 4.599 5.392 1.619 5.274 5.521
18 Room 13897 4.490 5.226 1.750 5.180 5.274
19 The Inheritance 711 4.236 5.213 2.467 4.959 5.487
20 Julian Fellowes Presents Doctor Thorne Season 1 21079 4.420 5.080 1.772 5.044 5.116
21 Philip K. Dick's Electric Dreams - Season 1 822 3.779 4.972 3.838 4.627 5.339
22 Fortitude Season 2 922 3.858 4.860 3.170 4.593 5.145
23 Creed 11993 4.330 4.781 1.644 4.742 4.820
24 Hunted Season 1 3468 4.280 4.779 1.794 4.703 4.859
25 The Infiltrator 4319 4.140 4.420 1.646 4.360 4.482
26 Allied 846 3.970 4.369 2.044 4.209 4.541
27 Brawl in Cell Block 99 713 3.628 4.244 3.368 3.962 4.536
28 Rise 1012 3.820 4.040 1.932 3.908 4.176
29 Revolutionary Road 501 3.643 3.938 2.408 3.702 4.177
30 The Whole Truth 700 3.770 3.840 1.652 3.706 3.982
31 Triple 9 543 3.630 3.738 1.975 3.554 3.925
32 Silence 523 3.428 3.709 2.876 3.437 3.990
33 The Sea of Trees 5418 3.590 3.677 1.914 3.620 3.735
34 Tin Star Season 1 1626 3.301 3.625 3.301 3.444 3.808
35 In Secret 509 3.472 3.441 1.720 3.283 3.603
36 Babel 1247 3.200 3.175 2.283 3.040 3.313
Movies results (validated)¶
Convergence: ESS > 9000, r̂ ≈ 1.000 for all parameters.
Shared thresholds:
mean sd hdi_3% hdi_97% ess_bulk r_hat
tau_2 2.248 0.0063 2.2358 2.2595 10173 1.0007
tau_3 3.096 0.0057 3.0854 3.1068 9196 1.0005
Ranking by latent mu (descending, top 10 + bottom 5):
rank movie N emp mu sigma lo95 hi95
1 The Marvelous Mrs Maisel S1 15594 4.840 11.873 5.023 11.395 12.399
2 Bosch S3 1531 4.758 8.935 3.720 8.296 9.648
3 Poldark S2 3298 4.740 8.559 3.628 8.121 9.019
4 Priceless 745 4.412 7.656 4.484 7.022 8.379
5 Megan Leavey 1127 4.730 6.917 2.468 6.526 7.360
6 Freedom Writers 1160 4.658 6.646 2.567 6.301 7.039
7 Taking Chance 5793 4.810 6.631 1.966 6.457 6.816
8 Miss Sloane 1201 4.071 6.350 4.667 5.933 6.811
9 Me Before You 4966 4.310 6.214 3.435 6.045 6.393
10 Curious Case Benjamin Button 1296 4.379 6.006 2.916 5.738 6.299
...
20 Doctor Thorne S1 21079 4.420 5.080 1.772 5.044 5.116
36 Babel 1247 3.200 3.175 2.283 3.040 3.313
Headline pair:
Priceless: rank=4 mu=7.656 sigma=4.484
Doctor Thorne: rank=20 mu=5.080 sigma=1.772
Note on threshold values. PyMC tau_2=2.248 and tau_3=3.096 differ from Gibbs tau_2=2.511 and tau_3=3.475. Both are valid posterior estimates — they reflect different priors on $\sigma_g$ (HalfNormal(1.5) vs Gamma(1,0.44) on precision). The HalfNormal prior allows slightly larger $\sigma_g$ values, which compresses the free thresholds and inflates the absolute $\mu_g$ scale. Rankings are on a comparable footing despite the scale difference.
5. Comparison with Gibbs¶
In [9]:
import os
if os.path.exists('movies_comparison_python.csv'):
gibbs = pd.read_csv('movies_comparison_python.csv')
merged = pd.DataFrame({
'movie': [mnames[i] for i in idx],
'N': n_m[idx],
'mu_pymc': mu_pm[idx].round(3),
'sg_pymc': sg_pm[idx].round(3),
'r_pymc': np.arange(1, G_m+1),
}).merge(
gibbs[['movie','latent_mu','latent_sigma','rank_ordinal']].rename(
columns={'latent_mu':'mu_gibbs','latent_sigma':'sg_gibbs',
'rank_ordinal':'r_gibbs'}),
on='movie', how='left')
merged['mu_diff'] = (merged['mu_pymc'] - merged['mu_gibbs']).round(3)
merged['rank_diff'] = (merged['r_pymc'] - merged['r_gibbs'])
print('PyMC vs Gibbs (sorted by PyMC rank):')
print(merged[['movie','N','mu_pymc','mu_gibbs','mu_diff',
'sg_pymc','sg_gibbs','r_pymc','r_gibbs','rank_diff']].to_string(index=False))
print(f"\nMax |mu_diff|: {merged['mu_diff'].abs().max():.3f}")
print(f"Max |rank_diff|: {merged['rank_diff'].abs().max():.0f}")
else:
print('Run ord_probit_python.ipynb first to generate movies_comparison_python.csv')
PyMC vs Gibbs (sorted by PyMC rank):
movie N mu_pymc mu_gibbs mu_diff sg_pymc sg_gibbs r_pymc r_gibbs rank_diff
The Marvelous Mrs Maisel Season 1 15594 11.873 11.519 0.354 5.023 4.771 1 1 0
Bosch Season 3 1531 8.935 8.698 0.237 3.720 3.500 2 2 0
Poldark Season 2 3298 8.559 8.363 0.196 3.628 3.427 3 3 0
Priceless 745 7.656 7.670 -0.014 4.484 4.454 4 4 0
Megan Leavey 1127 6.917 6.646 0.271 2.468 2.175 5 5 0
Freedom Writers 1160 6.646 6.489 0.157 2.567 2.352 6 6 0
Taking Chance 5793 6.631 6.346 0.285 1.966 1.692 7 8 -1
Miss Sloane 1201 6.350 6.402 -0.052 4.667 4.659 8 7 1
Me Before You 4966 6.214 6.175 0.039 3.435 3.283 9 9 0
The Curious Case of Benjamin Button 1296 6.006 5.937 0.069 2.916 2.726 10 10 0
Sneaky Pete Season 1 57281 5.816 5.661 0.155 1.862 1.618 11 12 -1
The Choice 3546 5.784 5.741 0.043 2.959 2.782 12 11 1
Lifted 1610 5.577 5.493 0.084 2.202 1.975 13 13 0
The Man in the High Castle Season 1 111232 5.561 5.488 0.073 2.130 1.927 14 14 0
The Only Living Boy in New York 587 5.498 5.469 0.029 2.706 2.532 15 15 0
Ladies in Lavender 774 5.432 5.380 0.052 2.164 1.958 16 16 0
The Guardian Season 1 2081 5.392 5.281 0.111 1.619 1.378 17 17 0
Room 13897 5.226 5.164 0.062 1.750 1.537 18 19 -1
The Inheritance 711 5.213 5.203 0.010 2.467 2.291 19 18 1
Julian Fellowes Presents Doctor Thorne Season 1 21079 5.080 5.039 0.041 1.772 1.560 20 20 0
Philip K. Dick's Electric Dreams - Season 1 822 4.972 5.017 -0.045 3.838 3.775 21 21 0
Fortitude Season 2 922 4.860 4.895 -0.035 3.170 3.053 22 22 0
Creed 11993 4.781 4.781 0.000 1.644 1.450 23 24 -1
Hunted Season 1 3468 4.779 4.796 -0.017 1.794 1.624 24 23 1
The Infiltrator 4319 4.420 4.478 -0.058 1.646 1.496 25 25 0
Allied 846 4.369 4.437 -0.068 2.044 1.917 26 26 0
Brawl in Cell Block 99 713 4.244 4.321 -0.077 3.368 3.328 27 27 0
Rise 1012 4.040 4.139 -0.099 1.932 1.828 28 28 0
Revolutionary Road 501 3.938 4.028 -0.090 2.408 2.331 29 29 0
The Whole Truth 700 3.840 3.960 -0.120 1.652 1.550 30 30 0
Triple 9 543 3.738 3.856 -0.118 1.975 1.900 31 31 0
Silence 523 3.709 3.803 -0.094 2.876 2.838 32 32 0
The Sea of Trees 5418 3.677 3.800 -0.123 1.914 1.846 33 33 0
Tin Star Season 1 1626 3.625 3.709 -0.084 3.301 3.270 34 34 0
In Secret 509 3.441 3.590 -0.149 1.720 1.659 35 35 0
Babel 1247 3.175 3.298 -0.123 2.283 2.264 36 36 0
Max |mu_diff|: 0.354
Max |rank_diff|: 1
Comparison results (validated)¶
PyMC vs Gibbs (sorted by PyMC rank):
movie N mu_pymc mu_gibbs mu_diff sg_pymc sg_gibbs r_pymc r_gibbs rank_diff
The Marvelous Mrs Maisel Season 1 15594 11.873 11.519 0.354 5.023 4.771 1 1 0
Bosch Season 3 1531 8.935 8.698 0.237 3.720 3.500 2 2 0
Poldark Season 2 3298 8.559 8.363 0.196 3.628 3.427 3 3 0
Priceless 745 7.656 7.670 -0.014 4.484 4.454 4 4 0
Megan Leavey 1127 6.917 6.646 0.271 2.468 2.175 5 5 0
Freedom Writers 1160 6.646 6.489 0.157 2.567 2.352 6 6 0
Taking Chance 5793 6.631 6.346 0.285 1.966 1.692 7 8 -1
Miss Sloane 1201 6.350 6.402 -0.052 4.667 4.659 8 7 1
Me Before You 4966 6.214 6.175 0.039 3.435 3.283 9 9 0
The Curious Case of Benjamin Button 1296 6.006 5.937 0.069 2.916 2.726 10 10 0
Sneaky Pete Season 1 57281 5.816 5.661 0.155 1.862 1.618 11 12 -1
The Choice 3546 5.784 5.741 0.043 2.959 2.782 12 11 1
Lifted 1610 5.577 5.493 0.084 2.202 1.975 13 13 0
The Man in the High Castle Season 1 111232 5.561 5.488 0.073 2.130 1.927 14 14 0
The Only Living Boy in New York 587 5.498 5.469 0.029 2.706 2.532 15 15 0
Ladies in Lavender 774 5.432 5.380 0.052 2.164 1.958 16 16 0
The Guardian Season 1 2081 5.392 5.281 0.111 1.619 1.378 17 17 0
Room 13897 5.226 5.164 0.062 1.750 1.537 18 19 -1
The Inheritance 711 5.213 5.203 0.010 2.467 2.291 19 18 1
Julian Fellowes Presents Doctor Thorne Season 1 21079 5.080 5.039 0.041 1.772 1.560 20 20 0
Philip K. Dick's Electric Dreams - Season 1 822 4.972 5.017 -0.045 3.838 3.775 21 21 0
Fortitude Season 2 922 4.860 4.895 -0.035 3.170 3.053 22 22 0
Creed 11993 4.781 4.781 0.000 1.644 1.450 23 24 -1
Hunted Season 1 3468 4.779 4.796 -0.017 1.794 1.624 24 23 1
The Infiltrator 4319 4.420 4.478 -0.058 1.646 1.496 25 25 0
Allied 846 4.369 4.437 -0.068 2.044 1.917 26 26 0
Brawl in Cell Block 99 713 4.244 4.321 -0.077 3.368 3.328 27 27 0
Rise 1012 4.040 4.139 -0.099 1.932 1.828 28 28 0
Revolutionary Road 501 3.938 4.028 -0.090 2.408 2.331 29 29 0
The Whole Truth 700 3.840 3.960 -0.120 1.652 1.550 30 30 0
Triple 9 543 3.738 3.856 -0.118 1.975 1.900 31 31 0
Silence 523 3.709 3.803 -0.094 2.876 2.838 32 32 0
The Sea of Trees 5418 3.677 3.800 -0.123 1.914 1.846 33 33 0
Tin Star Season 1 1626 3.625 3.709 -0.084 3.301 3.270 34 34 0
In Secret 509 3.441 3.590 -0.149 1.720 1.659 35 35 0
Babel 1247 3.175 3.298 -0.123 2.283 2.264 36 36 0
Max |mu_diff|: 0.354
Max |rank_diff|: 1
Findings. PyMC and Gibbs rank all 36 movies identically except for two adjacent pairs that swap one position: Taking Chance (PyMC 7 ↔ Gibbs 8) and Miss Sloane (PyMC 8 ↔ Gibbs 7); and Sneaky Pete / The Choice (11 ↔ 12), Room / The Inheritance (18 ↔ 19), Creed / Hunted (23 ↔ 24). In every case |rank_diff| ≤ 1 and the movies are back-to-back, so they are effectively tied.
Scale shift pattern. mu_diff is positive in the top half of the ranking and negative in the bottom half. This is a direct consequence of the threshold difference: PyMC tau_2=2.248 and tau_3=3.096 are lower than Gibbs tau_2=2.511 and tau_3=3.475, shifting the latent scale such that high-quality movies need slightly larger μ and low-quality movies end up with slightly smaller μ to reproduce the same observed star-rating proportions.
All three main findings replicated:
- Miss Sloane (high-polarisation outlier, rank 7–8) reproduced
- Priceless (rank 4, sigma=4.484) vs Doctor Thorne (rank 20, sigma=1.772) headline pair reproduced
- sigma values are marginally larger in PyMC (median sg_diff ≈ +0.20), consistent with the wider HalfNormal(1.5) prior vs the tighter Gamma(1,0.44) precision prior in Gibbs
5b. Ranking uncertainty — posterior overlap¶
The five rank-swap pairs between PyMC and Gibbs (|rank_diff|=1) are not a sign of disagreement between methods. They are adjacent movies whose posterior distributions of latent quality μ substantially overlap — the ranking between them is a near coin-flip in any single MCMC run.
Two complementary views:
- All 36 movies with 95% CrI — swap movies (red) have CrIs that visibly overlap their neighbours.
- KDE overlay for each swap pair — shows that P(movie A > movie B) ≈ 0.5 for all four pairs; the "swap" is sampling noise, not a model difference.
In [10]:
import matplotlib.pyplot as plt
from scipy.stats import gaussian_kde
from matplotlib.lines import Line2D
# ── Sorted data in PyMC rank order ────────────────────────────────────────────
# post_m['mu'].values: (36, 8000) — mnames order × posterior samples
mu_draws_all = post_m['mu'].values # (36, 8000)
mu_pm_s = mu_pm[idx] # posterior means, rank 1→36
mu_lo_s = np.quantile(mu_draws_all, 0.025, axis=1)[idx] # (36,)
mu_hi_s = np.quantile(mu_draws_all, 0.975, axis=1)[idx] # (36,)
movies_s = [mnames[i] for i in idx] # movie names, rank 1→36
def _s(name, n=25):
return name if len(name) <= n else name[:n - 1] + '…'
# Rows in merged (0-indexed, sorted by PyMC rank) where ranking swaps
swap_rows = sorted(merged.index[merged['rank_diff'] != 0].tolist())
# → [6, 7, 10, 11, 17, 18, 22, 23] (four adjacent pairs)
# ── Figure 1: All-movies CrI chart ────────────────────────────────────────────
fig, ax = plt.subplots(figsize=(10, 13))
y = np.arange(G_m - 1, -1, -1) # rank 1 at top
for i in range(G_m):
col = '#e74c3c' if i in swap_rows else '#2c7bb6'
ax.plot([mu_lo_s[i], mu_hi_s[i]], [y[i], y[i]], color=col, alpha=0.7, lw=1.6)
ax.scatter(mu_pm_s[i], y[i], color=col, s=22, zorder=3)
ax.set_yticks(y)
ax.set_yticklabels([_s(m) for m in movies_s], fontsize=8)
ax.set_xlabel('Posterior μ (95% CrI)', fontsize=10)
ax.set_title('Latent quality μ — posterior means and 95% credible intervals\n'
'(red = movies whose rank swaps between PyMC and Gibbs)', fontsize=10)
ax.legend(handles=[
Line2D([0],[0], color='#2c7bb6', marker='o', ms=5, lw=1.5, label='stable rank'),
Line2D([0],[0], color='#e74c3c', marker='o', ms=5, lw=1.5,
label='rank swap (|diff|=1, adjacent CrIs overlap)'),
], fontsize=8.5, loc='lower right')
plt.tight_layout()
plt.show()
# ── Figure 2: KDE overlay for the 4 swap pairs ───────────────────────────────
# Each pair: rows (r1, r2) in merged that are adjacent and have opposite rank_diff signs
swap_pair_rows = [(6, 7), (10, 11), (17, 18), (22, 23)]
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
for ax, (r1, r2) in zip(axes.flat, swap_pair_rows):
m1 = merged.iloc[r1]['movie']
m2 = merged.iloc[r2]['movie']
d1 = mu_draws_all[mnames.index(m1)] # (8000,) posterior draws for m1
d2 = mu_draws_all[mnames.index(m2)] # (8000,) posterior draws for m2
lo = min(d1.min(), d2.min()) - 0.4
hi = max(d1.max(), d2.max()) + 0.4
xs = np.linspace(lo, hi, 400)
y1 = gaussian_kde(d1)(xs)
y2 = gaussian_kde(d2)(xs)
ax.plot(xs, y1, '#2c7bb6', lw=2.2, label=f'{_s(m1, 22)} (rank {r1+1})')
ax.plot(xs, y2, '#e74c3c', lw=2.2, label=f'{_s(m2, 22)} (rank {r2+1})')
ax.fill_between(xs, np.minimum(y1, y2), alpha=0.22, color='purple', label='overlap')
ax.axvline(d1.mean(), color='#2c7bb6', ls='--', lw=1.1, alpha=0.6)
ax.axvline(d2.mean(), color='#e74c3c', ls='--', lw=1.1, alpha=0.6)
# Probability that rank r1 movie beats rank r2 movie
p12 = (d1 > d2).mean()
ax.text(0.97, 0.91, f'P(rank {r1+1} > rank {r2+1}) = {p12:.2f}',
transform=ax.transAxes, ha='right', fontsize=8.5,
bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.7))
ax.set_xlabel('Posterior μ')
ax.set_ylabel('Density')
ax.set_title(f'Swap pair: ranks {r1+1} ↔ {r2+1}', fontsize=9)
ax.legend(fontsize=7.5, loc='upper left')
fig.suptitle(
'Posterior μ distributions for the 4 rank-swap pairs\n'
'P ≈ 0.5 → near coin-flip; the PyMC/Gibbs disagreement is sampling noise, not a model difference',
fontsize=10)
plt.tight_layout()
plt.show()
Reading the two figures¶
Left — latent quality with credible intervals (all 36 movies). Ranked top to bottom by posterior mean $\mu$, each bar is a 95% CrI. The story is in where the bars separate. The top of the table (Mrs Maisel, Bosch, Poldark) and the bottom are cleanly resolved — those ranks are not in doubt. The red movies — the ones whose PyMC and Gibbs ranks differ — all sit in the crowded middle band, where a movie's CrI overlaps its neighbours' almost entirely. The ranking is precise at the extremes and genuinely uncertain in the middle, which is exactly where any two samplers will order near-ties differently.
Right — the four rank-swap pairs, up close. Each panel overlays the full posterior $\mu$ distributions of an adjacent pair. The two densities sit almost on top of one another, and the annotated $P(\text{higher-ranked} > \text{lower-ranked})$ hovers near 0.5 — a coin-flip. So a PyMC-vs-Gibbs rank difference of 1 is not a disagreement between the two methods: it is two movies the data cannot separate, and each independent MCMC run happens to break the tie its own way.
Takeaway. A rank computed from an ordinal model inherits the posterior uncertainty of the quantity it is built from — here, $\mu$. The honest object is the credible interval on the rank, not the point rank: the headline ordering (who is near the top, who is near the bottom) is robust, while mid-table positions should be read as overlapping bands rather than a strict 1-2-3. This is the ranking-level echo of the notebook's main point — that a single number (a mean, or here a point rank) hides structure the full distribution makes visible.
6. Save¶
In [11]:
comp.to_csv('movies_comparison_pymc.csv', index=False)
mu_draws = post_m['mu'].values.T
sg_draws = post_m['sigma'].values.T
pd.DataFrame(mu_draws, columns=mnames).to_csv('movies_mu_draws_pymc.csv', index=False)
pd.DataFrame(sg_draws, columns=mnames).to_csv('movies_sigma_draws_pymc.csv', index=False)
print('Saved: movies_comparison_pymc.csv')
print(f'Saved: posterior draws shape={mu_draws.shape} (draws x movies)')
Saved: movies_comparison_pymc.csv Saved: posterior draws shape=(8000, 36) (draws x movies)