""" Mixture-of-normals probit (Geweke & Keane 1997) — from-scratch Gibbs (NumPy + scipy.norm). Binary probit with a SCALE MIXTURE OF NORMALS disturbance (zero-mean components, mixing on precision — symmetric, leptokurtic; the variant GK found best for PSID LFP): y_t = 1[ x_t' beta + eps_t > 0 ], eps_t ~ sum_j p_j N(0, h_j^{-1}) Extends Albert-Chib (single-normal probit) and robit (t = a scale mixture) to a flexible, semiparametric symmetric error. Gibbs blocks: 1. latent y*_t -- truncated normal (Albert-Chib), sd = 1/sqrt(h_{s_t}) 2. component indicator s_t -- categorical, P(j) ∝ p_j sqrt(h_j) exp(-.5 h_j r_t^2) 3. beta -- weighted Gaussian (GLS with weights h_{s_t}) 4. h_j -- Gamma (conjugate precision per component) 5. p -- Dirichlet Identification: rescale each sweep so Var(eps)=sum_j p_j/h_j = 1 (the probit scale, as in conventional probit); components sorted by h to fix labeling. Choice probability has closed form: P(y=1|x) = sum_j p_j * Phi( sqrt(h_j) * x'beta ). """ import numpy as np from scipy.stats import norm def _rtruncnorm(mu, sd, d, rng): """Truncated normal: d=1 -> (0, inf), d=0 -> (-inf, 0). Inverse-CDF, vectorized.""" thr = norm.cdf(-mu / sd) # Phi at the 0-boundary (standardized) a = np.where(d == 1, thr, 0.0) b = np.where(d == 1, 1.0, thr) u = np.clip(a + (b - a) * rng.random(len(mu)), 1e-12, 1 - 1e-12) return mu + sd * norm.ppf(u) def simulate_mixprobit(n, beta, error='t', df=4, seed=0): """Synthetic probit data with a chosen (variance-1) error. error: 'normal','t','mix2'.""" rng = np.random.default_rng(seed) beta = np.asarray(beta, float); k = len(beta) X = np.column_stack([np.ones(n)] + [rng.normal(size=n) for _ in range(k - 1)]) if error == 'normal': eps = rng.standard_normal(n) elif error == 't': eps = rng.standard_t(df, n) / np.sqrt(df / (df - 2)) # standardized to var 1 elif error == 'mix2': # leptokurtic 2-component z = rng.random(n) < 0.85 raw = np.where(z, rng.normal(0, 0.6, n), rng.normal(0, 2.2, n)) eps = raw / raw.std() ystar = X @ beta + eps return X, (ystar > 0).astype(float), eps def mixprobit_gibbs(y, X, m=2, R=6000, burn=2000, seed=0, prior_sd_beta=10.0, a0=2.0, b0=1.0, alpha_dir=1.0): """Scale-mixture-of-normals probit via Gibbs. m = number of mixture components.""" rng = np.random.default_rng(seed) n, k = X.shape Hb = np.eye(k) / prior_sd_beta ** 2 beta = np.zeros(k) h = np.linspace(0.6, 1.8, m) if m > 1 else np.array([1.0]) p = np.ones(m) / m s = rng.integers(0, m, n) keep = R - burn B = np.zeros((keep, k)); P = np.zeros((keep, m)); H = np.zeros((keep, m)); LL = np.zeros((keep, n)) for g in range(R): mu = X @ beta # 1. latent utilities (Albert-Chib), per-obs sd from its component yl = _rtruncnorm(mu, 1.0 / np.sqrt(h[s]), y, rng) # 2. component indicators r = yl - mu logpj = np.log(p)[None, :] + 0.5 * np.log(h)[None, :] - 0.5 * h[None, :] * r[:, None] ** 2 logpj -= logpj.max(1, keepdims=True) pp = np.exp(logpj); pp /= pp.sum(1, keepdims=True) s = (rng.random(n)[:, None] > np.cumsum(pp, 1)).sum(1).clip(0, m - 1) # 3. beta | . (weighted Gaussian) w = h[s] V = np.linalg.inv((X * w[:, None]).T @ X + Hb) beta = V @ ((X * w[:, None]).T @ yl) + np.linalg.cholesky(V) @ rng.standard_normal(k) # 4. h_j | . (Gamma) r = yl - X @ beta for j in range(m): mk = s == j; nj = int(mk.sum()) h[j] = rng.gamma(a0 + nj / 2.0, 1.0 / (b0 + 0.5 * np.sum(r[mk] ** 2))) if nj > 0 \ else rng.gamma(a0, 1.0 / b0) # 5. p | . (Dirichlet) p = rng.dirichlet(alpha_dir + np.bincount(s, minlength=m)) # identification: fix Var(eps)=1, then sort components by precision v = np.sum(p / h); beta = beta / np.sqrt(v); h = h * v o = np.argsort(h); h = h[o]; p = p[o]; s = np.argsort(o)[s] if g >= burn: gg = g - burn; B[gg] = beta; P[gg] = p; H[gg] = h P1 = np.clip((p[None, :] * norm.cdf(np.sqrt(h)[None, :] * (X @ beta)[:, None])).sum(1), 1e-12, 1 - 1e-12) LL[gg] = np.where(y == 1, np.log(P1), np.log(1 - P1)) return dict(beta=B, p=P, h=H, ll=LL, m=m) def mixprobit_full_gibbs(y, X, m=2, R=8000, burn=3000, seed=0, prior_sd_beta=10.0, a0=4.0, kappa0=0.05, alpha_dir=2.0): """FULL (location-scale) mixture-of-normals probit: eps ~ sum_j p_j N(alpha_j, h_j^{-1}). Unlike the scale mixture above (symmetric), the free component means alpha_j let the disturbance be ASYMMETRIC -- a flexible/skewed link, the right tool for skewed dose-response curves (e.g. the Bliss beetle data, where probit/logit underfit but cloglog fits). Identification each sweep: recentre the mixture to mean 0 (location -> intercept) AND rescale to Var(eps)=sum_j p_j(alpha_j^2 + 1/h_j)=1 (the probit scale); components sorted by mean. Choice prob: P(y=1|x) = sum_j p_j * Phi( sqrt(h_j) * (x'beta + alpha_j) ).""" rng = np.random.default_rng(seed) n, k = X.shape Hb = np.eye(k) / prior_sd_beta ** 2 b0 = 0.5 * (a0 - 1) # prior mean component var ~ 0.5 beta = np.zeros(k) alpha = np.linspace(-0.5, 0.5, m) if m > 1 else np.zeros(1) h = np.ones(m); p = np.ones(m) / m s = rng.integers(0, m, n) keep = R - burn B = np.zeros((keep, k)); P = np.zeros((keep, m)); A = np.zeros((keep, m)) H = np.zeros((keep, m)); LL = np.zeros((keep, n)) for g in range(R): mu = X @ beta yl = _rtruncnorm(mu + alpha[s], 1.0 / np.sqrt(h[s]), y, rng) r = yl - mu logpj = np.log(p)[None, :] + 0.5 * np.log(h)[None, :] - 0.5 * h[None, :] * (r[:, None] - alpha[None, :]) ** 2 logpj -= logpj.max(1, keepdims=True) pp = np.exp(logpj); pp /= pp.sum(1, keepdims=True) s = (rng.random(n)[:, None] > np.cumsum(pp, 1)).sum(1).clip(0, m - 1) # beta | . (weighted GLS on yl - alpha_s) w = h[s] V = np.linalg.inv((X * w[:, None]).T @ X + Hb) beta = V @ ((X * w[:, None]).T @ (yl - alpha[s])) + np.linalg.cholesky(V) @ rng.standard_normal(k) # (alpha_j, h_j) | . (Normal-Inverse-Gamma) r = yl - X @ beta for j in range(m): mk = s == j; nj = int(mk.sum()); rj = r[mk] kn = kappa0 + nj; mn = rj.sum() / kn; an = a0 + nj / 2.0 bn = b0 + 0.5 * (np.sum((rj - rj.mean()) ** 2) if nj > 0 else 0.0) \ + 0.5 * (kappa0 * nj / kn) * (rj.mean() if nj > 0 else 0.0) ** 2 h[j] = rng.gamma(an, 1.0 / bn) alpha[j] = mn + np.sqrt(1.0 / (kn * h[j])) * rng.standard_normal() p = rng.dirichlet(alpha_dir + np.bincount(s, minlength=m)) # identification: recentre to mean 0, rescale to Var=1, sort by mean mbar = np.sum(p * alpha); beta[0] += mbar; alpha -= mbar sv = np.sqrt(np.sum(p * (alpha ** 2 + 1.0 / h))) beta = beta / sv; alpha = alpha / sv; h = h * sv ** 2 o = np.argsort(alpha); alpha = alpha[o]; h = h[o]; p = p[o]; s = np.argsort(o)[s] if g >= burn: gg = g - burn; B[gg] = beta; P[gg] = p; A[gg] = alpha; H[gg] = h P1 = np.clip((p[None, :] * norm.cdf(np.sqrt(h)[None, :] * ((X @ beta)[:, None] + alpha[None, :]))).sum(1), 1e-12, 1 - 1e-12) LL[gg] = np.where(y == 1, np.log(P1), np.log(1 - P1)) return dict(beta=B, p=P, alpha=A, h=H, ll=LL, m=m) def prob_curve(out, X): """Posterior-mean P(y=1|X) for a full-mixture-probit fit, evaluated at rows of X.""" b = out['beta'].mean(0); p = out['p'].mean(0); a = out['alpha'].mean(0); h = out['h'].mean(0) eta = X @ b return (p[None, :] * norm.cdf(np.sqrt(h)[None, :] * (eta[:, None] + a[None, :]))).sum(1) def waic(out): """WAIC from pointwise log-likelihood draws (lower = better).""" ll = out['ll']; mx = ll.max(0) lppd = (np.log(np.exp(ll - mx).mean(0)) + mx) pw = ll.var(0) return float(-2 * (lppd.sum() - pw.sum())) def error_density(out, grid): """Posterior-mean fitted disturbance density on `grid` (the estimated shock distribution).""" p = out['p'].mean(0); h = out['h'].mean(0) return sum(p[j] * norm.pdf(grid, 0, 1 / np.sqrt(h[j])) for j in range(out['m'])) def error_density_full(out, grid): """Fitted disturbance density for a FULL-mixture fit (components have free means alpha_j).""" p = out['p'].mean(0); a = out['alpha'].mean(0); h = out['h'].mean(0) return sum(p[j] * norm.pdf(grid, a[j], 1 / np.sqrt(h[j])) for j in range(out['m']))