Bayesian Multivariate Regression

Python · R  ·  Download conjugate sampler

Model

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.

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

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