Multinomial Probit

Definition

The multinomial probit (MNP) model extends binary probit to $m \geq 3$ unordered alternatives via a random utility framework. Each consumer $h$ at time $t$ assigns latent utility $y_{h,t,j}$ to alternative $j$ and selects the alternative with highest utility; the choice indicators are generated from an underlying multivariate Normal distribution on the latent utilities, making the MNP the multinomial analogue of probit.

Key Ideas

Random utility structure: $y_{h,t} = X_{h,t}\beta_h + \varepsilon_{h,t}$ , $\varepsilon_{h,t} \sim \mathcal{N}(0, \Lambda)$ , where $y_{h,t}$ is an $m$ -vector of latent utilities and $I_{h,t} = \arg\max_j y_{h,t,j}$ is the observed choice.
IIA and its avoidance: a scalar (diagonal) covariance $\Lambda = \text{diag}(1, \lambda_2, \ldots, \lambda_m)$ imposes the Independence of Irrelevant Alternatives (IIA) property; a full $\Lambda$ or a correlated random coefficient distribution on $\beta_h$ breaks IIA without requiring a full $m \times m$ unrestricted correlation structure.
Identification: one utility scale normalization is required ( $\lambda_1 = 1$ set to the first alternative), plus a location normalization (utility differences, not levels, are identified); McCulloch-Rossi (1994) work in the space of utility differences.
Hierarchical random coefficients: in the marketing panel setting (Allenby-Rossi 1999), $\beta_h \sim \mathcal{N}(\bar\beta, V_\beta)$ with conjugate priors on $(\bar\beta, V_\beta)$ ; the resulting hierarchical Bayes MNP (HB-MNP) delivers household-level posterior estimates $\pi(\beta_h|\text{data})$ as a by-product of the Gibbs chain.

How It Works

Data Augmentation (Albert-Chib 1993b, eqs. 19–21)

Augment with the full latent utility vector $\mathbf{Z}_{h,t} \sim \mathcal{N}_m(X_{h,t}\beta_h, \Lambda)$ . The observed choice imposes $I_{h,t} = j$ iff $Z_{h,t,j} > Z_{h,t,k}$ for all $k \neq j$ . The augmented Gibbs draws $\mathbf{Z}_{h,t}$ by accept-reject: propose from $\mathcal{N}_m(X_{h,t}\beta_h, \Lambda)$ and accept if the observed choice is the argmax. Efficiency degrades for large $m$ .

McCulloch-Rossi (1994) Gibbs for the Hierarchical MNP

The six-block conditional structure for the Allenby-Rossi (1999) model (eqs. 11–16) draws successively from:

$I_{h,t} | y_{h,t}$ — observed choice index (deterministic)
$y_{h,t} | X_{h,t}, \beta_h, \Lambda$ — latent utility draw (multivariate truncated Normal)
$\beta_h | \bar\beta, V_\beta$ — household coefficients (multivariate Normal, conjugate)
$\bar\beta | \bar{\bar\beta}, a V_\beta$ — population mean (Normal, conjugate)
$V_\beta | v_0, V_0$ — population covariance (Inverse-Wishart, conjugate)
$\Lambda | \nu, s$ — diagonal scale matrix (Inverse-Gamma per element, conjugate)

All blocks are conjugate; no Metropolis-Hastings (M-H) steps are needed for this specification.

Household-Level Inference

Individual posteriors emerge from marginalisation of the joint: $\pi(\beta_h | \text{data}) = \int \pi(\{\beta_i\}, \bar\beta, V_\beta | \text{data})\, d\beta_{-h}\, d\bar\beta\, dV_\beta$ computed by keeping Gibbs draws of $\beta_h$ and discarding the rest. Shrinkage toward the population mean $\bar\beta$ is automatic: sparse households (few purchases) lean heavily on the prior; data-rich households are identified from their own observations.

Demographic Regression Mean (Rossi-McCulloch-Allenby 1996)

A direct extension of the McCulloch-Rossi Gibbs adds observable demographic covariates $z_h$ to the prior mean of $\beta_h$ :

$\beta_h = \Delta z_h + v_h, \quad v_h \sim \mathcal{N}(0, V_\beta), \quad h = 1, \ldots, H$

where $z_h$ is a $d \times 1$ vector of demographics and $\Delta$ is a $k \times d$ matrix of regression coefficients. This allows brand preferences and price sensitivities to vary systematically with income, family size, retirement status, etc., in addition to unobserved heterogeneity captured by $V_\beta$ .

The Gibbs cycle gains one additional block: $\Delta | \{\beta_h\}, V_\beta, \{z_h\}$ — a standard normal posterior from the multivariate regression. The resulting $R^2$ decomposition reveals how much household variation is explained by demographics vs. unobservable sources ( $\rho^2 = 1 - \text{Var}(\epsilon)/\text{Var}(\beta)$ ). In the tuna dataset (Rossi-McCulloch-Allenby 1996), demographics explain only 7–33% of coefficient heterogeneity and just 7% of price sensitivity variation — the bulk is unobservable.

The practical implication for target marketing: demographic data alone yields only 12% incremental gain over blanket couponing, while even one purchase observation with causal variables yields 56%. See Rossi-McCulloch-Allenby (1996).

McCulloch-Polson-Rossi (2000) — ID Prior

The non-identified (NID) approach (McCulloch-Rossi 1994) places the prior on the full unidentified space and reports marginal posteriors of the identified parameters. The induced prior on the identified $\beta$ is a $\sqrt{\chi^2} \times \text{Normal}$ — heavy-tailed and awkward to assess informatively. The ID prior (McCulloch-Polson-Rossi 2000) places the prior directly on the identified space by fixing $\sigma_{11} = 1$ with probability one.

Reparameterisation: Partition $\varepsilon = (U, Z')'$ and define $\gamma = E(UZ)$ , $\Phi = \Sigma_Z - \gamma\gamma'$ (Schur complement). Then $\Sigma = \begin{pmatrix} 1 & \gamma' \\ \gamma & \Phi + \gamma\gamma' \end{pmatrix}$ and the one-to-one map allows: $\gamma \sim \mathcal{N}(\bar\gamma, B^{-1}), \quad \Phi^{-1} \sim W(\kappa, C)$

Gibbs sampler (four blocks, all conjugate, no tuning parameters): $\beta \sim$ Normal; $W_{ij} \sim$ truncated univariate Normal; $\gamma \mid \Phi, \beta, \{W_i\} \sim$ Normal (regression of $Z$ on $U$ ); $\Phi \mid \gamma, \beta, \{W_i\} \sim$ inverted Wishart.

Prior assessment: For $E(\Sigma) = I$ , the prior reduces to two scalars $(\kappa, \tau)$ . Recommended default: $\kappa = p+2$ , $\tau = 1/8$ . Improper $\Phi$ prior is dangerous — analytically it is highly informative about the smallest eigenvalue of $\Sigma$ (pushes it toward zero), causing sampler failure in high dimensions with little data.

Trade-off vs. NID: Informative priors on identified $\beta$ are straightforward; hierarchical extensions $\beta_j \sim \mathcal{N}(\bar\beta, V_\beta)$ apply naturally. Cost: higher chain autocorrelation. The Imai-van Dyk (2005) marginal augmentation approach later showed that the NID chain can be dramatically accelerated, partially closing this gap. See Imai-van Dyk (2005) and McCulloch-Polson-Rossi (2000).

Chib-Greenberg (1998) — Correlated Binary Responses (Multivariate Probit)

Distinct from MNP (one choice among $J$ alternatives), the multivariate probit (MVP) models $J$ simultaneous binary outcomes $Y_i=(Y_{i1},\ldots,Y_{iJ})'$ with $Z_i\sim \mathcal{N}_J(X_i\beta,\Sigma)$ and $Y_{ij}=I(z_{ij}>0)$ . $\Sigma$ must be in correlation form (unit diagonal) for identifiability.

Three-block Gibbs sampler (Chib-Greenberg 1998): (1) $Z_i|y_i,\beta,\Sigma$ — multivariate normal (MVN) truncated to $B_i=\prod_j B_{ij}$ , sampled component-wise via $J$ univariate truncated Normals (Geweke 1991); (2) $\beta|Z,\Sigma$ — conjugate Normal (eq. 6); (3) $\sigma|y,Z,\beta$ — M-H with tailored independence or reflection chain (eq. 7) using Newton-Raphson mode and inverse-Hessian covariance. For large $J$ , the correlation parameters are partitioned into blocks updated sequentially. Marginal likelihood via the Chib (1995) identity extended with kernel density estimation for the $\sigma$ ordinate.

Applications: Six Cities wheezing ( $n=537$ , $J=4$ , equi-correlated vs. unrestricted vs. independent; Bayes factor (BF) $_{2,3}=4.63\times10^{49}$ ); Panel Study of Income Dynamics (PSID) labour force participation ( $n=520$ , $J=7$ , 21 free $\sigma$ ). The MVP framework for binary longitudinal data is the direct parent of Chib-Carlin (1999) Algorithms 5 and 7. See Chib-Greenberg (1998) and Gibbs Sampler.

Marginal Data Augmentation (Imai-van Dyk 2005)

The identification constraint $\sigma_{11} = 1$ makes the standard Gibbs cycle awkward. The unidentifiable scale $\alpha$ (defined by $\tilde{W}_i = \alpha W_i$ for any $\alpha > 0$ ) creates a tension: fixing $\alpha$ slows mixing, but augmenting with it requires a prior that induces an interpretable distribution on the identifiable parameters.

The framework (Meng-van Dyk 1999): Define two strategies for using the working parameter $\alpha$ in a data-augmented sampler:

Conditional augmentation (Scheme 2): Fix $\alpha$ or condition on it when drawing $W$ — the augmented-data model $p(Y,W\mid\theta)$ is less diffuse, producing higher chain autocorrelation.
Marginal augmentation (Scheme 1): Average over a working prior for $\alpha$ before drawing $W$ — $p(Y,W\mid\theta)$ is more diffuse, improving the geometric rate of convergence. The Marginalization Strategy (Meng-van Dyk 1999) proves this can only improve convergence.

Under this framework:

McCulloch-Rossi (1994) uses Scheme 2 — that is why it is slower.
Nobile (1998) adds Step 4 that converts to Scheme 1 — that is why it helps; but it uses the same uninterpretable prior and has an error in the Metropolis acceptance probability $R$ (substitutes $(\tilde\beta, \tilde\Sigma)$ for $(\beta, \Sigma)$ ), distorting the stationary distribution by a factor of $c^{-(2-k-p(p-1)/2)}$ .
McCulloch et al. (2000) fixes $\alpha$ , eliminating the unidentifiable parameter from the chain — Scheme 2 without the benefit of a diffuse augmented-data model; slowest.

New prior specification (Imai-van Dyk 2005):

$\beta \sim \mathcal{N}(\beta_0, A^{-1})$ directly on the identifiable $\beta$ ( $A = 0$ allowed for flat prior)
$p(\Sigma) \propto |\Sigma|^{-(\nu+p)/2}\,[\text{trace}(S\Sigma^{-1})]^{-\nu(p-1)/2}$ subject to $\sigma_{11} = 1$ — the distribution of $\tilde\Sigma/\tilde\sigma_{11}$ where $\tilde\Sigma \sim \text{inv-Wishart}(\nu, \tilde S)$
Working prior: $\alpha^2 \mid \Sigma \sim \alpha_0^2\,\text{trace}(S\Sigma^{-1})/\chi^2_{\nu(p-1)}$ (closed-form conjugate draw — no MH needed)

Algorithm 1 (Scheme 1) — three-step Gibbs with marginal augmentation:

For each $i$ , draw $W_i$ components via univariate truncated normals; then draw $\alpha^2$ from the working prior and set $\tilde{W}_i = \alpha W_i$ (Scheme 1 marginalizes out $\alpha$ )
Draw $(\alpha^2, \tilde\beta) \mid \tilde{W}, \Sigma$ jointly (conjugate normal/chi-squared); set $\beta = \tilde\beta/\alpha$
Draw $\tilde\Sigma \mid \tilde\beta, \tilde{W} \sim \text{inv-Wishart}$ ; set $\Sigma = \tilde\Sigma/\tilde\sigma_{11}$ and $W_i = \tilde{W}_i/\sqrt{\tilde\sigma_{11}}$

Algorithm 2 handles $\beta_0 \neq 0$ via an alternative augmented-data model $\tilde{W}_i = \alpha(W_i - X_i\beta)$ ; Algorithm 1 Scheme 1 is preferred when $\beta_0 = 0$ .

Empirical performance: Detergent data (6 brands, 2657 households) — Algorithm 1 converges in ~1000 draws vs. 60,000–80,000 for McCulloch-Rossi and Nobile (corrected) from a bad starting value; Dutch parliamentary elections (4 choices, 20 covariates) — Algorithm 1 Scheme 1 has clearly lower lag-100 autocorrelations than all competitors.

See Imai and van Dyk (2005) and Data Augmentation.

Why It Matters

Avoids the IIA restriction of multinomial logit without requiring the $(m-1)T_h$ -dimensional integral of classical correlated probit maximum likelihood (eq. 19 in the paper).
Full Bayesian posterior over $\{\beta_h\}$ enables household-specific decision rules: optimal prices, targeted coupons, segment identification — all requiring nonlinear functions of parameters for which point estimates introduce overconfidence (Allenby-Rossi 1999, §4).
Continuous normal mixing distribution for $\beta_h$ outperforms finite mixture models: the finite mixture constrains individual posteriors to the convex hull of mass points, drastically understating tail heterogeneity; formal marginal log-likelihood (MLL) comparison on ketchup data ( $N=1401$ , 8191 purchases) gives $\Delta \ln p = 2070$ log-units in favour of HB.
Foundation for the broad class of hierarchical Bayes discrete-choice models now standard in quantitative marketing and applied microeconometrics.

Open Questions

The accept-reject draw for $\mathbf{Z}_{h,t}$ degrades for large $m$ (many alternatives). Later work develops more efficient multivariate truncated Normal samplers (Geweke 1991; Hajivassiliou-Ruud 1994 Geweke-Hajivassiliou-Keane (GHK) simulator).
The normal mixing distribution may be misspecified when heterogeneity is multimodal or heavy-tailed; McCulloch-Rossi (1996) introduce a mixture of normals for $\pi(\beta_h|\bar\beta, V_\beta)$ as a robustness check.
Scalability: with many alternatives and long panels, the draw of the latent utility matrix is computationally intensive; sparse data structures and parallel Gibbs chains are practical remedies.