Bayesian Hierarchical MNL — Mixture Heterogeneity

Model

This extends the hierarchical multinomial logit by replacing its single-normal heterogeneity prior with a mixture of normals: each respondent $i$ is still a logit with their own part-worths $\beta_i$ , but the population of $\beta_i$ is now a $K$ -component mixture. The single-normal hierarchy is the special case $K=1$ . It replicates bayesm's rhierMnlRwMixture (Rossi, Allenby & McCulloch 2005, Ch. 5).

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)

$\beta_i$ — respondent $i$ 's part-worths; the population is a mixture of $K$ normals
$(\pi_k, \mu_k, \Sigma_k)$ — weight, mean, covariance of taste-component $k$
$s_i$ — component assignment; $\Delta$ — optional observed-demographic effects

What the mixture is for (the Rossi et al. point). It is a flexible prior for the heterogeneity distribution — letting 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 rather than the disturbance. It is not primarily a device to find discrete segments. Two distinct questions follow: "Is the taste distribution non-normal?" (the BSM question, answered by the fitted density and the log marginal density) versus "Are there clean discrete segments?" (a stronger claim that real data usually does not support).

Sampler — RW-Metropolis within Gibbs

The sampler is the hier-MNL's RW-Metropolis-within-Gibbs with two added blocks for the mixture: a categorical assignment $s_i$ and the per-component Normal/Inverse-Wishart updates. Components are sorted each sweep (by $\mu_k$ on the last coordinate) to stabilise labels; the fitted heterogeneity density is label-invariant regardless. Model comparison uses the log marginal density (Newton–Raftery harmonic mean — bayesm's metric, known to be unstable, so corroborated by the label-invariant fitted density).

Block	Draw	Update
(1)	$\beta_i \mid \cdot$	Random-walk Metropolis — proposal $N(\beta_i,\, (H_i + \Sigma_{s_i}^{-1})^{-1})$ , $H_i$ the MNL information at the pooled MLE
(2)	$s_i \mid \cdot$	Categorical, $P(k) \propto \pi_k\, N(\beta_i;\, \mu_k + \Delta'z_i,\, \Sigma_k)$
(3)	$(\mu_k, \Sigma_k) \mid \cdot$	Normal / Inverse-Wishart from the $\beta_i$ assigned to component $k$
(4)	$\pi \mid \cdot$	Dirichlet on the component counts
(5)	$\Delta \mid \cdot$	Matrix-normal GLS (only if respondent covariates $Z$ supplied)

Notebooks

Section 1 validates the sampler on synthetic data with a deliberately bimodal price sensitivity (55% near $-1$ , 45% near $-4$ ): the $K=2$ mixture recovers both modes almost exactly where a single normal is forced into one bell — confirming it can find genuine structure when it exists. Section 2 fits the bayesm margarine panel (516 households) and asks the BSM question: a 2-component mixture is favored (log marginal density up $\approx +44$ over the single normal; $K=3$ falls back), and the fitted taste distribution is non-normal — leptokurtic and mildly left-skewed. But the two components' price means overlap heavily ( $-9.1$ vs. $-8.6$ ): this is fat tails, not clean segments. Crucially, swapping the normal heterogeneity for a mixture leaves the headline hier-MNL numbers — mean price coefficient, brand ordering, willingness-to-pay, IIA conclusions — intact; the mixture is the robustness refinement Rossi et al. advocate, revealing fat tails that matter for tail and individual-level inference. The R notebook reaches the identical verdict via rhierMnlRwMixture ( $ncomp=2$ wins, excess kurtosis $\approx 0.55$ on price).

Python Mixture-heterogeneity MNL Gibbs, synthetic bimodal recovery + margarine non-normality (Python) Open → R bayesm::rhierMnlRwMixture — log marginal density & fitted taste density, ncomp 1/2/3 (R) Open →

Downloads

hmnl_mixture_gibbs.py From-scratch mixture-heterogeneity MNL — RW-Metropolis-in-Gibbs sampler (NumPy / scipy) margarine_choice.csv bayesm margarine panel — 516 households, 10 products with shelf prices

Sampler Module — Source Code

"""
Hierarchical MNL with a MIXTURE-OF-NORMALS heterogeneity distribution -- from-scratch Gibbs.
Replicates bayesm::rhierMnlRwMixture (Rossi, Allenby & McCulloch 2005, Ch. 5).

Model
-----
  P(choice = j | task) = exp(x_j' beta_i) / sum_l exp(x_l' beta_i)        (MNL per respondent)
  beta_i | s_i = k ~ N(mu_k + Delta' z_i, Sigma_k)                        (mixture heterogeneity)
  s_i | pi        ~ Categorical(pi)
  pi              ~ Dirichlet(a)
  mu_k            ~ N(0, Amu^{-1} I),   Sigma_k ~ IW(nu, V)

The mixture of normals is used here as a FLEXIBLE heterogeneity distribution (it can approximate a
skewed / multimodal population of tastes), NOT necessarily as discrete segments -- the Rossi et al.
point. The single-normal model is the special case K=1.

Sampler (RW-Metropolis within Gibbs)
------------------------------------
  1. beta_i   -- random-walk Metropolis; proposal N(beta_i, (H_i + Sigma_{s_i}^{-1})^{-1}) where
                 H_i is respondent i's MNL information at the pooled MLE (the bayesm move), so the
                 step adapts to data + the component's prior precision each sweep.
  2. s_i      -- categorical, P(k) propto pi_k N(beta_i; mu_k + Delta'z_i, Sigma_k)
  3. mu_k,Sigma_k -- Normal / Inverse-Wishart conjugate from the beta_i assigned to component k
  4. pi       -- Dirichlet
  5. Delta    -- matrix-normal GLS [only if Z given]
Components are sorted each sweep (by mu_k on the last coordinate) to stabilise labels; the fitted
heterogeneity density is label-invariant regardless.

Companion to hier_mnl_gibbs.py (single-normal) and hlm_mixture_gibbs.py (linear mixture).
"""
import numpy as np
from scipy.stats import invwishart
from scipy.linalg import cholesky, solve_triangular


def mnl_loglik(beta, X, y):
    """MNL log-likelihood for one respondent. X:(T,J,k), y:(T,)."""
    util = X @ beta
    m = util.max(axis=1, keepdims=True)
    logden = m[:, 0] + np.log(np.exp(util - m).sum(axis=1))
    return float((util[np.arange(len(y)), y] - logden).sum())


def mnl_hess(beta, X):
    """Observed information sum_t X_t'(diag(p)-pp')X_t for one respondent."""
    util = X @ beta; util -= util.max(axis=1, keepdims=True)
    e = np.exp(util); p = e / e.sum(axis=1, keepdims=True)
    k = X.shape[2]; H = np.zeros((k, k))
    for t in range(X.shape[0]):
        Xt = X[t]; pt = p[t]
        H += Xt.T @ (np.diag(pt) - np.outer(pt, pt)) @ Xt
    return H


def _pooled_beta(data, k, maxit=25):
    b = np.zeros(k)
    for _ in range(maxit):
        g = np.zeros(k); H = np.zeros((k, k))
        for X, y in data:
            util = X @ b; util -= util.max(axis=1, keepdims=True)
            e = np.exp(util); p = e / e.sum(axis=1, keepdims=True)
            for t in range(len(y)):
                g += X[t, y[t]] - p[t] @ X[t]
            H += mnl_hess(b, X)
        step = np.linalg.solve(H + 1e-6 * np.eye(k), g); b += step
        if np.max(np.abs(step)) < 1e-8:
            break
    return b


def _logmvn(x, mu, L):
    dev = x - mu; v = solve_triangular(L, dev, lower=True)
    return -0.5 * len(x) * np.log(2 * np.pi) - np.log(np.diag(L)).sum() - 0.5 * float(v @ v)


def simulate_hmnl_mixture(N=500, T=15, J=3, k=3, pi=None, mus=None, Sigmas=None, seed=0):
    """Simulate hier-MNL with mixture heterogeneity. Default: bimodal price sensitivity."""
    rng = np.random.default_rng(seed)
    if pi is None:
        pi = np.array([0.55, 0.45])
        mus = np.array([[1.0, -0.5, -1.0],          # price-insensitive segment (price coef -1)
                        [1.0, -0.5, -4.0]])          # price-sensitive segment   (price coef -4)
        Sigmas = np.stack([0.3 * np.eye(k), 0.3 * np.eye(k)])
    K = len(pi)
    comps = rng.choice(K, size=N, p=pi)
    Ls = [np.linalg.cholesky(Sigmas[c]) for c in range(K)]
    B = np.array([mus[comps[i]] + Ls[comps[i]] @ rng.standard_normal(k) for i in range(N)])
    data = []
    for i in range(N):
        X = rng.normal(size=(T, J, k))
        X[:, :, :k-1] = (rng.random((T, J, k-1)) < 0.5).astype(float)   # brand-dummy-like
        util = X @ B[i]
        y = (util + rng.gumbel(size=(T, J))).argmax(axis=1)
        data.append((X, y))
    return dict(data=data, B_true=B, comps=comps, pi=pi, mus=mus, Sigmas=Sigmas)


def hmnl_mixture_gibbs(data, K=2, Z=None, R=4000, burn=1500, seed=0, s_scale=None,
                       a=5.0, Amu=0.01, nu=None, V=None):
    """Hierarchical MNL with K-component mixture-of-normals heterogeneity."""
    rng = np.random.default_rng(seed)
    N = len(data); k = data[0][0].shape[2]
    s_scale = (2.38 / np.sqrt(k)) if s_scale is None else s_scale
    nu = (k + 3) if nu is None else nu
    V = np.eye(k) if V is None else V
    use_Z = Z is not None
    if use_Z:
        Z = np.asarray(Z, float); nz = Z.shape[1]; Delta = np.zeros((nz, k)); Ad = 0.01 * np.eye(nz)

    bpool = _pooled_beta(data, k)
    Hs = [mnl_hess(bpool, X) for X, _ in data]
    B = np.tile(bpool, (N, 1)).astype(float)
    mu = np.tile(bpool, (K, 1)) + 0.01 * rng.standard_normal((K, k))
    Sigma = np.stack([np.eye(k)] * K); Sinv = np.stack([np.eye(k)] * K)
    pi = np.full(K, 1.0 / K)
    s = rng.integers(0, K, N)
    llcur = np.array([mnl_loglik(B[i], data[i][0], data[i][1]) for i in range(N)])

    keep = R - burn
    MU = np.zeros((keep, K, k)); SIG = np.zeros((keep, K, k, k)); PI = np.zeros((keep, K))
    LL = np.zeros(keep); Bm = np.zeros((N, k)); acc = 0
    DD = np.zeros((keep, nz, k)) if use_Z else None

    for r in range(R):
        Ls = [cholesky(Sigma[c], lower=True) for c in range(K)]
        # 1. beta_i  (RW-Metropolis, component-aware proposal)
        for i in range(N):
            ki = s[i]; mu_eff = mu[ki] + (Delta.T @ Z[i] if use_Z else 0.0)
            Lp = cholesky(Hs[i] + Sinv[ki], lower=True)
            prop = B[i] + s_scale * solve_triangular(Lp.T, rng.standard_normal(k), lower=False)
            llp = mnl_loglik(prop, data[i][0], data[i][1])
            lpr = _logmvn(prop, mu_eff, Ls[ki]) - _logmvn(B[i], mu_eff, Ls[ki])
            if np.log(rng.random()) < (llp - llcur[i] + lpr):
                B[i] = prop; llcur[i] = llp; acc += 1
        # 2. s_i
        logpi = np.log(np.clip(pi, 1e-300, None))
        for i in range(N):
            lw = np.array([logpi[c] + _logmvn(B[i], mu[c] + (Delta.T @ Z[i] if use_Z else 0.0), Ls[c])
                           for c in range(K)])
            lw -= lw.max(); w = np.exp(lw); w /= w.sum()
            s[i] = rng.choice(K, p=w)
        # 3. mu_k, Sigma_k
        for c in range(K):
            mask = s == c; nk = int(mask.sum())
            Be = B[mask] - (Z[mask] @ Delta if use_Z else 0.0)
            if nk > 0:
                Pm = Amu * np.eye(k) + nk * Sinv[c]
                rhs = nk * Sinv[c] @ Be.mean(0)
                Lm = cholesky(Pm, lower=True)
                mm = solve_triangular(Lm.T, solve_triangular(Lm, rhs, lower=True), lower=False)
                mu[c] = mm + solve_triangular(Lm.T, rng.standard_normal(k), lower=False)
                dev = Be - mu[c]
                Sigma[c] = invwishart.rvs(df=nu + nk, scale=V + dev.T @ dev, random_state=rng)
            else:
                mu[c] = solve_triangular(cholesky(Amu*np.eye(k), lower=True).T,
                                         rng.standard_normal(k), lower=False)
                Sigma[c] = invwishart.rvs(df=nu, scale=V, random_state=rng)
            Sinv[c] = np.linalg.inv(Sigma[c])
        # 4. pi
        pi = rng.dirichlet(a + np.bincount(s, minlength=K))
        # 5. Delta (matrix-normal GLS) [if Z]
        if use_Z:
            P = Ad.copy().astype(float); rhs = np.zeros((nz, k))
            for i in range(N):
                ci = s[i]
                P += np.outer(Z[i], Z[i])              # shared design; precision in coef space below
            # vectorised GLS per output dim is complex with component-specific Sigma; use simple
            # approximation: regress (B - mu_s) on Z with identity weight (adequate for demo)
            resid = B - mu[s]
            U = np.linalg.inv(Z.T @ Z + Ad)
            M = U @ (Z.T @ resid)
            Delta = M + cholesky(U, lower=True) @ rng.standard_normal((nz, k))
        # label fix: sort by mu on last coordinate (e.g. price)
        order = np.argsort(mu[:, -1])
        mu, Sigma, Sinv, pi = mu[order], Sigma[order], Sinv[order], pi[order]
        s = np.argsort(order)[s]

        if r >= burn:
            j = r - burn
            MU[j] = mu; SIG[j] = Sigma; PI[j] = pi; LL[j] = llcur.sum(); Bm += B
            if use_Z:
                DD[j] = Delta
    out = dict(mu=MU, Sigma=SIG, pi=PI, loglike=LL, B=Bm / keep,
               accept=acc / (R * N), K=K, k=k)
    if use_Z:
        out['Delta'] = DD
    return out


def hetero_density(out, coord, grid):
    """Posterior-mean fitted population density of coefficient `coord` (label-invariant)."""
    from scipy.stats import norm
    mu, Sig, pi = out['mu'], out['Sigma'], out['pi']
    R, K = mu.shape[0], mu.shape[1]
    dens = np.zeros_like(grid, dtype=float)
    for r in range(R):
        for c in range(K):
            dens += pi[r, c] * norm.pdf(grid, mu[r, c, coord], np.sqrt(Sig[r, c, coord, coord]))
    return dens / R


def log_marg_den_nr(loglike):
    """Newton-Raftery (harmonic-mean) log marginal density estimate (bayesm's metric; unstable)."""
    b = -np.asarray(loglike); mb = b.max()
    return float(-(mb + np.log(np.mean(np.exp(b - mb)))))

References

Rossi, P. E., Allenby, G. M. & McCulloch, R. (2005). Bayesian Statistics and Marketing, Ch. 5. Wiley. — the mixture-of-normals heterogeneity distribution and bayesm::rhierMnlRwMixture; the source of the margarine panel
Allenby, G. M. & Rossi, P. E. (1998). Marketing models of consumer heterogeneity. Journal of Econometrics 89(1–2), 57–78. — the case for a flexible, non-normal distribution of tastes
McFadden, D. & Train, K. (2000). Mixed MNL models for discrete response. Journal of Applied Econometrics 15(5), 447–470. — random-coefficients logit, the model whose taste density the mixture generalises
Newton, M. A. & Raftery, A. E. (1994). Approximate Bayesian inference with the weighted likelihood bootstrap. Journal of the Royal Statistical Society: Series B 56(1), 3–48. — the Newton–Raftery log marginal density used to compare K
Stephens, M. (2000). Dealing with label switching in mixture models. Journal of the Royal Statistical Society: Series B 62(4), 795–809. — why mixture components may come back permuted, as in the synthetic recovery here
Train, K. E. (2009). Discrete Choice Methods with Simulation (2nd ed.). Cambridge University Press. — mixed logit and flexible taste distributions