Hierarchical MNL with mixture-of-normals heterogeneity¶
Rossi, Allenby & McCulloch (Bayesian Statistics and Marketing, Ch. 5) — the mixture as a flexible taste distribution¶
References: Rossi, Allenby & McCulloch (2005), BSM Ch. 5 (rhierMnlRwMixture); Allenby & Rossi (1998).
| Section | Content |
|---|---|
| Model | Random-coefficients MNL whose coefficient distribution is a mixture of normals |
| Section 1 | Synthetic validation — recover a known bimodal taste distribution |
| Section 2 | bayesm margarine — the population of tastes is non-normal (fat-tailed), à la Rossi et al. |
From-scratch sampler: hmnl_mixture_gibbs.py. Companion to hier_mnl_gibbs.py (single-normal) and the linear hlm_mixture_gibbs.py.
Model — a flexible distribution of tastes¶
$$P(\text{choose }j\mid \text{task}) = \frac{e^{x_j'\beta_i}}{\sum_l e^{x_l'\beta_i}}, \qquad \beta_i \mid s_i{=}k \sim N(\mu_k + \Delta'z_i,\ \Sigma_k),\quad s_i\sim\text{Cat}(\pi).$$
Each respondent is a logit with their own part-worths $\beta_i$; the population of $\beta_i$ is a mixture of normals. The single-normal hierarchy is the special case $K=1$.
What the mixture is for (the Rossi et al. point). It is a flexible prior for the heterogeneity distribution — it lets the population of tastes be skewed or fat-tailed, which a single normal cannot. This is the same "flexible distribution" role as Geweke–Keane's mixture-of-normals error, applied here to the random-coefficient density instead of the disturbance. It is not primarily a device to find discrete segments. So there are two different questions:
- "Is the taste distribution non-normal?" — the BSM question, answered by the fitted density and the log marginal density.
- "Are there clean discrete segments?" — a stronger claim (stable cluster membership), usually not supported in real data.
Sampler (RW-Metropolis within Gibbs): (1) $\beta_i$ — random-walk Metropolis, proposal $N(\beta_i,(H_i+\Sigma_{s_i}^{-1})^{-1})$ ($H_i$ = MNL information at the pooled MLE); (2) $s_i$ — categorical; (3) $\mu_k,\Sigma_k$ — Normal/Inverse-Wishart; (4) $\pi$ — Dirichlet; (5) $\Delta$ — GLS [if Z].
Model comparison: the log marginal density (Newton–Raftery harmonic mean — bayesm's metric; known to be unstable, so we corroborate with the fitted heterogeneity density, which is label-invariant).
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy.stats import norm
from hmnl_mixture_gibbs import (simulate_hmnl_mixture, hmnl_mixture_gibbs,
hetero_density, log_marg_den_nr)
print('hmnl_mixture_gibbs loaded')
hmnl_mixture_gibbs loaded
1. Synthetic validation (bimodal tastes)¶
500 respondents, 15 tasks, 3 alternatives, $k=3$ (2 brand dummies + price). Price sensitivity is bimodal: 55% near $-1$ (insensitive) and 45% near $-4$ (sensitive). A single normal must average them into one bell; the mixture should recover both modes.
sim = simulate_hmnl_mixture(N=500, T=15, J=3, k=3, seed=1)
print('truth: pi', sim['pi'], ' price means -1 / -4')
sfit = {K: hmnl_mixture_gibbs(sim['data'], K=K, R=3000, burn=1000, seed=2) for K in (1, 2)}
for K, o in sfit.items():
extra = '' if K == 1 else f" pi={o['pi'].mean(0).round(2)} price means={o['mu'].mean(0)[:,2].round(2)}"
print(f"K={K}: logMargDen={log_marg_den_nr(o['loglike']):.1f} accept={o['accept']:.2f}{extra}")
truth: pi [0.55 0.45] price means -1 / -4 K=1: logMargDen=-4657.0 accept=0.35 K=2: logMargDen=-4645.6 accept=0.34 pi=[0.45 0.55] price means=[-4.01 -0.95]
grid = np.linspace(sim['B_true'][:,2].min()-1, sim['B_true'][:,2].max()+1, 300)
plt.figure(figsize=(8,4.6))
plt.hist(sim['B_true'][:,2], bins=40, density=True, color='lightgray', label='true price coefs')
plt.plot(grid, hetero_density(sfit[1],2,grid), 'b--', lw=2, label='single normal (K=1)')
plt.plot(grid, hetero_density(sfit[2],2,grid), 'r-', lw=2.2, label='mixture (K=2)')
plt.xlabel('price sensitivity'); plt.ylabel('density')
plt.title('Synthetic: mixture recovers the bimodal taste distribution')
plt.legend(); plt.tight_layout()
plt.savefig('hmnl_mixture_synth.png', dpi=120, bbox_inches='tight'); plt.show()
1. Results — synthetic¶
The K=2 mixture recovers the truth almost exactly — weights ≈ 0.45 / 0.55 and price means ≈ −4.0 / −1.0 (vs. truth −4/−1 at 0.45/0.55), and the fitted density is clearly bimodal where the single normal is forced into one bell. Log marginal density improves K=1→K=2. The sampler can find genuine structure when it exists — so any "no segments" finding on real data is about the data, not the method.
2. margarine: is the population of tastes non-normal?¶
The bayesm margarine household panel (516 households, ~8.7 purchases each, 10 products). Design = 9 brand intercepts (base = Parkay stick) + price ($k=10$). We fit $K=1,2,3$ and ask the BSM question: does a mixture reveal a non-normal taste distribution?
cp = pd.read_csv('margarine_choice.csv')
pricecols = ['PPk_Stk','PBB_Stk','PFl_Stk','PHse_Stk','PGen_Stk','PImp_Stk','PSS_Tub','PPk_Tub','PFl_Tub','PHse_Tub']
J, k = 10, 10
data = []
for h in cp.hhid.unique():
g = cp[cp.hhid == h]; T = len(g); pr = g[pricecols].values
X = np.zeros((T, J, k)); yv = g.choice.values.astype(int) - 1
for t in range(T):
for j in range(J):
if j >= 1: X[t, j, j-1] = 1.0 # 9 brand intercepts (alt 1 = Parkay stick = base)
X[t, j, 9] = pr[t, j] # own price
data.append((X, yv))
print(f'{len(data)} households, avg {np.mean([len(y) for _,y in data]):.1f} purchases each')
516 households, avg 8.7 purchases each
mfit = {K: hmnl_mixture_gibbs(data, K=K, R=4000, burn=1500, seed=3) for K in (1, 2, 3)}
print(f"{'K':>2}{'logMargDen':>12}{'accept':>8} price means per component")
for K, o in mfit.items():
print(f"{K:>2}{log_marg_den_nr(o['loglike']):>12.0f}{o['accept']:>8.2f} {o['mu'].mean(0)[:,9].round(2)}")
print(f"\nbest by log marginal density: K={max(mfit, key=lambda K: log_marg_den_nr(mfit[K]['loglike']))}")
K logMargDen accept price means per component 1 -4097 0.30 [-8.88] 2 -4053 0.30 [-9.08 -8.59] 3 -4098 0.29 [-9.83 -9.31 -8.91] best by log marginal density: K=2
# fitted heterogeneity density: single normal (K=1) vs mixture (K=2), for PRICE and a brand intercept
labels = ['int_BB','int_Fl','int_HseStk','int_Gen','int_Imp','int_SSTub','int_PkTub','int_FlTub','int_HseTub','PRICE']
fig, ax = plt.subplots(1, 2, figsize=(13, 4.6))
for axi, coord, name in [(ax[0], 9, 'price sensitivity'), (ax[1], 6, 'Parkay-tub part-worth')]:
bm = mfit[2]['B'][:, coord]
g = np.linspace(bm.min()-2, bm.max()+2, 300)
axi.hist(bm, bins=40, density=True, color='lightgray', label='household posterior means')
axi.plot(g, hetero_density(mfit[1], coord, g), 'b--', lw=2, label='single normal (K=1)')
axi.plot(g, hetero_density(mfit[2], coord, g), 'r-', lw=2.2, label='mixture (K=2)')
axi.set_xlabel(name); axi.set_ylabel('density'); axi.set_title(f'Heterogeneity: {name}'); axi.legend()
plt.tight_layout()
plt.savefig('hmnl_mixture_margarine.png', dpi=120, bbox_inches='tight'); plt.show()
# non-normality of the fitted mixture (posterior-mean skew / excess kurtosis per coordinate)
o = mfit[2]; mu, Sig, pi = o['mu'], o['Sigma'], o['pi']
print('K=2 fitted taste distribution -- non-normality by coordinate:')
print(f" {'coef':<11}{'skew':>8}{'exkurt':>9}")
for j in range(k):
m = (pi*mu[:,:,j]).sum(1); v = (pi*(Sig[:,:,j,j]+mu[:,:,j]**2)).sum(1)-m**2
m3 = (pi*(mu[:,:,j]**3+3*mu[:,:,j]*Sig[:,:,j,j])).sum(1)-3*m*v-m**3
m4 = (pi*(mu[:,:,j]**4+6*mu[:,:,j]**2*Sig[:,:,j,j]+3*Sig[:,:,j,j]**2)).sum(1)-4*m*m3-6*m**2*v-m**4
print(f" {labels[j]:<11}{np.mean(m3/v**1.5):>8.2f}{np.mean(m4/v**2-3):>9.2f}")
K=2 fitted taste distribution -- non-normality by coordinate: coef skew exkurt int_BB -0.08 0.24 int_Fl 0.09 0.47 int_HseStk -0.09 0.58 int_Gen -0.12 0.48 int_Imp -0.28 0.28 int_SSTub -0.26 0.90 int_PkTub -0.41 0.87 int_FlTub -0.04 0.14 int_HseTub -0.10 0.12 PRICE 0.08 0.51
2. Results — margarine¶
| K | log marg. density | price means |
|---|---|---|
| 1 | −4097 | [−8.88] |
| 2 | −4053 | [−9.08, −8.59] |
| 3 | −4098 | [−9.83, −9.31, −8.91] |
A 2-component mixture is favored (log marginal density up ≈ +44 over the single normal; K=3 falls back) — and the fitted taste distribution is non-normal: leptokurtic (excess kurtosis ≈ +0.5 for price, +0.9 for the Parkay-tub and stick-SS part-worths) and mildly left-skewed (skew ≈ −0.4 for Parkay-tub). The population is a tight core of typical households plus a heavier tail of atypical ones than a normal allows — exactly Rossi et al.'s point that a single normal is too restrictive for the heterogeneity distribution.
But note what it is not: the two components' price means are close (−9.1 vs −8.6) and overlap heavily — this is fat tails / skew, not clean discrete segments. So the mixture earns its keep here as a flexible distribution, not as a segmentation device — consistent with the linear-model finding (cheese/Lenk) that discrete segments don't survive, while flexible-density non-normality does.
The unifying theme. Mixtures of normals pay off in their flexible-distribution role — on the error (Geweke–Keane: earnings, beetle) and on the heterogeneity density (Rossi et al.: margarine) — far more reliably than as a clustering/segmentation device.
Does the mixture change the interpretation vs. the single-normal hierarchical MNL?¶
Compared with the single-normal hier-MNL on the same margarine data (the hier-MNL project), the answer is no substantive change — this is a robustness refinement, not a revision.
What stays the same (the economic story is robust):
| Quantity | Single-normal hier-MNL | Mixture (here) |
|---|---|---|
| Mean price coefficient | ≈ −8.9 [−9.4, −8.4] | ≈ −8.9 (essentially identical) |
| Brand preference ordering | Parkay stick/tub & Fleischmann top; generics & house-tub ~−5 below base | unchanged |
| Willingness-to-pay (brand part-worth / |price|) | generics need ~$0.56–0.59 discount, etc. | unchanged (a ratio of means) |
| IIA / substitution | relaxed across people by the hierarchy; MNP contrast | unchanged |
(The mixture's mean price coefficient is the π-weighted average across components: −9.08 at weight 0.62 and −8.59 at weight 0.38 give −8.89 — not the −9.1 of the first component alone. It lands essentially on top of the single-normal's −8.90, which is exactly the point: the mixture reshapes the *tails of the taste distribution, not its centre.)*
The population-mean part-worths, price sensitivity, WTP, and the IIA/substitution conclusions are all stable when we swap the normal heterogeneity for a mixture — the earlier hier-MNL results are not artifacts of the normality assumption.
What the mixture adds (a refinement of the heterogeneity shape): the taste distribution is non-normal — leptokurtic / fat-tailed (excess kurtosis ≈ 0.5 on price, ≈ 0.9–1.5 on some brand part-worths; mild skew). Consequences:
- a heavier tail of atypical households (very price-sensitive or strongly brand-loyal) than a single normal allows — the normal understates the extremes;
- individual-level $\beta_i$ for outlier households are over-shrunk toward the global mean under the normal; the fat-tailed prior shrinks them less → modestly better individual prediction/targeting at the tails and a higher marginal likelihood;
- but aggregate elasticity, average WTP, and brand equity are essentially unchanged because the mean is stable.
Two cautions against over-reading:
- No targetable segments — the two components overlap heavily (price means ≈ −9.1/−8.6); this is fat tails, not two customer types.
- This run used no demographics (Z) — the hier-MNL's demo effects (income → less price-sensitive +0.72; family size → more −0.70) are a separate observed-heterogeneity story that still stands; a mixture would add fat tails on top of those Δ effects, not replace them.
Bottom line: the mixture is the robustness check Rossi et al. advocate — it confirms the single-normal hier-MNL's headline numbers and reveals that the residual heterogeneity is fat-tailed, which matters for tail and individual-level inference but leaves mean part-worths, WTP, and IIA conclusions intact.