Bayesian Multivariate Probit Models

Model

Let $i = 1, \ldots, N$ index households and $j = 1, \ldots, J$ index brands. The binary multivariate probit model posits a latent utility structure:

• $\mathbf{Z}$ — $N \times J$ matrix of latent utilities
• $\mathbf{X}$ — $N \times K$ design matrix
• $\mathbf{B}$ — $J \times K$ matrix of regression coefficients
• $\boldsymbol{\Sigma}$ — $J \times J$ covariance matrix of utilities

The observed binary outcomes are:

Identification

$\boldsymbol{\Sigma}$ is only identified up to scale. We impose the standard normalisation at each draw:

so the stored draws are of the identified parameters $(\mathbf{B}^*, \mathbf{R})$.

Priors

where $\mathcal{IW}(\mathbf{V}, \nu)$ is the inverse-Wishart distribution with scale $\mathbf{V}$ and degrees of freedom $\nu$, so $\mathbb{E}[\boldsymbol{\Sigma}] = \mathbf{V}/(\nu - J - 1)$.

Gibbs Sampler

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

Block 1 — Latent utilities $\mathbf{W}$ (Albert–Chib)

For each brand $j$, sample from a conditional truncated normal:

with conditional mean and variance from the multivariate normal partition:

Block 2 — Regression coefficients $\mathbf{B}$

Let $\mathbf{C} = \operatorname{chol}(\boldsymbol{\Sigma}^{-1})$ (upper Cholesky). Posterior:

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

Residuals $\mathbf{E} = \mathbf{W} - \mathbf{X}\mathbf{B}^\top$. Posterior:

Normalise: compute $\mathbf{T}, \mathbf{R}, \mathbf{B}^*$ and store.

Notebooks

Gibbs Sampler — Source Code

Download mvprobit_gibbs.py

"""
Bayesian Multivariate Probit — Gibbs Sampler
=============================================
Python translation of B_mvprobit7b.m (Edwards & Allenby 2003).

Model
-----
    Z = X @ B.T + E,   E ~ N(0, Sigma)   (N x J latent utilities)
    y_ij = 1  if  z_ij > 0,  else 0

Identification
--------------
Sigma is not identified in scale — only the correlation matrix is.
At each draw we store the *normalized* parameters:
    B_norm   = TT @ B       where TT = diag(1/sqrt(diag(Sigma)))
    R        = TT @ Sigma @ TT    (correlation matrix)

Priors
------
    vec(B) ~ N(Bpmean, Bpcov)
    Sigma  ~ IW(V, v)          (inverse-Wishart, df = v, scale = V)

Gibbs blocks (per iteration)
-----------------------------
1. W | B, Sigma, y  — truncated normal, brand by brand (Albert-Chib)
2. B | W, Sigma     — multivariate normal
3. Sigma | W, B     — inverse-Wishart, then normalize
"""

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


# ---------------------------------------------------------------------------
# Helper: vectorized truncated-normal sampler
# ---------------------------------------------------------------------------

def _truncnorm_rnd(mu, var, lo, hi):
    """
    Sample x ~ TN(mu, var) restricted to (lo, hi).

    Parameters
    ----------
    mu  : (n,) array of means
    var : scalar variance (same for all observations)
    lo, hi : scalar bounds  (use ±np.inf for one-sided truncation)

    Returns
    -------
    (n,) array of draws
    """
    std = np.sqrt(var)
    a = (lo - mu) / std
    b = (hi - mu) / std
    return stats.truncnorm(a, b, loc=mu, scale=std).rvs()


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

def mvprobit_gibbs(y, x, prior, ndraws, b0, burn=0, verbose=100, A_age=None):
    """
    Gibbs sampler for the Bayesian multivariate probit model.

    Parameters
    ----------
    y      : (N, J) int array  — binary outcomes (0/1)
    x      : (N, K) float array — design matrix (include intercept)
    prior  : dict with keys
               'Bpmean' : (K*J,)   prior mean of vec(B)   [col-major: J varies first]
               'Bpcov'  : (K*J, K*J) prior covariance of vec(B)
               'v'      : scalar    inverse-Wishart df
               'V'      : (J, J)    inverse-Wishart scale matrix
    ndraws : int  — total number of Gibbs draws (including burn-in)
    b0     : dict with keys
               'Beta'  : (K*J,)  initial beta vector
               'Sigma' : (J, J)  initial covariance matrix
               'W'     : (N, J)  initial latent utilities
    burn   : int  — number of burn-in draws to discard (default 0)
    verbose: int  — print progress every this many draws (0 = silent)

    A_age  : (J, n_a) float array or None
               If provided, imposes B = A_age @ Gamma (constrained model).
               Gamma is (n_a, K); only n_a*K parameters are sampled instead of J*K.
               Example for Chib & Greenberg (1998): A_age = [[1,-2],[1,-1],[1,0],[1,1]]
               with age centred at 9. n_a = 2, K = 2 → 4 free parameters.

    Returns
    -------
    dict with:
      'B_draws'     : (ndraws-burn, K*J)   normalized beta draws
      'Sigma_draws' : (ndraws-burn, J*J)   lower-triangle of correlation matrix (col-major)
      'Gamma_draws' : (ndraws-burn, n_a*K) normalized Gamma draws (only when A_age given)
      'lastB'       : (K*J,)   final unidentified beta draw
      'lastSigma'   : (J, J)   final unidentified Sigma draw
      'lastW'       : (N, J)   final latent utility draw
    """
    N, J = y.shape
    _, K = x.shape
    KJ = K * J

    # ── Prior ──────────────────────────────────────────────────────────────
    Bpmean    = prior['Bpmean'].copy()           # (KJ,)
    Bpcov     = prior['Bpcov'].copy()            # (KJ, KJ)
    iBpcov    = np.linalg.inv(Bpcov)             # (KJ, KJ)
    iBpcovmean = iBpcov @ Bpmean                  # (KJ,)
    v         = prior['v']                        # scalar df
    V         = prior['V'].copy()                 # (J, J)

    # ── Initial values ─────────────────────────────────────────────────────
    beta0  = b0['Beta'].copy()    # (KJ,)  — unidentified
    Sigma0 = b0['Sigma'].copy()   # (J, J) — unidentified covariance
    W0     = b0['W'].copy()       # (N, J) — latent utilities

    # ── Pre-computed fixed quantities ───────────────────────────────────────
    # Index sets for each brand: who purchased (y1) and who didn't (y0)
    y1 = [np.where(y[:, j] == 1)[0] for j in range(J)]
    y0 = [np.where(y[:, j] == 0)[0] for j in range(J)]

    ds   = np.arange(J)          # brand index array [0, 1, ..., J-1]
    XtX  = x.T @ x               # (K, K)

    # ── Constrained-model setup (A_age given) ───────────────────────────────
    constrained = A_age is not None
    if constrained:
        n_a   = A_age.shape[1]           # number of age-basis columns (2 for C&G)
        KnA   = K * n_a                  # free parameters in Gamma
        iGcov = np.eye(KnA) / prior.get('beta_var', Bpcov[0, 0])

    # ── Storage ─────────────────────────────────────────────────────────────
    keep      = ndraws - burn
    B_out     = np.zeros((keep, KJ))
    Sigma_out = np.zeros((keep, J * J))
    Gamma_out = np.zeros((keep, KnA)) if constrained else None

    # ── Gibbs loop ──────────────────────────────────────────────────────────
    for it in range(ndraws):

        # ── Compute X @ B  (N x J) ─────────────────────────────────────────
        # MATLAB: xb = x * reshape(beta0, J, K)'
        # reshape uses column-major so that brand index varies fastest
        B_mat = beta0.reshape(J, K, order='F')   # (J, K)
        xb    = x @ B_mat.T                       # (N, J)

        iSigma0 = np.linalg.inv(Sigma0)

        # ── Block 1: Sample W (latent utilities) ────────────────────────────
        # For each brand j, sample z_j from its conditional distribution
        # given z_{-j}, restricted to the correct half-line.
        #
        # Conditional of z_j | z_{-j}:
        #   mean = xb[:,j] + F' @ (W[:,-j] - xb[:,-j])   (row-by-row)
        #   var  = Sigma[j,j] - Sigma[-j,j]' @ F
        # where F = Sigma[-j,-j] \ Sigma[-j,j]
        for j in range(J):
            s2 = ds[ds != j]                              # brands other than j

            # Conditional regression coefficients onto other brands
            F = np.linalg.solve(Sigma0[np.ix_(s2, s2)], Sigma0[s2, j])  # (J-1,)

            # Conditional variance (scalar)
            Wvarj = Sigma0[j, j] - Sigma0[s2, j] @ F

            # Conditional mean (N,)
            Wmeanj = F @ (W0[:, s2] - xb[:, s2]).T + xb[:, j]

            # Sample from truncated normal
            W0[y1[j], j] = _truncnorm_rnd(Wmeanj[y1[j]], Wvarj,  0,      np.inf)
            W0[y0[j], j] = _truncnorm_rnd(Wmeanj[y0[j]], Wvarj, -np.inf, 0)

        # ── Block 2: Sample beta ────────────────────────────────────────────
        C = np.linalg.cholesky(iSigma0).T         # (J, J) upper triangular

        if constrained:
            # Constrained update: B = A_age @ Gamma, Gamma is (n_a, K)
            # Posterior of vec(Gamma'): kron(A'Σ⁻¹A, X'X) + iGcov
            AtsigA  = A_age.T @ iSigma0 @ A_age           # (n_a, n_a)
            iBcov_G = np.kron(AtsigA, XtX) + iGcov        # (KnA, KnA)
            # Sufficient statistic: vec(X'W Σ⁻¹ A_age)
            Temp_G  = (x.T @ W0 @ iSigma0 @ A_age).ravel(order='F')  # (KnA,)
            Gmean   = np.linalg.solve(iBcov_G, Temp_G)
            gamma   = np.random.multivariate_normal(Gmean, np.linalg.inv(iBcov_G))
            Gamma   = gamma.reshape(K, n_a, order='F').T   # (n_a, K)
            B_mat   = A_age @ Gamma                        # (J, K)
            beta0   = B_mat.ravel(order='F')               # (KJ,)
        else:
            # Unconstrained update (original)
            iBcov = np.kron(XtX, C @ C.T) + iBpcov        # (KJ, KJ)
            Temp1 = C @ (W0 @ C).T @ x                    # (J, K)
            Temp2 = Temp1.reshape(KJ, order='F')           # (KJ,)
            Bmean = np.linalg.solve(iBcov, Temp2 + iBpcovmean)
            Bcov  = np.linalg.inv(iBcov)
            beta0 = np.random.multivariate_normal(Bmean, Bcov)

        # ── Block 3: Sample Sigma ───────────────────────────────────────────
        B_mat = beta0.reshape(J, K, order='F')     # (J, K)
        Es    = W0 - x @ B_mat.T                   # (N, J) residuals
        S     = Es.T @ Es                           # (J, J)

        # Posterior: iSigma ~ Wishart(inv(S + V), N + v)
        iSigma0 = wishart(df=N + v, scale=np.linalg.inv(S + V)).rvs()
        Sigma0  = np.linalg.inv(iSigma0)

        # ── Identification normalization ────────────────────────────────────
        # Scale so that diag(Sigma) = 1 (correlation matrix)
        TT = np.diag(1.0 / np.sqrt(np.diag(Sigma0)))   # (J, J)

        # ── Store (after burn-in) ───────────────────────────────────────────
        if it >= burn:
            idx = it - burn
            B_mat_norm = TT @ B_mat
            B_out[idx, :] = B_mat_norm.ravel(order='F')

            if constrained:
                # Normalized Gamma: TT @ A_age @ Gamma = B_norm → Gamma_norm = pinv(A_age) @ B_norm
                Gamma_norm = np.linalg.lstsq(A_age, B_mat_norm, rcond=None)[0]  # (n_a, K)
                Gamma_out[idx, :] = Gamma_norm.ravel(order='F')

            R = TT @ Sigma0 @ TT
            Sigma_out[idx, :] = np.tril(R).reshape(J * J, order='F')

        if verbose and (it + 1) % verbose == 0:
            print(f"  Iteration {it+1:>5} / {ndraws}", flush=True)

    out = {
        'B_draws'     : B_out,
        'Sigma_draws' : Sigma_out,
        'lastB'       : beta0,
        'lastSigma'   : Sigma0,
        'lastW'       : W0,
    }
    if constrained:
        out['Gamma_draws'] = Gamma_out
    return out


# ---------------------------------------------------------------------------
# Convenience: default flat prior and zero initialisation
# ---------------------------------------------------------------------------

def make_prior(J, K, v=None, V=None, beta_var=100.0):
    """
    Flat-ish default prior.
    v = J + 1 (minimal df for proper IW), V = identity.
    """
    KJ = K * J
    if v is None:
        v = J + 1
    if V is None:
        V = np.eye(J)
    return {
        'Bpmean': np.zeros(KJ),
        'Bpcov' : beta_var * np.eye(KJ),
        'v'     : v,
        'V'     : V,
    }


def make_init(y, x, J, K):
    """
    Simple initialisation: W from standard truncated normal, beta=0, Sigma=I.
    """
    N = y.shape[0]
    W0 = np.zeros((N, J))
    for j in range(J):
        W0[y[:, j] == 1, j] =  0.5
        W0[y[:, j] == 0, j] = -0.5
    return {
        'Beta' : np.zeros(K * J),
        'Sigma': np.eye(J),
        'W'    : W0,
    }

Data

The Scotch Whisky dataset (Corbit & Jain 1994, as used in Edwards & Allenby 2003) records whether each of 2,219 households purchased each of 21 Scotch brands — a binary outcome matrix of 2,219 × 21.

The Ohio Children's Wheeze data (Fitzmaurice & Laird 1993; the Steubenville cohort of the Harvard Six Cities study) records wheeze status for 537 children at ages 7–10 alongside maternal smoking — the primary example in Chib & Greenberg (1998). It is the geepack::ohio dataset, shipped here as a CSV so the notebook runs without a network connection.

References

• Chib, S. & Greenberg, E. (1998). Analysis of multivariate probit models. Biometrika 85(2), 347–361. — the foundational Bayesian MVP paper; the correlation and Gamma estimates are benchmarked against its Table 4
• 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 data-augmentation Gibbs step, extended here to J correlated equations
• Fitzmaurice, G. M. & Laird, N. M. (1993). A likelihood-based method for analysing longitudinal binary responses. Biometrika 80(1), 141–151. — source of the Ohio wheeze cohort
• Edwards, Y. D. & Allenby, G. M. (2003). Multivariate analysis of multiple response data. Journal of Marketing Research 40(3), 321–334. — the Scotch-whisky application and the bayesm lineage
• Rossi, P. E., Allenby, G. M. & McCulloch, R. (2005). Bayesian Statistics and Marketing. Wiley. — bayesm::rmvpGibbs, the R cross-check