MCMC算法基本原理

### MCMC算法基本原理

MCMC（Markov Chain Monte Carlo）是一种用于从复杂概率分布中抽样的强大工具。该方法结合了两种核心概念：蒙特卡罗模拟和马尔科夫链。

#### 蒙特卡罗方法简介

蒙特卡罗方法通过随机抽样来估计数值解，特别适用于难以解析求解的情况。这种方法依赖于大量样本的统计特性来逼近真实值[^1]。对于某些复杂的积分计算或高维空间上的概率密度函数评估，直接获得精确结果往往非常困难甚至不可能；此时利用蒙特卡罗技术可以有效地近似这些量。

#### 马尔科夫链基础

马尔科夫链描述了一种状态转移过程，在这种过程中下一个时刻的状态仅取决于当前所处的状态而与其他历史无关。具体来说就是给定一系列可能的状态集合S={s₁,s₂,...,sn}以及相应的转移矩阵P=(pij)，其中 pij 表示从 si 到 sj 的一步转移概率，则整个系统的演化遵循如下规则：

\[ P(X_{t+1}=s_j|X_t=s_i,X_{t-1},...,X_0)=P(X_{t+1}=s_j|X_t=s_i) \]

当满足特定条件时（如遍历性和周期性），经过足够长时间运行之后，无论初始位置在哪一点出发，最终都会收敛到平稳分布π(s)。

#### Metropolis-Hastings算法介绍

作为MCMC家族中最著名的成员之一，Metropolis-Hastings提供了一个通用框架用来构建能够渐进地接近目标分布π(x)的马尔可夫链。其工作流程大致如下：

1. 初始化参数θ₀；
2. 对每一轮迭代k=1,2,...
- 提议新候选点θ* ~ q(· | θ_(k−1))；
- 计算接受率α=min{1,r}，其中 \( r=\frac{\pi(\theta^{*})q(\theta^{(k-1)}|\theta^{*})}{\pi(\theta^{(k-1)})q(\theta^{*}|\theta^{(k-1)})}\) ；
- 抛硬币决定是否采纳建议：
- 若u~U[0,1]<α则令θ(k)=θ*，
- 否则保持不变θ(k)=θ(k−1)。

重复上述步骤直至收集足够的样本数量为止[^2]。

import numpy as np
from scipy.stats import norm

def metropolis_hastings(target_dist, proposal_func, n_samples, init_state):
    samples = []
    current_state = init_state
    
    for _ in range(n_samples):
        proposed_state = proposal_func(current_state)

        acceptance_ratio = min(
            target_dist(proposed_state)/target_dist(current_state), 
            1
        )
        
        u = np.random.uniform()
        if u < acceptance_ratio:
            current_state = proposed_state
            
        samples.append(current_state)
    
    return np.array(samples)


# Example usage with normal distribution
mu, sigma = 0., 1.
proposal_stddev = .5


def gaussian_target(x):
    return norm.pdf(x, loc=mu, scale=sigma)


proposals = lambda x: np.random.normal(loc=x, scale=proposal_stddev)

samples = metropolis_hastings(gaussian_target, proposals, 10000, mu)