采样实例

M-H采样实例

这里的例子是一个M-H采样的实例.

• 这里已经得到了目标分布, 即平稳分布, 是一个正态分布, 均值3, 标准差2.

• 马尔可夫链状态转移矩阵$Q$的条件转移概率是: 以$i$为均值, 标准差为1的正态分布, $j$所在位置的概率.

代码:

$\alpha(i,j) = min\{ \frac{\pi(j)Q(j,i)}{\pi(i)Q(i,j)},1\}$

由于这里的$Q$指的是正态分布, 因此$Q(i,j)=Q(j,i)$, 则有$\alpha(i,j) = min\{ \frac{\pi(j)}{\pi(i)},1\}$

import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt

def final_dist_prob(value):
    return norm.pdf(value, loc=3, scale=2)


num_samples = 1000
t = 0

pi = 0  # 初始状态值
samples = []

while t < num_samples:
    pi_star = norm.rvs(loc=pi, scale=1, size=1)[0]  # 从条件概率分布中采样
    alpha = min(1, final_dist_prob(pi_star) / (final_dist_prob(pi) + 10e-8))  # 计算接受率

    u = np.random.rand()
    if u < alpha:
        samples.append(pi_star)
        pi = pi_star
        t += 1

x_range = np.linspace(-3, 9, 200)
plt.plot(x_range, final_dist_prob(x_range))
plt.hist(samples, bins=50, normed=1, alpha=0.5)

plt.show()

最后采样的结果图:

image

可以看到采样得到的样本分布与已知的正态分布十分接近.

Gibbs二维采样实例

二维Gibbs采样实例