第19章 马尔可夫链蒙特卡罗法

习题19.1

  用蒙特卡罗积分法求

解答：

解答思路：

1. 给出蒙特卡罗积分计算方法。
2. 将被积分的函数分解，得到其数学期望表示形式。
3. 自编程实现，使用样本均值求近似计算积分。

解答步骤

第1步：蒙特卡罗积分计算方法

  根据书中第19.1.3节的关于蒙特卡罗积分法的描述：

  假设有一个函数 $h(x)$，目标是计算该函数的积分

$$$$

于是函数 $h(x)$ 的积分可以表示为函数 $f(x)$ 关于概率密度函数 $p(x)$ 的数学期望。实际上, 给定一个概率密度函数 $p(x)$, 只要取 $\displaystyle f(x)=\frac{h(x)}{p(x)}$，就可得上式。就是说, 任何一个函数的积分都可以表示为某一个函数的数学期望的形式。而函数的数学期望又可以通过函数的样本均值估计。于是，就可以利用样本均值来近似计算积分。

$$$$

f(x) = \frac{h(x)}{p(x)} = \frac{\displaystyle x^2 \exp \left( - \frac{x^2}{2} \right)}{\displaystyle \frac{1}{\sqrt{2 \pi}} \exp \left(- \frac{x^2}{2} \right)} = \sqrt{2 \pi} x^{2}

\int_{\mathcal{X}} x^2 \exp \left( - \frac{x^2}{2} \right) \mathrm{d} x = \int_{\mathcal{X}} f(x) p(x) \mathrm{d} x = E_{p(x)}[f(x)]

具有马尔可夫性的随机序列$X = \{ X_0, X_1, \cdots, X_t, \cdots \}$称为马尔可夫链，或马尔可夫过程。条件概率分布$P(X_t | X_{t-1})$称为马尔可夫链的转移概率分布。转移概率分布决定了马尔可夫链的特性。

  根据书中第19.2.2节的状态转移概率：

也就是说，时刻0从状态 $j$ 出发, 时刻 $t$ 到达状态 $i$ 的概率大于0，则称此马尔可夫链 $X$ 是不可约的（irreducible），否则称马尔可夫链是可约的（reducible）。

第3步：非周期的定义

  根据书中第19章的定义19.4的非周期的定义：

定义19.4（非周期） 设有马尔可夫链 $X = \{X_0, X_1, \cdots, X_t, \cdots \}$，状态空间为 $\mathcal{S}$，对于任意状态 $i \in \mathcal{S}$，如果时刻0从状态 $i$ 出发，$t$ 时刻返回状态的所有时间长 $\{ t: P( X_t = i | X_0 = i ) > 0 \}$ 的最大公约数是1，则称此马尔可夫链 $X$ 是非周期的（aperiodic），否则称马尔可夫链是周期的（periodic）。

第4步：证明如果马尔可夫链是不可约的，则所有状态周期都相同

  如果该马尔可夫链是不可约的，则对于 $\forall i, j \in S$，必然存在时刻 $m, n (m > 0, n > 0)$ 满足

根据状态转移概率的定义，有 $p_{ij}^{m} = P ( X_m = i | X_0 = j), p_{ji}^{n} = P ( X_n = j | X_0 = i)$

  可得

  则必然存在 $s$ 使得 $P_{j j}^{s}>0$，且满足

  假设状态 $i, j$ 的周期分别是 $d(i), d(j)$，由上两式可知

  所以 $s \bmod d(i) = 0$。

  又因为 $P_{j j}^{s} > 0$，所以 $s \bmod d(j) = 0$。因此 $d(i) \bmod d(j) = 0$。同理可得 $d(j) \bmod d(i) = 0$

  故 $d(i)==d(j)$，可得如果马尔可夫链是不可约的，所有状态周期都是相同的。

第5步：证明原命题

  在不可约的马尔可夫链中如果有一个状态的是非周期的，即其周期为1，则其他所有状态的周期也为1，所以该马尔可夫链是非周期的。原命题得证。

习题19.3

  验证具有以下转移概率矩阵的马尔可夫链是可约的，但是非周期的。

解答：

解答思路：

1. 给出平稳分布的定义
2. 计算该马尔可夫链的平稳分布
3. 验证该马尔可夫链是可约的
4. 验证该马尔可夫链是非周期的

解答步骤：

第1步：平稳分布的定义

  根据书中第19章的定义19.2的平稳分布的定义：

  定义19.2（平稳分布） 设有马尔可夫链 $X=\{ X_0, X_1, \cdots, X_t, \cdots \}$，其状态空间为 $\mathcal{S}$，转移概率矩阵为 $P = (p_{i j})$，如果存在状态空间 $\mathcal{S}$ 上的一个分布

$$$$

\pi_1 \ \pi_2 \ \vdots \end{array} \right]

则称 $\pi$ 为马尔可夫链 $X = \{X_0, X_1, \cdots, X_t, \cdots \}$ 的平稳分布。

  根据书中第19章的引理19.1：

  引理19.1 给定一个马尔可夫链 $X = \{X_0, X_1, \cdots, X_t, \cdots \}$，状态空间为 $\mathcal{S}$，转移概率矩阵为 $P = (p_{i j})$，则分布 $\pi = (\pi_1, \pi_2, \cdots )^T$ 为 $X$ 的平稳分布的充分必要条件是 $\pi = (\pi_1, \pi_2, \cdots )^T$ 是下列方程组的解：

$$$$

\displaystyle x_i = \sum_j p_{i j} x_j, \quad i = 1,2, \cdots \ x_i \geqslant 0, \quad i = 1, 2, \cdots \ \displaystyle \sum_i x_i = 1 \end{array}

\begin{array}{l} \displaystyle x_1 = \frac{1}{2} x_1 + \frac{1}{2} x_2 \ \displaystyle x_2 = \frac{1}{2} x_1 + \frac{1}{2} x_3 \ \displaystyle x_3 = \frac{1}{2} x_2 \ \displaystyle x_4 = \frac{1}{2} x_3 + 1 x_4 \ x_1 + x_2 + x_3 + x_4 = 1 \ x_i \geqslant 0, \quad i = 1, 2, 3, 4 \end{array}

\pi = (0, 0, 0, 1)^T

P=\left[\begin{array}{cccc} 0 & 1 / 2 & 0 & 0 \ 1 & 0 & 1 / 2 & 0 \ 0 & 1 / 2 & 0 & 1 \ 0 & 0 & 1 / 2 & 0 \end{array}\right]

\begin{array}{l} \displaystyle x_1 = \frac{1}{2} x_2 \ \displaystyle x_2 = 1 x_1 + \frac{1}{2} x_3 \ \displaystyle x_3 = \frac{1}{2} x_2 + 1 x_4 \ \displaystyle x_4 = \frac{1}{2} x_3 \ x_1 + x_2 + x_3 + x_4 = 1 \ x_i \geqslant 0, \quad i = 1, 2, 3, 4 \end{array}

\pi=(\frac{1}{6}, \frac{1}{3}, \frac{1}{3}, \frac{1}{6})^T

或简写为

$$$$

P = \left [ \begin{array}{cccc} 1 & 0 & \cdots & 0 \ 0 & 1 & \cdots & 0 \ \vdots & \vdots & \ddots & \vdots \ 0 & 0 & 0 & 1 \end{array} \right ]

  （c）从区间$(0,1)$中按均匀分布随机抽取一个数$u$。
     若 $u \leqslant \alpha (x, x')$，则状态 $x_i = x'$；否则，状态 $x_i = x$。
（3）得到样本集合 $\{ x_{m+1}, x_{m+2}, \cdots, x_n \}$
计算

$$$$

其中 $x_{j}$ 表示随机变量 $x$ 的第 $j$ 个分量，$j=1,2, \cdots, k$，而 $x^{(i)}$ 表示马尔可夫链在时刻 $i$ 的状态

$$$$

其中分量 $1, 2, \cdots, j-1$ 已经更新。然后，按照接受概率

$$$$

{p \left( x{'}j^{(i)} | x{-j}^{(i)} \right) q \left( x{'}j^{(i)}, x_j^{(i-1)} | x{-j}^{(i)} \right)} {p \left( x_j^{(i-1)} | x_{-j}^{(i)} \right) q \left( x_j^{(i-1)}, x{'}j^{(i)} | x{-j}^{(i)} \right)} \right }

3. 推导单分量 Metropolis-Hastings 算法

单分量 Metropolis-Hastings 算法
输入：抽样的目标分布的密度函数 $p(x)$，函数 $f(x)$，其中 $x = \left( x_1, x_2, \cdots, x_k \right)^T$；
输出：$p(x)$ 的随机样本 $\{ x^{(m+1)}, x^{(m+2)}, \cdots, x^{(n)} \}$，函数样本均值 $f_{m n}$；
参数：收敛步数 $m$, 迭代步数 $n$。
（1）任意选择一个初始值 $x_0$
（2）对 $i = 1, 2, \cdots, n$ 循环执行
  设在第 $(i-1)$ 次迭代结束时分量 $x_{j}$ 的取值为 $x_{j}^{(i-1)}$
  （a）对 $j=1,2, \cdots, k$ 循环执行
    （i）在第 $i$ 次迭代的第 $j$ 步，由建议分布 $q (x_j^{(i-1)}, x_j | x_{-j}^{(i)})$ 抽样产生分量 $x_j$ 的候选值 $x{'}_j^{(i)}$，
        这里 $x_{-j}^{(i)}$ 表示在第 $i$ 次迭代的第 $(j-1)$ 步后的 $x^{(i)}$ 除去 $x_j^{(i-1)}$ 的所有值，即

    （ii）计算接受概率

    （iii）从区间 $(0,1)$ 中按均匀分布随机抽取一个数 $u$。
        若 $u \leqslant \alpha \left( x_j^{(i-1)}, x{'}_j^{(i)} | x_{-j}^{(i)} \right )$，则状态 $x_j^{(i)} = x{'}_j^{(i)}$；否则，状态 $x_j^{(i)} = x_j^{(i-1)}$。
（3）得到样本集合 $\{ x^{(m+1)}, x^{(m+2)}, \cdots, x^{(n)} \}$
计算

习题19.7

  假设进行伯努利实验，后验概率为 $P(\theta | y)$，其中变量 $y \in \{ 0, 1 \}$ 表示实验可能的结果，变量 $\theta$ 表示结果为 1 的概率。再假设先验概率 $P(\theta)$ 遵循 Beta 分布 $B(\alpha, \beta)$，其中 $\alpha = 1, \beta = 1$；似然函数 $P(y | \theta)$ 遵循二项分布$\text{Bin}(n, k, \theta)$，其中 $n = 10, k = 4$，即实验进行 10 次其中结果为 1 的次数为 4。试用 Metropolis-Hastings 算法求后验概率分布 $P(\theta | y) \propto P(\theta) P(y | \theta)$ 的均值和方差。（提示：可采用 Metropolis 选择，即假设建议分布是对称的。）

解答：

解题思路

1. 给出 Metropolis-Hastings 算法
2. 写出目标分布和建议分布
3. 自编程实现 Metropolis-Hastings 算法求解

解答步骤

第1步：Metropolis-Hastings 算法

  Metropolis-Hastings 算法的步骤详见书中第373~374页算法19.2。

第2步：写出目标分布和建议分布

  根据题意，可知后验概率 $P(\theta) \propto B(1, 1)$，后验函数 $P(y | \theta) \propto \text{Bin}(10, 4, \theta)$

  则后验概率分布

  可取建议分布为 $B(1, 1)$，接受分布为 $\text{Bin}(10, 4, \theta)$

第3步：自编程实现 Metropolis-Hastings 算法求解

import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import beta, binom

class MetropolisHastings:
    def __init__(self, proposal_dist, accepted_dist, m=1e4, n=1e5):
        """
        Metropolis Hastings

        :param proposal_dist: 建议分布
        :param accepted_dist: 接受分布
        :param m: 收敛步数
        :param n: 迭代步数
        """
        self.proposal_dist = proposal_dist
        self.accepted_dist = accepted_dist
        self.m = m
        self.n = n

    @staticmethod
    def __calc_acceptance_ratio(q, p, x, x_prime):
        """
        计算接受概率

        :param q: 建议分布
        :param p: 接受分布
        :param x: 上一状态
        :param x_prime: 候选状态
        """
        prob_1 = p.prob(x_prime) * q.joint_prob(x_prime, x)
        prob_2 = p.prob(x) * q.joint_prob(x, x_prime)
        alpha = np.min((1., prob_1 / prob_2))
        return alpha

    def solve(self):
        """
        Metropolis Hastings 算法求解
        """
        all_samples = np.zeros(self.n)
        # (1) 任意选择一个初始值
        x_0 = np.random.random()
        for i in range(int(self.n)):
            x = x_0 if i == 0 else all_samples[i - 1]
            # (2.a) 从建议分布中抽样选取
            x_prime = self.proposal_dist.sample()
            # (2.b) 计算接受概率
            alpha = self.__calc_acceptance_ratio(self.proposal_dist, self.accepted_dist, x, x_prime)
            # (2.c) 从区间 (0,1) 中按均匀分布随机抽取一个数 u
            u = np.random.uniform(0, 1)
            # 根据 u <= alpha，选择 x 或 x_prime 进行赋值
            if u <= alpha:
                all_samples[i] = x_prime
            else:
                all_samples[i] = x

        # (3) 随机样本集合
        samples = all_samples[self.m:]
        # 函数样本均值
        dist_mean = samples.mean()
        # 函数样本方差
        dist_var = samples.var()
        return samples[self.m:], dist_mean, dist_var

    @staticmethod
    def visualize(samples, bins=50):
        """
        可视化展示
        :param samples: 抽取的随机样本集合
        :param bins: 频率直方图的分组个数
        """
        fig, ax = plt.subplots()
        ax.set_title('Metropolis Hastings')
        ax.hist(samples, bins, alpha=0.7, label='Samples Distribution')
        ax.set_xlim(0, 1)
        ax.legend()
        plt.show()

class ProposalDistribution:
    """
    建议分布
    """

    @staticmethod
    def sample():
        """
        从建议分布中抽取一个样本
        """
        # B(1,1)
        return beta.rvs(1, 1, size=1)

    @staticmethod
    def prob(x):
        """
        P(X = x) 的概率
        """
        return beta.pdf(x, 1, 1)

    def joint_prob(self, x_1, x_2):
        """
        P(X = x_1, Y = x_2) 的联合概率
        """
        return self.prob(x_1) * self.prob(x_2)

class AcceptedDistribution:
    """
    接受分布
    """

    @staticmethod
    def prob(x):
        """
        P(X = x) 的概率
        """
        # Bin(4, 10)
        return binom.pmf(4, 10, x)

import warnings
warnings.filterwarnings("ignore")

# 收敛步数
m = 1000
# 迭代步数
n = 10000

# 建议分布
proposal_dist = ProposalDistribution()
# 接受分布
accepted_dist = AcceptedDistribution()

metropolis_hastings = MetropolisHastings(proposal_dist, accepted_dist, m, n)

# 使用 Metropolis-Hastings 算法进行求解
samples, dist_mean, dist_var = metropolis_hastings.solve()
print("均值:", dist_mean)
print("方差:", dist_var)

# 对结果进行可视化
metropolis_hastings.visualize(samples, bins=20)

均值: 0.4214185307283342
方差: 0.019466091658555083

习题19.8

  设某试验可能有五种结果，其出现的概率分别为

模型含有两个参数 $\theta$ 和 $\eta$，都介于 0 和 1 之间。现有 22 次试验结果的观测值为

其中 $y_i$ 表示 22 次试验中第 $i$ 个结果出现的次数，$i = 1,2,\cdots,5$。试用吉布斯抽样估计参数 $\theta$ 和 $\eta$ 的均值和方差。

解答

解题思路

1. 给出吉布斯抽样算法
2. 定义目标概率的密度函数和各参数的满条件分布
3. 自编程实现吉布斯抽样算法进行求解

解答步骤

第1步：吉布斯抽样算法

  根据书中第19.5.2节的吉布斯抽样算法的描述：

算法19.3(吉布斯抽样)
输入：目标概率分布的密度函数 $p(x)$，函数 $f(x)$；
输出：$p(x)$ 的随机样本 $x_{m+1}, x_{m+2}, \cdots, x_{n}$，函数样本均值 $f_{m n}$；
参数：收敛步数 $m$，迭代步数 $n$。
（1）初始化。给出初始样本 $x^{(0)} = \left( x_1^{(0)}, x_2^{(0)}, \cdots, x_k^{(0)} \right)^T$。
（2）对 $i$ 循环执行
设第 $(i-1)$ 次迭代结束时的样本为 $x^{(i-1)} = \left( x_1^{(i-1)}, x_2^{(i-1)}, \cdots, x_k^{(i-1)} \right)^T$，则第 $i$ 次迭代进行如下几步操作：

$$$$

\begin{array}{l} \text {(1) 由满条件分布 } p \left( x_1 | x_2^{(i-1)}, \cdots, x_k^{(i-1)} \right) \text { 抽取 } x_1^{(i)} \ \quad \quad \vdots \ \text {(j) 由满条件分布 } p \left( x_j | x_1^{(i)}, \cdots, x_{j-1}^{(i)}, x_{j+1}^{(i-1)}, \cdots, x_k^{(i-1)} \right) \text{ 抽取 } x_j^{(i)} \ \quad \quad \vdots \ \text {(k) 由满条件分布 } p \left( x_k | x_1^{(i)}, \cdots, x_{k-1}^{(i)} \right) \text { 抽取 } x_k^{(i)} \end{array} \right.

（4）计算

$$$$

\begin{aligned} P &= (p_1, p_2, p_3, p_4, p_5) \ &= (\frac{\theta}{4}+\frac{1}{8},\ \frac{\theta}{4}, \ \frac{\eta}{4}, \ \frac{\eta}{4}+\frac{3}{8}, \ \frac{1}{2}(1-\theta-\eta)) \end{aligned}

\begin{array}{l} \displaystyle \frac{\theta}{4}+\frac{1}{8} > 0 \ \displaystyle \frac{\theta}{4} > 0 \ \displaystyle \frac{\eta}{4} > 0 \ \displaystyle \frac{\eta}{4}+\frac{3}{8} > 0 \ \displaystyle \frac{1}{2}(1-\theta-\eta) > 0 \ \displaystyle \frac{\theta}{4}+\frac{1}{8} + \frac{\theta}{4} + \frac{\eta}{4} + \frac{\eta}{4}+\frac{3}{8} + \frac{1}{2}(1-\theta-\eta) = 1 \ 0 < \theta < 1 \ 0 < \eta < 1 \ \end{array}

\begin{array}{l} \theta + \eta < 1 \ 0 < \theta < 1 \ 0 < \eta < 1 \ \end{array}

P (\theta | y, \eta) = \left( \left( \frac{\theta}{4}+\frac{1}{8} \right)^{14}, \ \frac{\theta}{4}, \ \frac{\eta}{4}, \ \frac{\eta}{4}+\frac{3}{8}, \ \left( \frac{1}{2} (1-\theta-\eta) \right)^5 \right), \quad \theta \in (0, 1-\eta) \ P (\eta | y, \theta) = \left( \left( \frac{\theta}{4}+\frac{1}{8} \right)^{14}, \ \frac{\theta}{4}, \ \frac{\eta}{4}, \ \frac{\eta}{4}+\frac{3}{8}, \ \left( \frac{1}{2} (1-\theta-\eta) \right)^5 \right), \quad \eta \in (0, 1-\theta)