""" Hierarchical MNL with a MIXTURE-OF-NORMALS heterogeneity distribution -- from-scratch Gibbs. Replicates bayesm::rhierMnlRwMixture (Rossi, Allenby & McCulloch 2005, Ch. 5). Model ----- P(choice = j | task) = exp(x_j' beta_i) / sum_l exp(x_l' beta_i) (MNL per respondent) beta_i | s_i = k ~ N(mu_k + Delta' z_i, Sigma_k) (mixture heterogeneity) s_i | pi ~ Categorical(pi) pi ~ Dirichlet(a) mu_k ~ N(0, Amu^{-1} I), Sigma_k ~ IW(nu, V) The mixture of normals is used here as a FLEXIBLE heterogeneity distribution (it can approximate a skewed / multimodal population of tastes), NOT necessarily as discrete segments -- the Rossi et al. point. The single-normal model is the special case K=1. Sampler (RW-Metropolis within Gibbs) ------------------------------------ 1. beta_i -- random-walk Metropolis; proposal N(beta_i, (H_i + Sigma_{s_i}^{-1})^{-1}) where H_i is respondent i's MNL information at the pooled MLE (the bayesm move), so the step adapts to data + the component's prior precision each sweep. 2. s_i -- categorical, P(k) propto pi_k N(beta_i; mu_k + Delta'z_i, Sigma_k) 3. mu_k,Sigma_k -- Normal / Inverse-Wishart conjugate from the beta_i assigned to component k 4. pi -- Dirichlet 5. Delta -- matrix-normal GLS [only if Z given] Components are sorted each sweep (by mu_k on the last coordinate) to stabilise labels; the fitted heterogeneity density is label-invariant regardless. Companion to hier_mnl_gibbs.py (single-normal) and hlm_mixture_gibbs.py (linear mixture). """ import numpy as np from scipy.stats import invwishart from scipy.linalg import cholesky, solve_triangular def mnl_loglik(beta, X, y): """MNL log-likelihood for one respondent. X:(T,J,k), y:(T,).""" util = X @ beta m = util.max(axis=1, keepdims=True) logden = m[:, 0] + np.log(np.exp(util - m).sum(axis=1)) return float((util[np.arange(len(y)), y] - logden).sum()) def mnl_hess(beta, X): """Observed information sum_t X_t'(diag(p)-pp')X_t for one respondent.""" util = X @ beta; util -= util.max(axis=1, keepdims=True) e = np.exp(util); p = e / e.sum(axis=1, keepdims=True) k = X.shape[2]; H = np.zeros((k, k)) for t in range(X.shape[0]): Xt = X[t]; pt = p[t] H += Xt.T @ (np.diag(pt) - np.outer(pt, pt)) @ Xt return H def _pooled_beta(data, k, maxit=25): b = np.zeros(k) for _ in range(maxit): g = np.zeros(k); H = np.zeros((k, k)) for X, y in data: util = X @ b; util -= util.max(axis=1, keepdims=True) e = np.exp(util); p = e / e.sum(axis=1, keepdims=True) for t in range(len(y)): g += X[t, y[t]] - p[t] @ X[t] H += mnl_hess(b, X) step = np.linalg.solve(H + 1e-6 * np.eye(k), g); b += step if np.max(np.abs(step)) < 1e-8: break return b def _logmvn(x, mu, L): dev = x - mu; v = solve_triangular(L, dev, lower=True) return -0.5 * len(x) * np.log(2 * np.pi) - np.log(np.diag(L)).sum() - 0.5 * float(v @ v) def simulate_hmnl_mixture(N=500, T=15, J=3, k=3, pi=None, mus=None, Sigmas=None, seed=0): """Simulate hier-MNL with mixture heterogeneity. Default: bimodal price sensitivity.""" rng = np.random.default_rng(seed) if pi is None: pi = np.array([0.55, 0.45]) mus = np.array([[1.0, -0.5, -1.0], # price-insensitive segment (price coef -1) [1.0, -0.5, -4.0]]) # price-sensitive segment (price coef -4) Sigmas = np.stack([0.3 * np.eye(k), 0.3 * np.eye(k)]) K = len(pi) comps = rng.choice(K, size=N, p=pi) Ls = [np.linalg.cholesky(Sigmas[c]) for c in range(K)] B = np.array([mus[comps[i]] + Ls[comps[i]] @ rng.standard_normal(k) for i in range(N)]) data = [] for i in range(N): X = rng.normal(size=(T, J, k)) X[:, :, :k-1] = (rng.random((T, J, k-1)) < 0.5).astype(float) # brand-dummy-like util = X @ B[i] y = (util + rng.gumbel(size=(T, J))).argmax(axis=1) data.append((X, y)) return dict(data=data, B_true=B, comps=comps, pi=pi, mus=mus, Sigmas=Sigmas) def hmnl_mixture_gibbs(data, K=2, Z=None, R=4000, burn=1500, seed=0, s_scale=None, a=5.0, Amu=0.01, nu=None, V=None): """Hierarchical MNL with K-component mixture-of-normals heterogeneity.""" rng = np.random.default_rng(seed) N = len(data); k = data[0][0].shape[2] s_scale = (2.38 / np.sqrt(k)) if s_scale is None else s_scale nu = (k + 3) if nu is None else nu V = np.eye(k) if V is None else V use_Z = Z is not None if use_Z: Z = np.asarray(Z, float); nz = Z.shape[1]; Delta = np.zeros((nz, k)); Ad = 0.01 * np.eye(nz) bpool = _pooled_beta(data, k) Hs = [mnl_hess(bpool, X) for X, _ in data] B = np.tile(bpool, (N, 1)).astype(float) mu = np.tile(bpool, (K, 1)) + 0.01 * rng.standard_normal((K, k)) Sigma = np.stack([np.eye(k)] * K); Sinv = np.stack([np.eye(k)] * K) pi = np.full(K, 1.0 / K) s = rng.integers(0, K, N) llcur = np.array([mnl_loglik(B[i], data[i][0], data[i][1]) for i in range(N)]) keep = R - burn MU = np.zeros((keep, K, k)); SIG = np.zeros((keep, K, k, k)); PI = np.zeros((keep, K)) LL = np.zeros(keep); Bm = np.zeros((N, k)); acc = 0 DD = np.zeros((keep, nz, k)) if use_Z else None for r in range(R): Ls = [cholesky(Sigma[c], lower=True) for c in range(K)] # 1. beta_i (RW-Metropolis, component-aware proposal) for i in range(N): ki = s[i]; mu_eff = mu[ki] + (Delta.T @ Z[i] if use_Z else 0.0) Lp = cholesky(Hs[i] + Sinv[ki], lower=True) prop = B[i] + s_scale * solve_triangular(Lp.T, rng.standard_normal(k), lower=False) llp = mnl_loglik(prop, data[i][0], data[i][1]) lpr = _logmvn(prop, mu_eff, Ls[ki]) - _logmvn(B[i], mu_eff, Ls[ki]) if np.log(rng.random()) < (llp - llcur[i] + lpr): B[i] = prop; llcur[i] = llp; acc += 1 # 2. s_i logpi = np.log(np.clip(pi, 1e-300, None)) for i in range(N): lw = np.array([logpi[c] + _logmvn(B[i], mu[c] + (Delta.T @ Z[i] if use_Z else 0.0), Ls[c]) for c in range(K)]) lw -= lw.max(); w = np.exp(lw); w /= w.sum() s[i] = rng.choice(K, p=w) # 3. mu_k, Sigma_k for c in range(K): mask = s == c; nk = int(mask.sum()) Be = B[mask] - (Z[mask] @ Delta if use_Z else 0.0) if nk > 0: Pm = Amu * np.eye(k) + nk * Sinv[c] rhs = nk * Sinv[c] @ Be.mean(0) Lm = cholesky(Pm, lower=True) mm = solve_triangular(Lm.T, solve_triangular(Lm, rhs, lower=True), lower=False) mu[c] = mm + solve_triangular(Lm.T, rng.standard_normal(k), lower=False) dev = Be - mu[c] Sigma[c] = invwishart.rvs(df=nu + nk, scale=V + dev.T @ dev, random_state=rng) else: mu[c] = solve_triangular(cholesky(Amu*np.eye(k), lower=True).T, rng.standard_normal(k), lower=False) Sigma[c] = invwishart.rvs(df=nu, scale=V, random_state=rng) Sinv[c] = np.linalg.inv(Sigma[c]) # 4. pi pi = rng.dirichlet(a + np.bincount(s, minlength=K)) # 5. Delta (matrix-normal GLS) [if Z] if use_Z: P = Ad.copy().astype(float); rhs = np.zeros((nz, k)) for i in range(N): ci = s[i] P += np.outer(Z[i], Z[i]) # shared design; precision in coef space below # vectorised GLS per output dim is complex with component-specific Sigma; use simple # approximation: regress (B - mu_s) on Z with identity weight (adequate for demo) resid = B - mu[s] U = np.linalg.inv(Z.T @ Z + Ad) M = U @ (Z.T @ resid) Delta = M + cholesky(U, lower=True) @ rng.standard_normal((nz, k)) # label fix: sort by mu on last coordinate (e.g. price) order = np.argsort(mu[:, -1]) mu, Sigma, Sinv, pi = mu[order], Sigma[order], Sinv[order], pi[order] s = np.argsort(order)[s] if r >= burn: j = r - burn MU[j] = mu; SIG[j] = Sigma; PI[j] = pi; LL[j] = llcur.sum(); Bm += B if use_Z: DD[j] = Delta out = dict(mu=MU, Sigma=SIG, pi=PI, loglike=LL, B=Bm / keep, accept=acc / (R * N), K=K, k=k) if use_Z: out['Delta'] = DD return out def hetero_density(out, coord, grid): """Posterior-mean fitted population density of coefficient `coord` (label-invariant).""" from scipy.stats import norm mu, Sig, pi = out['mu'], out['Sigma'], out['pi'] R, K = mu.shape[0], mu.shape[1] dens = np.zeros_like(grid, dtype=float) for r in range(R): for c in range(K): dens += pi[r, c] * norm.pdf(grid, mu[r, c, coord], np.sqrt(Sig[r, c, coord, coord])) return dens / R def log_marg_den_nr(loglike): """Newton-Raftery (harmonic-mean) log marginal density estimate (bayesm's metric; unstable).""" b = -np.asarray(loglike); mb = b.max() return float(-(mb + np.log(np.mean(np.exp(b - mb)))))