蒙特卡罗与MCMC

简述

蒙特卡罗方法是一种用多次取样来逐渐逼近样本真实分布的方法。

另外蒙特卡罗并不是人的名字，而是位于摩纳哥的一个赌场的名字。

简单的例子

最经典的例子应该是几何概型中的算$\pi$的数值。也就是高中课本上的撒豆子问题。

通过多次随机取样，取样点会随机分布在一个正方形的区域里，然后在该区域中画圆，就会是样本分成两部分。一部分在圆内，另一部分在圆外。然后通过算圆内的样本数占总样本数的比例可以近似算的圆的面积，最后通过圆面积公式能算得$\pi$的取值。

实操

$$\frac{1}{4}S =\frac{1}{4}\pi R^{2} = \frac{n}{N} \\ n:圆内区域的样本数 \\ N:样本总数$$

import random

sampleNum = 10000
counter = 0
plt.figure(figsize=(8, 8))

x1 = list()
y1 = list()
x2 = list()
y2 = list()
for i in range(sampleNum):
    x = random.random()
    y = random.random()
    if x**2 + y**2 < 1:
        counter += 1
        x1.append(x)
        y1.append(y)
    else:
        x2.append(x)
        y2.append(y)
   
print(4*counter/sampleNum)
plt.scatter(x2, y2, color='b')
plt.scatter(x1, y1, color='r')

plt.show()

随着$sampleNum$的增大，最终结果应该会越来越接近3.14。

蒙特卡罗求积分

对于一个函数$f(x)$,我们想在求这个函数在$x\in [a, b]$上的积分。我们会发现它的积分可能会比较难显式地求出来。这个时候我们可以通过蒙特卡罗的方法将他取样出来。

比如说我们现在打算求$f(x)$的积分，姑且假定在区间$x \in [a, b]$上，$f(x) \ge 0$。此时，我们可以设定一个上界$U \ge \mathrm{sup}_{x \in [a, b]} f(x)$。 然后随机从均匀分布$X \sim \mathrm{Unif}(a, b)$和$Y \sim \mathrm{Unif}(0, U)$上采样$\{x_{i}, x_{j}\}$。当样本数量足够大的时候，

$$\int_{a}^{b}f(x)dx = \frac{\sum_{(x_{i}, y_{i})}\delta_{(y_{i} \le f(x_{i}))}}{N_{tol}}$$

上面的过程中，我们需要对两个变量进行采样，那么是否可以只对一个变量进行采样的方法呢？

回想一下黎曼积分的思想，我们可以将积分看作在$x$轴上进行均匀地切割，然后把对应的函数值加在一起。基于类似的思想，我们可以只对$X\sim \mathrm{Unif}(a, b)$进行采样，采样结果记作$\{x_{i}\}_{N}$。 此时，我们可以

$$\int_{a}^{b} f(x) dx = \sum_{i}^{N}f(x_{i}) \frac{b-a}{N},$$

其中$\frac{b-a}{N}$可以看作均匀切割后的每一段的长度，其实可以对应积分中的$dx$。 通过这种方法，我们不再需要对$Y$进行采样。

蒙特卡罗求期望

对于一个函数$f(x)$,我们想在求这个函数在$x$上的期望$\mathbb{E}_{x\sim\pi(x)}[f(x)]$，其中$\pi(x)$是随机变量$X$服从的分布的概率密度函数或者分布律。 根据期望的定义我们有

$$\mathbb{E}_{x\sim \pi(x)}[f(x)] = \int_{x}\pi(x)f(x) dx.$$

假设我们对$X$进行采样，样本记作$\{x_{i}\}_{N}$。 根据大数定律，则有

$$\mathbb{E}_{x\sim \pi(x)}[f(x)] = \frac{1}{N}\sum_{i}^{N}f(x_{i}).$$

举例来讲，如果现在$X \sim \mathrm{Unif}(a, b)$， 那么

$$\mathbb{E}_{x\sim \pi(x)}[f(x)] = \int_{x}\frac{1}{b-a}f(x) dx = \frac{1}{N}\sum_{i}^{N}f(x_{i}).$$

重要性采样

有的时候直接从目标分布（注意，这个分布是已知的）采样非常困难，代价高昂，或者一定样本的到的精度不高时。 我们可以引入一个容易采样的建议分布（Proposal Distribution），并赋予每个样本相应的权重，来无偏地估计目标分布的期望值。

假设我们选择的易于采样的建议分布的pdf为$q(x)$，

$$\begin{aligned} \mathbb{E}_{x\sim \pi(x)}[f(x)] &= \int_{x}\pi(x)f(x) dx\\ &= \int_{x}q(x)\frac{\pi(x)}{q(x)}f(x) dx \\ &= \mathrm{E}_{x \sim q(x)}\left[\frac{\pi(x)}{q(x)}f(x)\right]. \end{aligned}$$

如果我们把$\frac{\pi(x)}{q(x)}$定义为$w(x)$，那么

$$\mathrm{E}_{x \sim q(x)}\left[\frac{\pi(x)}{q(x)}f(x)\right] = \mathrm{E}_{x \sim q(x)}\left[w(x)f(x)\right]$$

可以看作我们在建议分布上采样的时候，目标函数$f(x)$的值要乘上一个权重$w(x)$。 对于不同的$x$，权重不同，说明对应的$f(x)$的重要性不同，这个便是该小节标题中“重要性”的来源。

对$X$进行采样之后，再次根据大数定律，

$$\mathrm{E}_{x \sim q(x)}\left[w(x)f(x)\right] = \frac{1}{N}\sum_{i}^{N}w(x_{i})f(x_{i}).$$

与上面的直观理解相似，“重要性”体现在，对于不同$x_{i}$，权重$w(x_{i})$不同，$f(x_{i})$的”贡献“自然不同。 对上面的式子中，我们需要注意的是，目标分布$\pi(x)$，我们选取的建议分布$q(x)$和目标函数$f(x)$都是已知的。我们只是想通过采样之后得到$f(x)$的期望。

拒绝采样

重要性采样是对目标函数$f(x)$添加权重，但我们的得到的样本$\{x_{i}\}_{N}$是服从建议分布的，而不是服从目标分布。 那我们是否有办法真正的得到从目标分布中采样的来的样本呢？ 答案是肯定的。

因为目标分布的直接采样困难， 在拒绝采样的过程中，我们还是要先选取一个易于采样的建议分布$q(x)$。 与重要性采样不同的是，我们还要选取一个常数$M$，使得

$$M\cdot q(x) \ge \pi(x).$$

自然，我们接下来应该做的是直接从$q(x)$中采样。 在重要性采样中，我们对于没一个采到的样本都保留下来，进行利用。 而在拒绝采样中，对于采到的样本$x_{i}$，我们以$\alpha_{\mathrm{accept}}$的概率进行保留。自然以$\alpha_{\mathrm{reject}} = 1 - \alpha_{\mathrm{accept}}$的概率丢弃。

对于$x_{i}$接受概率的计算：

$$\alpha_{\mathrm{accept}} = \frac{\pi(x_{i})}{M\cdot q(x_{i})}.$$

为了实现概率接受和概率丢弃，我们可以随机生成一个服从$\mathrm{Unif}(0, 1)$的随机数$u$，如果$u$落在$[0, \alpha_{\mathrm{accept}}]$，$x_{i}$就会被接受，反之，如果$u \in (\alpha_{\mathrm{accept}}, 1]$，$x_{i}$便会被丢弃。

拒绝采样的直觉理解是很简单的，

图中，$\tilde{p}(z)$便是目标分布（即，本文中的$\pi(x)$）。$k$和$q(z)$为选取的$M$和$q(x)$。 对与样本$z_{0}$而言，$z_0$在$\tilde{p}(z)$中被采集到的可能性比上在$kq(z)$中被采集到的可能性应该为两个pdf在$z_{0}$处的高度比。

MCMC

拒绝采样适用于低维简单分布，而 MCMC 专为高维复杂分布设计，通过构建马尔可夫链使其平稳分布收敛于目标分布。 MCMC采样很神奇的点在于，每次采样的时候不是完全独立，而是上一个采样的样本会对下一采样的样本有一定约束。

马尔可夫链

对于一组随机变量$\boldsymbol{X} = \{X_{t}\}_{t \ge 0}$，其中$t$代表时刻，如果有

$$P(X_{t} \mid \{X_{t}\}_{0 \le t \le t-1}) = P(X_{t} \mid X_{t-1}),$$

那么，称$\boldsymbol{X} = \{X_{t}\}_{t \ge 0}$为马尔可夫链。 进一步，$P(X_{t}\mid X_{t-1})$被称为$X_{t-1}$到$X_{t}$的转移概率分布。

例如，我们在这里假设$\boldsymbol{X}_{\le t} = \{X_{t}\}_{\le t} \in [1, 2, 3]^{(t+1)}$， 我们定义

$$\begin{aligned} \boldsymbol{P} = \left[ \begin{array}{ccc} p_{11} & p_{12} & p_{13} \\ p_{21} & p_{22} & p_{23} \\ p_{31} & p_{32} & p_{33} \end{array} \right], \end{aligned}$$

其中，$P_{ji} = P(X_{t}=j \mid X_{t-1}=i)$。 我们称上述的矩阵为转移概率矩阵。

我们进一步定义

$$\pi(t) = [\pi_{1}(t), \pi_{2}(t), \pi_{3}(t)]^{\top},$$

为$t$时刻的状态分布， 其中$\pi_{i}(t) = P(X_{t}=i)$。 $\pi(t) = [\pi_{1}(t), \pi_{2}(t), \pi_{3}(t)]$即为$t$时刻目标分布的分布律（这里考虑的是离散的情况）。 同时，$\pi_{t}$还可以写作

$$\pi(t)^{\top} = \pi(0)^{\top}P^{t}.$$

我们都知道，对于一个转移过程，若

$$\pi(t)^{\top} = \pi(t)^{\top}P，$$

则称$\pi(t)$是马尔可夫链的平稳分布。

一个马尔可夫链拥有唯一平稳分布通常需要具备一下三点性质：

• 不可约性（Irreducible）：系统中的任何一个状态都可以通过有限步转移到达其他所有状态。
• 非周期性（Aperiodic）：系统从某状态出发，返回该状态的可能步数不具有公约数，避免了系统出现死循环振荡。
• 正常返性（Positive Recurrent）：系统平均返回到某一状态的时间是有限的。

最后，我们定义一种特殊的马尔可夫链。

对于任意时刻$t$，如果一个马尔可夫链满足细致平衡方程（detailed balance equation）:

$$\forall i \forall j, \qquad p_{ji}\pi_{i} = p_{ij}\pi_{j},$$

则称该马尔可夫链为可逆马尔可夫链。

定理 1
满足细致平衡方程的状态分布$\pi$是该马尔可夫链的平稳分布。

MCMC采样

如定理1所言，我们只要能找到一个满足细致平衡方程的$\boldsymbol{P}$，那么，马尔可夫转移每次选择的状态（在本文中为1 or 2 or 3）可以作为采的样本，并且该样本符合目标分布$\pi$。

接下来的问题变成如何去寻找满足条件的$\boldsymbol{P}$。

我们可以先任取一个转移矩阵$\boldsymbol{Q}$。此时我们一般有

$$\pi_{i}q_{ji} \ne \pi_{j}q_{ij}.$$

但我们可以在不等式两边进行适当缩放，然后不等式变为等式

$$\pi_{i}q_{ji}\alpha_{ji} = \pi_{j}q_{ij}\alpha_{ij},$$

其中$\alpha_{ji}$和$\alpha_{ij}$为缩放系数。

进一步，要是等式恒成立，我们可以简单地令

$$\alpha_{ji} = \pi_{j}q_{ij}, \qquad \alpha_{ij} = \pi_{i}q_{ji}.$$

此时，等式两边都为$\pi_{j}\cdot q_{ij}\cdot \pi_{i}\cdot q_{ji}$。

最后，如果我们把 $\alpha_{i}q_{ji}$看作$\tilde{q}_{ji}$，$\alpha_{j}q_{ij}$看作$\tilde{q}_{ij}$，便得到了细致平衡方程

$$\pi_{i}\tilde{q}_{ji} = \pi_{j}\tilde{q}_{ij}.$$

注意这里的平稳分布是$\pi_{i}$和$\pi_{j}$，但是$\tilde{q}_{ji}$和$\tilde{q}_{ij}$并未进行归一化，所以应该不能是一个合理的转移矩阵。

那么，我们是否还可以用$\widetilde{\boldsymbol{Q}}$来进行采样呢？ 答案是可以的。 实际上，我们可以这么定义$p_{ji}$和$\tilde{q}_{ji}$:

$$\begin{aligned} p_{ji} = \left\{ \begin{array}{ll} \tilde{q}_{ji}, & i \ne j, \\ \tilde{q}_{ji} + (1 - \sum_{j \in \mathrm{supp}(\boldsymbol{X})}\tilde{q}_{ji}), & j = i. \end{array} \right. \end{aligned}$$

等价的

$$p_{ji} = \tilde{q}_{ji} + (1 - \sum_{j \in \mathrm{supp}(\boldsymbol{X})}\tilde{q}_{ji})\delta_{i=j}.$$

直觉上，我们可以认为当两次采样的状态（值）不一样的时候，直接用$\tilde{q}_{ji}$代替$p_{ji}$，而当$i=j$的时候，我们在用$\tilde{q}_{ii}$代替$p_{ii}$的基础上，还加入了让$\sum_{j \in \mathrm{supp}(\boldsymbol{X})}\tilde{q}_{ji}$与$1$的差值。 这样子我们便得到了一个归一化的$p_{ji}$了。

至于这样的$p_{ji}$是否满足细致平稳呢？ 我们发现，对于任意$\pi_{i}$，$\pi_{j}$，$f_{i, j}$，和$f_{j,i}$都有

$$f_{i, j}\delta_{i=j}\pi_{i} = f_{j, i}\delta_{i=j}\pi_{j}.$$

具体来说，当$i \ne j$时，上式必然成立。当$i = j$的时候，$f_{i, i}\pi_{i} = f_{j, j}\pi_{j}$也是必然成立的。 所以，当$\widetilde{\boldsymbol{Q}}$满足细致平衡条件的时候，多出来的$i=j$的项并不影响$\boldsymbol{P}$也满足细致平衡条件。

一言蔽之就是说，通过调整转移矩阵的对角线元素可以完成归一化，并且不影响细致平稳条件的成立。

Metropolis算法

接下来写一下到目前为止的算法。 我们的目标实在满足细致平衡条件的$\widetilde{\boldsymbol{Q}}$上采样，如果$\widetilde{\boldsymbol{Q}}$比较难以采样，我们便在容易采样的$\boldsymbol{Q}$上采样。 注意，$\widetilde{q}_{ji}$和${q}_{ji}$差一个$\alpha_{ji}$的倍数。 当我们通过${q}_{ji}$从状态$i$选择状态$j$的时候，如果$i \ne j$，我们会发现我们实际上需要的选择的概率为$\tilde{q}_{ji}$。 与拒绝采样相似，$\alpha_{ji}$可以被看作是一个接受概率。 换言之，当我们处于$i$状态，而下一个被选择的状态是与$i$不同的$j$的时候，我们以概率$\alpha_{ji}$接受这个选择，以$1-\alpha_{ji}$的概率拒绝这个选择。 当我们拒绝选择状态$j$的时候，我们会留在状态$i$，即，有采样了一次$i$。

需要注意的是，我们由状态$i$到状态$i$时的接受率并不是$\alpha_{ii}$而是$1-\alpha_{ji}$。 这就是之前所讲的，通过调整转移矩阵的对角线元素可以完成归一化，实现最终以$p_{ji}$进行转移。

最后，Metropolis算法总结如下

Metropolis采样

初始状态的$x_{0}$, t时刻状态为$x_{t}$。
通过如下流程采样样本$x_{t+1}$:

1. 通过$q_{t+1, t}$采样一个样本$x^{*}_{t+1}$。
2. 从$\mathrm{Unif}(0, 1)$中，生成一个随机数$u$。
3. 如果$u \le \alpha_{t+1, t}$，$x_{t+1} \gets \alpha_{t+1, t} \cdot x^{*}_{t+1}$；否则，$x_{t+1} \gets x_{t}$。

这个算法与拒绝采样相似，但不完全一直。 拒绝采样中如果样本被拒绝那就是拒绝了，但Metropolis算法的话，会把原地踏步把目前的状态计入时刻$t$的样本。

Metropolis-Hastings 采样

通过$\alpha_{ji}$来接受拒绝下一个状态的选择其实会面临一个问题。 如果$\alpha_{ji}$太小，那采样会一直原地踏步，无法对其他状态进行探索。

为了解决这个问题，我们先回到细致平稳方程

$$\pi_{i}q_{ji}\alpha_{ji} = \pi_{j}q_{ij}\alpha_{ij}.$$

正如我前面所说，$\alpha_{ji}$和$\alpha_{ij}$是一个缩放，那么我们能否在这个式子的基础上，左右两边再同乘上一个常数$c_{i}$，让$c_{i}\cdot \alpha_{ji}$来作为接受率呢？ 显然是可以的，但是接受率一定要是一个介于0到1之间的数，我们再乘上$c_{i}$之后，$c_{i}\cdot \alpha_{ji}$和$c_{i}\cdot \alpha_{ij}$都不能超多一。 对于$c_{i}\cdot \alpha_{ji}$，我们何时可以将它缩放到最大呢？ 当$c_{i}\cdot \alpha_{ji} = 1$或者$c_{i}\cdot \alpha_{ij}=1$的时候。

令$\mathcal{A}_{ji}:=c_{i}\cdot \alpha_{ji}$，则有

$$\begin{aligned} \mathcal{A}_{ji} = \left\{ \begin{array}{ll} 1, & \alpha_{ji} \ge \alpha_{ij},\\ \frac{\pi_{j}q_{ij}}{\pi_{i}q_{ji}}, & \alpha_{ji} < \alpha_{ij}. \end{array} \right. \end{aligned}$$

最后，我们只需要让新的接受率

$$\mathcal{A}_{ji} = \mathrm{min}(1, \frac{\pi_{j}q_{ij}}{\pi_{i}q_{ji}})$$

就好。

如果，$q_{ij} = q_{ji}$， 则

$$\mathcal{A}_{ji} = \mathrm{min}(1, \frac{\pi_{j}}{\pi_{i}})$$

Metropolis-Hastings 采样

初始状态的$x_{0}$, t时刻状态为$x_{t}$。
通过如下流程采样样本$x_{t+1}$:

1. 通过$q_{t+1, t}$采样一个样本$x^{*}_{t+1}$。
2. 从$\mathrm{Unif}(0, 1)$中，生成一个随机数$u$。
3. 如果$u \le \mathcal{A}_{t+1, t}=\mathrm{min}(1, \frac{\pi_{t+1}q_{t, t+1}}{\pi_{t}q_{t+1, t}})$，$x_{t+1} \gets \alpha_{t+1, t} \cdot x^{*}_{t+1}$；否则，$x_{t+1} \gets x_{t}$。

一段M-H采样高斯分布的python 代码

import random
import math
from scipy.stats import norm
import matplotlib.pyplot as plt

def norm_dist_prob(theta):
    y = norm.pdf(theta, loc=10, scale=5)
    return y

T = 5000
pi = [0 for i in range(T)]
sigma = 1
t = 0
while t < T-1:
    t = t + 1
    pi_star = norm.rvs(loc=pi[t - 1], scale=sigma, size=1, random_state=None)   #状态转移进行随机抽样
    alpha = min(1, (norm_dist_prob(pi_star[0]) / norm_dist_prob(pi[t - 1])))   #alpha值

    u = random.uniform(0, 1)
    if u < alpha:
        pi[t] = pi_star[0]
    else:
        pi[t] = pi[t - 1]


plt.scatter(pi, norm.pdf(pi, loc=10, scale=5),label='Target Distribution', c= 'red')
num_bins = 50
plt.hist(pi, num_bins, density=1, facecolor='green', alpha=0.7,label='Samples Distribution')
plt.legend()
plt.show()

吉布斯采样（Gibbs Sampling）

接下来，我们看看更现实的应用场景，高维数据的情况。

假设$\boldsymbol{X}_{t}= \{X_{t, 1}, X_{t, 2}, \cdots, X_{t, d}\}$是一个$d$维的随机变量。 Gibbs采样每个时刻$t$，只对$d$维中的一维进行采样。 比如说我们接下来要对第$h$维采样，我们将时刻$t+1$的样本记作$\boldsymbol{X}_{t+1} = \boldsymbol{X}_{t, [X_{t, h}=X_{t+1, h}]}$ ，我们定义转移概率（这里的多维变量，我们用最稳妥的写法）为

$$q(\boldsymbol{X}_{t+1} \mid \boldsymbol{X}_{t}) = P(X_{t+1, h} \mid \boldsymbol{X}_{t, -h})= \frac{P(\boldsymbol{X}_{t, [X_{t, h}=X_{t+1, h}]})}{P(\boldsymbol{X}_{t, -h})},$$

其中$\boldsymbol{X}_{t, -h} = \boldsymbol{X}_{t} \setminus \{X_{t, h}\}$。 这个式子意味着，我们每次固定第$h$维以外的维度的取值，然后只更新第$h$维的取值。 然后，我们定义接受率

$$\begin{aligned} \mathcal{A}_{t+1, t} &:= \mathrm{min}\left(1, \frac{P(\boldsymbol{X}_{t+1})q(\boldsymbol{X}_{t}\mid \boldsymbol{X}_{t+1})}{P(\boldsymbol{X}_{t})q(\boldsymbol{X}_{t+1}\mid \boldsymbol{X}_{t})}\right) \\ &= \mathrm{min}\left(1, \frac{P(\boldsymbol{X}_{t, [X_{t, h}=X_{t+1, h}]})P(X_{t, h} \mid \boldsymbol{X}_{t, -h})}{P(\boldsymbol{X}_{t})P(X_{t, h} \mid \boldsymbol{X}_{t, -h})}\right) \\ &= \mathrm{min}\left(1, \frac{P(\boldsymbol{X}_{t, [X_{t, h}=X_{t+1, h}]})P(X_{t, h} \mid \boldsymbol{X}_{t, -h})P(\boldsymbol{X}_{t, -h})}{P(\boldsymbol{X}_{t})P(\boldsymbol{X}_{t, [X_{t, h}=X_{t+1, h}]})}\right) \\ &= \mathrm{min}\left(1, \frac{P(\boldsymbol{X}_{t, [X_{t, h}=X_{t+1, h}]})P(X_{t})}{P(\boldsymbol{X}_{t})P(\boldsymbol{X}_{t, [X_{t, h}=X_{t+1, h}]})}\right) \\ &= \mathrm{min}(1, 1) \\ &= 1. \end{aligned}$$

所以，对于每一个时刻$t$所作的选择，Gibbs采样都不做拒绝地将其作为样本保留下来。

Gibbs采样

初始状态为$\boldsymbol{x}_{0}$, t时刻状态的分布为$\boldsymbol{x}_{t}$。
通过如下流程采样样本$\boldsymbol{x}_{t+1}$:

1. 均匀的从$[d]$中选取一个维度$h$。
2. 计算$q(\boldsymbol{X}_{t+1} \mid \boldsymbol{X}_{t})= \frac{P(\boldsymbol{X}_{t, [X_{t, h}=X_{t, h}]})}{P(\boldsymbol{X}_{t, -h})}$
3. 通过$q(\boldsymbol{X}_{t+1} \mid \boldsymbol{X}_{t})$采样第$h$维的一个状态$x_{t+1, h}$。
4. $x_{t+1} \gets \boldsymbol x_{t, [x_{t, h}=x_{t+1, h}]}$。

最后一步，是将第$d$维的数据替换成新采样的状态。

import random
import math
from scipy.stats import beta
from scipy.stats import norm
import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D
from scipy.stats import multivariate_normal

samplesource = multivariate_normal(mean=[5,-1], cov=[[1,0.5],[0.5,2]])

def p_ygivenx(x, m1, m2, s1, s2):
    return (random.normalvariate(m2 + rho * s2 / s1 * (x - m1), math.sqrt(1 - rho ** 2) * s2))

def p_xgiveny(y, m1, m2, s1, s2):
    return (random.normalvariate(m1 + rho * s1 / s2 * (y - m2), math.sqrt(1 - rho ** 2) * s1))

N = 5000
K = 20
x_res = []
y_res = []
z_res = []
m1 = 5
m2 = -1
s1 = 1
s2 = 2

rho = 0.5
y = m2

for i in range(N):
    for j in range(K):
        x = p_xgiveny(y, m1, m2, s1, s2)   #y给定得到x的采样
        y = p_ygivenx(x, m1, m2, s1, s2)   #x给定得到y的采样
        z = samplesource.pdf([x,y])
        x_res.append(x)
        y_res.append(y)
        z_res.append(z)

num_bins = 50
plt.scatter(x_res, norm.pdf(x_res, loc=5, scale=1),label='Target Distribution x', c= 'green')
plt.scatter(y_res, norm.pdf(y_res, loc=-1, scale=2),label='Target Distribution y', c= 'orange')
plt.hist(x_res, num_bins, density=1, facecolor='Cyan', alpha=0.5,label='x')
plt.hist(y_res, num_bins, density=1, facecolor='magenta', alpha=0.5,label='y')
plt.title('Histogram')
plt.legend()
plt.show()

参考

1. https://spaces.ac.cn/archives/8084
2. https://zhuanlan.zhihu.com/p/253784711
3. Chatgpt