Bayesian Multivariate Regression

Python · R  ·  Download conjugate sampler

Model

Seven brands of canned tuna share a shelf, and whatever moves one week's sales — a holiday, a heatwave, a rival's promotion — tends to move all seven together. Running seven separate regressions throws that fact away. Fitting them jointly keeps it, and in the particular case where every equation uses the same regressors, the joint answer is available in closed form: there is a formula for the posterior, so nothing has to be simulated at all.

For NN observations and MM dependent variables, the model is:

Y=XB+E,Et∼iidNM(0, Σ)\mathbf{Y} = \mathbf{X}\mathbf{B} + \mathbf{E}, \qquad \mathbf{E}_t \overset{\text{iid}}{\sim} \mathcal{N}_M(\mathbf{0},\, \boldsymbol{\Sigma})

All MM equations share the same N×KN \times K design matrix X\mathbf{X}. When equations have different regressors, use the SUR framework instead.

Prior — Normal–Inverse-Wishart

B∣Σ∼MN(Bˉ,  Σ⊗A−1)\mathbf{B} \mid \boldsymbol{\Sigma} \sim \mathcal{MN}(\bar{\mathbf{B}},\; \boldsymbol{\Sigma} \otimes \mathbf{A}^{-1})
Σ∼IW(ν0,  V0)\boldsymbol{\Sigma} \sim \mathcal{IW}(\nu_0,\; \mathbf{V}_0)

where IW(ν0,V0)\mathcal{IW}(\nu_0, \mathbf{V}_0) is the inverse-Wishart with E[Σ]=V0/(ν0−M−1)\mathbb{E}[\boldsymbol{\Sigma}] = \mathbf{V}_0 / (\nu_0 - M - 1). The Normal-IW is conjugate: the posterior has the same form.

Posterior

The posterior hyperparameters are available in closed form:

A~=A+X⊤X,B~=A~−1 ⁣(ABˉ+X⊤Y)\tilde{\mathbf{A}} = \mathbf{A} + \mathbf{X}^\top\mathbf{X}, \qquad \tilde{\mathbf{B}} = \tilde{\mathbf{A}}^{-1}\!\bigl(\mathbf{A}\bar{\mathbf{B}} + \mathbf{X}^\top\mathbf{Y}\bigr)
V~=V0+Y⊤Y+Bˉ⊤ABˉ−B~⊤A~B~,ν~=ν0+N\tilde{\mathbf{V}} = \mathbf{V}_0 + \mathbf{Y}^\top\mathbf{Y} + \bar{\mathbf{B}}^\top\mathbf{A}\bar{\mathbf{B}} - \tilde{\mathbf{B}}^\top\tilde{\mathbf{A}}\tilde{\mathbf{B}}, \qquad \tilde{\nu} = \nu_0 + N
(B,Σ)∣Y,X  ∼  MN ⁣-IW ⁣(B~,  A~−1,  ν~,  V~)(\mathbf{B},\boldsymbol{\Sigma}) \mid \mathbf{Y},\mathbf{X} \;\sim\; \mathcal{MN}\!\text{-}\mathcal{IW}\!\left(\tilde{\mathbf{B}},\; \tilde{\mathbf{A}}^{-1},\; \tilde{\nu},\; \tilde{\mathbf{V}}\right)

Because the posterior is exact, each call to rmultireg\texttt{rmultireg} (R) or rmultireg_np_batch\texttt{rmultireg\_np\_batch} (Python) produces independent draws — no MCMC, no burn-in, no convergence diagnostics needed. That exemption is worth spelling out, because R-hat, effective sample size and burn-in are the standard checks applied to almost every other posterior in this collection. Those checks exist because a Markov chain only approaches its target: it begins wherever it was initialised, needs a stretch of wandering before its draws are representative, and every draw is correlated with the one before it. Here the conjugate prior leaves the posterior as a Normal–Inverse-Wishart — the standard conjugate pairing for a multivariate Normal, a Normal over the coefficients and an Inverse-Wishart over the covariance matrix, chosen precisely because it returns a posterior of the same family as the prior — with known parameters, so a draw is generated from that distribution directly — representative from the very first one, and statistically independent of all the others. There is no chain to converge, so there is nothing to diagnose.

Notebooks

The from-scratch notebook first validates the conjugate sampler on simulated data, then applies it to a tuna demand system — M=7M=7 brands' log-movement regressed on log-prices and promotion, so B\mathbf{B} is an elasticity matrix and Σ\boldsymbol{\Sigma} captures the cross-brand error correlation. Three engines agree: the from-scratch conjugate i.i.d. sampler, R's bayesm rmultireg, and a PyMC/NUTS fit with an LKJ prior on Σ\boldsymbol{\Sigma} — a prior over correlation matrices with one shape parameter, tunable from uniform across every valid correlation matrix toward the identity, and the usual choice when a covariance needs a prior and conjugacy is unavailable. All recover the same elasticity matrix (negative own-price for 6 of 7 brands), the same strong cross-equation error correlation (with a credible interval excluding zero), and the same collinearity artifact on one brand — validating the conjugate result against HMC on a different Σ\boldsymbol{\Sigma} parameterisation. A closing figure makes the exactness claim concrete rather than asserted: one elasticity is sampled twice from the same known target, once by the conjugate route and once by a random-walk Metropolis sampler. The conjugate draws are uncorrelated and representative immediately; the Metropolis chain spends its first ≈38\approx 38 draws walking in from its starting value and thereafter needs roughly 36 draws to carry the information of one conjugate draw. Burn-in, R^\hat{R} and effective sample size exist to manage exactly those two problems — which is why a conjugate posterior reports none of them.

Python Conjugate Normal–Inverse-Wishart sampler — synthetic validation + tuna elasticity system (Python) Open → Python PyMC cross-check — NUTS with an LKJ prior on Σ, same tuna system (Python) Open → R R cross-check — bayesm::rmultireg on the tuna system (R) Open →

Downloads

mvr_conjugate.py Conjugate Normal–Inverse-Wishart multivariate-regression sampler (NumPy) tuna_lmove.csv Tuna demand system — response (log-movement / sales for 7 brands) tuna_X.csv Tuna demand system — design matrix (log-prices + promotion)

Conjugate Sampler — Source Code

Download mvr_conjugate.py 

"""
mvr_conjugate.py
----------------
Conjugate posterior sampler for the multivariate regression model

    Y = X B + E,   E_t ~ N(0, Sigma)

Prior:  B | Sigma ~ MN(Bbar, Sigma x A^{-1})
        Sigma ~ IW(nu, V)

rmultireg_np() mirrors R bayesm::rmultireg exactly (column-major B storage).
"""

import numpy as np
import pandas as pd
from scipy.stats import invwishart


def rmultireg_np(Y, X, Bbar, A, nu, V, rng):
    """
    One IID draw from the conjugate Normal-IW posterior.

    Parameters
    ----------
    Y    : (N, M) array of dependent variables
    X    : (N, K) design matrix  (same for all equations)
    Bbar : (K, M) prior mean of B
    A    : (K, K) prior precision on B  (small = diffuse)
    nu   : scalar IW degrees of freedom  (must be > M-1)
    V    : (M, M) IW scale matrix
    rng  : numpy Generator  (e.g. np.random.default_rng(42))

    Returns
    -------
    B     : (K, M) draw from posterior
    Sigma : (M, M) draw from posterior
    """
    N, M = Y.shape
    K    = X.shape[1]

    XtX = X.T @ X
    XtY = X.T @ Y

    # Posterior hyperparameters
    A_tilde = A + XtX
    A_tilde_inv = np.linalg.inv(A_tilde)
    B_tilde = A_tilde_inv @ (A @ Bbar + XtY)

    V_tilde = (V
               + Y.T @ Y
               + Bbar.T @ A @ Bbar
               - B_tilde.T @ A_tilde @ B_tilde)
    nu_tilde = nu + N

    # Draw Sigma ~ IW(nu_tilde, V_tilde)
    Sigma = invwishart.rvs(df=nu_tilde, scale=V_tilde, random_state=rng)

    # Draw B | Sigma ~ MN(B_tilde, Sigma x A_tilde^{-1})
    L_Ainv = np.linalg.cholesky(A_tilde_inv)
    L_Sig  = np.linalg.cholesky(Sigma)
    Z      = rng.standard_normal((K, M))
    B      = B_tilde + L_Ainv @ Z @ L_Sig.T

    return B, Sigma


def rmultireg_np_batch(Y, X, Bbar, A, nu, V, R, seed=42):
    """
    Draw R IID samples from the conjugate posterior.

    Returns
    -------
    B_draws     : (R, K*M)  column-major, matching R as.vector(B)
    Sigma_draws : (R, M*M)  column-major
    """
    N, M = Y.shape
    K    = X.shape[1]
    rng  = np.random.default_rng(seed)

    # Precompute posterior hyperparameters once
    XtX     = X.T @ X
    XtY     = X.T @ Y
    A_tilde = A + XtX
    A_tilde_inv = np.linalg.inv(A_tilde)
    B_tilde = A_tilde_inv @ (A @ Bbar + XtY)
    V_tilde = (V
               + Y.T @ Y
               + Bbar.T @ A @ Bbar
               - B_tilde.T @ A_tilde @ B_tilde)
    nu_tilde = nu + N

    L_Ainv = np.linalg.cholesky(A_tilde_inv)

    B_draws     = np.empty((R, K * M))
    Sigma_draws = np.empty((R, M * M))

    for r in range(R):
        Sigma = invwishart.rvs(df=nu_tilde, scale=V_tilde, random_state=rng)
        L_Sig = np.linalg.cholesky(Sigma)
        Z     = rng.standard_normal((K, M))
        B     = B_tilde + L_Ainv @ Z @ L_Sig.T
        B_draws[r]     = B.ravel(order='F')   # column-major = R's as.vector
        Sigma_draws[r] = Sigma.ravel(order='F')

    return B_draws, Sigma_draws


def posterior_summary(B_draws, K, M, eq_names=None, param_names=None,
                      B_true=None, ols_B=None):
    """
    Build a tidy DataFrame summarising B posterior draws.

    B_draws : (R, K*M) column-major array (output of rmultireg_np_batch)
    """
    if eq_names    is None: eq_names    = [f'y{m+1}' for m in range(M)]
    if param_names is None: param_names = [f'x{k}'   for k in range(K)]

    rows = []
    for m in range(M):
        for k in range(K):
            d = B_draws[:, m * K + k]
            row = dict(
                equation  = eq_names[m],
                regressor = param_names[k],
                bayes_mean= d.mean(),
                bayes_sd  = d.std(),
                CI_lo     = np.quantile(d, 0.025),
                CI_hi     = np.quantile(d, 0.975),
            )
            if B_true is not None:
                row['truth'] = B_true[k, m]
            if ols_B is not None:
                row['ols'] = ols_B[k, m]
            rows.append(row)

    df = pd.DataFrame(rows)
    if B_true is not None:
        df['covers'] = (df['CI_lo'] <= df['truth']) & (df['truth'] <= df['CI_hi'])
    return df

References