import numpy as np, pandas as pd, matplotlib.pyplot as plt
from mvr_conjugate import rmultireg_np_batch, posterior_summary
np.set_printoptions(suppress=True)

Ys = pd.read_csv('mvr_Y.csv').values.astype(float); Xs = pd.read_csv('mvr_X.csv').values.astype(float)
Ns, Ms = Ys.shape; Ks = Xs.shape[1]
trueB = pd.read_csv('mvr_true_B.csv', index_col=0).values
if trueB.shape != (Ks, Ms): trueB = trueB.T
ols = np.linalg.lstsq(Xs, Ys, rcond=None)[0]
Bd, Sd = rmultireg_np_batch(Ys, Xs, np.zeros((Ks,Ms)), 0.01*np.eye(Ks), Ms+3, (Ms+2)*np.eye(Ms), R=10000, seed=1)
summ = posterior_summary(Bd, Ks, Ms, B_true=trueB, ols_B=ols)
print('B: max |posterior mean - OLS| = %.4f   (diffuse prior -> should match OLS)' % (summ['bayes_mean'].values - summ['ols'].values).__abs__().max())
print('95%% credible-interval coverage of the true B: %d / %d' % (summ['covers'].sum(), len(summ)))
print('posterior-mean Sigma vs truth:')
print(np.round(Sd.mean(0).reshape(Ms,Ms,order='F'),3)); print('truth:'); print(pd.read_csv('mvr_true_Sigma.csv', index_col=0).values.round(3))

B: max |posterior mean - OLS| = 0.0010   (diffuse prior -> should match OLS)
95% credible-interval coverage of the true B: 12 / 12
posterior-mean Sigma vs truth:
[[1.357 0.477 0.067]
 [0.477 0.809 0.329]
 [0.067 0.329 0.557]]
truth:
[[1.45 0.63 0.21]
 [0.63 0.89 0.35]
 [0.21 0.35 0.49]]

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 store-weeks, M=%d brands, K=%d regressors' % (N, M, K))

Bd, Sd = rmultireg_np_batch(Y, X, np.zeros((K,M)), 0.01*np.eye(K), M+3, (M+2)*np.eye(M), R=20000, seed=42)
Bmat = Bd.mean(0).reshape((K, M), order='F')                  # (K, M) posterior-mean coefficients
E_cross = Bmat[1:1+M, :].T                                    # (sales eq j, price brand k): rows lp1..lp7
Sig = np.array([s.reshape(M, M, order='F') for s in Sd])      # (R, M, M)
dd  = np.sqrt(np.diagonal(Sig, axis1=1, axis2=2))             # (R, M)
Corr = Sig / (dd[:, :, None] * dd[:, None, :])               # (R, M, M) per-draw correlation
Corr_mean = Corr.mean(0)
print('own-price elasticities (diagonal):', np.round(np.diag(E_cross), 2))

Tuna: N=338 store-weeks, M=7 brands, K=15 regressors

own-price elasticities (diagonal): [-4.29 -4.66 -5.36 -4.87 -4.02  3.11 -2.36]

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])>0.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\n(off-diagonal ≠ 0 ⇒ joint modeling pays)'); 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*0.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\n(diagonal = own-price, off-diagonal = cross-price)'); plt.colorbar(im1, ax=ax[1], fraction=.046)
plt.tight_layout(); plt.savefig('mvr_heatmaps.png', dpi=120, bbox_inches='tight'); plt.show()

fig, ax = plt.subplots(1, 2, figsize=(13.5, 5))
own = np.array([Bd[:, j*K + (1+j)] for j in range(M)])        # own-price draws per brand: B[lp_j, eq_j]
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)'); ax[0].legend(fontsize=8)
# strongest off-diagonal correlation
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\n95%% CrI [%.2f, %.2f] — excludes 0'%(np.percentile(cd,2.5), np.percentile(cd,97.5))); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig('mvr_elasticities.png', dpi=120, bbox_inches='tight'); plt.show()

# ── Same target, two samplers: exact conjugate draws vs random-walk Metropolis ──
# Target: the conditional posterior of one coefficient, B[1,0] (StarKist own-price
# elasticity) given Sigma at its posterior mean. Under the conjugate model this is exactly
# Normal, so both routes are aiming at an identical, known distribution.
A_pri  = 0.01 * np.eye(K)
XtXA   = X.T @ X + A_pri
Bbar   = np.linalg.solve(XtXA, X.T @ Y)                 # posterior mean of B given Sigma
Sig_pm = Sig.mean(0)
m_t    = Bbar[1, 0]
v_t    = Sig_pm[0, 0] * np.linalg.inv(XtXA)[1, 1]       # exact conditional variance
s_t    = np.sqrt(v_t)
R_CMP  = 4000

rng_c = np.random.default_rng(11)
iid   = m_t + s_t * rng_c.standard_normal(R_CMP)        # what the conjugate sampler does

step  = 0.35 * s_t                                       # a reasonable, not-tuned proposal
rng_m = np.random.default_rng(12)
mh    = np.empty(R_CMP)
cur   = m_t - 6 * s_t                                    # started away from the mode
acc   = 0
for t in range(R_CMP):
    prop = cur + step * rng_m.standard_normal()
    if np.log(rng_m.uniform()) < -0.5*((prop-m_t)**2 - (cur-m_t)**2)/v_t:
        cur = prop; acc += 1
    mh[t] = cur
acc /= R_CMP

def acf(x, L=40):
    z = x - x.mean(); d = z @ z
    return np.array([1.0] + [(z[k:] @ z[:-k]) / d for k in range(1, L+1)])

def ess(x):                                              # initial-positive-sequence estimator
    r = acf(x, min(400, len(x)//4)); s = 0.0
    for k in range(1, len(r)-1, 2):
        p = r[k] + r[k+1]
        if p <= 0: break
        s += p
    return len(x) / (1 + 2*s)

ess_i, ess_m = ess(iid), ess(mh)
GREY = '#cbd5e0'
fig, ax = plt.subplots(1, 3, figsize=(16.5, 4.4))

sl = slice(0, 700)
ax[0].axhline(m_t, color=GREY, lw=1.8, zorder=1)
ax[0].plot(np.arange(700), iid[sl], color='seagreen',  lw=.6, zorder=3, label='Conjugate (i.i.d.)')
ax[0].plot(np.arange(700), mh[sl],  color='indianred', lw=.6, zorder=3, label='Metropolis (MCMC)')
ax[0].set_xlabel('draw'); ax[0].set_ylabel('own-price elasticity, StarKist')
ax[0].legend(fontsize=8); ax[0].set_title('First 700 draws', fontsize=10)

L = 40
ax[1].axhline(0, color='k', lw=.8)
ax[1].vlines(np.arange(L+1) - .12, 0, acf(iid, L), color='seagreen',  lw=2, label='Conjugate')
ax[1].vlines(np.arange(L+1) + .12, 0, acf(mh,  L), color='indianred', lw=2, label='Metropolis')
ax[1].set_xlabel('lag'); ax[1].set_ylabel('autocorrelation')
ax[1].legend(fontsize=8); ax[1].set_title('Are consecutive draws related?', fontsize=10)

nn = np.arange(1, R_CMP+1)
ax[2].axhline(m_t, color=GREY, lw=1.8, zorder=1, label='true posterior mean')
ax[2].plot(nn, np.cumsum(iid)/nn, color='seagreen',  lw=1.1, zorder=3, label='Conjugate')
ax[2].plot(nn, np.cumsum(mh)/nn,  color='indianred', lw=1.1, zorder=3, label='Metropolis')
ax[2].set_xscale('log'); ax[2].set_xlabel('draws so far (log scale)')
ax[2].set_ylabel('running mean'); ax[2].legend(fontsize=8)
ax[2].set_title('Converging on the answer', fontsize=10)

fig.suptitle('Graph 5 — the same posterior, sampled two ways', y=1.02)
plt.tight_layout(); plt.show()

print(f'Target (known exactly):  mean {m_t:.4f}   sd {s_t:.4f}')
print(f'{"":26s} {"mean":>9s} {"sd":>8s} {"ESS":>9s} {"ESS / draw":>11s}')
print('-' * 68)
print(f'{"Conjugate (i.i.d.)":26s} {iid.mean():9.4f} {iid.std():8.4f} {ess_i:9.0f} {ess_i/R_CMP:11.2f}')
print(f'{"Metropolis (MCMC)":26s} {mh.mean():9.4f} {mh.std():8.4f} {ess_m:9.0f} {ess_m/R_CMP:11.2f}')
print(f'\n(The conjugate row scores below {R_CMP} only because ESS is itself estimated from a finite')
print(f'sample -- with genuinely independent draws the true value is the draw count.)')
print(f'\nMetropolis acceptance rate: {acc:.0%};  it needs about {ess_i/ess_m:.0f} draws to carry the')
print(f'information of one conjugate draw, and its first {int(np.argmax(np.abs(mh - m_t) < 2*s_t))} draws are')
print(f'still walking in from the starting value -- the burn-in that has to be thrown away.')
print(f'\nThe left panel is the whole argument. The green series has no memory: each draw is generated')
print(f'from the target itself, so it is already representative and the first draw is as good as the')
print(f'thousandth. The red series is a walk -- it starts wherever it was put, takes time to arrive,')
print(f'and thereafter each draw sits close to the one before it. R-hat exists to detect chains that')
print(f'have not arrived; effective sample size exists to price the correlation between neighbouring')
print(f'draws; burn-in exists to discard the walking-in. With exact draws all three are answered by')
print(f'construction, which is why this notebook reports none of them.')

Target (known exactly):  mean -4.2941   sd 0.2791
                                mean       sd       ESS  ESS / draw
--------------------------------------------------------------------
Conjugate (i.i.d.)           -4.2908   0.2790      3402        0.85
Metropolis (MCMC)            -4.2843   0.2892        96        0.02

(The conjugate row scores below 4000 only because ESS is itself estimated from a finite
sample -- with genuinely independent draws the true value is the draw count.)

Metropolis acceptance rate: 89%;  it needs about 36 draws to carry the
information of one conjugate draw, and its first 38 draws are
still walking in from the starting value -- the burn-in that has to be thrown away.

The left panel is the whole argument. The green series has no memory: each draw is generated
from the target itself, so it is already representative and the first draw is as good as the
thousandth. The red series is a walk -- it starts wherever it was put, takes time to arrive,
and thereafter each draw sits close to the one before it. R-hat exists to detect chains that
have not arrived; effective sample size exists to price the correlation between neighbouring
draws; burn-in exists to discard the walking-in. With exact draws all three are answered by
construction, which is why this notebook reports none of them.

#	brand	note
1	StarKist
2	Chicken of the Sea (`ChickSea`)
3	Bumble Bee Solid (`BB Solid`)	← Bumble Bee
4	Bumble Bee Chunk (`BB Chunk`)	← Bumble Bee
5	Geisha
6	Bumble Bee Large (`BB Large`)	← Bumble Bee
7	H-H Chunk (`HH Chunk`)	store/other

Bayesian Multivariate Regression — from-scratch conjugate sampler¶

`mvr_conjugate.py` · the tuna demand system¶

1. Synthetic validation¶

2. Tuna demand system (real data)¶

Graphs 1 & 2 — cross-equation error correlation and the elasticity matrix¶

Graphs 3 & 4 — own-price elasticities and the strongest cross-equation correlation¶

Graph 5 — what independent draws buys, against the same target¶

Analysis of the findings¶

Bayesian Multivariate Regression — from-scratch conjugate sampler¶

mvr_conjugate.py · the tuna demand system¶

1. Synthetic validation¶

2. Tuna demand system (real data)¶

Graphs 1 & 2 — cross-equation error correlation and the elasticity matrix¶

Graphs 3 & 4 — own-price elasticities and the strongest cross-equation correlation¶

Graph 5 — what independent draws buys, against the same target¶

Analysis of the findings¶

`mvr_conjugate.py` · the tuna demand system¶