贝叶斯实证推断：先验、似然、后验与 MCMC

*==============================================================*
* 贝叶斯 OLS + 贝叶斯 Probit：bayesmh 完整流程 + 收敛诊断
*==============================================================*
clear all
set more off
set seed 20260911

* --- 1. 模拟数据：Y = 1 + 0.5*X + eps ---
set obs 500
gen X = rnormal()
gen eps = rnormal()
gen Y = 1 + 0.5*X + eps
gen D = (1 + 0.5*X + eps > 0)        // Probit 的二值结果

* --- 2. 贝叶斯 OLS：bayesmh ---
* 语法：bayesmh 似然, likelihood(...) prior(...)
* 无信息正态先验 beta ~ N(0, 100)；sigma ~ Half-Normal
bayesmh Y, likelihood(regress({y} {beta:X _cons}, var({var}))) ///
    prior({beta:X _cons}:, flat)            ///  无信息先验
    prior({var}:, igamma(1e-3,1e-3))        ///  方差逆伽马
    chains(4)                                ///  4 条链
    montecarlo(saving("bayes_ols.ster", replace) seed(20260911)) ///
    nchains(4) burnin(1000) mcmcsize(2000)
* 收敛诊断：bayesstats 内置
bayesstats summary {beta:X _cons}
bayesstats ic                            // DIC / WAIC 模型比较
* trace plot
bayesgraph diagnostics {beta:X}

* --- 3. 贝叶斯 Probit：Albert-Chib 数据增广 ---
* bayesmh 自动处理 probit 似然
bayesmh D, likelihood(probit({y} {b:X _cons})) ///
    prior({b:X _cons}:, flat) ///
    chains(4) burnin(1000) mcmcsize(2000)
bayesstats summary {b:X}
* 解释：后验 P(b_X > 0 | data) 直接可读，无需 p 值

* --- 4. 收敛诊断（必做） ---
* R-hat < 1.01；ESS > 400；trace plot 平稳
* bayesstats ess  输出 effective sample size
* bayesstats grubin  Gelman-Rubin R-hat
bayesstats ess {beta:X _cons}
bayesstats grubin, estimates("bayes_ols.ster")

# ==============================================================
# 贝叶斯 Probit：PyMC + NUTS + ArviZ 收敛诊断
# ==============================================================
import numpy as np
import pymc as pm
import arviz as az
import matplotlib.pyplot as plt

rng = np.random.default_rng(20260911)

# --- 1. 模拟数据：与 Stata 一致 ---
n = 500
X = rng.normal(size=n)
eps = rng.normal(size=n)
Y_latent = 1 + 0.5*X + eps
D = (Y_latent > 0).astype(int)

# --- 2. 定义贝叶斯 Probit 模型 ---
with pm.Model() as probit_model:
    # 先验：无信息正态
    beta_0 = pm.Normal('beta_0', mu=0, sigma=10)
    beta_X = pm.Normal('beta_X', mu=0, sigma=10)
    # 潜变量与观测（PyMC 内置 Bernoulli + probit link）
    p = pm.math.invprobit(beta_0 + beta_X*X)
    Y_obs = pm.Bernoulli('Y_obs', p=p, observed=D)

# --- 3. NUTS 采样：4 链，2000 迭代，1000 warmup ---
with probit_model:
    trace = pm.sample(draws=2000, tune=1000, chains=4,
                      random_seed=20260911, target_accept=0.9)

# --- 4. 收敛诊断 ---
print(az.summary(trace, var_names=['beta_0', 'beta_X']))
# 关键指标：
#   r_hat < 1.01  ✓
#   ess_bulk > 400 ✓
#   ess_tail > 400 ✓
az.plot_trace(trace, var_names=['beta_0', 'beta_X']); plt.tight_layout()
az.plot_posterior(trace, var_names=['beta_X'], ref_val=0); plt.show()

# --- 5. 后验推断：P(beta_X > 0 | data) ---
post_beta_X = trace.posterior['beta_X'].values.flatten()
p_positive = (post_beta_X > 0).mean()
print(f"P(beta_X > 0 | data) = {p_positive:.3f}  (true = 0.5)")
print(f"95% 可信区间: {np.percentile(post_beta_X, [2.5, 97.5])}")

# ==============================================================
# 贝叶斯层次模型：组级随机效应
# ==============================================================
import numpy as np, pymc as pm, arviz as az
rng = np.random.default_rng(20260911)

J, n_j = 8, 50
alpha_true = rng.normal(0, 1.5, J)                # 组级截距
X = rng.normal(size=J*n_j)
group = np.repeat(np.arange(J), n_j)
Y = alpha_true[group] + 0.5*X + rng.normal(scale=0.5, size=J*n_j)

with pm.Model() as hier:
    mu_a   = pm.Normal('mu_a', 0, 5)
    sigma_a = pm.HalfNormal('sigma_a', 1)
    alpha_j = pm.Normal('alpha_j', mu=mu_a, sigma=sigma_a, shape=J)
    beta   = pm.Normal('beta', 0, 5)
    sigma  = pm.HalfNormal('sigma', 1)
    pm.Normal('Y_obs', mu=alpha_j[group] + beta*X, sigma=sigma, observed=Y)

with hier:
    trace_h = pm.sample(draws=2000, tune=1000, chains=4,
                        random_seed=20260911, target_accept=0.9)
print(az.summary(trace_h, var_names=['mu_a','sigma_a','beta']))