蒙特卡洛积分法在贝叶斯推断中的后验期望估计

题目描述
在贝叶斯统计中，后验分布的期望是关键参数（如后验均值）。设参数θ的后验分布为π(θ|X)，其中X是观测数据。目标计算后验期望E[g(θ)|X] = ∫ g(θ)π(θ|X)dθ，其中g(θ)是任意函数（如θ²表示后验二阶矩）。若后验分布形式复杂（如非共轭先验），解析积分不可行。需通过蒙特卡洛积分法，从后验分布中抽取样本{θ₁, θ₂, ..., θ_N}，用样本均值近似积分。

解题过程

1. 问题转化

   • 后验期望E[g(θ)|X]本质是函数g(θ)关于后验分布的积分。蒙特卡洛方法的核心思想是：若能从π(θ|X)中独立抽样，则样本均值收敛于期望值。
   • 关键步骤：
     （1）生成后验分布样本：使用MCMC（如Metropolis-Hastings）或直接抽样法。
     （2）计算g(θ)在样本上的取值并求平均。

2. 抽样方法选择

   • 若后验分布有标准形式（如共轭先验），可直接抽样（如Beta分布用numpy.random.beta）。
   • 对于复杂后验，常用Metropolis-Hastings算法：
     • 提议分布q(θ*|θ)选择对称分布（如正态分布），简化接受率计算：
       α = min(1, π(θ*|X)/π(θ|X))。
     • 迭代生成样本链，去除燃烧期（burn-in）后保留有效样本。

3. 蒙特卡洛估计实现

   • 设抽取N个独立样本{θ₁, ..., θ_N}，估计量为：
     Ê[g(θ)|X] = (1/N) Σᵢ g(θᵢ)。
   • 由大数定律，当N→∞时，估计值依概率收敛到真实期望。
   • 误差估计：计算样本标准差s = √[Σᵢ (g(θᵢ) - Ê)²/(N-1)]，则标准误为s/√N。

4. 实际案例演示

   • 例：估计正态分布均值μ的后验期望。假设先验μ~N(0,1)，观测数据X={x₁,...,xₙ}~N(μ,1)，则后验π(μ|X)为N(∑xᵢ/(n+1), 1/(n+1))。
     • 直接抽样：从N(∑xᵢ/(n+1), 1/(n+1))生成样本。
     • 若g(μ)=μ²，则计算样本的平方均值。
   • 复杂案例：用MCMC抽样逻辑回归后验，计算系数期望。

5. 收敛性验证

   • 绘制样本路径图检查平稳性。
   • 计算多链的Gelman-Rubin统计量（若<1.05认为收敛）。
   • 增加样本量N直至标准误小于预设容差（如0.001）。

6. 优化技巧

   • 方差缩减：若后验抽样效率低，可采用重要性采样，引入辅助分布q(θ)近似π(θ|X)，则估计量为(1/N) Σᵢ g(θᵢ)π(θᵢ|X)/q(θᵢ)。
   • 针对高维后验，使用Hamiltonian Monte Carlo减少自相关性。

通过以上步骤，蒙特卡洛积分将复杂后验期望转化为可计算的样本均值，成为贝叶斯计算的核心工具。