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

Two applications, both classics of the literature, and the sampler is judged against a published benchmark on the first.

Ohio Children's Wheeze — 537 children observed at ages 7, 8, 9 and 10, with maternal smoking as the covariate. This is Chib & Greenberg's (1998) own primary example, fitted here in their constrained form ($\mathbf{B} = \mathbf{A}_J\boldsymbol{\Gamma}$, four free coefficients plus six correlations), which makes it a benchmark rather than an illustration. The posterior means reproduce their Table 4 maximum-likelihood fit almost exactly: the six within-child correlations come back $0.585,\ 0.529,\ 0.688,\ 0.566,\ 0.565,\ 0.632$ against the paper's $0.584,\ 0.521,\ 0.688,\ 0.586,\ 0.562,\ 0.631$ — four of the six agree to within 0.003 and the worst discrepancy is 0.020 — with the intercept at $-1.114$ against $-1.118$ and the age slope at $-0.079$ against $-0.079$. Substantively, wheeze declines with age and maternal smoking raises it ($+0.150$ at age 9), though that interval covers zero (SD 0.101). The correlations are the real finding: every pair sits between 0.53 and 0.69, so a child who wheezes at 7 is far more likely to wheeze at 10, and any analysis treating the four visits as four independent observations would overstate its own precision.

Scotch Whisky — 2,219 households by 21 brands, which means 46,599 latent utilities, 21 coefficients and 231 covariance parameters drawn every sweep. The fit reproduces every brand's purchase rate to within 0.002 (Chivas Regal 0.363 observed against 0.363 implied). What makes the model worth the cost is the correlation matrix. Ordinary Pearson correlations computed on the raw 0/1 indicators span only $[-0.133,\ 0.360]$ — binarising a continuous preference crushes the association it is meant to measure. The tetrachoric correlations the Gibbs sampler estimates, which are correlations of the underlying latent utilities, span $[-0.323,\ 0.675]$. Decomposing that estimated matrix into principal components recovers 5 components with eigenvalue above 1, explaining 58.4% of the variance — against the 5 components and 57.0% reported by Edwards & Allenby (2003) — and the first component correlates $-0.932$ with brand rating, matching the paper's $-0.932$ to three decimals.

What it costs. The coefficients mix acceptably (ESS 3,200–8,700 out of 55,000 draws) but the correlations do not: the worst pair, $R[\text{the Singleton},\ \text{Chivas Regal}]$, returns an ESS of 210 — roughly one usable draw in every 260 — which is why the run is 60,000 draws rather than the few thousand the coefficients alone would need. R's bayesm rmvpGibbs reproduces the same brand rates and coefficients on the Scotch data.

Between them the two Python notebooks carry twelve figures. Four are worth opening first, all on the Scotch side: the $21\times21$ posterior-mean correlation heat-map; the Attenuation: Pearson vs tetrachoric plot, which shows directly how much association a plain binary correlation throws away; its companion comparing pairwise-MLE against the joint Gibbs tetrachorics; and the $\text{PC1}\times\text{PC2}$ loadings plot, which separates single malts from blends without ever having been told which is which. The Ohio notebook is deliberately leaner — a prevalence-by-age-and-smoking chart, trace plots for the worst-mixing parameters, and its own correlation heat-map; the comparison against Chib & Greenberg's Table 4 is reported there as a table rather than a figure.

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