Bayesian Multinomial Probit Models

Model

For $i=1,\ldots,n$ observations and $p$ mutually exclusive alternatives, define utility differences relative to the base (alternative 1):

The observed choice is $y_i = \arg\max_j\, z_{ij}$. Modeling utility differences removes the unidentified mean: the latent structure is

Identification

Only one scale restriction is needed: $\sigma_{11} = 1$. After each inverse-Wishart draw the draws are normalised:

so the stored draws are of the identified parameters $(\boldsymbol{\beta}^*, \boldsymbol{\Sigma}^*)$.

Priors

where $\mathcal{IW}(\mathbf{V},\nu)$ is the inverse-Wishart with scale $\mathbf{V}$ and degrees of freedom $\nu$.

Gibbs Sampler

Initialise $\boldsymbol{\beta}^{(0)}, \boldsymbol{\Sigma}^{(0)}, \mathbf{W}^{(0)}$. At each iteration:

Block 1 — Latent utility differences $\mathbf{W}$ (truncated normal)

For each non-base alternative $j = 1,\ldots,p-1$, draw from a conditional truncated normal whose bounds depend on the observed choice $y_i$:

Let $c = y_i$ (1-indexed). The truncation bounds for $w_{ij}$ are:

Block 2 — Regression coefficients $\boldsymbol{\beta}$

Given $\mathbf{W}$ and $\boldsymbol{\Sigma}$, the posterior is multivariate normal with precision and mean:

Block 3 — Covariance matrix $\boldsymbol{\Sigma}$

Residuals $\mathbf{E} = \mathbf{W} - \tilde{\mathbf{X}}\boldsymbol{\beta}^\top$ ($n \times (p-1)$). Posterior:

MCMC Diagnostics

Block 1 (latent utility sampling) is the binding constraint on mixing. Each component $w_{ij}$ is drawn from a univariate truncated normal whose conditional variance shrinks as the off-diagonal correlations in $\boldsymbol{\Sigma}$ grow — producing slow exploration when alternatives are highly substitutable. The n_sweeps parameter repeats the full Block 1 sweep multiple times per iteration, improving ESS in absolute terms. However, a structural limit applies: when brand intercepts are absent (M1), $\boldsymbol{\beta}$ and $\boldsymbol{\Sigma}$ are jointly weakly identified — scaling one and adjusting the other leaves the likelihood nearly unchanged — creating a slow drift direction that additional sweeps cannot eliminate. Adding brand intercepts (M2, M3) breaks this confounding and restores near-geometric mixing.

Notebooks

Applied to the IRI margarine scanner-panel dataset (516 households, 4,470 purchase occasions, 10 brands) from Rossi, Allenby & McCulloch (2005), distributed with the bayesm R package. Three models are estimated: price-only (M1), brand intercepts + price for the top 4 brands (M2), and the full 10-brand model (M3).

Downloads

Gibbs Sampler — Source Code

"""
Bayesian Multinomial Probit — Gibbs Sampler
===========================================
McCulloch & Rossi (1994) / McCulloch, Polson & Rossi (1999).

Model
-----
For n observations, p mutually exclusive alternatives (1-indexed).
Define utility differences relative to the base (alternative 1):

    w_ij = z_ij - z_i1,   j = 2, ..., p
    w_i  ~ N(X_i @ beta, Sigma),   X_i : (p-1) x k
    y_i  = argmax_j z_ij   (observed choice, values 1..p)

Identification
--------------
Only one scale restriction: Sigma[0, 0] = 1.
After each IW draw:
    s      = sqrt(Sigma[0, 0])
    Sigma <- Sigma / s^2
    beta  <- beta  / s

Priors
------
    beta  ~ N(betabar, A^{-1})     A : k x k precision matrix
    Sigma ~ IW(nu, V)              V : (p-1) x (p-1)

Gibbs blocks
------------
1. w_i | y_i, beta, Sigma  — sequential univariate truncated Normal
   Bounds depend on which alternative was chosen:
     y = 1 (base)  : all w_j in (-inf, 0)
     y = c > 1, j = c-1 (chosen) : w_j in (max(0, max other w), +inf)
     y = c > 1, j != c-1         : w_j in (-inf, w_{c-1})
2. beta | W, Sigma  — multivariate Normal
3. Sigma | W, beta  — Inverse-Wishart, then normalize Sigma[0,0] = 1
"""

import numpy as np
from scipy import stats
from scipy.stats import wishart


# ---------------------------------------------------------------------------
# Block 1: sample latent utility differences W
# ---------------------------------------------------------------------------

def _draw_W(W, X_cube, MU, y, Sigma, pm1, n_sweeps=1):
    """
    Update W (n x pm1) via sequential univariate truncated normals.
    Vectorised over observations for each component j.

    n_sweeps > 1 repeats the full j-sweep to improve mixing when Sigma has
    high correlations (low conditional variance → slow single-sweep mixing).
    The Sigma-dependent parameters (F, cv) are precomputed once per call.

    Parameters
    ----------
    W        : (n, pm1)     current latent utilities (modified in-place)
    X_cube   : (n, pm1, k)  reshaped design matrix
    MU       : (n, pm1)     conditional means  X_cube @ beta
    y        : (n,) int     choices 1..p
    Sigma    : (pm1, pm1)   current covariance (Sigma[0,0] = 1)
    pm1      : int          p - 1
    n_sweeps : int          number of full Gibbs sweeps per call (default 1)
    """
    n  = len(y)
    ds = np.arange(pm1)

    # 0-indexed position of the chosen non-base alt in w (-1 for base choice)
    chosen_j = np.where(y > 1, y - 2, -1).astype(int)

    # Precompute Sigma-dependent conditional parameters (constant across sweeps)
    cond_params = []
    for j in range(pm1):
        s2 = ds[ds != j]
        if len(s2) > 0:
            Sig_s2  = Sigma[np.ix_(s2, s2)]
            sig_js2 = Sigma[j, s2]
            F  = np.linalg.solve(Sig_s2, sig_js2)
            cv = max(Sigma[j, j] - sig_js2 @ F, 1e-12)
        else:
            F, cv = np.zeros(0), max(Sigma[j, j], 1e-12)
        cond_params.append((s2, F, cv))

    # Masks that don't change across sweeps
    m1 = (y == 1)                           # base chosen
    m2 = [y == j + 2 for j in range(pm1)]  # j is the chosen non-base
    m3 = [(y > 1) & (y != j + 2) for j in range(pm1)]  # another non-base chosen

    for _ in range(n_sweeps):
        for j in range(pm1):
            s2, F, cv = cond_params[j]
            std_val = np.sqrt(cv)

            # Conditional mean: mu_j + F @ (w_{s2} - mu_{s2})
            if len(s2) > 0:
                cond_means = MU[:, j] + (W[:, s2] - MU[:, s2]) @ F
            else:
                cond_means = MU[:, j].copy()

            lo = np.full(n, -np.inf)
            hi = np.full(n,  np.inf)

            # base chosen → all w < 0
            hi[m1] = 0.0

            # j is the chosen non-base → w_j > max(0, other w)
            if m2[j].any():
                if len(s2) > 0:
                    lo[m2[j]] = np.maximum(
                        0.0,
                        W[np.ix_(np.where(m2[j])[0], s2)].max(axis=1)
                    )
                else:
                    lo[m2[j]] = 0.0

            # another non-base chosen → w_j < w_{chosen}
            if m3[j].any():
                idx3 = np.where(m3[j])[0]
                hi[m3[j]] = W[idx3, chosen_j[idx3]]

            a = (lo - cond_means) / std_val
            b = (hi - cond_means) / std_val
            W[:, j] = stats.truncnorm.rvs(a, b, loc=cond_means, scale=std_val)

    return W


# ---------------------------------------------------------------------------
# Block 2: sample beta
# ---------------------------------------------------------------------------

def _draw_beta(W, X_cube, Sigma, betabar, A, brand_price=False):
    """
    beta | W, Sigma ~ N(b_hat, B_hat)

    B_hat = (A + sum_i X_i' Sigma^{-1} X_i)^{-1}
    b_hat = B_hat @ (A @ betabar + sum_i X_i' Sigma^{-1} w_i)

    brand_price : bool
        Fast path for X_i = [I_pm1 | d_i] (brand dummies + one price column).
        Reduces complexity from O(N*k*pm1^2) to O(N*pm1^2) by exploiting the
        block structure: sum_i X_i'S X_i = [[N*S, S*sum(d)], [sum(d)'*S, sum d_i'Sd_i]]
    """
    Sig_inv = np.linalg.inv(Sigma)
    n, pm1, k = X_cube.shape

    if brand_price and k == pm1 + 1:
        # X_i = [I_pm1 | d_i];  last column of X_cube is the price diff vector
        D = X_cube[:, :, -1]             # (n, pm1) price diffs per obs

        SigD  = D @ Sig_inv              # (n, pm1)   d_i' Sig_inv, transposed
        SigW  = W @ Sig_inv              # (n, pm1)   w_i' Sig_inv, transposed

        W_sum = W.sum(axis=0)            # (pm1,)
        D_sum = D.sum(axis=0)            # (pm1,)
        DSD   = np.einsum('ni,ni->', D, SigD)    # scalar  sum_i d_i'Sd_i
        DSW   = np.einsum('ni,ni->', D, SigW)    # scalar  sum_i d_i'Sw_i

        XtSiX = np.empty((k, k))
        XtSiX[:pm1, :pm1] = n * Sig_inv          # N * S
        XtSiX[:pm1, -1]   = Sig_inv @ D_sum      # S * sum(d)
        XtSiX[-1,  :pm1]  = XtSiX[:pm1, -1]     # symmetric
        XtSiX[-1,  -1]    = DSD
        XtSiX += A

        XtSiW = np.empty(k)
        XtSiW[:pm1] = Sig_inv @ W_sum
        XtSiW[-1]   = DSW
        XtSiW += A @ betabar

    else:
        Xt    = X_cube.transpose(0, 2, 1)          # (n, k, pm1)
        XtSiX = A + np.einsum('nkp,pq,nqm->km', Xt, Sig_inv, X_cube)
        XtSiW = A @ betabar + np.einsum('nkp,pq,nq->k', Xt, Sig_inv, W)

    B_hat = np.linalg.inv(XtSiX)
    b_hat = B_hat @ XtSiW
    return np.random.multivariate_normal(b_hat, B_hat)


# ---------------------------------------------------------------------------
# Block 3: sample Sigma + normalize
# ---------------------------------------------------------------------------

def _draw_Sigma(W, X_cube, beta, nu, V):
    """
    Sigma | W, beta ~ IW(V + S, nu + n), then normalize Sigma[0,0] = 1.
    Returns (Sigma_norm, beta_norm).
    """
    n   = len(W)
    MU  = np.einsum('nij,j->ni', X_cube, beta)   # (n, pm1)
    E   = W - MU                                   # (n, pm1) residuals
    S   = E.T @ E                                  # (pm1, pm1)

    nu_post = nu + n
    V_post  = V + S

    iSig  = wishart(df=nu_post, scale=np.linalg.inv(V_post)).rvs()
    Sigma = np.linalg.inv(iSig)

    s     = np.sqrt(Sigma[0, 0])
    return Sigma / s**2, beta / s


# ---------------------------------------------------------------------------
# Main sampler
# ---------------------------------------------------------------------------

def mnprobit_gibbs(Data, Prior, Mcmc):
    """
    Gibbs sampler for the Bayesian multinomial probit.

    Parameters
    ----------
    Data : dict
        'p'  int            number of alternatives (including base)
        'y'  (n,) int       observed choices, values 1..p
        'X'  (n*(p-1), k)   stacked design matrix (utility diffs vs base)
    Prior : dict
        'betabar'  (k,)       prior mean
        'A'        (k, k)     prior precision  (inverse of prior covariance)
        'nu'       int        IW degrees of freedom
        'V'        (p-1,p-1)  IW scale matrix
    Mcmc : dict
        'R'        int        total draws
        'keep'     int        store every keep-th draw  (default 1)
        'nprint'   int        print progress every nprint draws (0 = silent)
        'n_sweeps'    int   Gibbs sweeps over w per iteration (default 1;
                           use 3-5 when ESS is low due to high correlations)
        'brand_price' bool  Fast Block-2 path when X=[I_pm1|price] (default False;
                           set True for brand-dummy + price models — ~5x speedup)
        'beta0'    (k,)       optional starting beta     (default zeros)
        'Sigma0'   (p-1,p-1)  optional starting Sigma   (default identity)

    Returns
    -------
    dict
        'betadraw'   (R//keep, k)         posterior draws of identified beta
        'sigmadraw'  (R//keep, (p-1)^2)   posterior draws of identified Sigma
                                           stored row-major: Sigma.ravel()
    """
    p   = int(Data['p'])
    y   = np.asarray(Data['y'],  dtype=int)
    X   = np.asarray(Data['X'],  dtype=float)
    n   = len(y)
    pm1 = p - 1
    k   = X.shape[1]

    if X.shape[0] != n * pm1:
        raise ValueError(f"X must have n*(p-1) = {n*pm1} rows; got {X.shape[0]}")
    if set(np.unique(y)) != set(range(1, p + 1)):
        raise ValueError(f"y must contain all values 1..{p}")

    betabar = np.asarray(Prior['betabar'], float)
    A       = np.asarray(Prior['A'],       float)
    nu      = int(Prior['nu'])
    V       = np.asarray(Prior['V'],       float)

    R           = int(Mcmc['R'])
    keep        = int(Mcmc.get('keep',        1))
    nprint      = int(Mcmc.get('nprint',      500))
    n_sweeps    = int(Mcmc.get('n_sweeps',    1))
    brand_price = bool(Mcmc.get('brand_price', False))

    beta  = np.asarray(Mcmc.get('beta0',  np.zeros(k)),  float)
    Sigma = np.asarray(Mcmc.get('Sigma0', np.eye(pm1)),  float)

    # Reshape X: (n, pm1, k) for einsum ops
    X_cube = X.reshape(n, pm1, k)

    # Initialise W to satisfy choice constraints
    W = np.zeros((n, pm1))
    for i in range(n):
        c = int(y[i])
        W[i, :] = -0.5
        if c > 1:
            W[i, c - 2] = 0.5

    if nprint > 0:
        print(f"MNP Gibbs: n={n}, p={p}, k={k}, R={R}, keep={keep}, "
              f"n_sweeps={n_sweeps}, brand_price={brand_price}", flush=True)

    n_store     = R // keep
    beta_draws  = np.zeros((n_store, k))
    sigma_draws = np.zeros((n_store, pm1 * pm1))
    store_idx   = 0

    for r in range(1, R + 1):

        # Block 1 — latent utilities (n_sweeps full passes)
        MU = np.einsum('nij,j->ni', X_cube, beta)    # (n, pm1)
        W  = _draw_W(W, X_cube, MU, y, Sigma, pm1, n_sweeps)

        # Block 2 — beta
        beta = _draw_beta(W, X_cube, Sigma, betabar, A, brand_price)

        # Block 3 — Sigma  (+ normalize sigma_11 = 1)
        Sigma, beta = _draw_Sigma(W, X_cube, beta, nu, V)

        if r % keep == 0:
            beta_draws[store_idx]  = beta
            sigma_draws[store_idx] = Sigma.ravel()   # row-major
            store_idx += 1

        if nprint > 0 and r % nprint == 0:
            print(f"  Iteration {r:6d} / {R}", flush=True)

    return {'betadraw': beta_draws, 'sigmadraw': sigma_draws}


# ---------------------------------------------------------------------------
# Convenience helpers
# ---------------------------------------------------------------------------

def make_prior(pm1, k, nu=None, V=None, beta_var=100.0):
    """
    Default weakly-informative prior matching rmnpGibbs default.
    nu = pm1 + 3,  V = nu * I_{pm1}.
    """
    if nu is None:
        nu = pm1 + 3
    if V is None:
        V = float(nu) * np.eye(pm1)   # matches rmnpGibbs default V = nu * I
    return {
        'betabar': np.zeros(k),
        'A'      : (1.0 / beta_var) * np.eye(k),
        'nu'     : nu,
        'V'      : V,
    }


def make_init(p, k):
    """Zero beta, identity Sigma."""
    return {'beta0': np.zeros(k), 'Sigma0': np.eye(p - 1)}

References

• McCulloch, R. & Rossi, P. E. (1994). An exact likelihood analysis of the multinomial probit model. Journal of Econometrics 64(1–2), 207–240. — the data-augmentation Gibbs sampler implemented here, and the weak identification of Σ
• McCulloch, R., Polson, N. G. & Rossi, P. E. (2000). A Bayesian analysis of the multinomial probit model with fully identified parameters. Journal of Econometrics 99(1), 173–193. — the fully-identified re-parameterisation of Σ
• Albert, J. H. & Chib, S. (1993). Bayesian analysis of binary and polychotomous response data. Journal of the American Statistical Association 88(422), 669–679. — the latent-utility truncated-normal augmentation at the core of every block
• Rossi, P. E., Allenby, G. M. & McCulloch, R. (2005). Bayesian Statistics and Marketing. Wiley. — the bayesm::rmnpGibbs cross-check and the margarine scanner-panel data
• Imai, K. & van Dyk, D. A. (2005). A Bayesian analysis of the multinomial probit model using marginal data augmentation. Journal of Econometrics 124(2), 311–334. — marginal augmentation that improves mixing of the correlation parameters