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 that confounding — but it does not buy better mixing, and the runs are explicit about this. M2's effective sample sizes are comparable to M1's, and M3, with nine correlated equations, is lower still: excluding the pinned $\sigma_{11}=1$, no monitored parameter clears 0.6 % of its draws. The sampler here is dependable but expensive, and the posteriors below rest on long runs — M3 is 50,000 draws over 41 minutes — rather than on efficient ones.

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).

Synthetic data first, where the truth is known ($N=2{,}000$, four alternatives, price only). The price coefficient is recovered: true $\beta^* = -2.000$, posterior mean $-1.851$, 95% interval $[-2.065,\ -1.648]$ covering it. $\boldsymbol{\Sigma}^*$ is another matter — $\Sigma_{22}$ comes back $0.742$ against a true $1.000$, and $\Sigma_{23}$ at $0.230$ against $0.400$, a largest error of $0.258$, with the true values sitting out in the tails of their own posteriors — the notebook's opening figure puts the three posteriors side by side against the values that generated them, and the contrast between the tight, well-centred $\beta$ and the two wide, displaced $\boldsymbol{\Sigma}$ panels is the clearest statement of the point on the page. The decisive detail is that this is not an implementation fault: R's bayesm ($\beta^* = -1.847$, $\Sigma_{22} = 0.739$) and PyMC's NUTS ($-1.881$, $0.744$) land in the same wrong place. Three independent samplers, two languages and two entirely different algorithms agreeing on the same biased $\boldsymbol{\Sigma}$ is what makes the point a property of the data rather than of the code: choice data identify $\boldsymbol{\beta}$ sharply and $\boldsymbol{\Sigma}$ weakly — McCulloch & Rossi's (1994) weak identification, shown rather than cited.

Margarine, and what happens when brand intercepts are missing. M1 (price only, top four brands, $N=3{,}377$) returns $\gamma^* = -1.518$, 95% interval $[-1.701,\ -1.329]$, with $P(\gamma^*<0)=1.000$ — but the implied error correlations all come back around $0.92\text{–}0.93$, which is not credible for three distinct margarines and is the tell that $\boldsymbol{\Sigma}$ is absorbing structure the mean function should be carrying. M2 adds brand intercepts and moves both halves at once: the price coefficient triples to $\gamma^* = -4.483$ and the correlations collapse to $0.36$, $0.47$ and $0.15$. bayesm agrees closely ($\gamma^* = -4.470$, intercepts within 0.03). The brand ordering can then be read straight off: House Stick sits $-1.170$ below Parkay and Blue Bonnet $-0.557$, while Shedd Tub at $-0.084$ is statistically indistinguishable from Parkay, its interval $[-0.405,\ +0.188]$ straddling zero.

Willingness to pay. Dividing a brand intercept by the price coefficient converts utility into dollars, and here the answer is unusually concrete: $-\alpha_{\text{Shedd}}/\gamma^* = 0.084/4.483 \approx \$0.02$. The average household would pay only two cents more for Parkay than for Shedd Tub, so a Shedd price cut beyond two cents would in principle flip the average buyer. Because it is a ratio, this figure is free of the utility scale that $\sigma_{11}=1$ fixes — which is what makes it comparable across models, and the hierarchical multinomial logit reaches the same number from a completely different error structure.

The full ten-brand model (M3, $N=4{,}470$, $k=10$, 50,000 draws — a 41-minute run) gives $\gamma^* = -4.122$, $[-4.587,\ -3.670]$, with intercepts spanning Blue Bonnet at $-0.619$ down to House Tub at $-3.087$. Its posterior-predictive shares at observed mean prices track the real ones closely (Parkay 0.428 against 0.395 observed, Blue Bonnet 0.127 against 0.156, Generic 0.069 against 0.070); forcing every price equal scrambles them entirely (Parkay falls to 0.120, Shedd Tub rises to 0.251), which is the model's way of saying that price, not brand loyalty, does most of the work in this market. The notebook carries twelve figures. Three carry most of the argument: Choice data identify $\beta$ sharply, $\boldsymbol{\Sigma}$ weakly on the synthetic run; Omitted-variable attenuation, visualised, which is the M1-to-M2 story in one panel; and the posterior-predictive market-share charts, repeated for M1, M2 and M3, which put the model's implied shares next to the observed ones under several price scenarios at once. The diagnostics table above is backed by a grid of traces, running means and autocorrelation plots for every monitored parameter.

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