Scotch Whisky — Binary Multivariate Probit via Gibbs Sampler¶
Reference: Edwards, Y. D. & Allenby, G. M. (2003). Multivariate Analysis of Multiple Response Data. Journal of Marketing Research, 40(3), 321–334.
Model¶
Let $i = 1, \ldots, N$ index households and $j = 1, \ldots, J$ index brands. The binary multivariate probit model posits a latent utility structure:
$$ \mathbf{Z} = \mathbf{X}\mathbf{B}^\top + \mathbf{E}, \qquad \mathbf{E} \sim \mathcal{MN}\bigl(\mathbf{0},\, \mathbf{I}_N \otimes \boldsymbol{\Sigma}\bigr) $$
where
- $\mathbf{Z}$ is the $N \times J$ matrix of latent utilities
- $\mathbf{X}$ is the $N \times K$ design matrix (here intercept-only, $K=1$)
- $\mathbf{B}$ is the $J \times K$ matrix of regression coefficients (brand thresholds)
- $\boldsymbol{\Sigma}$ is the $J \times J$ covariance matrix of utilities
The observed binary outcomes are
$$ y_{ij} = \mathbf{1}\{z_{ij} > 0\}, \qquad i=1,\ldots,N,\; j=1,\ldots,J $$
Identification¶
$\boldsymbol{\Sigma}$ is only identified up to scale. We impose the standard normalisation at each draw:
$$ \mathbf{T} = \operatorname{diag}\!\left(\sigma_{11}^{-1/2},\ldots,\sigma_{JJ}^{-1/2}\right), \qquad \mathbf{R} = \mathbf{T}\boldsymbol{\Sigma}\mathbf{T}, \qquad \mathbf{B}^* = \mathbf{T}\mathbf{B} $$
so the stored draws are of the identified parameters $(\mathbf{B}^*, \mathbf{R})$.
Priors¶
$$ \operatorname{vec}(\mathbf{B}) \sim \mathcal{N}\bigl(\mathbf{b}_0,\, \mathbf{B}_0\bigr) $$
$$ \boldsymbol{\Sigma} \sim \mathcal{IW}(\mathbf{V},\, \nu) $$
where $\mathcal{IW}(\mathbf{V}, \nu)$ denotes the inverse-Wishart distribution with scale matrix $\mathbf{V}$ and degrees of freedom $\nu$, so $\mathbb{E}[\boldsymbol{\Sigma}] = \mathbf{V}/(\nu - J - 1)$.
| Parameter | Value used | Interpretation |
|---|---|---|
| $\mathbf{b}_0$ | $\mathbf{0}_{KJ}$ | flat centre on brand thresholds |
| $\mathbf{B}_0$ | $100\cdot\mathbf{I}_{KJ}$ | diffuse (near-flat) |
| $\nu$ | $J+1$ | minimal proper df for IW |
| $\mathbf{V}$ | $\mathbf{I}_J$ | prior mean $\boldsymbol{\Sigma} \propto \mathbf{I}$ |
Gibbs Sampler Algorithm¶
Initialise $\mathbf{B}^{(0)},\boldsymbol{\Sigma}^{(0)},\mathbf{W}^{(0)}$. At each iteration $t$:
Block 1 — Latent utilities $\mathbf{W}$ (Albert–Chib data augmentation)¶
For each brand $j = 1,\ldots,J$, sample $w_{ij}$ from its conditional truncated normal:
$$ w_{ij} \mid \mathbf{w}_{i,-j},\mathbf{B},\boldsymbol{\Sigma} \;\sim\; \mathcal{TN}\!\left(\mu_{ij\mid -j},\;\sigma^2_{j\mid -j};\; \begin{cases}(0,+\infty) & y_{ij}=1 \\ (-\infty,0) & y_{ij}=0\end{cases}\right) $$
with conditional mean and variance from the multivariate normal partition:
$$ \boldsymbol{f}_j = \boldsymbol{\Sigma}_{-j,-j}^{-1}\boldsymbol{\Sigma}_{-j,j}, \qquad \sigma^2_{j\mid -j} = \Sigma_{jj} - \boldsymbol{\Sigma}_{-j,j}^\top\boldsymbol{f}_j $$
$$ \mu_{ij\mid -j} = (\mathbf{x}_i\mathbf{B}^\top)_j + \boldsymbol{f}_j^\top\bigl(\mathbf{w}_{i,-j} - (\mathbf{x}_i\mathbf{B}^\top)_{-j}\bigr) $$
Block 2 — Regression coefficients $\mathbf{B}$¶
Let $\mathbf{C} = \operatorname{chol}(\boldsymbol{\Sigma}^{-1})$ (upper Cholesky). Posterior:
$$ \operatorname{vec}(\mathbf{B}) \mid \mathbf{W},\boldsymbol{\Sigma} \;\sim\; \mathcal{N}\!\left(\boldsymbol{\mu}_B,\;\boldsymbol{\Omega}_B^{-1}\right) $$
$$ \boldsymbol{\Omega}_B = \mathbf{X}^\top\mathbf{X}\otimes\mathbf{C}\mathbf{C}^\top + \mathbf{B}_0^{-1}, \qquad \boldsymbol{\mu}_B = \boldsymbol{\Omega}_B^{-1} \bigl(\operatorname{vec}(\mathbf{C}(\mathbf{W}\mathbf{C})^\top\mathbf{X}) + \mathbf{B}_0^{-1}\mathbf{b}_0\bigr) $$
Block 3 — Covariance matrix $\boldsymbol{\Sigma}$¶
Residuals $\mathbf{E} = \mathbf{W} - \mathbf{X}\mathbf{B}^\top$. Posterior:
$$ \boldsymbol{\Sigma}^{-1} \mid \mathbf{W},\mathbf{B} \;\sim\; \mathcal{W}\!\left((\mathbf{E}^\top\mathbf{E}+\mathbf{V})^{-1},\;N+\nu\right), \qquad \boldsymbol{\Sigma} = \bigl(\boldsymbol{\Sigma}^{-1}\bigr)^{-1} $$
Normalise: compute $\mathbf{T},\mathbf{R},\mathbf{B}^*$ and store.
import os, sys, time
_conda_lib = r"C:\Users\user\anaconda3\envs\pymc-env\Library\bin"
if _conda_lib not in os.environ.get("PATH", ""):
os.environ["PATH"] = _conda_lib + os.pathsep + os.environ.get("PATH", "")
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.special import ndtr
import importlib, mvprobit_gibbs
importlib.reload(mvprobit_gibbs)
from mvprobit_gibbs import mvprobit_gibbs as gibbs_sampler, make_prior, make_init
print("Imports OK")
Imports OK
1. Load data¶
raw = pd.read_excel("Scotch.xls", sheet_name="Scotch", header=4)
raw = raw.dropna(how="all").dropna(axis=1, how="all")
raw = raw.rename(columns={"Unnamed: 2": "household_id"})
BRANDS = [c for c in raw.columns if c != "household_id"]
Y = raw[BRANDS].fillna(0).values.astype(int)
N, J = Y.shape
BRAND_TYPE = {
"Chivas.Regal": "Blended",
"Dewar.s.White.Label": "Blended",
"Johnnie.Walker.Black.Label": "Blended",
"J...B": "Blended",
"Johnnie.Walker.Red.Label": "Blended",
"Other.Brands": "Other",
"Glenlivet": "Single Malt",
"Cutty.Sark": "Blended",
"Glenfiddich": "Single Malt",
"Pinch..Haig.": "Blended",
"Clan.MacGregor": "Blended",
"Ballantine": "Blended",
"Macallan": "Single Malt",
"Passport": "Blended",
"Black...White": "Blended",
"Scoresby.Rare": "Blended",
"Grants": "Blended",
"Ushers": "Blended",
"White.Horse": "Blended",
"Knockando": "Single Malt",
"the.Singleton": "Single Malt",
}
BRAND_PRICE = {
"Chivas.Regal": 21.99,
"Dewar.s.White.Label": 17.99,
"Johnnie.Walker.Black.Label": 22.99,
"J...B": 18.99,
"Johnnie.Walker.Red.Label": 18.99,
"Other.Brands": np.nan,
"Glenlivet": 22.99,
"Cutty.Sark": 15.99,
"Glenfiddich": 39.99,
"Pinch..Haig.": 24.99,
"Clan.MacGregor": 10.00,
"Ballantine": 14.90,
"Macallan": 32.99,
"Passport": 10.90,
"Black...White": 12.10,
"Scoresby.Rare": 10.60,
"Grants": 12.50,
"Ushers": 13.56,
"White.Horse": 16.99,
"Knockando": 33.99,
"the.Singleton": 28.99,
}
BRAND_ORIGIN = {
"Chivas.Regal": "Foreign",
"Dewar.s.White.Label": "Foreign",
"Johnnie.Walker.Black.Label": "Foreign",
"J...B": "Foreign",
"Johnnie.Walker.Red.Label": "Foreign",
"Other.Brands": "Other",
"Glenlivet": "Foreign",
"Cutty.Sark": "Foreign",
"Glenfiddich": "Foreign",
"Pinch..Haig.": "Foreign",
"Clan.MacGregor": "Domestic",
"Ballantine": "Foreign",
"Macallan": "Foreign",
"Passport": "Domestic",
"Black...White": "Foreign",
"Scoresby.Rare": "Domestic",
"Grants": "Foreign",
"Ushers": "Foreign",
"White.Horse": "Foreign",
"Knockando": "Foreign",
"the.Singleton": "Foreign",
}
BRAND_SYM = {
"Chivas.Regal": "CHR", "Dewar.s.White.Label": "DWL",
"Johnnie.Walker.Black.Label": "JWB", "J...B": "JaB",
"Johnnie.Walker.Red.Label": "JWR", "Other.Brands": "OTH",
"Glenlivet": "GLT", "Cutty.Sark": "CTY", "Glenfiddich": "GFH",
"Pinch..Haig.": "PCH", "Clan.MacGregor": "MCG", "Ballantine": "BAL",
"Macallan": "MCL", "Passport": "PAS", "Black...White": "BaW",
"Scoresby.Rare": "SCY", "Grants": "GRT", "Ushers": "USH",
"White.Horse": "WHT", "Knockando": "KND", "the.Singleton": "SGT",
}
print(f"N households : {N}")
print(f"J brands : {J}")
print(f"Avg brands/HH: {Y.sum(axis=1).mean():.2f}")
rates = Y.mean(axis=0)
purchase_rates = pd.Series(rates, index=BRANDS).sort_values(ascending=False)
df = purchase_rates.round(3).to_frame("rate")
df["type"] = [BRAND_TYPE.get(b, "Other") for b in df.index]
df["price"] = [BRAND_PRICE.get(b, np.nan) for b in df.index]
df["origin"] = [BRAND_ORIGIN.get(b, "Other") for b in df.index]
df
N households : 2219 J brands : 21 Avg brands/HH: 2.29
| rate | type | |
|---|---|---|
| Chivas.Regal | 0.363 | Blended |
| Dewar.s.White.Label | 0.233 | Blended |
| Johnnie.Walker.Black.Label | 0.226 | Blended |
| J...B | 0.206 | Blended |
| Johnnie.Walker.Red.Label | 0.191 | Blended |
| Other.Brands | 0.187 | Other |
| Glenlivet | 0.160 | Single Malt |
| Cutty.Sark | 0.153 | Blended |
| Glenfiddich | 0.151 | Single Malt |
| Pinch..Haig. | 0.053 | Blended |
| Clan.MacGregor | 0.046 | Blended |
| Ballantine | 0.045 | Blended |
| Macallan | 0.043 | Single Malt |
| Passport | 0.037 | Blended |
| Black...White | 0.037 | Blended |
| Scoresby.Rare | 0.036 | Blended |
| Grants | 0.033 | Blended |
| Ushers | 0.030 | Blended |
| White.Horse | 0.028 | Blended |
| Knockando | 0.021 | Single Malt |
| the.Singleton | 0.014 | Single Malt |
fig, ax = plt.subplots(figsize=(11, 4))
purchase_rates.plot(kind="bar", ax=ax, color="steelblue", edgecolor="white")
ax.set_ylabel("Purchase rate")
ax.set_title("Scotch whisky brand purchase rates (N=2,219 households)")
plt.tight_layout()
plt.show()
2. Set up model¶
X = np.ones((N, 1)) # intercept-only, K=1
K = 1
prior = make_prior(J, K, v=J+1, V=np.eye(J), beta_var=100.0)
b0 = make_init(Y, X, J, K)
print(f"W parameters : {N*J:,}")
print(f"Beta params : {K*J}")
print(f"Sigma params : {J*(J+1)//2}")
W parameters : 46,599 Beta params : 21 Sigma params : 231
3. Run Gibbs sampler¶
NDRAWS = 60000
BURN = 5000
print(f"Running Gibbs sampler: {NDRAWS} draws, {BURN} burn-in ...")
t0 = time.time()
results = gibbs_sampler(Y, X, prior, ndraws=NDRAWS, b0=b0, burn=BURN, verbose=500)
elapsed = time.time() - t0
print(f"Done in {elapsed/60:.1f} minutes.")
Running Gibbs sampler: 60000 draws, 5000 burn-in ... Iteration 500 / 60000 Iteration 1000 / 60000 Iteration 1500 / 60000 Iteration 2000 / 60000 Iteration 2500 / 60000 Iteration 3000 / 60000 Iteration 3500 / 60000 Iteration 4000 / 60000 Iteration 4500 / 60000 Iteration 5000 / 60000 Iteration 5500 / 60000 Iteration 6000 / 60000 Iteration 6500 / 60000 Iteration 7000 / 60000 Iteration 7500 / 60000 Iteration 8000 / 60000 Iteration 8500 / 60000 Iteration 9000 / 60000 Iteration 9500 / 60000 Iteration 10000 / 60000 Iteration 10500 / 60000 Iteration 11000 / 60000 Iteration 11500 / 60000 Iteration 12000 / 60000 Iteration 12500 / 60000 Iteration 13000 / 60000 Iteration 13500 / 60000 Iteration 14000 / 60000 Iteration 14500 / 60000 Iteration 15000 / 60000 Iteration 15500 / 60000 Iteration 16000 / 60000 Iteration 16500 / 60000 Iteration 17000 / 60000 Iteration 17500 / 60000 Iteration 18000 / 60000 Iteration 18500 / 60000 Iteration 19000 / 60000 Iteration 19500 / 60000 Iteration 20000 / 60000 Iteration 20500 / 60000 Iteration 21000 / 60000 Iteration 21500 / 60000 Iteration 22000 / 60000 Iteration 22500 / 60000 Iteration 23000 / 60000 Iteration 23500 / 60000 Iteration 24000 / 60000 Iteration 24500 / 60000 Iteration 25000 / 60000 Iteration 25500 / 60000 Iteration 26000 / 60000 Iteration 26500 / 60000 Iteration 27000 / 60000 Iteration 27500 / 60000 Iteration 28000 / 60000 Iteration 28500 / 60000 Iteration 29000 / 60000 Iteration 29500 / 60000 Iteration 30000 / 60000 Iteration 30500 / 60000 Iteration 31000 / 60000 Iteration 31500 / 60000 Iteration 32000 / 60000 Iteration 32500 / 60000 Iteration 33000 / 60000 Iteration 33500 / 60000 Iteration 34000 / 60000 Iteration 34500 / 60000 Iteration 35000 / 60000 Iteration 35500 / 60000 Iteration 36000 / 60000 Iteration 36500 / 60000 Iteration 37000 / 60000 Iteration 37500 / 60000 Iteration 38000 / 60000 Iteration 38500 / 60000 Iteration 39000 / 60000 Iteration 39500 / 60000 Iteration 40000 / 60000 Iteration 40500 / 60000 Iteration 41000 / 60000 Iteration 41500 / 60000 Iteration 42000 / 60000 Iteration 42500 / 60000 Iteration 43000 / 60000 Iteration 43500 / 60000 Iteration 44000 / 60000 Iteration 44500 / 60000 Iteration 45000 / 60000 Iteration 45500 / 60000 Iteration 46000 / 60000 Iteration 46500 / 60000 Iteration 47000 / 60000 Iteration 47500 / 60000 Iteration 48000 / 60000 Iteration 48500 / 60000 Iteration 49000 / 60000 Iteration 49500 / 60000 Iteration 50000 / 60000 Iteration 50500 / 60000 Iteration 51000 / 60000 Iteration 51500 / 60000 Iteration 52000 / 60000 Iteration 52500 / 60000 Iteration 53000 / 60000 Iteration 53500 / 60000 Iteration 54000 / 60000 Iteration 54500 / 60000 Iteration 55000 / 60000 Iteration 55500 / 60000 Iteration 56000 / 60000 Iteration 56500 / 60000 Iteration 57000 / 60000 Iteration 57500 / 60000 Iteration 58000 / 60000 Iteration 58500 / 60000 Iteration 59000 / 60000 Iteration 59500 / 60000 Iteration 60000 / 60000 Done in 36.4 minutes.
4. Extract posterior draws¶
B_draws = results['B_draws'] # (keep, K*J)
Sigma_draws = results['Sigma_draws'] # (keep, J*J)
keep = B_draws.shape[0]
B_mean = B_draws.mean(axis=0).reshape(J, K, order='F') # (J, K)
R_lower = Sigma_draws.reshape(keep, J, J, order='F')
R_draws = R_lower + np.triu(R_lower.transpose(0, 2, 1), k=1)
R_mean = R_draws.mean(axis=0)
R_std = R_draws.std(axis=0)
print(f"Stored draws : {keep:,}")
print(f"B_draws shape: {B_draws.shape}")
print(f"R_mean shape : {R_mean.shape}")
Stored draws : 55,000 B_draws shape: (55000, 21) R_mean shape : (21, 21)
5. MCMC Diagnostics¶
Check before interpreting results:
- Trace plots — visual mixing and stationarity
- ACF — autocorrelation to determine thinning
- ESS — effective independent draws
- Running mean — burn-in adequacy
from statsmodels.tsa.stattools import acf as compute_acf
N_MONITOR = 4 # worst-mixing parameters to show from each group
def ess(chain, max_lag=500):
n = len(chain)
acf_v = compute_acf(chain, nlags=min(max_lag, n // 2), fft=True)
ci = 1.96 / np.sqrt(n)
rho = 0.0
for k in range(1, len(acf_v)):
if abs(acf_v[k]) < ci:
break
rho += acf_v[k]
return n / max(1.0, 1 + 2 * rho)
# ── ESS for every parameter ───────────────────────────────────────────────────
print("Computing ESS for all parameters ...", end=" ", flush=True)
beta_ess = np.array([ess(B_draws[:, j]) for j in range(J)])
rows_tri, cols_tri = np.tril_indices(J, k=-1)
corr_ess = np.array([ess(Sigma_draws[:, r + c * J])
for r, c in zip(rows_tri, cols_tri)])
print("done.")
# ── Select N_MONITOR worst-mixing from each group ─────────────────────────────
worst_beta = np.argsort(beta_ess)[:N_MONITOR]
worst_corr = np.argsort(corr_ess)[:N_MONITOR]
monitor_beta_idx = list(worst_beta)
monitor_beta_lbl = [BRANDS[i] for i in monitor_beta_idx]
monitor_pairs = [(int(rows_tri[i]), int(cols_tri[i])) for i in worst_corr]
monitor_corr_lbl = [f"R[{BRANDS[r][:8]},{BRANDS[c][:8]}]"
for r, c in monitor_pairs]
print(f"\nWorst-mixing beta parameters (N={N_MONITOR}):")
for i in worst_beta:
print(f" {BRANDS[i]:<35s} ESS = {beta_ess[i]:6.0f}")
print(f"\nWorst-mixing correlation pairs (N={N_MONITOR}):")
for idx, (r, c) in zip(worst_corr, monitor_pairs):
print(f" R[{BRANDS[r]}, {BRANDS[c]}] ESS = {corr_ess[idx]:6.0f}")
# ── Chain arrays used by all downstream diagnostic cells ─────────────────────
beta_chains = B_draws[:, monitor_beta_idx]
corr_chains = np.column_stack([
Sigma_draws[:, r + c * J] for r, c in monitor_pairs
])
all_chains = np.hstack([beta_chains, corr_chains])
all_labels = monitor_beta_lbl + monitor_corr_lbl
n_params = all_chains.shape[1]
iters = np.arange(1, keep + 1)
print(f"\nMonitoring {n_params} parameters over {keep:,} draws")
Computing ESS for all parameters ... done. Worst-mixing beta parameters (N=4): the.Singleton ESS = 3221 Knockando ESS = 5727 Ushers ESS = 8050 White.Horse ESS = 8706 Worst-mixing correlation pairs (N=4): R[the.Singleton, Chivas.Regal] ESS = 210 R[the.Singleton, Ushers] ESS = 253 R[the.Singleton, Pinch..Haig.] ESS = 305 R[the.Singleton, Scoresby.Rare] ESS = 315 Monitoring 8 parameters over 55,000 draws
# ── Trace plots ───────────────────────────────────────────────────────────────
fig, axes = plt.subplots(n_params, 1, figsize=(12, 2.2 * n_params), sharex=True)
for ax, chain, lbl in zip(axes, all_chains.T, all_labels):
ax.plot(iters, chain, lw=0.3, color="steelblue", alpha=0.8)
ax.axhline(chain.mean(), color="red", lw=1, linestyle="--")
ax.set_ylabel(lbl, fontsize=8, rotation=0, ha="right", labelpad=65)
axes[-1].set_xlabel("Post-burn-in draw")
fig.suptitle("Trace plots", fontsize=11, y=1.01)
plt.tight_layout()
plt.show()
# ── ACF plots ─────────────────────────────────────────────────────────────────
NLAGS = 80
ci = 1.96 / np.sqrt(keep)
fig, axes = plt.subplots(n_params, 1, figsize=(12, 2.2 * n_params), sharex=True)
for ax, chain, lbl in zip(axes, all_chains.T, all_labels):
acf_vals = compute_acf(chain, nlags=NLAGS, fft=True)[1:]
ax.bar(np.arange(1, NLAGS+1), acf_vals, color="steelblue", width=0.8, alpha=0.7)
ax.axhline( ci, color="red", lw=1, linestyle="--")
ax.axhline(-ci, color="red", lw=1, linestyle="--")
ax.axhline(0, color="black", lw=0.5)
ax.set_ylim(-0.3, 1.0)
ax.set_ylabel(lbl, fontsize=8, rotation=0, ha="right", labelpad=65)
axes[-1].set_xlabel("Lag")
fig.suptitle("ACF — bars inside red lines = no thinning needed", fontsize=11, y=1.01)
plt.tight_layout()
plt.show()
# ── Effective Sample Size (monitored parameters) ─────────────────────────────
ess_vals = [ess(c) for c in all_chains.T]
ess_df = pd.DataFrame({
"parameter": all_labels,
"draws" : keep,
"ESS" : [round(e) for e in ess_vals],
"ESS/draw" : [f"{e/keep:.3f}" for e in ess_vals],
"thin_by" : [f"~{max(1, round(keep/e))}" for e in ess_vals],
}).set_index("parameter")
print(ess_df.to_string())
print(f"\nMin ESS = {min(ess_vals):.0f} ({min(ess_vals)/keep:.1%} efficiency)")
print(f"Recommended thinning: keep every {max(1, round(keep/min(ess_vals)))}th draw")
draws ESS ESS/draw thin_by parameter the.Singleton 55000 3221 0.059 ~17 Knockando 55000 5727 0.104 ~10 Ushers 55000 8050 0.146 ~7 White.Horse 55000 8706 0.158 ~6 R[the.Sing,Chivas.R] 55000 210 0.004 ~262 R[the.Sing,Ushers] 55000 253 0.005 ~218 R[the.Sing,Pinch..H] 55000 305 0.006 ~181 R[the.Sing,Scoresby] 55000 315 0.006 ~175 Min ESS = 210 (0.4% efficiency) Recommended thinning: keep every 262th draw
# ── Running mean ──────────────────────────────────────────────────────────────
fig, axes = plt.subplots(n_params, 1, figsize=(12, 2.2 * n_params), sharex=True)
for ax, chain, lbl in zip(axes, all_chains.T, all_labels):
rm = np.cumsum(chain) / iters
ax.plot(iters, rm, lw=1.2, color="steelblue")
ax.axhline(chain.mean(), color="red", lw=1, linestyle="--")
ax.set_ylabel(lbl, fontsize=8, rotation=0, ha="right", labelpad=65)
axes[-1].set_xlabel("Post-burn-in draw")
fig.suptitle("Running mean — flat = burn-in adequate", fontsize=11, y=1.01)
plt.tight_layout()
plt.show()
6. Brand intercepts vs observed purchase rates¶
With the Gibbs sampler, $\Phi(\beta_j)$ should match the observed rate (no sequential GHK bias).
beta_int = B_mean[:, 0]
impl_rates = ndtr(beta_int)
comparison = pd.DataFrame({
"brand" : BRANDS,
"observed" : rates.round(3),
"model_implied" : impl_rates.round(3),
"beta" : beta_int.round(3),
}).set_index("brand").sort_values("observed", ascending=False)
comparison
| observed | model_implied | beta | |
|---|---|---|---|
| brand | |||
| Chivas.Regal | 0.363 | 0.363 | -0.351 |
| Dewar.s.White.Label | 0.233 | 0.234 | -0.726 |
| Johnnie.Walker.Black.Label | 0.226 | 0.227 | -0.749 |
| J...B | 0.206 | 0.207 | -0.817 |
| Johnnie.Walker.Red.Label | 0.191 | 0.193 | -0.867 |
| Other.Brands | 0.187 | 0.187 | -0.889 |
| Glenlivet | 0.160 | 0.160 | -0.995 |
| Cutty.Sark | 0.153 | 0.154 | -1.020 |
| Glenfiddich | 0.151 | 0.151 | -1.031 |
| Pinch..Haig. | 0.053 | 0.052 | -1.621 |
| Clan.MacGregor | 0.046 | 0.047 | -1.674 |
| Ballantine | 0.045 | 0.045 | -1.696 |
| Macallan | 0.043 | 0.044 | -1.710 |
| Passport | 0.037 | 0.038 | -1.780 |
| Black...White | 0.037 | 0.037 | -1.791 |
| Scoresby.Rare | 0.036 | 0.036 | -1.795 |
| Grants | 0.033 | 0.034 | -1.828 |
| Ushers | 0.030 | 0.031 | -1.869 |
| White.Horse | 0.028 | 0.028 | -1.918 |
| Knockando | 0.021 | 0.021 | -2.030 |
| the.Singleton | 0.014 | 0.013 | -2.224 |
fig, ax = plt.subplots(figsize=(6, 6))
ax.scatter(comparison["observed"], comparison["model_implied"], s=60, color="steelblue")
lims = [0, 0.42]
ax.plot(lims, lims, "k--", lw=1, label="45° line")
for _, row in comparison.iterrows():
ax.annotate(row.name.split(".")[0], (row["observed"], row["model_implied"]),
fontsize=7, xytext=(4, 2), textcoords="offset points")
ax.set_xlabel("Observed purchase rate")
ax.set_ylabel("Model-implied $\\Phi(\\beta)$")
ax.set_title("Gibbs: observed vs model-implied rates")
ax.legend()
plt.tight_layout()
plt.show()
Note — MCMC mixing and purchase rate.
The Singleton has the most extreme intercept ($\beta = -2.22$, purchase rate 1.4%), and it is also the worst-mixing parameter in the sampler (Section 5: beta ESS = 2,967; correlations involving the Singleton have ESS as low as 226 out of 55,000 draws).
The reason is structural: with so few households purchasing the Singleton, the likelihood provides almost no information to anchor its latent utility. The sampler must explore a wide range of $z_{i,\text{Singleton}}$ values at each iteration, producing a highly autocorrelated chain. The same mechanism affects other rare brands (Knockando, Macallan) and any correlation pair that involves them.
7. Top correlations¶
rows, cols = np.tril_indices(J, k=-1)
corr_df = pd.DataFrame({
"brand_1" : [BRANDS[r] for r in rows],
"brand_2" : [BRANDS[c] for c in cols],
"correlation": R_mean[rows, cols].round(3),
"post_std" : R_std[rows, cols].round(3),
}).sort_values("correlation", key=abs, ascending=False)
corr_df.head(20)
| brand_1 | brand_2 | correlation | post_std | |
|---|---|---|---|---|
| 183 | Knockando | Macallan | 0.675 | 0.053 |
| 34 | Glenfiddich | Glenlivet | 0.598 | 0.032 |
| 207 | the.Singleton | Ushers | 0.556 | 0.082 |
| 205 | the.Singleton | Scoresby.Rare | 0.516 | 0.078 |
| 151 | Ushers | Scoresby.Rare | 0.496 | 0.073 |
| 203 | the.Singleton | Passport | 0.494 | 0.083 |
| 202 | the.Singleton | Macallan | 0.492 | 0.083 |
| 179 | Knockando | Glenfiddich | 0.484 | 0.062 |
| 177 | Knockando | Glenlivet | 0.471 | 0.060 |
| 133 | Grants | Passport | 0.468 | 0.069 |
| 208 | the.Singleton | White.Horse | 0.467 | 0.086 |
| 152 | Ushers | Grants | 0.453 | 0.080 |
| 72 | Macallan | Glenlivet | 0.442 | 0.050 |
| 149 | Ushers | Passport | 0.433 | 0.079 |
| 164 | White.Horse | Ballantine | 0.430 | 0.076 |
| 102 | Black...White | Ballantine | 0.428 | 0.071 |
| 167 | White.Horse | Black...White | 0.425 | 0.075 |
| 88 | Passport | Clan.MacGregor | 0.408 | 0.068 |
| 166 | White.Horse | Passport | 0.406 | 0.080 |
| 209 | the.Singleton | Knockando | 0.406 | 0.092 |
8. Correlation heatmap¶
short = [b.replace(".", " ").replace("Johnnie Walker ", "JW ") for b in BRANDS]
fig, ax = plt.subplots(figsize=(14, 11))
mask = np.triu(np.ones_like(R_mean, dtype=bool), k=1)
sns.heatmap(
R_mean, mask=mask,
annot=True, fmt=".2f", center=0,
cmap="RdBu_r", vmin=-0.5, vmax=0.5,
xticklabels=short, yticklabels=short,
ax=ax, linewidths=0.3, annot_kws={"size": 7},
)
ax.set_title(
f"Posterior mean purchase correlation matrix\n"
f"Edwards-Allenby Gibbs | {keep:,} draws | N={N}, J={J}"
)
plt.xticks(rotation=45, ha="right", fontsize=8)
plt.yticks(rotation=0, fontsize=8)
plt.show()
Correlation findings¶
1. Single-malt enthusiast cluster (correlations 0.42–0.68)¶
The four core single malts — Knockando, Macallan, Glenfiddich, Glenlivet — dominate the top of the ranking. The single strongest correlation in the entire matrix is Knockando–Macallan (0.675), followed by Glenfiddich–Glenlivet (0.598). All four are Speyside malts. This reflects a clear category-enthusiast segment: consumers who buy single malts buy several of them, exploring across brands rather than substituting.
The tetrachoric framing matters here: these brands all have low purchase rates (4–16%), so their Pearson/$\phi$ coefficients are severely attenuated. The latent-utility correlations reveal that the underlying preference intensity for these brands co-moves strongly, even though few households actually buy any of them.
2. Value-blend cluster (correlations 0.41–0.50)¶
A second tight cluster emerges among the budget blends: Ushers, Grants, Passport, Scoresby Rare, White Horse, Black & White, Ballantine, Clan MacGregor. Pairs within this group are densely represented in the top 20, with Ushers–Scoresby Rare (0.496) and Grants–Passport (0.468) as the strongest. This is a price-sensitive switcher segment — households buying one value blend tend to buy others, likely driven by promotions rather than brand loyalty.
The Singleton also appears in this cluster (correlations with Ushers 0.56, Scoresby Rare 0.52, Passport 0.49), sitting above several within-cluster pairs. This is likely a small-sample artefact — with a purchase rate of only 1.4%, posterior SDs for Singleton pairs are 0.08–0.09 and its chain mixes poorly (see Section 5 and the note in Section 6). The Singleton is better treated as an outlier than a genuine cluster member.
3. Mainstream premium blends are absent from the top¶
Chivas Regal, Dewar's, Johnnie Walker Black/Red, J&B, Cutty Sark — the five most-purchased brands — produce no pairs in the top 20. Two explanations:
- High penetration dilutes the co-purchase signal across heterogeneous household types.
- These brands likely substitute within-tier: a Dewar's household is not a Johnnie Walker household.
The JW Red–JW Black correlation falls below the top-20 threshold despite brand-family overlap, suggesting consumers trade up from Red to Black rather than buy both.
4. No cross-cluster spillover¶
Single malts do not correlate strongly with mainstream blends, and value blends do not correlate with single malts. The two enthusiast groups shop entirely within their own segment.
Summary: three latent consumer types emerge — single-malt collectors, value-blend switchers, and mainstream-blend loyalists — with almost no between-segment sharing. The multivariate probit recovers this preference structure from binary data; marginal purchase rates alone would not reveal it.
9. PCA of the Posterior Correlation Matrix¶
Principal components analysis of R decomposes the joint purchase structure into orthogonal dimensions ordered by the variance they explain (Edwards & Allenby 2003, Figure 3). The eigenvalues of a $J\times J$ correlation matrix sum to $J=21$; components with eigenvalue $>1$ explain more than a single brand's worth of variance (Kaiser criterion).
The paper reports 5 components $>1$, accounting for 57% of variance, with eigenvalues 4.68, 2.85, 1.87, 1.36, 1.12, and interprets them as:
| Component | Associated with |
|---|---|
| PC1 | Brand market share (purchase rate) |
| PC2 | Type of Scotch (single malt vs blend) |
| PC3 | Price |
| PC4–5 | Residual structure |
from scipy.stats import pearsonr
from matplotlib.lines import Line2D
# ── Eigendecomposition ────────────────────────────────────────────────────────
eigvals, eigvecs = np.linalg.eigh(R_mean)
eigvals = eigvals[::-1] # descending order
eigvecs = eigvecs[:, ::-1]
total_var = eigvals.sum() # = J = 21 for a correlation matrix
pct_var = eigvals / total_var * 100
cum_var = np.cumsum(pct_var)
n_gt1 = int((eigvals > 1).sum())
print(f"Eigenvalues > 1 : {n_gt1} ({cum_var[n_gt1-1]:.1f}% of variance)")
print(f"Paper reports : 5 (57.0% of variance)")
print(f"\nTop-10 eigenvalues:")
print(f" {'PC':<5} {'lambda':>8} {'% var':>7} {'cum %':>7}")
for i in range(10):
mark = " <--" if eigvals[i] > 1 else ""
print(f" PC{i+1:<3} {eigvals[i]:8.3f} {pct_var[i]:7.2f} {cum_var[i]:7.2f}{mark}")
# ── Scree plot ────────────────────────────────────────────────────────────────
fig, ax = plt.subplots(figsize=(9, 4))
colors_scree = ["steelblue" if v > 1 else "lightsteelblue" for v in eigvals]
ax.bar(range(1, J + 1), eigvals, color=colors_scree, edgecolor="white")
ax.axhline(1, color="red", lw=1.2, linestyle="--", label="Kaiser criterion (\u03bb = 1)")
ax.set_xlabel("Principal component")
ax.set_ylabel("Eigenvalue")
ax.set_title("Scree plot — PCA of posterior mean correlation matrix R")
ax.set_xticks(range(1, J + 1))
ax.legend()
plt.tight_layout()
plt.show()
# ── Biplot PC1 x PC2 ─────────────────────────────────────────────────────────
pc1 = eigvecs[:, 0]
pc2 = eigvecs[:, 1]
type_colors = {"Single Malt": "darkorange", "Blended": "steelblue", "Other": "gray"}
pt_colors = [type_colors.get(BRAND_TYPE.get(b, "Other")) for b in BRANDS]
sym_labels = [BRAND_SYM.get(b, b[:4]) for b in BRANDS]
fig, ax = plt.subplots(figsize=(9, 7))
for x, y, lbl, col in zip(pc1, pc2, sym_labels, pt_colors):
ax.scatter(x, y, color=col, s=70, zorder=3, edgecolors="white", linewidths=0.5)
ax.annotate(lbl, (x, y), fontsize=8, xytext=(5, 3), textcoords="offset points")
ax.axhline(0, color="k", lw=0.5)
ax.axvline(0, color="k", lw=0.5)
ax.set_xlabel(f"PC1 ({pct_var[0]:.1f}% of variance)")
ax.set_ylabel(f"PC2 ({pct_var[1]:.1f}% of variance)")
ax.set_title("Component loadings — PC1 \u00d7 PC2\n(orange = Single Malt, blue = Blended)")
handles = [Line2D([0],[0], marker="o", color="w", markerfacecolor=v,
markersize=9, label=k)
for k, v in type_colors.items() if k != "Other"]
ax.legend(handles=handles, fontsize=9)
plt.tight_layout()
plt.show()
# ── Attribute correlations with PC1–PC3 ───────────────────────────────────────
attr_rate = rates
attr_type = np.array([1 if BRAND_TYPE.get(b) == "Single Malt" else 0
for b in BRANDS], dtype=float)
attr_price = np.array([BRAND_PRICE.get(b, np.nan) for b in BRANDS])
attr_origin = np.array([1 if BRAND_ORIGIN.get(b) == "Domestic" else 0
for b in BRANDS], dtype=float)
price_mask = ~np.isnan(attr_price) # exclude Other.Brands
print(f"\nCorrelation of PC loadings with brand attributes:")
print(f" {'Attribute':<10} {'PC1':>8} {'PC2':>8} {'PC3':>8}")
print(" " + "-" * 38)
for name, attr in [("Rate", attr_rate),
("Type", attr_type),
("Price", attr_price),
("Origin", attr_origin)]:
row = []
for k in range(3):
pc = eigvecs[:, k]
mask = price_mask if name == "Price" else np.ones(J, dtype=bool)
r, p = pearsonr(pc[mask], attr[mask])
sig = "*" if p < 0.05 else " "
row.append(f"{r:+.3f}{sig}")
print(f" {name:<10} {' '.join(row)}")
print(" * p < 0.05")
print(f"\n Paper: PC1 vs Rate r = -0.932")
Eigenvalues > 1 : 5 (58.4% of variance) Paper reports : 5 (57.0% of variance) Top-10 eigenvalues: PC lambda % var cum % PC1 4.791 22.82 22.82 <-- PC2 2.925 13.93 36.74 <-- PC3 1.954 9.30 46.05 <-- PC4 1.449 6.90 52.94 <-- PC5 1.151 5.48 58.43 <-- PC6 0.957 4.56 62.98 PC7 0.896 4.27 67.25 PC8 0.870 4.14 71.39 PC9 0.786 3.74 75.13 PC10 0.710 3.38 78.52
Correlation of PC loadings with brand attributes: Attribute PC1 PC2 PC3 -------------------------------------- Rate -0.932* +0.223 +0.546* Type +0.220 -0.736* +0.246 Price -0.131 -0.593* +0.544* Origin +0.303 -0.033 -0.416 * p < 0.05 Paper: PC1 vs Rate r = -0.932
Interpretation¶
The five components with eigenvalue $>1$ account for 58.4% of the total variance in purchase correlations (paper reports 57%), with individual eigenvalues 4.79, 2.93, 1.95, 1.45, 1.15 — closely replicating Table 5 of Edwards & Allenby (2003).
PC1 — Market share ($r = -0.932^*$ with purchase rate). The dominant dimension separates popular brands (negative loadings: Chivas Regal, Dewar's, Johnnie Walker) from rare ones (positive loadings: Knockando, the Singleton, Macallan). This is the Ehrenberg (1995) double jeopardy pattern: brands with low market share are bought by fewer consumers who also buy them less often.
PC2 — Type of Scotch ($r = -0.736^*$ with Single Malt indicator; $r = -0.593^*$ with price). Single malts load negatively (top of the biplot per the paper) and blended Scotches positively. Price loads alongside type because single malts are systematically more expensive ($23–$40) than budget blends ($10–$17); once type is controlled, price adds little independent information on this axis.
PC3 — Price, conditional on type ($r = +0.544^*$ with price; $r = +0.546^*$ with rate). After type is absorbed by PC2, PC3 separates the mainstream premium blends (Chivas $21.99, JW Black $22.99 — expensive and widely purchased) from the cheap domestic blends (Clan MacGregor $10.00, Scoresby Rare $10.60, Passport $10.90). The simultaneous correlation with rate reflects that the premium mainstream blends occupy a distinct position: high share and high price.
Origin (Domestic vs Foreign) is not independently significant on any component, as all three domestic brands (Clan MacGregor, Scoresby Rare, Passport) fall at the cheap end of the price range and are already captured by PC2–PC3.
10. Tetrachoric Correlations — Three Estimators¶
| Estimator | Method | Notes |
|---|---|---|
| Pearson / phi | np.corrcoef(Y.T) |
Biased toward 0 — bounded by marginal rates |
| Pairwise MLE tetrachoric | Bivariate probit MLE per pair | Corrects attenuation; pairs treated independently |
Gibbs tetrachoric R_mean |
Edwards-Allenby (2003) joint posterior mean | All pairs estimated simultaneously; positive-definite by construction; full posterior uncertainty |
Pearson correlations on binary data are bounded by $[\phi_{\min}, \phi_{\max}]$ determined by the marginal purchase rates, so they systematically understate latent associations. Both tetrachoric estimators correct for this. The Gibbs version is preferred because the full $J\times J$ correlation matrix is estimated jointly, guaranteeing positive-definiteness and producing posterior intervals.
from scipy.stats import multivariate_normal, norm as sp_norm
from scipy.optimize import minimize_scalar
def _bvn_cdf(h, k, rho):
rho = float(np.clip(rho, -0.9999, 0.9999))
return multivariate_normal.cdf([h, k], mean=[0., 0.], cov=[[1., rho], [rho, 1.]])
def _tetrachoric_pair(yj, yk):
pj, pk = yj.mean(), yk.mean()
# edge case: no variation
if pj in (0., 1.) or pk in (0., 1.):
return np.nan
tau_j, tau_k = sp_norm.ppf(pj), sp_norm.ppf(pk)
n11 = int(((yj == 1) & (yk == 1)).sum())
n10 = int(((yj == 1) & (yk == 0)).sum())
n01 = int(((yj == 0) & (yk == 1)).sum())
n00 = int(((yj == 0) & (yk == 0)).sum())
def neg_ll(rho):
p11 = _bvn_cdf(tau_j, tau_k, rho)
p10 = max(pj - p11, 1e-12)
p01 = max(pk - p11, 1e-12)
p00 = max(1 - pj - pk + p11, 1e-12)
p11 = max(p11, 1e-12)
return -(n11*np.log(p11) + n10*np.log(p10) + n01*np.log(p01) + n00*np.log(p00))
res = minimize_scalar(neg_ll, bounds=(-0.9999, 0.9999), method="bounded")
return res.x
# ── Pairwise MLE tetrachoric ─────────────────────────────────────────────────
print("Computing pairwise MLE tetrachoric correlations (210 pairs) …", end=" ")
R_tetra = np.eye(J)
for j in range(J):
for k in range(j + 1, J):
rho = _tetrachoric_pair(Y[:, j], Y[:, k])
R_tetra[j, k] = rho
R_tetra[k, j] = rho
print("done.")
print(f" min eigenvalue: {np.linalg.eigvalsh(R_tetra).min():.4f} "
f"(>0 → positive-definite)")
# ── Naive Pearson on binary Y ─────────────────────────────────────────────────
R_pearson = np.corrcoef(Y.T)
print(f" Pearson range: [{R_pearson[np.tril_indices(J,-1)].min():.3f}, "
f"{R_pearson[np.tril_indices(J,-1)].max():.3f}]")
Computing pairwise MLE tetrachoric correlations (210 pairs) … done. min eigenvalue: -0.0025 (>0 → positive-definite) Pearson range: [-0.133, 0.360]
# ── Three-way scatter comparison ────────────────────────────────────────────
rows, cols = np.tril_indices(J, k=-1)
gibbs_vals = R_mean[rows, cols]
tetra_vals = R_tetra[rows, cols]
pearson_vals = R_pearson[rows, cols]
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
ax = axes[0]
ax.scatter(pearson_vals, gibbs_vals, s=20, alpha=0.55, color="steelblue", edgecolors="none")
lim = [min(pearson_vals.min(), gibbs_vals.min()) - 0.03,
max(pearson_vals.max(), gibbs_vals.max()) + 0.03]
ax.plot(lim, lim, "k--", lw=1, label="y = x")
ax.set_xlabel("Pearson / phi (binary Y)")
ax.set_ylabel("Gibbs tetrachoric $\\bar{R}$")
ax.set_title("Attenuation: Pearson vs Tetrachoric")
ax.set_xlim(lim); ax.set_ylim(lim)
ax.legend(fontsize=8)
ax = axes[1]
ax.scatter(tetra_vals, gibbs_vals, s=20, alpha=0.55, color="darkorange", edgecolors="none")
lim2 = [min(tetra_vals.min(), gibbs_vals.min()) - 0.03,
max(tetra_vals.max(), gibbs_vals.max()) + 0.03]
ax.plot(lim2, lim2, "k--", lw=1, label="y = x")
ax.set_xlabel("Pairwise MLE tetrachoric")
ax.set_ylabel("Gibbs tetrachoric $\\bar{R}$")
ax.set_title("Pairwise MLE vs Gibbs (joint) tetrachoric")
ax.set_xlim(lim2); ax.set_ylim(lim2)
ax.legend(fontsize=8)
plt.tight_layout()
plt.show()
# ── Summary statistics ───────────────────────────────────────────────────────
print(f"{'Method':<26} {'mean':>6} {'std':>6} {'min':>6} {'max':>6}")
print("-" * 56)
for name, vals in [("Pearson (binary)", pearson_vals),
("Pairwise MLE tetrachoric", tetra_vals),
("Gibbs tetrachoric", gibbs_vals)]:
print(f"{name:<26} {vals.mean():6.3f} {vals.std():6.3f} "
f"{vals.min():6.3f} {vals.max():6.3f}")
Method mean std min max -------------------------------------------------------- Pearson (binary) 0.102 0.082 -0.133 0.360 Pairwise MLE tetrachoric 0.276 0.201 -0.258 0.772 Gibbs tetrachoric 0.158 0.175 -0.323 0.675
11. Save draws¶
np.save("scotch_B_draws.npy", B_draws)
np.save("scotch_Sigma_draws.npy", Sigma_draws)
np.save("scotch_R_mean.npy", R_mean)
print("Saved: scotch_B_draws.npy, scotch_Sigma_draws.npy, scotch_R_mean.npy")
Saved: scotch_B_draws.npy, scotch_Sigma_draws.npy, scotch_R_mean.npy