Bayesian Multivariate Regression — PyMC companion

Companion to mvr_python.ipynb (from-scratch conjugate). Same tuna demand system, same four graphs — here via PyMC / NUTS: each row $\mathbf y_t\sim N_M(\mathbf B'\mathbf x_t,\Sigma)$, with $B_{kj}\sim N(0,10)$ and an LKJ prior on $\Sigma$ (correlation $\times$ scales). NUTS samples $B$ and $\Sigma$ jointly; we read off the cross-equation correlation and the elasticity matrix from the posterior.

import os, sys
conda_lib = os.path.join(os.environ.get('CONDA_PREFIX',''), 'Library', 'bin')
if os.path.isdir(conda_lib) and conda_lib not in os.environ.get('PATH',''):
    os.environ['PATH'] = conda_lib + ';' + os.environ.get('PATH','')      # Windows DLL fix
import numpy as np, pandas as pd, matplotlib.pyplot as plt, pymc as pm, pytensor.tensor as pt, arviz as az
Y = pd.read_csv('tuna_lmove.csv').values.astype(float); X = pd.read_csv('tuna_X.csv').values.astype(float)
N, M = Y.shape; K = X.shape[1]; brands = ['StarKist','ChickSea','BB Solid','BB Chunk','Geisha','BB Large','HH Chunk']
print('Tuna: N=%d, M=%d brands, K=%d' % (N, M, K))
with pm.Model() as mvr_model:
    B = pm.Normal('B', 0.0, 10.0, shape=(K, M))
    chol, corr, stds = pm.LKJCholeskyCov('cov', n=M, eta=1.5, sd_dist=pm.HalfNormal.dist(1.0, shape=M), compute_corr=True)
    Sigma = pm.Deterministic('Sigma', chol @ chol.T)
    pm.MvNormal('Y', mu=pt.dot(X, B), chol=chol, observed=Y)
    idata = pm.sample(1000, tune=1000, chains=2, target_accept=0.9, random_seed=42, progressbar=False)
print('max r_hat %.3f' % float(az.summary(idata, var_names=['B'])['r_hat'].max()))

Tuna: N=338, M=7 brands, K=15

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (2 chains in 2 jobs)

NUTS: [B, cov]

Sampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 148 seconds.

We recommend running at least 4 chains for robust computation of convergence diagnostics

max r_hat 1.000

Bs = idata.posterior['B'].values.reshape(-1, K, M)              # (S, K, M)
Bmat = Bs.mean(0); E_cross = Bmat[1:1+M, :].T                   # rows lp1..lp7 -> (sales eq, price)
Sig = idata.posterior['Sigma'].values.reshape(-1, M, M)
dd = np.sqrt(np.diagonal(Sig, axis1=1, axis2=2)); Corr = Sig/(dd[:,:,None]*dd[:,None,:]); Corr_mean = Corr.mean(0)
print('own-price elasticities (diagonal):', np.round(np.diag(E_cross), 2))

fig, ax = plt.subplots(1, 2, figsize=(13.5, 5.5))
im0 = ax[0].imshow(Corr_mean, cmap='RdBu_r', vmin=-1, vmax=1)
for i in range(M):
    for j in range(M): ax[0].text(j,i,f'{Corr_mean[i,j]:.2f}',ha='center',va='center',fontsize=7,color='white' if abs(Corr_mean[i,j])>.5 else 'black')
ax[0].set_xticks(range(M)); ax[0].set_xticklabels(brands, rotation=45, ha='right'); ax[0].set_yticks(range(M)); ax[0].set_yticklabels(brands)
ax[0].set_title('Σ error correlation across brands (PyMC)'); plt.colorbar(im0, ax=ax[0], fraction=.046)
vmax=np.abs(E_cross).max(); im1=ax[1].imshow(E_cross, cmap='RdBu_r', vmin=-vmax, vmax=vmax)
for i in range(M):
    for j in range(M): ax[1].text(j,i,f'{E_cross[i,j]:.2f}',ha='center',va='center',fontsize=7,color='white' if abs(E_cross[i,j])>vmax*.5 else 'black')
ax[1].set_xticks(range(M)); ax[1].set_xticklabels([f'P:{b}' for b in brands],rotation=45,ha='right'); ax[1].set_yticks(range(M)); ax[1].set_yticklabels([f'sales:{b}' for b in brands])
ax[1].set_title('Price-elasticity matrix (PyMC)'); plt.colorbar(im1, ax=ax[1], fraction=.046)
plt.tight_layout(); plt.savefig('mvr_pymc_heatmaps.png', dpi=120, bbox_inches='tight'); plt.show()

own-price elasticities (diagonal): [-4.31 -4.66 -5.46 -4.88 -4.16  2.96 -2.39]

Left — the posterior-mean error correlation across the seven brands. The off-diagonals are far from zero (the standout being the Bumble Bee Solid–Large pair), which is exactly why the seven demand equations are worth estimating jointly rather than one at a time. Right — the price-elasticity matrix: a negative diagonal (own-price) and positive off-diagonals (cross-price ⇒ the brands are substitutes). Both panels reproduce the from-scratch conjugate sampler's picture, here from NUTS under an LKJ prior on Σ.

fig, ax = plt.subplots(1, 2, figsize=(13.5, 5))
own = np.array([Bs[:, 1+j, j] for j in range(M)]); mean=own.mean(1); lo=np.percentile(own,2.5,1); hi=np.percentile(own,97.5,1); yy=np.arange(M)
ax[0].errorbar(mean, yy, xerr=[mean-lo, hi-mean], fmt='o', color='firebrick', capsize=4)
ax[0].axvline(0,color='gray',ls=':'); ax[0].axvline(-1,color='steelblue',ls='--',lw=1,label='unit elastic (−1)')
ax[0].set_yticks(yy); ax[0].set_yticklabels(brands); ax[0].invert_yaxis(); ax[0].set_xlabel('own-price elasticity'); ax[0].set_title('Own-price elasticities (95% CrI, PyMC)'); ax[0].legend(fontsize=8)
iu, ju = np.triu_indices(M, 1); k=np.argmax(np.abs(Corr_mean[iu,ju])); i_,j_=iu[k],ju[k]; cd=Corr[:,i_,j_]
ax[1].hist(cd, bins=50, color='steelblue', edgecolor='white', density=True); ax[1].axvline(cd.mean(),color='firebrick',lw=2,label='mean %.2f'%cd.mean()); ax[1].axvline(0,color='black',ls=':')
ax[1].set_xlabel('error correlation corr(%s, %s)'%(brands[i_],brands[j_])); ax[1].set_yticks([])
ax[1].set_title('Strongest cross-equation correlation (PyMC)\n95%% CrI [%.2f, %.2f]'%(np.percentile(cd,2.5),np.percentile(cd,97.5))); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig('mvr_pymc_elasticities.png', dpi=120, bbox_inches='tight'); plt.show()

Left — each brand's own-price elasticity with its 95% credible interval. Six of the seven sit far below the unit-elastic line (≈ −4 to −5): canned-tuna buyers switch brands readily. The lone positive value is the known Bumble Bee collinearity artefact — the same one the conjugate sampler produces. Right — the full posterior of the largest off-diagonal error correlation; its 95% CrI excludes zero, so the cross-equation link is real signal, not noise.

Results

• PyMC's NUTS posterior reproduces the conjugate results: the same error-correlation pattern (the strong $B3$–$B6$ link), the same elasticity matrix (negative own-price for 6 of 7 brands, the same $B6$ collinearity artifact), and the same strongest cross-equation correlation with a CrI excluding 0.
• Two engines agree — the diffuse-prior conjugate i.i.d. sampler and HMC on the LKJ parameterisation land in the same place, validating both.