用python实现 MCMC的实例分析

MCMC（Markov Chain Monte Carlo）是一种用于模拟复杂系统的方法。它通过构造一个马尔可夫链，来实现从目标分布中抽样的过程。

下面我们通过一个简单的例子来演示如何用Python实现MCMC。

假设我们有一个未知的概率分布 $p(x)$，我们希望从中抽取样本。我们可以使用MCMC方法来构造一个马尔可夫链，从而实现从该分布中抽样。这里我们选择Metropolis-Hastings算法来实现MCMC。

具体实现过程如下：

首先，我们需要定义一个目标概率分布 $p(x)$。在这个例子中，我们选择一个标准的正态分布：

import numpy as np

def target_distribution(x):
    return np.exp(-0.5 * x**2) / np.sqrt(2*np.pi)

接着，我们需要构造一个马尔可夫链。在Metropolis-Hastings算法中，我们需要定义一个转移概率分布 $q(x'|x)$，该分布用于从当前状态 $x$ 转移到新状态 $x'$。这里我们选择一个标准正态分布：

def proposal_distribution(x):
    return np.random.normal(x, 1)

然后，我们需要实现Metropolis-Hastings算法。具体实现过程如下：

def metropolis_hastings(n_iterations, initial_value):
    chain = [initial_value]
    for i in range(n_iterations):
        current_value = chain[-1]
        proposed_value = proposal_distribution(current_value)
        acceptance_probability = min(1, target_distribution(proposed_value) / target_distribution(current_value))
        if np.random.uniform(0, 1) < acceptance_probability:
            chain.append(proposed_value)
        else:
            chain.append(current_value)
    return chain

最后，我们可以运行MCMC算法并绘制抽样结果：

import matplotlib.pyplot as plt

chain = metropolis_hastings(10000, 0)
plt.hist(chain, bins=100, density=True)
x = np.linspace(-5, 5, 1000)
plt.plot(x, target_distribution(x))
plt.show()

该程序将生成一个直方图，其中包含从目标分布中抽取的样本，同时还绘制了目标分布的曲线。可以看到，MCMC算法可以很好地抽取目标分布中的样本。

完整代码如下：