Hierarchical Linear Model Mixture -- Python Gibbs Sampler¶
Module: hlm_mixture_gibbs.py
Replicates: bayesm::rhierLinearMixture (R)
Reference: Rossi, Allenby & McCulloch (2005), Bayesian Statistics and Marketing, Ch. 5.
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, §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)$.
In [1]:
import os, sys, subprocess, re, time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import gaussian_kde
# Windows DLL fix for conda envs
if sys.platform == 'win32':
import ctypes
ctypes.CDLL('vcruntime140_1.dll', mode=ctypes.RTLD_GLOBAL)
sys.path.insert(0, os.getcwd())
from hlm_mixture_gibbs import (rhier_linear_mixture, simulate_hlm_mixture,
posterior_summary)
np.set_printoptions(precision=3, suppress=True)
pd.set_option('display.float_format', '{:.4f}'.format)
print('hlm_mixture_gibbs loaded OK')
hlm_mixture_gibbs loaded OK
Part 1 -- Synthetic Data: Parameter Recovery¶
Simulate from a known K=2 mixture with well-separated components and verify the Gibbs sampler recovers component means, mixing weights, and store betas.
In [2]:
K_syn = 2
p_syn = 3 # intercept + 2 slopes
true_pi = np.array([0.40, 0.60])
true_mu = np.array([[5.0, -2.0, 1.0],
[9.0, -1.0, 0.5]])
true_Sigma = np.array([0.5 * np.eye(p_syn),
0.5 * np.eye(p_syn)])
regdata_syn, betas_syn, comps_syn = simulate_hlm_mixture(
nreg=120, p=p_syn, K=K_syn,
true_pi=true_pi, true_mu=true_mu, true_Sigma=true_Sigma,
true_sigma2=1.0, ni=30, seed=0)
print(f'Simulated: 120 stores x 30 obs, K={K_syn}, p={p_syn}')
print(f'True component assignments: {np.bincount(comps_syn)}')
print(f'True mu[0] (component 0): {true_mu[0]}')
print(f'True mu[1] (component 1): {true_mu[1]}')
Simulated: 120 stores x 30 obs, K=2, p=3 True component assignments: [41 79] True mu[0] (component 0): [ 5. -2. 1.] True mu[1] (component 1): [ 9. -1. 0.5]
In [3]:
Prior_syn = {
'a': 5.,
'mubar': np.zeros(p_syn),
'Amu': 0.1,
'nu': p_syn + 3,
'V': float(p_syn + 2) * np.eye(p_syn),
'nu_e': 3.0,
'ssq': 1.0,
}
Mcmc_syn = {'R': 15000, 'burn': 5000, 'keep': 1, 'nprint': 5000}
t0 = time.perf_counter()
out_syn = rhier_linear_mixture(regdata_syn, K=K_syn,
Prior=Prior_syn, Mcmc=Mcmc_syn, seed=42)
print(f'Done in {time.perf_counter()-t0:.1f}s')
Iter 5000/15000 (40s)
Iter 10000/15000 (81s)
Iter 15000/15000 (121s) Done in 2.0 min | kept: 10000 Done in 121.1s
In [4]:
# Component means (sorted by intercept ascending -- matches our label-switching rule)
mu_pm = out_syn['mudraw'].mean(0) # (K, p)
pi_pm = out_syn['pidraw'].mean(0) # (K,)
print('Component parameter recovery (sorted by intercept ascending):')
print(f' True pi: {true_pi}')
print(f' Est pi: {pi_pm.round(3)}')
print()
for k in range(K_syn):
print(f' Component {k}: true mu = {true_mu[k]} | est mu = {mu_pm[k].round(3)}')
# Store beta recovery
beta_pm_syn = out_syn['betadraw'].mean(axis=2) # (nreg, p)
ols_syn = np.array([
np.linalg.lstsq(d['X'], d['y'], rcond=None)[0]
for d in regdata_syn
])
for j, name in enumerate(['Intercept', 'Slope x1', 'Slope x2']):
corr_ols = np.corrcoef(betas_syn[:, j], ols_syn[:, j])[0, 1]
corr_bayes = np.corrcoef(betas_syn[:, j], beta_pm_syn[:, j])[0, 1]
print(f' {name}: r(true,OLS)={corr_ols:.3f} r(true,Bayes)={corr_bayes:.3f}')
# Component assignment accuracy
s_pm = np.round(out_syn['sdraw'].mean(0)).astype(int)
acc = (s_pm == comps_syn).mean()
print(f'\nMode component assignment accuracy: {acc:.1%}')
Component parameter recovery (sorted by intercept ascending): True pi: [0.4 0.6] Est pi: [0.356 0.644] Component 0: true mu = [ 5. -2. 1.] | est mu = [ 5.068 -2.052 0.866] Component 1: true mu = [ 9. -1. 0.5] | est mu = [ 9.027 -1.104 0.502] Intercept: r(true,OLS)=0.995 r(true,Bayes)=0.995 Slope x1: r(true,OLS)=0.976 r(true,Bayes)=0.976 Slope x2: r(true,OLS)=0.970 r(true,Bayes)=0.970 Mode component assignment accuracy: 100.0%
Part 2 -- Cheese Data¶
Replicate the rhierLinearMixture analysis from hlm_mixture_bayesm.ipynb
using the Python Gibbs sampler.
Target values from bayesm (K=2, no Z, 25K draws, 5K burn):
| Component | Weight | Intercept | log-price | Display |
|---|---|---|---|---|
| 1 (lower $\mu$) | 0.317 | 8.326 | -1.876 | 1.787 |
| 2 (higher $\mu$) | 0.683 | 9.765 | -2.056 | 1.144 |
In [5]:
R_EXE = r'C:\Program Files\R\R-4.6.0\bin\x64\Rscript.exe'
if not os.path.exists('cheese_raw.csv'):
r_code = '''
userLib <- file.path(Sys.getenv('USERPROFILE'), 'R', 'win-library', '4.6')
.libPaths(c(userLib, .libPaths()))
library(bayesm)
data(cheese)
write.csv(cheese, 'cheese_raw.csv', row.names=FALSE)
cat('Exported cheese_raw.csv\n')
'''
with open('_tmp_export.R', 'w') as f:
f.write(r_code)
res = subprocess.run([R_EXE, '_tmp_export.R'], capture_output=True, text=True)
os.remove('_tmp_export.R')
print(res.stdout.strip())
else:
print('cheese_raw.csv already exists')
cheese = pd.read_csv('cheese_raw.csv')
retailers = sorted(cheese['RETAILER'].unique())
nreg_c = len(retailers) # 88
nvar_c = 3 # intercept, log(price), disp
regdata_c = []
for ret in retailers:
sub = cheese[cheese['RETAILER'] == ret]
y = np.log(sub['VOLUME'].values)
X = np.column_stack([np.ones(len(sub)),
np.log(sub['PRICE'].values),
sub['DISP'].values])
regdata_c.append({'y': y, 'X': X})
# OLS residual variances for ssq prior
ssq_c = np.array([
max(float(np.sum((d['y'] - d['X'] @ np.linalg.lstsq(d['X'], d['y'], rcond=None)[0])**2)
/ max(len(d['y']) - nvar_c, 1)), 1e-4)
for d in regdata_c
])
# Large-market flag (centred -- matches bayesm Z_c)
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])
Z_c = (large_raw - large_raw.mean())[:, None] # (88, 1) centred, no intercept
print(f'Retailers: {nreg_c} Large market: {large_raw.sum()}')
print(f'ssq range: [{ssq_c.min():.4f}, {ssq_c.max():.4f}]')
print(f'Z (centred large-market) mean: {Z_c.mean():.6f}')
cheese_raw.csv already exists Retailers: 88 Large market: 24 ssq range: [0.0052, 0.3450] Z (centred large-market) mean: 0.000000
§A K=2 Mixture -- No Z¶
Pure mixture with two components. Prior matches the bayesm notebook exactly:
Amu = 0.1(precision for $\mu_k$, so SD $\approx$ 3.16)nu = 6,V = 5I(IW for $\Sigma_k$)a = 5(symmetric Dirichlet)
In [6]:
Prior_K2 = {
'a': 5.,
'mubar': np.zeros(nvar_c),
'Amu': 0.1,
'nu': nvar_c + 3, # = 6
'V': float(nvar_c + 2) * np.eye(nvar_c), # = 5*I
'nu_e': 3.0,
'ssq': ssq_c,
}
Mcmc_K2 = {'R': 25000, 'burn': 5000, 'keep': 1, 'nprint': 5000}
print('Running K=2 mixture (no Z)...')
t0 = time.perf_counter()
out_K2 = rhier_linear_mixture(regdata_c, K=2,
Prior=Prior_K2, Mcmc=Mcmc_K2, seed=11)
print(f'Done in {time.perf_counter()-t0:.1f}s')
Running K=2 mixture (no Z)...
Iter 5000/25000 (31s)
Iter 10000/25000 (62s)
Iter 15000/25000 (97s)
Iter 20000/25000 (132s)
Iter 25000/25000 (168s) Done in 2.8 min | kept: 20000 Done in 167.5s
In [7]:
mu_K2 = out_K2['mudraw'].mean(0) # (2, 3)
pi_K2 = out_K2['pidraw'].mean(0) # (2,)
# bayesm reference (from hlm_mixture_bayesm.ipynb)
# Note: bayesm reports component 1 = higher intercept (9.765), component 2 = lower (8.326)
# Our sampler sorts ascending, so component 0 = lower, component 1 = higher
ref_K2 = {
'mu': np.array([[8.326, -1.876, 1.787],
[9.765, -2.056, 1.144]]),
'pi': np.array([0.317, 0.683]),
}
print('K=2 mixture (no Z) -- posterior means vs bayesm:')
print(f'\n{"":25s} {"Component 0 (low)":>20s} {"Component 1 (high)":>20s}')
print(f'{"Mixing weight":25s} {pi_K2[0]:>20.3f} {pi_K2[1]:>20.3f}')
print(f'{"bayesm weight":25s} {ref_K2["pi"][0]:>20.3f} {ref_K2["pi"][1]:>20.3f}')
print()
for j, name in enumerate(['Intercept', 'log-price', 'Display']):
print(f' {name:10s} Python: {mu_K2[0,j]:6.3f} / {mu_K2[1,j]:6.3f} '
f' bayesm: {ref_K2["mu"][0,j]:6.3f} / {ref_K2["mu"][1,j]:6.3f}')
K=2 mixture (no Z) -- posterior means vs bayesm:
Component 0 (low) Component 1 (high)
Mixing weight 0.511 0.489
bayesm weight 0.317 0.683
Intercept Python: 9.458 / 10.522 bayesm: 8.326 / 9.765
log-price Python: -1.964 / -2.365 bayesm: -1.876 / -2.056
Display Python: 2.454 / 2.799 bayesm: 1.787 / 1.144
Label switching: what is going on, and how to fix it¶
What is going on. A finite-mixture likelihood is invariant to permuting the component labels — relabelling "component 0" and "component 1" leaves $p(y\mid\theta)$ unchanged. When the components are well separated the sampler still effectively stays in one labelling, but here they overlap heavily, so the Gibbs chain wanders freely between the $K! = 2$ equivalent labellings. On top of that a component can briefly empty out ($\pi_k\to0$): with almost no stores assigned, its mean is sampled from the diffuse prior $\mathcal N(0,A_\mu^{-1})$ and spikes to extreme values until it refills.
The first figure (raw chains) shows both symptoms: the two $\mu_k$ traces trade places, the weights $\pi_0,\pi_1$ each sweep all of $[0,1]$ as mirror images, and synchronized spikes appear whenever a component vanishes. This is non-identification, not a convergence failure (mixing is otherwise fast) — and it is itself evidence that there are no well-separated store segments. It also explains why different engines report different per-component means: each just froze a different labelling.
Plausible fixes (roughly in order of effort):
- Post-hoc relabelling by an ordering rule — within each draw, sort components by one coordinate (here $\mu^{\text{int}}$) before summarising. Cheap, and the default bayesm uses; this is what the second figure does. Weakness: a single-coordinate sort can mislabel when components differ on several coefficients.
- Identifiability / ordering constraint in the sampler — impose $\mu^{\text{int}}_0<\mu^{\text{int}}_1$ (reorder or reject violating draws). Stops switching during sampling, but can bias estimates when the true means are close (as they are here).
- Loss-based relabelling algorithms — Stephens (2000) KL-relabelling, the ECR algorithm, or pivotal reordering: relabel each draw to minimise a loss against a reference labelling. More robust than a coordinate sort for multi-dimensional components (R:
label.switching). - Report only label-invariant quantities — the predictive density, WAIC/LOO, the number of occupied components, and per-store membership entropy do not depend on labels. These are what we trust for the substantive conclusion.
- Tame the empty-component spikes — a tighter prior on $\mu_k$ (larger $A_\mu$) keeps an emptied component's mean from running to $\pm10$; or use a deliberately sparse / overfitted-mixture prior (small Dirichlet $a$) that empties spare components and simply count the occupied ones.
- Accept the verdict / reduce $K$ — since no store commits (membership clusters at 0.5), $K=1$ (single normal) is the honest model; the mixture is a flexible-density refinement, not a segmentation device.
In [8]:
# RAW trace plots -- illustrates LABEL SWITCHING (components swap; weights mirror across [0,1]; spikes = empty component)
fig, axes = plt.subplots(2, 3, figsize=(14, 5))
param_labels = ['Intercept', 'log-price', 'Display']
for j, ax in enumerate(axes[0]):
ax.plot(out_K2['mudraw'][:, 0, j], lw=0.5, color='steelblue', label='Comp 0')
ax.plot(out_K2['mudraw'][:, 1, j], lw=0.5, color='firebrick', label='Comp 1')
ax.axhline(ref_K2['mu'][0, j], color='steelblue', ls='--', lw=1.2)
ax.axhline(ref_K2['mu'][1, j], color='firebrick', ls='--', lw=1.2)
ax.set_title(f'mu_k -- {param_labels[j]}'); ax.legend(fontsize=8)
for j, ax in enumerate(axes[1]):
if j < 2:
ax.plot(out_K2['pidraw'][:, j], lw=0.5, color=['steelblue', 'firebrick'][j])
ax.axhline(ref_K2['pi'][j], color='k', ls='--', lw=1.2)
ax.set_title(f'pi[{j}]')
else:
ax.hist(out_K2['sdraw'].mean(0), bins=30, color='orchid', edgecolor='white')
ax.set_title('Store membership probability (mean draw)')
ax.set_xlabel('P(component 1)')
plt.suptitle('K=2 (no Z) -- RAW trace plots: label switching [dashed = bayesm reference]')
plt.tight_layout(); plt.savefig('hlm_mixture_trace_K2_raw.png', dpi=120, bbox_inches='tight'); plt.show()
In [9]:
# Component MEANS relabelled by intercept (order-by-intercept). Only the means change;
# the weights pi and store membership are label-invariant -> see the raw figure above.
mu = out_K2['mudraw']
order = np.argsort(mu[:, :, 0], axis=1) # per draw: order components by intercept
mu_rl = np.take_along_axis(mu, order[:, :, None], axis=1)
ro = np.argsort(ref_K2['mu'][:, 0]); ref_mu = ref_K2['mu'][ro]
fig, axes = plt.subplots(1, 3, figsize=(14, 3.4))
param_labels = ['Intercept', 'log-price', 'Display']
lab = ['Comp 0 (low intercept)', 'Comp 1 (high intercept)']
for j, ax in enumerate(axes):
ax.plot(mu_rl[:, 0, j], lw=0.5, color='steelblue', label=lab[0])
ax.plot(mu_rl[:, 1, j], lw=0.5, color='firebrick', label=lab[1])
ax.axhline(ref_mu[0, j], color='steelblue', ls='--', lw=1.2)
ax.axhline(ref_mu[1, j], color='firebrick', ls='--', lw=1.2)
ax.set_title(f'mu_k -- {param_labels[j]}'); ax.legend(fontsize=8)
plt.suptitle('K=2 (no Z) -- component MEANS relabelled by intercept (now in stable bands) [dashed = bayesm reference]')
plt.tight_layout(); plt.savefig('hlm_mixture_trace_K2.png', dpi=120, bbox_inches='tight'); plt.show()
What the relabelling shows (and its limit). Ordering by intercept cleanly separates the $\mu^{\text{int}}$ traces (red high-intercept component stays above blue), so the mean chains are now readable. But the weight traces $\pi_0,\pi_1$ in the raw figure above still sweep $[0,1]$ (so they are not re-plotted here): because the two components overlap on the intercept as well, which component counts as low-intercept flips on negligible differences and is not tied to which one carries more weight. A one-coordinate sort therefore only partly identifies the model here — a loss-based relabelling (fix #3) would do better — and the residual ambiguity is one more sign that there are no genuinely distinct segments.
In [10]:
# Store-level posterior mean betas
beta_pm_K2 = out_K2['betadraw'].mean(axis=2) # (88, 3)
# Hard component assignment (modal s across draws)
seg_K2 = np.array([np.bincount(out_K2['sdraw'][:, i]).argmax()
for i in range(nreg_c)])
print(f'Store count per component: {np.bincount(seg_K2)}')
# Price elasticity by segment
for k in range(2):
idx = seg_K2 == k
pe = beta_pm_K2[idx, 1]
print(f' Component {k} ({idx.sum()} stores): '
f'mean={pe.mean():.3f} sd={pe.std():.3f}')
# Density overlay: price elasticities by component
fig, ax = plt.subplots(figsize=(8, 4))
xs = np.linspace(-5, 1, 300)
colors = ['steelblue', 'firebrick']
for k in range(2):
vals = beta_pm_K2[seg_K2 == k, 1]
if len(vals) > 2:
kd = gaussian_kde(vals, bw_method='silverman')
ax.plot(xs, kd(xs), color=colors[k], lw=2,
label=f'Component {k} (n={len(vals)})')
ax.axvline(mu_K2[k, 1], color=colors[k], ls='--', lw=1.2)
kd_all = gaussian_kde(beta_pm_K2[:, 1], bw_method='silverman')
ax.plot(xs, kd_all(xs), 'k-', lw=1.5, label='All stores')
ax.set_xlabel('log-price elasticity')
ax.set_title('K=2 mixture -- price elasticities by component')
ax.legend()
plt.tight_layout()
plt.show()
Store count per component: [66 22] Component 0 (66 stores): mean=-1.905 sd=0.625 Component 1 (22 stores): mean=-2.871 sd=0.709
§B K=2 Mixture -- With Centred Z¶
Add the centred large-market dummy as a unit covariate. $\Delta$ captures a systematic mean shift of all $\beta_i$ attributable to market size, layered on top of the mixture structure.
bayesm target (K=2+Z):
| Component 0 (low $\mu$) | Component 1 (high $\mu$) | |
|---|---|---|
| Weight | 0.102 | 0.898 |
| Intercept | 3.310 | 10.188 |
| log-price | -0.799 | -2.088 |
| Display | 0.983 | 0.995 |
Delta (large-market shift): intercept +0.808, log-price -0.030, display +0.280
In [11]:
Prior_K2Z = {
'a': 5.,
'mubar': np.zeros(nvar_c),
'Amu': 0.1,
'nu': nvar_c + 3,
'V': float(nvar_c + 2) * np.eye(nvar_c),
'nu_e': 3.0,
'ssq': ssq_c,
'Ad': 0.01 * np.eye(nvar_c), # (p, p) precision for delta (nz=1)
'deltabar': np.zeros(nvar_c),
}
Mcmc_K2Z = {'R': 25000, 'burn': 5000, 'keep': 1, 'nprint': 5000}
print('Running K=2 mixture (with Z)...')
t0 = time.perf_counter()
out_K2Z = rhier_linear_mixture(regdata_c, K=2, Z=Z_c,
Prior=Prior_K2Z, Mcmc=Mcmc_K2Z, seed=22)
print(f'Done in {time.perf_counter()-t0:.1f}s')
Running K=2 mixture (with Z)...
Iter 5000/25000 (46s)
Iter 10000/25000 (91s)
Iter 15000/25000 (136s)
Iter 20000/25000 (178s)
Iter 25000/25000 (220s) Done in 3.7 min | kept: 20000 Done in 220.3s
In [12]:
mu_K2Z = out_K2Z['mudraw'].mean(0) # (2, 3)
pi_K2Z = out_K2Z['pidraw'].mean(0) # (2,)
delta_pm = out_K2Z['Deltadraw'].mean(0) # (p*nz,) = (3,) for nz=1
ref_K2Z = {
'mu': np.array([[3.310, -0.799, 0.983],
[10.188, -2.088, 0.995]]),
'pi': np.array([0.102, 0.898]),
'delta': np.array([0.808, -0.030, 0.280]),
}
print('K=2+Z -- component means and weights:')
print(f'\n Weights: Python {pi_K2Z.round(3)} bayesm {ref_K2Z["pi"]}')
print()
for j, name in enumerate(['Intercept', 'log-price', 'Display']):
print(f' {name:10s} Python: {mu_K2Z[0,j]:7.3f} / {mu_K2Z[1,j]:7.3f} '
f' bayesm: {ref_K2Z["mu"][0,j]:7.3f} / {ref_K2Z["mu"][1,j]:7.3f}')
print('\nDelta (large-market shift on beta means):')
for j, name in enumerate(['Intercept', 'log-price', 'Display']):
print(f' {name:10s} Python: {delta_pm[j]:+.3f} bayesm: {ref_K2Z["delta"][j]:+.3f}')
K=2+Z -- component means and weights: Weights: Python [0.427 0.573] bayesm [0.102 0.898] Intercept Python: 9.471 / 10.465 bayesm: 3.310 / 10.188 log-price Python: -1.982 / -2.354 bayesm: -0.799 / -2.088 Display Python: 3.012 / 2.317 bayesm: 0.983 / 0.995 Delta (large-market shift on beta means): Intercept Python: +0.735 bayesm: +0.808 log-price Python: +0.009 bayesm: -0.030 Display Python: +0.282 bayesm: +0.280
§C K=3 Mixture -- No Z¶
Try three components. bayesm finds all three non-trivial: weights approximately [0.374, 0.443, 0.182] (reported as sorted by bayesm's ordering).
In [13]:
Prior_K3 = {
'a': 5.,
'mubar': np.zeros(nvar_c),
'Amu': 0.1,
'nu': nvar_c + 3,
'V': float(nvar_c + 2) * np.eye(nvar_c),
'nu_e': 3.0,
'ssq': ssq_c,
}
Mcmc_K3 = {'R': 25000, 'burn': 5000, 'keep': 1, 'nprint': 5000}
print('Running K=3 mixture (no Z)...')
t0 = time.perf_counter()
out_K3 = rhier_linear_mixture(regdata_c, K=3,
Prior=Prior_K3, Mcmc=Mcmc_K3, seed=33)
print(f'Done in {time.perf_counter()-t0:.1f}s')
Running K=3 mixture (no Z)...
Iter 5000/25000 (38s)
Iter 10000/25000 (77s)
Iter 15000/25000 (114s)
Iter 20000/25000 (155s)
Iter 25000/25000 (194s) Done in 3.2 min | kept: 20000 Done in 194.1s
In [14]:
mu_K3 = out_K3['mudraw'].mean(0) # (3, 3)
pi_K3 = out_K3['pidraw'].mean(0) # (3,)
print('K=3 mixture -- posterior means (sorted by intercept):')
print(f'\n Mixing weights: {pi_K3.round(3)}')
print(f' Approx store counts: {(pi_K3 * nreg_c).round(1)}')
print()
for k in range(3):
print(f' Component {k}: '
f'intercept={mu_K3[k,0]:.3f} '
f'log-price={mu_K3[k,1]:.3f} '
f'display={mu_K3[k,2]:.3f}')
if min(pi_K3) < 0.05:
print('\nWarning: near-empty component -- K=3 may be over-specified.')
else:
print('\nAll three components have non-trivial weight.')
K=3 mixture -- posterior means (sorted by intercept): Mixing weights: [0.074 0.494 0.432] Approx store counts: [ 6.5 43.5 38. ] Component 0: intercept=1.368 log-price=-0.248 display=0.446 Component 1: intercept=9.906 log-price=-2.058 display=2.599 Component 2: intercept=10.560 log-price=-2.389 display=2.881 All three components have non-trivial weight.
§D Comparison: Single Normal vs Mixture Models¶
Compare the distribution of store-level price elasticities across
the three estimation approaches: K=1 (HLM from hlm_python.ipynb
if available), K=2, and K=3. Mixture models provide store segmentation
as a by-product.
In [15]:
beta_pm_K3 = out_K3['betadraw'].mean(axis=2) # (88, 3)
ols_c = np.array([
np.linalg.lstsq(d['X'], d['y'], rcond=None)[0]
for d in regdata_c
])
datasets = [
('OLS', ols_c, 'grey'),
('K=2 mixture', beta_pm_K2, 'steelblue'),
('K=3 mixture', beta_pm_K3, 'firebrick'),
]
# Try to load K=1 (HLM) betas if available
if os.path.exists('cheese_hlm_python_betas.csv'):
hlm_py = pd.read_csv('cheese_hlm_python_betas.csv')
b_hlm = hlm_py[hlm_py['model']=='py_A'].sort_values('retailer')
beta_hlm = b_hlm[['intercept','log_price','disp']].values
datasets.insert(1, ('HLM (K=1)', beta_hlm, 'darkorange'))
print('Cross-store SD of log-price elasticity:')
for label, betas, _ in datasets:
print(f' {label:15s}: {betas[:, 1].std():.3f}')
fig, axes = plt.subplots(1, len(datasets), figsize=(4*len(datasets), 4), sharey=False)
xs = np.linspace(-5, 1, 300)
breaks = np.arange(-5, 1.2, 0.4)
for ax, (label, betas, color) in zip(axes, datasets):
ax.hist(betas[:, 1], bins=breaks, color=color, alpha=0.7, edgecolor='white')
ax.set_title(label)
ax.set_xlabel('log-price elasticity')
ax.set_xlim(-5, 1)
plt.suptitle('Distribution of store-level price elasticities')
plt.tight_layout()
plt.show()
Cross-store SD of log-price elasticity: OLS : 1.765 HLM (K=1) : 0.764 K=2 mixture : 0.770 K=3 mixture : 0.772
Save Results¶
In [16]:
# K=2 store betas + segment
df_K2 = pd.DataFrame({
'retailer': retailers,
'large_market': large_raw,
'intercept': beta_pm_K2[:, 0],
'log_price': beta_pm_K2[:, 1],
'disp': beta_pm_K2[:, 2],
'component': seg_K2,
'model': 'py_K2',
})
df_K2.to_csv('cheese_hlm_mixture_K2_python.csv', index=False)
print('Saved: cheese_hlm_mixture_K2_python.csv')
# Component summary
comp_K2 = pd.DataFrame({
'component': [0, 1],
'weight': pi_K2.round(3),
'mu_int': mu_K2[:, 0].round(3),
'mu_price': mu_K2[:, 1].round(3),
'mu_disp': mu_K2[:, 2].round(3),
})
comp_K2.to_csv('cheese_mixture_K2_components_python.csv', index=False)
print('Saved: cheese_mixture_K2_components_python.csv')
print()
print('K=2 component summary:')
print(comp_K2.to_string(index=False))
Saved: cheese_hlm_mixture_K2_python.csv
Saved: cheese_mixture_K2_components_python.csv
K=2 component summary:
component weight mu_int mu_price mu_disp
0 0.5110 9.4580 -1.9640 2.4540
1 0.4890 10.5220 -2.3650 2.7990
Model selection and segment membership¶
In [17]:
# Model selection (WAIC vs K) + segment-membership uncertainty
XtX = np.stack([d['X'].T @ d['X'] for d in regdata_c]); Xty = np.stack([d['X'].T @ d['y'] for d in regdata_c])
yty = np.array([d['y'] @ d['y'] for d in regdata_c]); Tn = np.array([len(d['y']) for d in regdata_c], float)
p = regdata_c[0]['X'].shape[1]
def waic(o, ss=10):
mud, Sigd, pid, s2d = o['mudraw'], o['Sigmadraw'], o['pidraw'], o['sigmasqdraw']
Rk, K = mud.shape[0], mud.shape[1]; idx = np.arange(0, Rk, ss); LL = np.empty((len(idx), nreg_c))
for ri, r in enumerate(idx):
s2 = s2d[r]; llk = np.empty((nreg_c, K))
for k in range(K):
Sk = Sigd[r,k]; Si = np.linalg.inv(Sk); mu = np.tile(mud[r,k], (nreg_c,1))
Xte = Xty - np.einsum('ipq,iq->ip', XtX, mu)
ete = yty - 2*np.sum(mu*Xty,1) + np.einsum('ip,ipq,iq->i', mu, XtX, mu)
M = np.linalg.inv(Si[None] + XtX/s2[:,None,None]); quad = ete/s2 - np.einsum('ip,ipq,iq->i', Xte, M, Xte)/s2**2
B = np.eye(p)[None] + np.einsum('pq,iqr->ipr', Sk, XtX)/s2[:,None,None]
llk[:,k] = np.log(pid[r,k]+1e-300) - 0.5*(Tn*np.log(2*np.pi) + Tn*np.log(s2) + np.linalg.slogdet(B)[1] + quad)
m = llk.max(1); LL[ri] = m + np.log(np.exp(llk-m[:,None]).sum(1))
mm = LL.max(0); return -2*((np.log(np.exp(LL-mm).mean(0))+mm).sum() - LL.var(0).sum())
out_K1 = rhier_linear_mixture(regdata_c, K=1, Prior=Prior_K2, Mcmc=dict(R=10000, burn=2000, keep=1, nprint=0), seed=11)
waics = {1: waic(out_K1), 2: waic(out_K2), 3: waic(out_K3)}
print('WAIC:', {k: round(v,1) for k,v in waics.items()})
# membership probabilities (K=2): fraction of draws each store is assigned to each component
S = out_K2['sdraw']; memb = np.stack([np.bincount(S[:,i], minlength=2)/S.shape[0] for i in range(nreg_c)])
maxp = memb.max(1)
fig, ax = plt.subplots(1, 2, figsize=(13, 4.4))
ks = sorted(waics); ax[0].plot(ks, [waics[k] for k in ks], 'o-', color='steelblue', lw=2)
best = min(waics, key=waics.get); ax[0].scatter([best],[waics[best]], color='firebrick', s=90, zorder=5, label=f'best: K={best}')
ax[0].set_xticks(ks); ax[0].set_xlabel('number of components K'); ax[0].set_ylabel('WAIC (lower = better)')
ax[0].set_title('Model selection: K=2 mildly best, then worse'); ax[0].legend()
ax[1].hist(maxp, bins=20, color='steelblue', edgecolor='white')
ax[1].axvline(0.9, color='firebrick', ls='--', label='strong commitment (0.9)')
ax[1].set_xlabel('max posterior membership probability per store'); ax[1].set_ylabel('# stores'); ax[1].set_xlim(0.5, 1.0)
ax[1].set_title('Segment membership (K=2): stores cluster near 0.5\n-> components overlap, no clean segments'); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig('hlm_mixture_select.png', dpi=120, bbox_inches='tight'); plt.show()
Done in 0.9 min | kept: 8000
WAIC: {1: np.float64(-838.9), 2: np.float64(-843.3), 3: np.float64(-840.0)}
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; see the model-selection figure). 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.5) 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.