""" 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, }