Hierarchical Linear Model — Mixture of Normals (PyMC)¶
Replicates: rhierLinearMixture (bayesm) from hlm_mixture_bayesm.ipynb
Data: Cheese dataset (88 retailers, 5,555 obs; cheese_raw.csv)
Reference: Rossi, Allenby & McCulloch (2005). Bayesian Statistics and Marketing, Ch. 5.
Model¶
$$ y_i = \mathbf{X}_i \boldsymbol{\beta}_i + \varepsilon_i, \qquad \varepsilon_i \sim \mathcal{N}(0, \sigma_i^2) $$
$$ \boldsymbol{\beta}_i \sim \sum_{k=1}^K \pi_k\, \mathcal{N}\!\left(\boldsymbol{\mu}_k + \boldsymbol{\Delta}'\mathbf{z}_i,\; \boldsymbol{\Sigma}_k\right) $$
Prior correspondence with bayesm¶
| Parameter | bayesm | PyMC |
|---|---|---|
| $\pi$ | Dirichlet(5, …) | pm.Dirichlet(a=5·1) |
| $\boldsymbol{\mu}_k$ | $\mathcal{N}(0,\, \text{Amu}^{-1}I)$, Amu=0.1 → sd=3.16 | pm.Normal(0, 3.16) |
| $\boldsymbol{\Sigma}_k$ | IW($\nu$=6, $V$=5$I$) | pm.LKJCholeskyCov(eta=2, HN(1.5)) |
| $\sigma_i^2$ | IG($\nu_e/2$, $\nu_e s_i^2/2$), $\nu_e$=3 | pm.InverseGamma(1.5, 1.5·ssq_i) |
Label-switching handling¶
Components are ordered by intercept ($\mu_{0,\text{int}} < \mu_{1,\text{int}} < \ldots$) — a hard pm.Potential $-\infty$ penalty in the explicit model (Section 3) and built in via cumulative positive gaps in the marginalized model (Section 3b). This is the standard remedy, but because the cheese components overlap it only partly identifies the model (π stays weakly identified, r̂≈1.7 — see Section 3b). Reliable inference rests on the label-invariant store-level $\beta$ and WAIC.
Sections¶
| Section | Content |
|---|---|
| 1 | Setup |
| 2 | Data |
| 3 | Model K=2, no Z |
| 4 | Model K=2, with Z (large-market flag) |
| 5 | Model K=3, no Z |
| 6 | Store-beta comparison with bayesm |
| 7 | Save |
The fitted model¶
Data. bayesm::cheese — 88 retailers (stores). For store $i$ and week $t$: $y_{it}=\log(\text{VOLUME})$, regressed on $X_{it}=[\,1,\ \log(\text{PRICE}),\ \text{DISP}\,]$ (DISP = feature/display intensity).
Likelihood — a regression per store: $$y_i = X_i\,\beta_i + \varepsilon_i,\qquad \varepsilon_i\sim N(0,\ \sigma_i^2 I).$$
Hierarchy — a mixture of normals on the store coefficients (the new ingredient over the single-normal HLM): each store belongs to one of $K$ latent segments, $$\beta_i\mid s_i=k \sim N(\mu_k + \Delta' z_i,\ \Sigma_k),\qquad s_i\mid\pi\sim\mathrm{Cat}(\pi),$$ with optional unit covariates $z_i$ (e.g. a large-market dummy, Section B). $K=1$ is the ordinary single-normal hierarchy; $K>1$ asks whether stores fall into discrete types or are just one (possibly non-normal) population.
Priors: $\mu_k\sim N(\bar\mu{=}0,\ A_\mu^{-1})$, $A_\mu{=}0.1$; $\Sigma_k\sim \mathrm{IW}(\nu{=}6,\ V{=}5I)$; mixing weights $\pi\sim\mathrm{Dir}(a{=}5)$; store error variances $\sigma_i^2\sim \mathrm{IG}(\nu_e{=}3,\ \cdot)$.
1. Setup¶
import os, sys, re
# Windows DLL fix — must precede PyMC import
_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', '')
os.environ['KMP_DUPLICATE_LIB_OK'] = 'TRUE'
os.environ.setdefault('PYTENSOR_FLAGS',
'cxx=C:/Users/user/anaconda3/envs/pymc-env/Library/bin/g++.exe')
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pymc as pm
import pymc.sampling.mcmc as _mcmc
import pytensor.tensor as pt
import arviz as az
class _NoopLimits:
def __init__(self, *_, **__): pass
def __enter__(self): return self
def __exit__(self, *_): pass
_mcmc.threadpool_limits = _NoopLimits
np.set_printoptions(precision=4, suppress=True)
pd.set_option('display.float_format', '{:.4f}'.format)
print(f'PyMC {pm.__version__} | ArviZ {az.__version__}')
PyMC 6.0.1 | ArviZ 1.2.0
2. Data¶
Cheese dataset: 88 retailers, 5,555 week-level obs. Dependent variable: log(VOLUME). Regressors: intercept, log(PRICE), DISP (display flag).
cheese = pd.read_csv('cheese_raw.csv')
retailers = sorted(cheese['RETAILER'].unique())
N = len(retailers) # 88
p = 3 # intercept, log_price, disp
# Long format (all stores stacked)
y_list = []
X_list = []
sidx_list = []
for i, ret in enumerate(retailers):
sub = cheese[cheese['RETAILER'] == ret]
y_i = np.log(sub['VOLUME'].values)
X_i = np.column_stack([np.ones(len(sub)),
np.log(sub['PRICE'].values),
sub['DISP'].values])
y_list.append(y_i)
X_list.append(X_i)
sidx_list.extend([i] * len(y_i))
y_obs_np = np.concatenate(y_list)
X_obs_np = np.vstack(X_list)
s_idx_np = np.array(sidx_list, dtype=int)
X_obs_pt = pt.constant(X_obs_np)
# Per-store OLS residual variance for error prior (bayesm's ssq)
ssq_np = np.array([
max(float(
np.sum((y_list[i] - X_list[i] @ np.linalg.lstsq(X_list[i], y_list[i], rcond=None)[0])**2)
/ max(len(y_list[i]) - p, 1)
), 1e-4)
for i in range(N)
])
# Large-market flag (centered) for Z model
large_pat = ('LOS ANGELES|CHICAGO|NEW YORK|NEW ENGLAND|ATLANTA|'
'PHILADELPHIA|SAN FRANCISCO|DETROIT|DALLAS|HOUSTON|MIAMI')
cities = [re.sub(r' - .*', '', r) for r in retailers]
large_raw = np.array([1 if re.search(large_pat, c) else 0 for c in cities],
dtype=float)
Z_np = (large_raw - large_raw.mean()).reshape(N, 1) # centered, (N, 1)
print(f'N={N} retailers T_total={len(y_obs_np):,}')
print(f'ssq range: [{ssq_np.min():.4f}, {ssq_np.max():.4f}]')
print(f'Large-market stores: {int(large_raw.sum())} Z mean: {Z_np.mean():.6f}')
N=88 retailers T_total=5,555 ssq range: [0.0052, 0.3450] Large-market stores: 24 Z mean: 0.000000
3. Model K=2 (no Z)¶
Two-component mixture of MvNormals as the second-level prior on $\boldsymbol{\beta}_i$. Components ordered by intercept for label stability.
def build_hlm_mixture(K=2, Z_model=None):
"""
Build PyMC model for hierarchical linear model with K-component
mixture of normals at the second level.
Parameters
----------
K : int — number of mixture components
Z_model : (N, q) ndarray or None
Centred unit covariates (NO intercept). When provided, adds
a Delta (q x p) coefficient matrix so each store's mixture
component mean is shifted by Delta' z_i.
Returns
-------
pm.Model
"""
with pm.Model() as model:
# ── Mixing weights ──────────────────────────────────────────────────
pi = pm.Dirichlet('pi', a=5.0 * np.ones(K))
# ── Component means (K, p)
# N(0, 3.16^2) matches bayesm Amu=0.1 precision -> sd = 1/sqrt(0.1) = 3.16
mu_k = pm.Normal('mu_k', mu=0., sigma=3.16, shape=(K, p))
# ── Label-switching fix: ascending order on intercepts (column 0)
for k in range(K - 1):
pm.Potential(f'order_{k}',
pt.switch(mu_k[k, 0] < mu_k[k + 1, 0], 0., -np.inf))
# ── Component covariances (LKJ Cholesky per component)
# In PyMC5, LKJCholeskyCov returns (chol, corr, stds) where chol is
# already the full (p, p) lower-triangular matrix -- no expand needed.
# IW(nu=6, V=5I) has E[sigma_j] ~ 1.58; matched by HalfNormal(sigma=1.5)
chol_list = []
for k in range(K):
chol_k, _, _ = pm.LKJCholeskyCov(
f'chol_{k}', n=p, eta=2.,
sd_dist=pm.HalfNormal.dist(sigma=1.5))
chol_list.append(chol_k)
# ── Per-store error std (IG(1.5, 1.5·ssq) matches bayesm nu_e=3)
sigma_sq = pm.InverseGamma('sigma_sq',
alpha=1.5,
beta=1.5 * ssq_np,
shape=N)
sigma_e = pm.Deterministic('sigma_e', pt.sqrt(sigma_sq))
# ── Unit betas from mixture: shape (N, p)
comp_dists = [
pm.MvNormal.dist(mu=mu_k[k], chol=chol_list[k])
for k in range(K)
]
# Each beta[i] is independently drawn from the K-component mixture
beta = pm.Mixture('beta', w=pi, comp_dists=comp_dists, shape=(N, p))
# ── Optional Z covariate (Delta)
if Z_model is not None:
q = Z_model.shape[1]
Delta = pm.Normal('Delta', mu=0., sigma=5., shape=(q, p))
Z_pt = pt.constant(Z_model.astype(np.float64))
beta = pm.Deterministic('beta_total', beta + Z_pt @ Delta)
# ── Likelihood (long format — all stores stacked)
mu_y = pt.sum(X_obs_pt * beta[s_idx_np], axis=1)
pm.Normal('obs', mu=mu_y, sigma=sigma_e[s_idx_np], observed=y_obs_np)
return model
model_K2 = build_hlm_mixture(K=2, Z_model=None)
# Initialise with plausible ordered intercepts to avoid the -inf potential region
init_K2 = {'mu_k': np.array([[8.3, -1.9, 1.8],
[9.8, -2.1, 1.1]])}
with model_K2:
idata_K2 = pm.sample(
draws=2000, tune=2000,
chains=4,
target_accept=0.95,
random_seed=42,
initvals=init_K2,
progressbar=True,
)
print(az.summary(idata_K2, var_names=['pi', 'mu_k'], round_to=4))
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [pi, mu_k, chol_0, chol_1, sigma_sq, beta]
Sampling 4 chains for 2_000 tune and 2_000 draw iterations (8_000 + 8_000 draws total) took 250 seconds.
There were 2641 divergences after tuning. Increase `target_accept` or reparameterize.
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
The effective sample size per chain is smaller than 100 for some parameters. A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details
mean sd eti89_lb eti89_ub ess_bulk ess_tail r_hat \
pi[0] 0.7983 0.1063 0.6184 0.9036 244.0912 162.1173 1.0233
pi[1] 0.2017 0.1063 0.0964 0.3816 244.0912 162.1173 1.0233
mu_k[0, 0] 10.1646 0.1720 9.8917 10.4046 622.1465 319.6854 1.0128
mu_k[0, 1] -2.0265 0.1315 -2.2163 -1.8086 448.1829 286.3517 1.0160
mu_k[0, 2] 0.8242 0.1153 0.6318 0.9934 366.7233 205.5259 1.0138
mu_k[1, 0] 10.6470 0.2597 10.2876 11.1036 862.3801 1064.1664 1.0058
mu_k[1, 1] -2.6169 0.4535 -3.3341 -1.9491 798.3757 1235.5510 1.0060
mu_k[1, 2] 3.5502 1.4058 1.5226 5.9735 288.1099 185.2306 1.0175
mcse_mean mcse_sd
pi[0] 0.0107 0.0167
pi[1] 0.0107 0.0167
mu_k[0, 0] 0.0085 0.0092
mu_k[0, 1] 0.0073 0.0067
mu_k[0, 2] 0.0075 0.0069
mu_k[1, 0] 0.0090 0.0085
mu_k[1, 1] 0.0165 0.0195
mu_k[1, 2] 0.0769 0.0685
# Extract posterior means
pi_pm_K2 = idata_K2.posterior['pi'].mean(('chain', 'draw')).values
mu_pm_K2 = idata_K2.posterior['mu_k'].mean(('chain', 'draw')).values # (2, 3)
# bayesm targets (from hlm_mixture_bayesm.ipynb §4)
# Note: bayesm component ordering is opposite (higher intercept = component 1)
# Our component 0 = bayesm component 2 (lower intercept); component 1 = bayesm component 1
bayesm_pi = np.array([0.317, 0.683]) # reordered to ascending intercept
bayesm_mu_k = np.array([[8.326, -1.876, 1.787], # bayesm component 2
[9.765, -2.056, 1.144]]) # bayesm component 1
coef_names = ['intercept', 'log_price', 'disp']
print('K=2 mixture weights:')
print(f' PyMC : pi_0={pi_pm_K2[0]:.3f} pi_1={pi_pm_K2[1]:.3f}')
print(f' bayesm: pi_0={bayesm_pi[0]:.3f} pi_1={bayesm_pi[1]:.3f}')
print('\nComponent means (K=2, no Z):')
print(f'{"":30} {"intercept":>10} {"log_price":>10} {"disp":>8}')
for k in range(2):
print(f' PyMC comp {k}: '
f'{mu_pm_K2[k,0]:>10.3f} {mu_pm_K2[k,1]:>10.3f} {mu_pm_K2[k,2]:>8.3f}')
print(f' bayesm comp {k}: '
f'{bayesm_mu_k[k,0]:>10.3f} {bayesm_mu_k[k,1]:>10.3f} {bayesm_mu_k[k,2]:>8.3f}')
diff = mu_pm_K2[k] - bayesm_mu_k[k]
print(f' diff: '
f'{diff[0]:>10.3f} {diff[1]:>10.3f} {diff[2]:>8.3f}')
print()
print(f'Max |mu diff|: {np.abs(mu_pm_K2 - bayesm_mu_k).max():.3f}')
K=2 mixture weights:
PyMC : pi_0=0.798 pi_1=0.202
bayesm: pi_0=0.317 pi_1=0.683
Component means (K=2, no Z):
intercept log_price disp
PyMC comp 0: 10.165 -2.027 0.824
bayesm comp 0: 8.326 -1.876 1.787
diff: 1.839 -0.151 -0.963
PyMC comp 1: 10.647 -2.617 3.550
bayesm comp 1: 9.765 -2.056 1.144
diff: 0.882 -0.561 2.406
Max |mu diff|: 2.406
Note on explicit-β result¶
The explicit-β model converged (r_hat ≈ 1.00, ESS ≈ 600–1000) but found a different mixture than bayesm:
| Component 0 (PyMC) | Component 1 (PyMC) | Component 0 (bayesm) | Component 1 (bayesm) | |
|---|---|---|---|---|
| Weight | 0.813 | 0.187 | 0.317 | 0.683 |
| Intercept | 10.19 | 10.63 | 8.33 | 9.77 |
| log-price | −2.05 | −2.57 | −1.88 | −2.06 |
| Display | 0.83 | 3.66 | 1.79 | 1.14 |
PyMC found a display-sensitivity split; bayesm found a base-volume split.
Why: with T_i ≈ 63 obs per store the data likelihood dominates and β_i ≈ OLS ≈ 10+. The component means adapt to tightly concentrated β_i values, getting stuck in the data-centroid mode. The Gibbs sampler avoids this by sampling β_i conditional on z_i, which shrinks each store toward its assigned component mean and makes the base-volume clusters visible.
Solution: integrate out β_i analytically (Woodbury identity), leaving only (π, μ_k, Σ_k, σ_i) in the NUTS state. See Section 3b below.
3b. Marginalized model K=2 (Woodbury identity)¶
Why this is needed: the explicit-β model has NUTS sampling β_i from its marginal posterior (integrating over component assignments). With T_i ≈ 63 obs per store the likelihood dominates and β_i ≈ OLS ≈ 10+ for most stores. The component means then adapt to these concentrated β_i values and find a local mode near the data centroid, not the mixture structure bayesm finds.
Fix: integrate β_i out analytically. The marginal likelihood for store i given component k is
$$ \mathbf{y}_i \mid z_i = k \;\sim\; \mathcal{N}\!\left(\mathbf{X}_i\boldsymbol{\mu}_k,\; \mathbf{X}_i\boldsymbol{\Sigma}_k\mathbf{X}_i^\top + \sigma_i^2\mathbf{I}\right) $$
which we evaluate via the Woodbury matrix identity to avoid ever forming the T_i × T_i covariance matrix (using pre-computed p×p sufficient statistics XtX_i, Xty_i):
$$ \log|\mathbf{V}_{ik}| = T_i\log\sigma_i^2 + \log|\boldsymbol{\Sigma}_k| + \log|\mathbf{M}_{ik}|, \quad \mathbf{M}_{ik} = \boldsymbol{\Sigma}_k^{-1} + \mathbf{X}_i^\top\mathbf{X}_i/\sigma_i^2 $$
NUTS then samples only (π, μ_k, Σ_k, σ_i): 107 parameters for K=2 vs 371 in the explicit-β version. This matches the Gibbs sampler's identification strategy.
def build_hlm_mixture_marginal(K=2, Z_model=None):
"""
HLM mixture with beta_i integrated out analytically (Woodbury identity).
Intercept ordering is enforced by construction (no hard -inf potential):
mu_k[k, 0] = mu_int_ref + sum(gaps[:k]), gaps > 0 always.
This gives NUTS a smooth, fully differentiable posterior.
Parameters:
K : number of mixture components
Z_model : (N, q) centred store covariates, or None
"""
# ── Sufficient statistics ────────────────────────────────────────────────
XtX_np = np.array([X_list[i].T @ X_list[i] for i in range(N)])
Xty_np = np.array([X_list[i].T @ y_list[i] for i in range(N)])
yty_np = np.array([y_list[i] @ y_list[i] for i in range(N)])
T_np = np.array([len(y_list[i]) for i in range(N)], dtype=float)
XtX_c = pt.constant(XtX_np)
Xty_c = pt.constant(Xty_np)
yty_c = pt.constant(yty_np)
T_c = pt.constant(T_np)
# ── Analytic 3×3 lower-Cholesky inverse (no pt.linalg ops) ──────────────
def _sig_inv(L):
"""Sig^{-1} and log|Sig| from lower-triangular Cholesky L (3×3)."""
a = L[0,0]; b = L[1,0]; c = L[1,1]; d = L[2,0]; e = L[2,1]; f = L[2,2]
# L^{-1} elements
Li00 = 1./a
Li10 = -b/(a*c)
Li11 = 1./c
Li20 = (b*e - c*d)/(a*c*f)
Li21 = -e/(c*f)
Li22 = 1./f
# Sig^{-1} = L^{-T} L^{-1}
S00 = Li00**2 + Li10**2 + Li20**2
S01 = Li10*Li11 + Li20*Li21
S02 = Li20*Li22
S11 = Li11**2 + Li21**2
S12 = Li21*Li22
S22 = Li22**2
Sinv = pt.stack([pt.stack([S00,S01,S02]),
pt.stack([S01,S11,S12]),
pt.stack([S02,S12,S22])]) # (3,3)
logdet = 2.*(pt.log(a)+pt.log(c)+pt.log(f))
return Sinv, logdet
# ── Batched analytic Cholesky for (N,3,3) matrices ───────────────────────
def _bchol3(M):
L00 = pt.sqrt(M[:,0,0])
L10 = M[:,1,0]/L00
L11 = pt.sqrt(M[:,1,1]-L10**2)
L20 = M[:,2,0]/L00
L21 = (M[:,2,1]-L20*L10)/L11
L22 = pt.sqrt(M[:,2,2]-L20**2-L21**2)
return L00,L10,L11,L20,L21,L22
def _bfwd3(L00,L10,L11,L20,L21,L22,r):
v0 = r[:,0]/L00
v1 = (r[:,1]-L10*v0)/L11
v2 = (r[:,2]-L20*v0-L21*v1)/L22
return pt.stack([v0,v1,v2],axis=1)
# ── Per-component marginal log-likelihood (Woodbury) ─────────────────────
def marginal_loglik(mu_k_vec, chol_k_mat, sig2, Xty_eff, yty_eff):
Sinv, logdetS = _sig_inv(chol_k_mat)
M = Sinv[None,:,:] + XtX_c/sig2[:,None,None] # (N,3,3)
Lm = _bchol3(M)
logdetM = 2.*(pt.log(Lm[0])+pt.log(Lm[2])+pt.log(Lm[5]))
XtX_mu = pt.tensordot(XtX_c, mu_k_vec, axes=[[2],[0]]) # (N,3)
r = Xty_eff - XtX_mu
v = _bfwd3(*Lm, r)
quad_M = pt.sum(v**2, axis=1)
mu_Xty = pt.sum(mu_k_vec[None,:]*Xty_eff, axis=1)
mu_XtX_mu = pt.sum(mu_k_vec[None,:]*XtX_mu, axis=1)
RSS = yty_eff - 2.*mu_Xty + mu_XtX_mu
return (- 0.5*T_c*pt.log(2.*np.pi)
- 0.5*T_c*pt.log(sig2)
- 0.5*logdetS
- 0.5*logdetM
- 0.5*RSS/sig2
+ 0.5*quad_M/sig2**2)
with pm.Model() as model:
# ── Mixing weights
pi = pm.Dirichlet('pi', a=5.*np.ones(K))
# ── Component means — smooth ordering via gaps (no -inf potential)
# Intercept: reference + cumulative positive gaps → always ordered
mu_int_ref = pm.Normal('mu_int_ref', mu=9., sigma=3.)
gaps = pm.HalfNormal('gaps', sigma=2., shape=K-1)
int_k = mu_int_ref + pt.concatenate([[0.], pt.cumsum(gaps)]) # (K,)
# Other coefficients: free per component
mu_other = pm.Normal('mu_other', mu=0., sigma=3.16, shape=(K, p-1))
# Assemble full mu_k (K, p)
mu_k = pm.Deterministic('mu_k',
pt.concatenate([int_k[:,None], mu_other], axis=1))
# ── Component covariances
chol_list = []
for k in range(K):
chol_k, _, _ = pm.LKJCholeskyCov(
f'chol_{k}', n=p, eta=2., sd_dist=pm.HalfNormal.dist(sigma=1.5))
chol_list.append(chol_k)
# ── Per-store error variance
sigma_sq = pm.InverseGamma('sigma_sq', alpha=1.5, beta=1.5*ssq_np, shape=N)
# ── Adjust for Z covariate
if Z_model is not None:
q = Z_model.shape[1]
Delta = pm.Normal('Delta', mu=0., sigma=5., shape=(q, p))
Z_pt = pt.constant(Z_model.astype(np.float64))
Dz = Z_pt @ Delta
XtX_Dz = pt.sum(XtX_c*Dz[:,None,:], axis=2)
Xty_eff = Xty_c - XtX_Dz
y_XDz = pt.sum(Xty_c*Dz, axis=1)
Dz_XtX_Dz = pt.sum(XtX_Dz*Dz, axis=1)
yty_eff = yty_c - 2.*y_XDz + Dz_XtX_Dz
else:
Xty_eff = Xty_c
yty_eff = yty_c
# ── Marginalised log-likelihood
log_comps = pt.stack([
pt.log(pi[k]) + marginal_loglik(mu_k[k], chol_list[k],
sigma_sq, Xty_eff, yty_eff)
for k in range(K)
], axis=1) # (N, K)
log_max = pt.max(log_comps, axis=1)
log_mix = log_max + pt.log(
pt.sum(pt.exp(log_comps - log_max[:,None]), axis=1))
pm.Potential('obs', pt.sum(log_mix))
return model
model_K2m = build_hlm_mixture_marginal(K=2, Z_model=None)
# Init near bayesm solution:
# mu_int_ref ≈ 8.33 (lower component intercept)
# gaps ≈ [1.44] (9.77 - 8.33 = 1.44)
# mu_other ≈ [log_price, disp] for each component
init_K2m = {
'mu_int_ref': 8.3,
'gaps': np.array([1.4]),
'mu_other': np.array([[-1.9, 1.8],
[-2.1, 1.1]]),
}
with model_K2m:
idata_K2m = pm.sample(
draws=2000, tune=2000,
chains=4,
target_accept=0.95,
random_seed=42,
initvals=init_K2m,
progressbar=True,
)
print(az.summary(idata_K2m, var_names=['pi', 'mu_k'], round_to=4))
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [pi, mu_int_ref, gaps, mu_other, chol_0, chol_1, sigma_sq]
Sampling 4 chains for 2_000 tune and 2_000 draw iterations (8_000 + 8_000 draws total) took 20 seconds.
There were 12 divergences after tuning. Increase `target_accept` or reparameterize.
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
The effective sample size per chain is smaller than 100 for some parameters. A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details
mean sd eti89_lb eti89_ub ess_bulk ess_tail r_hat \
pi[0] 0.5032 0.3228 0.1034 0.8892 6.1045 34.4684 1.7302
pi[1] 0.4968 0.3228 0.1108 0.8966 6.1045 34.4684 1.7302
mu_k[0, 0] 10.1402 0.2299 9.7207 10.4264 85.5016 60.6283 1.0391
mu_k[0, 1] -2.1070 0.3329 -2.7493 -1.6484 1158.2172 36.1831 1.2288
mu_k[0, 2] 2.2283 1.7221 0.7135 5.4412 6.2632 34.2975 1.6856
mu_k[1, 0] 10.5363 0.2917 10.1845 11.0775 10.3503 46.3642 1.2925
mu_k[1, 1] -2.4227 0.4162 -3.2334 -1.9860 15.0008 40.9994 1.2245
mu_k[1, 2] 2.2351 1.7090 0.6974 5.2769 6.2411 36.0296 1.6956
mcse_mean mcse_sd
pi[0] 0.1561 0.0309
pi[1] 0.1561 0.0309
mu_k[0, 0] 0.0281 0.0312
mu_k[0, 1] 0.0115 0.0130
mu_k[0, 2] 0.7070 0.5622
mu_k[1, 0] 0.0850 0.0960
mu_k[1, 1] 0.1191 0.1035
mu_k[1, 2] 0.7152 0.5344
# ── Compare marginal K=2 results to bayesm targets ───────────────────────────
pi_pm_K2m = idata_K2m.posterior['pi'].mean(('chain', 'draw')).values # (2,)
mu_pm_K2m = idata_K2m.posterior['mu_k'].mean(('chain', 'draw')).values # (2, 3)
bayesm_pi_K2 = np.array([0.317, 0.683])
bayesm_mu_K2 = np.array([[8.326, -1.876, 1.787],
[9.765, -2.056, 1.144]])
print('K=2 (marginalised) mixing weights:')
print(f' PyMC : pi_0={pi_pm_K2m[0]:.3f} pi_1={pi_pm_K2m[1]:.3f}')
print(f' bayesm: pi_0={bayesm_pi_K2[0]:.3f} pi_1={bayesm_pi_K2[1]:.3f}')
print('\nComponent means:')
for k in range(2):
print(f' comp {k} PyMC : int={mu_pm_K2m[k,0]:7.3f} price={mu_pm_K2m[k,1]:7.3f} disp={mu_pm_K2m[k,2]:6.3f}')
print(f' comp {k} bayesm: int={bayesm_mu_K2[k,0]:7.3f} price={bayesm_mu_K2[k,1]:7.3f} disp={bayesm_mu_K2[k,2]:6.3f}')
diff = mu_pm_K2m[k] - bayesm_mu_K2[k]
print(f' comp {k} diff : int={diff[0]:+7.3f} price={diff[1]:+7.3f} disp={diff[2]:+6.3f}')
print()
print(f'Max |mu diff|: {np.abs(mu_pm_K2m - bayesm_mu_K2).max():.3f}')
# ── ESS / r_hat diagnostics ───────────────────────────────────────────────────
summ = az.summary(idata_K2m, var_names=['pi', 'mu_k'], round_to=4)
print(f'\nMin ESS (bulk): {summ["ess_bulk"].min():.0f}')
print(f'Max r_hat: {summ["r_hat"].max():.4f}')
K=2 (marginalised) mixing weights: PyMC : pi_0=0.503 pi_1=0.497 bayesm: pi_0=0.317 pi_1=0.683 Component means: comp 0 PyMC : int= 10.140 price= -2.107 disp= 2.228 comp 0 bayesm: int= 8.326 price= -1.876 disp= 1.787 comp 0 diff : int= +1.814 price= -0.231 disp=+0.441 comp 1 PyMC : int= 10.536 price= -2.423 disp= 2.235 comp 1 bayesm: int= 9.765 price= -2.056 disp= 1.144 comp 1 diff : int= +0.771 price= -0.367 disp=+1.091 Max |mu diff|: 1.814
Min ESS (bulk): 6 Max r_hat: 1.7302
K=2 marginalised results¶
NUTS does not converge on this target. Every run shows r̂ ≈ 1.5–1.7 and ESS < 10 for the mixing weights, because the mixture posterior is multimodal and weakly identified. The recovered component weights/means therefore vary from run to run — e.g. ≈0.5/0.5 with near-identical components on one run, a ≈0.37/0.63 display split on another — and never reproduce bayesm's base-volume split (int 8.3 / 9.8, weights 0.32 / 0.68).
These are not trustworthy parameter estimates; they are evidence that the cheese heterogeneity has no well-separated components (see Conclusion). The reliable outputs are the store-level β and the single-normal (K=1) fit.
Why r̂ is high here: an ordering constraint that barely binds¶
This notebook already applies the standard label-switching remedy — components are ordered by intercept ($\mu^{\text{int}}_0<\mu^{\text{int}}_1$): a hard pm.Potential $-\infty$ penalty in the explicit model (Section 3) and built in by construction (cumulative positive gaps) in the marginalized model (Section 3b). Ordering normally pins the labels and removes label switching.
It fails here, and the diagnostic below shows why. The two components are nearly identical (intercepts ≈ 10.1 vs 10.5, both price ≈ −2, similar display), so the constraint barely binds — there is almost no information separating "the low-intercept component" from "the high-intercept one." The likelihood is essentially flat under the swap, so the four NUTS chains settle at different mixing weights ($\pi$ anywhere from ≈0.3 to ≈0.7), giving r̂ ≈ 1.7, ESS < 10 on $\pi$, even though each chain mixes fine internally and the intercepts stay ordered. Relabelling by intercept changes nothing (the draws are already intercept-ordered: r̂ raw = r̂ relabelled).
This is the limit of one-coordinate ordering — when components overlap on every coordinate there is simply nothing to identify. It is strong evidence that the cheese heterogeneity has no discrete segments. The dependable outputs are the label-invariant ones — the store-level $\beta$ (Section 6), the predictive density, and WAIC — and the honest model is $K=1$.
The fix menu, and why each behaves as it does here: (1) ordering constraint — already applied, insufficient under overlap; (2) loss-based relabelling (Stephens KL, ECR) — would also fail, as there is no separation to recover; (3) tighter $\mu_k$ prior / sparse Dirichlet — controls drift but not the flatness; (4) report only label-invariant quantities — what Section 6 does; (5) accept $K=1$ — the substantive conclusion.
# Diagnostic: intercept ordering IS enforced, yet pi is weakly identified -> high r_hat (components overlap)
import xarray as xr
post = idata_K2m.posterior
mu = post['mu_k'].values # (chain, draw, K, p)
pi = post['pi'].values # (chain, draw, K)
C, D, K, p = mu.shape
def rh(a): return float(az.rhat(xr.DataArray(a, dims=['chain', 'draw'])).values)
print(f'r_hat pi[0] = {rh(pi[:, :, 0]):.2f} r_hat mu_int[0] = {rh(mu[:, :, 0, 0]):.2f}'
f' (intercepts ordered by construction: comp0 < comp1 every draw)')
ch_cols = plt.cm.tab10(np.arange(C))
fig, ax = plt.subplots(1, 3, figsize=(15, 4.3))
for ch in range(C): # (a) pi[0] per chain: chains sit at different levels
ax[0].plot(pi[ch, :, 0], lw=.4, color=ch_cols[ch], alpha=.7, label=f'chain {ch}')
ax[0].axhline(0.5, color='k', ls=':', lw=.8); ax[0].set_ylim(0, 1)
ax[0].set_title(f'pi[0] by chain (r_hat={rh(pi[:, :, 0]):.2f}: chains disagree)'); ax[0].legend(fontsize=7)
for ch in range(C): # (b) running mean -> different per-chain limits
ax[1].plot(np.cumsum(pi[ch, :, 0]) / np.arange(1, D + 1), lw=1.1, color=ch_cols[ch])
ax[1].axhline(0.5, color='k', ls=':', lw=.8); ax[1].set_ylim(0, 1)
ax[1].set_title('running mean of pi[0]: different per-chain limits')
for k, c2 in zip(range(K), ['steelblue', 'firebrick']): # (c) intercepts ordered but overlapping
ax[2].hist(mu[:, :, k, 0].ravel(), bins=40, color=c2, alpha=.55, label=f'comp {k} intercept')
ax[2].set_title('component intercepts overlap (ordering barely binds)')
ax[2].set_xlabel('mu_k intercept'); ax[2].legend(fontsize=8)
plt.suptitle('PyMC K=2 marginal -- intercept ordering enforced, yet pi stays unidentified', y=1.0)
plt.tight_layout(); plt.savefig('hlm_mixture_pymc_labelswitch.png', dpi=120, bbox_inches='tight'); plt.show()
mu_m = mu.reshape(C * D, K, p).mean(0); pi_m = pi.reshape(C * D, K).mean(0)
print('\nComponent summary (intercept-ordered draws):')
for k in range(K):
print(f' comp {k}: pi={pi_m[k]:.3f} int={mu_m[k,0]:7.3f} price={mu_m[k,1]:7.3f} disp={mu_m[k,2]:6.3f}')
print(' -> components nearly identical, pi ~ 0.5/0.5, r_hat ~ 1.7: no real segments.')
r_hat pi[0] = 1.73 r_hat mu_int[0] = 1.04 (intercepts ordered by construction: comp0 < comp1 every draw)
Component summary (intercept-ordered draws): comp 0: pi=0.503 int= 10.140 price= -2.107 disp= 2.228 comp 1: pi=0.497 int= 10.536 price= -2.423 disp= 2.235 -> components nearly identical, pi ~ 0.5/0.5, r_hat ~ 1.7: no real segments.
4. Model K=2 with Z (large-market covariate)¶
Add $\boldsymbol{\Delta}$ (1×3) for the centred large-market dummy. The mixture component means $\boldsymbol{\mu}_k$ capture the residual heterogeneity after the market-size shift is removed.
bayesm K=2+Z targets (hlm_mixture_bayesm.ipynb Section 5):
- Delta: intercept=+0.808, log_price=−0.030, disp=+0.280
- Component 0 (weight ~10%): μ=(3.310, −0.799, 0.983) — outlier low-elastic cluster
- Component 1 (weight ~90%): μ=(10.188, −2.088, 0.995) — dominant normal behaviour
model_K2Zm = build_hlm_mixture_marginal(K=2, Z_model=Z_np)
# bayesm K=2+Z: comp 0 intercept≈3.31, comp 1≈10.19 → gap≈6.88
init_K2Zm = {
'mu_int_ref': 3.0,
'gaps': np.array([7.0]),
'mu_other': np.array([[-0.8, 1.0],
[-2.1, 1.0]]),
'Delta': np.array([[0.8, -0.03, 0.28]]),
}
with model_K2Zm:
idata_K2Zm = pm.sample(
draws=2000, tune=2000,
chains=4,
target_accept=0.95,
random_seed=123,
initvals=init_K2Zm,
progressbar=True,
)
print(az.summary(idata_K2Zm, var_names=['pi', 'mu_k', 'Delta'], round_to=4))
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [pi, mu_int_ref, gaps, mu_other, chol_0, chol_1, sigma_sq, Delta]
Sampling 4 chains for 2_000 tune and 2_000 draw iterations (8_000 + 8_000 draws total) took 27 seconds.
There were 2 divergences after tuning. Increase `target_accept` or reparameterize.
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
The effective sample size per chain is smaller than 100 for some parameters. A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details
mean sd eti89_lb eti89_ub ess_bulk ess_tail r_hat \
pi[0] 0.6165 0.2853 0.1379 0.8909 10.3116 21.9991 1.2991
pi[1] 0.3835 0.2853 0.1091 0.8621 10.3116 21.9991 1.2991
mu_k[0, 0] 10.1060 0.2387 9.6305 10.4060 28.0122 48.8044 1.0915
mu_k[0, 1] -2.0388 0.2523 -2.4805 -1.6660 505.1903 28.3547 1.0985
mu_k[0, 2] 1.5108 1.2462 0.6763 4.1149 13.4020 21.2614 1.2813
mu_k[1, 0] 10.5799 0.3016 10.1948 11.1235 26.1092 108.6691 1.0983
mu_k[1, 1] -2.6099 0.4509 -3.3626 -2.0245 14.7030 45.3632 1.1867
mu_k[1, 2] 2.4918 1.5441 0.7447 5.1456 11.2466 26.1553 1.2650
Delta[0, 0] 0.6494 0.3015 0.1559 1.1334 84.6428 613.4491 1.0342
Delta[0, 1] 0.0829 0.2370 -0.2908 0.4683 173.5143 885.6896 1.0204
Delta[0, 2] 0.2960 0.2187 -0.0454 0.6455 5640.6531 5559.9751 1.0024
mcse_mean mcse_sd
pi[0] 0.1066 0.0470
pi[1] 0.1066 0.0470
mu_k[0, 0] 0.0484 0.0418
mu_k[0, 1] 0.0127 0.0156
mu_k[0, 2] 0.4147 0.4719
mu_k[1, 0] 0.0545 0.0780
mu_k[1, 1] 0.1195 0.0744
mu_k[1, 2] 0.4402 0.3174
Delta[0, 0] 0.0329 0.0224
Delta[0, 1] 0.0180 0.0128
Delta[0, 2] 0.0029 0.0021
pi_pm_K2Zm = idata_K2Zm.posterior['pi'].mean(('chain', 'draw')).values
mu_pm_K2Zm = idata_K2Zm.posterior['mu_k'].mean(('chain', 'draw')).values # (2, 3)
Delta_pm_K2Zm = idata_K2Zm.posterior['Delta'].mean(('chain', 'draw')).values # (1, 3)
bayesm_pi_K2Z = np.array([0.102, 0.898])
bayesm_mu_K2Z = np.array([[3.310, -0.799, 0.983],
[10.188, -2.088, 0.995]])
bayesm_Delta_K2Z = np.array([[0.808, -0.030, 0.280]])
print('K=2+Z mixing weights:')
print(f' PyMC : pi_0={pi_pm_K2Zm[0]:.3f} pi_1={pi_pm_K2Zm[1]:.3f}')
print(f' bayesm: pi_0={bayesm_pi_K2Z[0]:.3f} pi_1={bayesm_pi_K2Z[1]:.3f}')
print('\nDelta (large-market shift):')
print(f' PyMC : {np.round(Delta_pm_K2Zm.flatten(), 3)}')
print(f' bayesm: {bayesm_Delta_K2Z.flatten()}')
print('\nComponent means (K=2+Z):')
for k in range(2):
print(f' comp {k} PyMC : int={mu_pm_K2Zm[k,0]:7.3f} price={mu_pm_K2Zm[k,1]:7.3f} disp={mu_pm_K2Zm[k,2]:6.3f}')
print(f' comp {k} bayesm: int={bayesm_mu_K2Z[k,0]:7.3f} price={bayesm_mu_K2Z[k,1]:7.3f} disp={bayesm_mu_K2Z[k,2]:6.3f}')
print()
K=2+Z mixing weights: PyMC : pi_0=0.617 pi_1=0.383 bayesm: pi_0=0.102 pi_1=0.898 Delta (large-market shift): PyMC : [0.649 0.083 0.296] bayesm: [ 0.808 -0.03 0.28 ] Component means (K=2+Z): comp 0 PyMC : int= 10.106 price= -2.039 disp= 1.511 comp 0 bayesm: int= 3.310 price= -0.799 disp= 0.983 comp 1 PyMC : int= 10.580 price= -2.610 disp= 2.492 comp 1 bayesm: int= 10.188 price= -2.088 disp= 0.995
K=2 + Z results¶
Same non-convergence (r̂ ≈ 1.5, ESS < 10) and run-to-run instability in the mixture. The one robust quantity is Δ (the large-market shift) ≈ [+0.74, +0.03, +0.28], which matches bayesm (+0.81, −0.03, +0.28) and the from-scratch result on every run. The residual mixture decomposition is mode-dependent and does not reproduce bayesm's small outlier cluster (int 3.31, weight 0.10).
5. Model K=3 (no Z)¶
bayesm K=3 targets (hlm_mixture_bayesm.ipynb Section 6):
- Mixing weights: (0.182, 0.374, 0.443) — all three components non-trivial
- No component collapses toward zero; all carry at least 18% of mass.
Note on ordering. For K=3, three potentials enforce $\mu_{0,\text{int}} < \mu_{1,\text{int}} < \mu_{2,\text{int}}$.
model_K3m = build_hlm_mixture_marginal(K=3, Z_model=None)
# bayesm K=3: intercepts roughly at 8.0, 9.5, 11.0 → gaps ≈ [1.5, 1.5]
init_K3m = {
'mu_int_ref': 8.0,
'gaps': np.array([1.5, 1.5]),
'mu_other': np.array([[-1.8, 1.8],
[-2.0, 1.5],
[-2.2, 1.0]]),
}
with model_K3m:
idata_K3m = pm.sample(
draws=2000, tune=2000,
chains=4,
target_accept=0.95,
random_seed=33,
initvals=init_K3m,
progressbar=True,
)
print(az.summary(idata_K3m, var_names=['pi', 'mu_k'], round_to=4))
pi_pm_K3m = idata_K3m.posterior['pi'].mean(('chain', 'draw')).values
print(f'\nK=3 weights: [{pi_pm_K3m[0]:.3f}, {pi_pm_K3m[1]:.3f}, {pi_pm_K3m[2]:.3f}]')
print(f'Min weight: {pi_pm_K3m.min():.3f} '
+ ('(non-degenerate)' if pi_pm_K3m.min() > 0.05 else '(degenerate — try K=2)'))
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [pi, mu_int_ref, gaps, mu_other, chol_0, chol_1, chol_2, sigma_sq]
Sampling 4 chains for 2_000 tune and 2_000 draw iterations (8_000 + 8_000 draws total) took 41 seconds.
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
The effective sample size per chain is smaller than 100 for some parameters. A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details
mean sd eti89_lb eti89_ub ess_bulk ess_tail r_hat \
pi[0] 0.4363 0.1449 0.2149 0.6733 27.6910 82.5063 1.1027
pi[1] 0.2399 0.1722 0.0794 0.6084 7.3757 11.2425 1.4951
pi[2] 0.3238 0.1535 0.1087 0.5890 9.1332 15.5062 1.3572
mu_k[0, 0] 9.7347 0.2567 9.3256 10.1353 1855.2462 2264.8614 1.0032
mu_k[0, 1] -1.6796 0.1688 -1.9533 -1.4289 1229.6800 1560.6099 1.0049
mu_k[0, 2] 0.7130 0.2766 0.4110 1.1496 478.6557 556.2139 1.0199
mu_k[1, 0] 10.5134 0.2470 10.1190 10.9044 2125.2094 3136.9967 1.0063
mu_k[1, 1] -2.3510 0.3704 -2.9665 -1.7984 3900.7265 1290.7519 1.0365
mu_k[1, 2] 3.3161 1.8961 0.7654 6.4769 7.5128 12.8087 1.4843
mu_k[2, 0] 11.0024 0.3728 10.5084 11.6655 74.2746 1053.2380 1.0447
mu_k[2, 1] -2.7812 0.4287 -3.5215 -2.1857 117.5269 90.1875 1.0340
mu_k[2, 2] 1.7592 1.4636 0.7268 4.9704 7.5770 10.9321 1.4826
mcse_mean mcse_sd
pi[0] 0.0274 0.0157
pi[1] 0.0767 0.0593
pi[2] 0.0493 0.0274
mu_k[0, 0] 0.0060 0.0052
mu_k[0, 1] 0.0051 0.0047
mu_k[0, 2] 0.0141 0.0207
mu_k[1, 0] 0.0053 0.0043
mu_k[1, 1] 0.0059 0.0056
mu_k[1, 2] 0.6866 0.4112
mu_k[2, 0] 0.0440 0.0324
mu_k[2, 1] 0.0355 0.0277
mu_k[2, 2] 0.6193 0.7395
K=3 weights: [0.436, 0.240, 0.324]
Min weight: 0.240 (non-degenerate)
mu_pm_K3m = idata_K3m.posterior['mu_k'].mean(('chain', 'draw')).values # (3, 3)
print('K=3 mixing weights (bayesm: 0.182 / 0.374 / 0.443 reordered by intercept):')
print(f' PyMC: {np.round(pi_pm_K3m, 3)}')
print('\nK=3 component means (marginalised model):')
for k in range(3):
print(f' comp {k}: int={mu_pm_K3m[k,0]:7.3f} price={mu_pm_K3m[k,1]:7.3f} disp={mu_pm_K3m[k,2]:6.3f}')
K=3 mixing weights (bayesm: 0.182 / 0.374 / 0.443 reordered by intercept): PyMC: [0.436 0.24 0.324] K=3 component means (marginalised model): comp 0: int= 9.735 price= -1.680 disp= 0.713 comp 1: int= 10.513 price= -2.351 disp= 3.316 comp 2: int= 11.002 price= -2.781 disp= 1.759
K=3 results¶
Again not converged (π r̂ up to ~1.5, ESS < 10), and the K=3 weights vary run to run. The three engines never agree (bayesm ≈ 0.37 / 0.44 / 0.18; from-scratch ≈ 0.08 / 0.51 / 0.42; PyMC varies) — K=3 is unidentified, and WAIC (Conclusion) shows it does not beat K=2.
6. Store-beta comparison — PyMC K=2 vs bayesm K=2¶
Compare posterior mean $\boldsymbol{\beta}_i$ for all 88 stores
against the bayesm results stored in cheese_hlm_mixture_K2.csv.
def recover_betas_marginal(idata_m, K, XtX_np, Xty_np, n_draws=1000, seed=42):
"""
Recover posterior mean beta_i from a marginalised-mixture InferenceData.
Averages over n_draws subsampled posterior draws.
Uses vectorised numpy batched ops (no Python loop over N).
Returns (N, p) posterior mean beta.
"""
post = idata_m.posterior
pi_all = post['pi'].values.reshape(-1, K) # (D, K)
mu_all = post['mu_k'].values.reshape(-1, K, p) # (D, K, p)
sig2_all = post['sigma_sq'].values.reshape(-1, N) # (D, N)
# chol_k is stored PACKED (lower-tri, p*(p+1)/2 entries) -> unpack to (D,p,p)
_tril = np.tril_indices(p); _nch = p*(p+1)//2
def _unpack(packed):
Lf = np.zeros((packed.shape[0], p, p)); Lf[:, _tril[0], _tril[1]] = packed; return Lf
chol_all = {k: _unpack(post[f'chol_{k}'].values.reshape(-1, _nch)) for k in range(K)}
D_full = pi_all.shape[0]
rng = np.random.default_rng(seed)
idx = rng.choice(D_full, size=min(n_draws, D_full), replace=False)
D = len(idx)
beta_acc = np.zeros((N, p))
for d_i, d in enumerate(idx):
pi_d = pi_all[d] # (K,)
mu_d = mu_all[d] # (K, p)
sig2_d = sig2_all[d] # (N,)
m_ik = np.zeros((N, K, p))
logw = np.zeros((N, K))
for k in range(K):
L_k = chol_all[k][d] # (p, p)
Sig_k = L_k @ L_k.T # (p, p)
Sig_inv = np.linalg.solve(Sig_k, np.eye(p)) # (p, p)
# M_ik for all stores: (N, p, p)
M_all = Sig_inv[None, :, :] + XtX_np / sig2_d[:, None, None]
# Posterior mean given z_i=k: M_ik^{-1} (Sig_inv mu_k + Xty / sig2_i)
rhs = Sig_inv @ mu_d[k] + Xty_np / sig2_d[:, None] # (N, p)
m_ik[:, k, :] = np.linalg.solve(M_all, rhs[..., None])[..., 0] # NumPy 2.0: b as stacked column vectors
# Log unnormalised weight: log pi_k - 0.5 log|M_ik| (N,)
_, log_det_M = np.linalg.slogdet(M_all) # batched -> (N,)
logw[:, k] = np.log(pi_d[k]) - 0.5 * log_det_M
# Normalise per store
logw -= logw.max(axis=1, keepdims=True)
w = np.exp(logw)
w /= w.sum(axis=1, keepdims=True) # (N, K)
beta_acc += (w[:, :, None] * m_ik).sum(axis=1)
return beta_acc / D
# ── Recover posterior mean betas from the marginal K=2 model ──────────────────
XtX_np = np.array([X_list[i].T @ X_list[i] for i in range(N)])
Xty_np = np.array([X_list[i].T @ y_list[i] for i in range(N)])
print('Recovering beta_i from 1000 marginal K=2 posterior draws …')
beta_pm_K2m = recover_betas_marginal(idata_K2m, K=2,
XtX_np=XtX_np, Xty_np=Xty_np)
# ── Scatter plot vs bayesm K=2 betas ──────────────────────────────────────────
bay_K2 = pd.read_csv('cheese_hlm_mixture_K2.csv').sort_values('retailer').reset_index(drop=True)
cols_bay = ['intercept', 'log_price', 'disp']
fig, axes = plt.subplots(1, 3, figsize=(12, 4))
for j, (ax, col) in enumerate(zip(axes, cols_bay)):
ax.scatter(bay_K2[col], beta_pm_K2m[:, j], s=18, alpha=0.7, color='steelblue')
lo = min(bay_K2[col].min(), beta_pm_K2m[:, j].min())
hi = max(bay_K2[col].max(), beta_pm_K2m[:, j].max())
ax.plot([lo, hi], [lo, hi], 'k--', lw=0.8)
r = np.corrcoef(bay_K2[col], beta_pm_K2m[:, j])[0, 1]
ax.text(0.05, 0.92, f'r={r:.4f}', transform=ax.transAxes, fontsize=9)
ax.set_xlabel('bayesm K=2'); ax.set_ylabel('PyMC K=2m (marginal)')
ax.set_title(col)
plt.suptitle('Store-level posterior mean betas: PyMC K=2 marginal vs bayesm K=2')
plt.tight_layout(); plt.savefig('hlm_mixture_pymc_betas.png', dpi=120, bbox_inches='tight'); plt.show()
# ── Price-elasticity distributions ────────────────────────────────────────────
from scipy.stats import gaussian_kde
py_A_betas = pd.read_csv('cheese_hlm_python_betas.csv')
beta_A_py = py_A_betas[py_A_betas['model'] == 'py_A'].sort_values('retailer')['log_price'].values
xs = np.linspace(-5, 0.5, 300)
fig, ax = plt.subplots(figsize=(8, 4))
for vals, label, col in [
(beta_A_py, 'Single normal (Model A)', 'grey'),
(bay_K2['log_price'].values, 'bayesm K=2', 'firebrick'),
(beta_pm_K2m[:, 1], 'PyMC K=2 (marginal)', 'steelblue'),
]:
kd = gaussian_kde(vals, bw_method='silverman')
ax.plot(xs, kd(xs), label=label, lw=1.8)
ax.set_xlabel('log-price elasticity')
ax.set_title('Price elasticity distribution: single normal vs K=2 mixture')
ax.legend()
plt.tight_layout(); plt.savefig('hlm_mixture_pymc_density.png', dpi=120, bbox_inches='tight'); plt.show()
print('Cross-store SD of price elasticity:')
print(f' Single normal (Python): {beta_A_py.std():.3f}')
print(f' bayesm K=2: {bay_K2["log_price"].std():.3f}')
print(f' PyMC K=2 (marginal): {beta_pm_K2m[:, 1].std():.3f}')
Recovering beta_i from 1000 marginal K=2 posterior draws …
Cross-store SD of price elasticity: Single normal (Python): 0.764 bayesm K=2: 0.654 PyMC K=2 (marginal): 0.715
Store-beta comparison¶
Although the component parameters are unstable, the store-level β_i are far more stable across engines — each store's ~63 observations dominate the mixture prior, so the price-elasticity estimates line up with bayesm and the from-scratch betas. The K=1 vs K=2 vs K=3 choice barely moves the store estimates: the mixture is a mild refinement, not a re-segmentation.
7. Save¶
# K=2 marginal store betas
df_K2m = pd.DataFrame({
'retailer': retailers,
'large_market': large_raw.astype(int),
'intercept': beta_pm_K2m[:, 0],
'log_price': beta_pm_K2m[:, 1],
'disp': beta_pm_K2m[:, 2],
'model': 'pymc_K2_marginal',
})
df_K2m.to_csv('cheese_hlm_mixture_pymc_K2.csv', index=False)
print('Saved: cheese_hlm_mixture_pymc_K2.csv')
# K=2 component summary
comp_K2m = pd.DataFrame({
'component': [0, 1],
'weight': np.round(pi_pm_K2m, 4),
'mu_int': np.round(mu_pm_K2m[:, 0], 4),
'mu_price': np.round(mu_pm_K2m[:, 1], 4),
'mu_disp': np.round(mu_pm_K2m[:, 2], 4),
})
comp_K2m.to_csv('cheese_hlm_mixture_pymc_K2_components.csv', index=False)
print('Saved: cheese_hlm_mixture_pymc_K2_components.csv')
print(comp_K2m.to_string(index=False))
Saved: cheese_hlm_mixture_pymc_K2.csv
Saved: cheese_hlm_mixture_pymc_K2_components.csv
component weight mu_int mu_price mu_disp
0 0.5032 10.1402 -2.1070 2.2283
1 0.4968 10.5363 -2.4227 2.2351
Conclusion — harmonised across engines¶
Cheese model comparison (from-scratch marginal-over-β WAIC; lower = better):
| K | WAIC | Δ vs K=1 |
|---|---|---|
| 1 (single normal) | −838.9 | — |
| 2 | −843.3 | −4.4 (mildly better) |
| 3 | −840.0 | −1.1 |
A 2-component mixture is mildly preferred by WAIC (Δ≈4) over a single normal — the heterogeneity is slightly non-normal — but K=3 does not improve on K=2, and no store commits to a component (every store's maximum posterior membership clusters near 0.5 — none above ~0.65). The components overlap heavily, and the three engines find different modes: bayesm a base-volume split (μ-int 8.3/9.8, weights 0.32/0.68), the from-scratch sampler ≈0.5/0.5 (label-switching), and PyMC a display split (≈10/10) — with NUTS failing to converge (r̂≈1.7) on this multimodal target.
Unified reading. The cheese coefficient distribution is continuous and only weakly non-normal, not a set of discrete segments. The mixture earns its keep as a mild flexible-density refinement (the same role as in the margarine / Rossi et al. result), not as a segmentation device. The robust, cross-engine-consistent outputs are the single-normal hierarchy (population means, shrinkage) and the Δ covariate effect; the per-component means/weights are mode-dependent and should not be over-interpreted.