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 — the probability the model assigns to the data after averaging over every value its parameters could take, so added flexibility is charged for automatically. It is estimated here by the Newton–Raftery harmonic mean, which reweights ordinary posterior draws into an estimate of that quantity. That estimator is bayesm's metric and is notoriously unstable — a single draw in the low-likelihood tail can dominate the whole average — so the verdict here is corroborated by the label-invariant fitted density rather than rested on it.

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