# Python贝叶斯统计 - 完整代码示例# 贝叶斯统计通过先验分布和数据更新得到后验分布import numpy as npimport matplotlib.pyplot as pltfrom scipy import stats# 第一部分贝叶斯推断基础 # 使用 PyMC 进行 MCMC 采样注意PyMC 可能需要安装# 由于环境可能没有 PyMC我们使用 scipy 和手动 MCMC 实现贝叶斯推断# 首先使用解析方法共轭先验然后演示 MCMC 采样# 1. 贝塔-二项分布模型共轭先验# 场景估计一枚硬币的正面概率# 先验Beta(α2, β2)假设硬币大致均匀# 数据10 次投掷中 7 次正面alpha_prior, beta_prior 2, 2 # Beta 先验参数n_trials, n_heads 10, 7 # 观测数据# 后验分布Beta(α 正面数, β 反面数)alpha_post alpha_prior n_headsbeta_post beta_prior (n_trials - n_heads)print(贝叶斯推断 - 硬币正面概率估计:)print(f先验分布: Beta({alpha_prior}, {beta_prior}))print(f观测数据: {n_trials} 次投掷, {n_heads} 次正面)print(f后验分布: Beta({alpha_post}, {beta_post}))# 后验均值和 95% 可信区间post_mean alpha_post / (alpha_post beta_post)post_std np.sqrt(alpha_post * beta_post / ((alpha_post beta_post)**2 * (alpha_post beta_post 1)))ci_lower stats.beta.ppf(0.025, alpha_post, beta_post)ci_upper stats.beta.ppf(0.975, alpha_post, beta_post)print(f后验均值: {post_mean:.4f})print(f后验标准差: {post_std:.4f})print(f95% 可信区间: [{ci_lower:.4f}, {ci_upper:.4f}])# 2. 网格近似法计算后验当共轭先验不可用时# 对于非共轭模型可以使用网格近似theta_grid np.linspace(0, 1, 1000) # 参数网格# 先验Beta(2, 2)prior_grid stats.beta.pdf(theta_grid, 2, 2)# 似然二项分布likelihood_grid stats.binom.pmf(n_heads, n_trials, theta_grid)# 后验未归一化posterior_grid prior_grid * likelihood_grid# 归一化posterior_grid / np.trapz(posterior_grid, theta_grid)# 从网格近似中抽样grid_samples np.random.choice(theta_grid, size10000, pposterior_grid/np.sum(posterior_grid))print(f\n网格近似后验均值: {np.mean(grid_samples):.4f})print(f网格近似后验 95% CI: [{np.percentile(grid_samples, 2.5):.4f}, {np.percentile(grid_samples, 97.5):.4f}])# 3. 马尔科夫链蒙特卡洛 (MCMC) 采样 - Metropolis-Hastings 算法def metropolis_hastings(log_posterior, init_val, n_samples10000, proposal_std0.1):Metropolis-Hastings MCMC 采样器samples [init_val]current init_valn_accepted 0for i in range(n_samples):# 从提议分布生成候选值正态分布candidate current np.random.normal(0, proposal_std)# 计算接受概率log_accept_ratio log_posterior(candidate) - log_posterior(current)# 接受或拒绝if np.log(np.random.random()) log_accept_ratio:current candidaten_accepted 1samples.append(current)return np.array(samples), n_accepted / n_samples# 定义对数后验非归一化def log_posterior_bernoulli(theta, data):Bernoulli 模型的对数后验Beta 先验if theta 0 or theta 1:return -np.inf# Beta(2,2) 先验log_prior stats.beta.logpdf(theta, 2, 2)# Bernoulli 似然n_heads np.sum(data)n_tails len(data) - n_headslog_lik n_heads * np.log(theta) n_tails * np.log(1 - theta)return log_prior log_lik# 生成数据coin_data np.array([1]*7 [0]*3) # 1正面, 0反面np.random.seed(42)mcmc_samples, accept_rate metropolis_hastings(lambda t: log_posterior_bernoulli(t, coin_data),init_val0.5, n_samples15000, proposal_std0.15)# 丢弃前 5000 个样本burn-inburn_in 5000mcmc_chain mcmc_samples[burn_in:]print(f\nMCMC (Metropolis-Hastings) 采样结果:)print(f接受率: {accept_rate:.3f} (理想范围 0.2-0.5))print(f后验均值: {np.mean(mcmc_chain):.4f})print(f后验中位数: {np.median(mcmc_chain):.4f})print(f95% 可信区间: [{np.percentile(mcmc_chain, 2.5):.4f}, {np.percentile(mcmc_chain, 97.5):.4f}])# 4. MCMC 诊断迹图和平稳性# 计算自相关以评估 MCMC 链的混合效果def autocorrelation(x, lag1):计算给定滞后的自相关return np.corrcoef(x[:-lag], x[lag:])[0, 1]lag1_autocorr autocorrelation(mcmc_chain)print(f滞后 1 自相关: {lag1_autocorr:.4f} (越小越好, 0.9 表示混合差))# 有效样本量估计 (ESS)def effective_sample_size(chain):估计 MCMC 链的有效样本量n len(chain)# 计算自相关和直到自相关变为负的滞后rho np.array([autocorrelation(chain, lagk) for k in range(1, min(50, n//2))])# 找到第一个负自相关的位置neg_idx np.where(rho 0)[0]cutoff neg_idx[0] if len(neg_idx) 0 else len(rho)ess n / (1 2 * np.sum(rho[:cutoff]))return essess effective_sample_size(mcmc_chain)print(f有效样本量: {ess:.0f} (实际样本: {len(mcmc_chain)}, 效率: {ess/len(mcmc_chain):.2f}))# 5. 贝叶斯线性回归# y α βx ε, ε ~ N(0, σ²)# 使用解析后验共轭先验np.random.seed(123)n_reg 50x_reg np.random.uniform(-3, 3, n_reg)true_alpha, true_beta, true_sigma 1.5, 2.0, 1.0y_reg true_alpha true_beta * x_reg np.random.normal(0, true_sigma, n_reg)# 手动最小二乘解X_design np.column_stack([np.ones(n_reg), x_reg])beta_ols np.linalg.inv(X_design.T X_design) X_design.T y_regy_pred_ols X_design beta_olssigma2_ols np.sum((y_reg - y_pred_ols)**2) / (n_reg - 2)# 贝叶斯后验无信息先验# 后验均值 OLS 估计, 后验方差 σ²(XX)⁻¹post_cov sigma2_ols * np.linalg.inv(X_design.T X_design)post_std_params np.sqrt(np.diag(post_cov))print(f\n贝叶斯线性回归 (无信息先验):)print(f真实参数: α{true_alpha}, β{true_beta})print(f后验均值: α{beta_ols[0]:.4f}, β{beta_ols[1]:.4f})print(f后验标准差: α{post_std_params[0]:.4f}, β{post_std_params[1]:.4f})# 后验预测分布x_new np.array([-2, 0, 2])X_new np.column_stack([np.ones(3), x_new])y_pred_new X_new beta_ols# 预测方差 σ² x_new * Cov(β) * x_newpred_var sigma2_ols np.diag(X_new post_cov X_new.T)pred_std np.sqrt(pred_var)print(f预测 (x{x_new}):)for i, x in enumerate(x_new):print(f x{x:.0f}: 预测 y{y_pred_new[i]:.3f} ± {1.96*pred_std[i]:.3f} (95% 预测区间))# 6. 贝叶斯模型比较WAIC 近似# 计算后验预测对数似然log_lik_matrix np.zeros((len(mcmc_chain), n_reg))for i in range(len(mcmc_chain)):alpha_i beta_ols[0] # 简化使用点估计beta_i beta_ols[1]y_pred_i alpha_i beta_i * x_reglog_lik_matrix[i] stats.norm.logpdf(y_reg, y_pred_i, np.sqrt(sigma2_ols))# WAIC 计算lppd np.sum(np.log(np.mean(np.exp(log_lik_matrix), axis0))) # 对数点预测密度p_waic np.sum(np.var(log_lik_matrix, axis0)) # 有效参数数waic -2 * (lppd - p_waic)print(f\n模型 WAIC: {waic:.2f} (越小越好))print(\n贝叶斯统计总结先验 数据 后验)print(MCMC 采样是处理复杂贝叶斯模型的核心计算工具)