Rejection sampling
拒绝采样法是生成一个随机变量样本的方法。
如果随机变量 $X$ 取值于 $[a, b]$，且密度函数 $p(x)$ 有上界 $M$
- 生成 $U_1, U_2 \sim U[0, 1]$
- 令 $X = a + (b-a) \cdot U_1$
- 如果 $U_2 < \displaystyle\frac{p(X)}{M}$，则输出 $X$；否则，重新开始执行 step 1
上述算法可以生成随机变量 $X$ 的样本。
或者对于 $X$ 选取试投密度 $g$，并找到一个尽量小的常数 $c$ 使得 $\displaystyle\frac{p(x)}{g(x)} \leq c$
- 生成随机样本 $X \sim g, \; U \sim U[0, 1]$
- 如果 $U < \displaystyle \frac{p(X)}{cg(X)}$，则输出 $X$；否则重新开始执行 step1
$$ \begin{aligned} \mathbb{P}(X \text { accepted }) &=\mathbb{P}\left(U \leq \frac{f(X)}{c g(X)}\right) \\ &=\int \mathbb{P}\left(U \leq \frac{f(x)}{c g(x)} \mid X=x\right) g(x) \mathrm{d} x \\ &=\int \frac{f(x)}{c g(x)} g(x) \mathrm{d} x \\ &=\frac{1}{c} \end{aligned} $$
接受的概率是 $1/c$ ，因此 $c$ 越小，算法的效率越高。
算法正确性的证明如下： $$ \begin{aligned} \mathbb{P}(X \leq t \mid X \text { accepted }) &=\frac{\mathbb{P}(X \leq t, X \text { accepted })}{\mathbb{P}(X \text { accepted })} \\ &=\frac{\mathbb{P}(X \leq t, X \text { accepted })}{1 / c} \\ &=c \mathbb{E}_{g} \mathbb{E}\left[\mathbf{1}\{x \leq t\} \mathbf{1}\left\{U \leq \frac{f(x)}{c g(x)}\right\} \mid X=x\right] \\ &=c \mathbb{E}_{g}\left[\mathbf{1}\{X \leq t\} \mathbb{E}\left[\mathbf{1}\left\{U \leq \frac{f(x)}{c g(x)}\right\} \mid X=x\right]\right] \\ &=c \mathbb{E}_{g}\left[\mathbf{1}\{X \leq t\} \frac{f(X)}{c g(X)}\right] \\ &=\int_{-\infty}^{\infty} \mathbf{1}\{x \leq t\} \frac{f(x)}{g(x)} g(x) d x \\ &=\int_{-\infty}^{t} f(x) d x \\ &=F(t) \end{aligned} $$
示例：生成密度函数为 $p(x)=\displaystyle\frac{1}{0.000336} x(1-x)^{3}, \quad 0.8<x<1$ 的随机数
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28import numpy as np import matplotlib.pyplot as plt from random import random p = lambda x: 1 / 0.000336 * x * (1 - x) ** 3 g = lambda x: 50 - 50 * x def LT(a=0, b=1): # 密度函数以 [a, b] 为底，单调下降的三角形 standard = min(random(), random()) return a + (b - a) * standard def rand_example(): c = 2 while True: x = LT(0.8, 1) y = random() if y < p(x) / (c * g(x)): return x r = [rand_example() for _ in range(10000)] plt.hist(r, bins=30, density=True) x = np.arange(0.8, 1, 0.01) plt.plot(x, p(x)) # 画出真实的概率密度曲线 plt.xticks([0.8, 0.9, 1]) plt.show()
Empirical Supremum Rejection Sampling
拒绝采样法的效率严重依赖于常数 $c$，当 $c = \displaystyle\sup_x \displaystyle \frac{p(x)}{g(x)}$ 时效率最高，一种改进的拒绝采样方法允许我们同时估计最优的常数 $c$ 和生成所需的随机样本，也就是说，我们没有必要事先对 $c$ 进行估计。
- 初始化一个对常数 $c$ 的估计 $\hat{c} > 1$
- 生成随机样本 $X \sim g, \; U \sim U[0, 1]$
- 如果 $U < \displaystyle \frac{p(X)}{cg(X)}$，则输出 $X$
- 更新 $\hat{c} = \max \left\{\hat{c}, \displaystyle\frac{p(X)}{g(X)}\right\}$
详细的分析可见：
- https://bookdown.org/rdpeng/advstatcomp/rejection-sampling.html#empirical-supremum-rejection-sampling
- https://academic.oup.com/biomet/article-abstract/89/4/745/242234
示例
1 2 3 4 5 6 7 8 9 10 11c = 1.1 def rand_example(): # Empirical Supremum Rejection Sampling global c while True: x = LT(0.8, 1) y = random() if y < p(x) / (c * g(x)): c = max(c, p(x) / g(x)) return x c = max(c, p(x) / g(x))
为计算 $I = \displaystyle\int_a^b h(x) \mathrm{d} x$，取 $U \sim U[a, b]$，则 $$ \mathrm{E}[h(U)] =\int_{a}^{b} h(u) \frac{1}{b-a} d u=\frac{I}{b-a} \;\Longrightarrow\; I =(b-a) \cdot \mathrm{E}[h(U)] $$ 于是可用 $\hat{I} = \displaystyle\frac{b - a}{N} \sum_{i=1}^N h(U_i)$ 来近似 $I$ .
只要 $h(U)$ 的期望存在且有限，根据大数定律，就有 $$ \hat{I} \to _{a.s.} I \quad(N \to \infty) $$ 根据中心极限定理 $$ \left . \left(\sum_{i=1}^N h(U_i) - N \cdot \mathrm{E}[h(U)]\right) \middle/ \sqrt{N \cdot\text{ Var}[h(U)]} \right. \;\to_d\, N(0, 1) \quad (N \to \infty) $$ 即有 $$ \text{error} = \hat{I} - I \sim N\left(0, \frac{(b-a)^2}{N} \text{ Var}[h(U)]\right) $$ 这时 $\text{Var}(\hat{I}) = \displaystyle\frac{(b-a)^2}{N} \text{Var}[h(U)]$ .
Importance Sampling
考虑以 $X \sim g(x)$ 为密度函数进行非均匀抽样，则 $$ \mathrm{E}_g \left[\frac{h(X)}{g(X)}\right] = \int \frac{h(x)}{g(x)} g(x) \mathrm{d} x = I $$ 当 $X$ 服从均匀分布时即为均匀抽样。
可用 $\hat{I} = \displaystyle\frac{1}{N}\sum_{i=1}^N \frac{h(X_i)}{g(X_i)}$来近似 $I$ ， 此时 $\text{Var}(\hat{I}) = \displaystyle\frac{1}{N} \text{Var}\left[\frac{h(X)}{g(X)}\right]$
因为 $$ \begin{aligned} \operatorname{Var}\left[\frac{h( {X})}{g( {X})}\right] &=\mathrm{E}\left(\frac{h^{2}( {X})}{g^{2}( {X})}\right)-\left(\mathrm{E}\frac{h( {X})}{g( {X})}\right)=\mathrm{E}\left(\frac{h^{2}( {X})}{g^{2}( {X})}\right)-\left(\int_{C} h( {x}) d {x}\right)^{2} \\ & \geq\left(\mathrm{E}\frac{|h( {X})|}{g( {X})}\right)^{2}-I^{2} \quad \text { (Jensen's inequality) } \\ &=\left(\int_{C}|h( {x})| d {x}\right)^{2}-I^{2} \end{aligned} $$ 当 $g(x)$ 与 $h(x)$ 形状接近时，$\hat{I}$ 的方差可以得到减小。
Markov Chain Monte Carlo
MCMC 是一种对高维随机向量抽样的方法，此方法模拟一个 Markov 链，使 Markov 链的平稳分布为目标分布，由此产生大量的近似服从目标分布的样本，但样本之间不是相互独立的。
给定一个分布 $\pi$，只要找到一个 Markov 链，使得 $\pi$ 是这个 Markov 链的平稳分布，就可以通过模拟 Markov 链的轨道来生成分布的随机样本。
MCMC 方法的核心就在于如何找到一个状态转移过程使得平稳分布是我们要模拟的分布。
一个充分条件是细致平稳条件 (detailed balance equation) $$ \pi(i) P(i, j)=\pi(j) P(j, i) $$
Metropolis-Hastings
找到一个恰当的转移概率可能是很难的，但是，从任意一个转移概率 $q(y \mid x)$ 出发，我们可以适当构造一个新的转移概率，来得到目标分布 $\pi(x)$ 的样本。Metropolis-Hastings 采样法是 MCMC 的一种实现。
首先，如果 $q(y \mid x)$ 不满足细致平稳条件，即 $$ \pi(x) q (y \mid x) \neq \pi(y) q(x \mid y) $$ 我们可以设计一个新的转移概率 $\alpha(y \mid x)$ 使得 $$ \pi(x) \cdot q (y \mid x) \alpha (y \mid x) \neq \pi(y) \cdot q(x \mid y) \alpha(x \mid y) $$ 取 $$ \alpha(y \mid x)=\min \left\{\frac{\pi(y) q(x \mid y)}{\pi(x) q(y \mid x)}, 1\right\} $$ 可使两边相等。实际上这时候的转移概率是 $$ K(y \mid x)=\alpha(y \mid x) q(y \mid x)+\mathbf{1}\{y=x\}\left[1-\int \alpha(s \mid x) q(s \mid x) \mathrm{d} s\right] $$ MH 算法可以描述为
- 给定初始状态 $X$
- 生成样本 $Y \sim q(y \mid X), \; U \sim U[0, 1]$
- 如果 $U \leq \alpha (Y\mid X)$ 则转移到状态 $X$，否则下一个状态还是 $X$
- 继续 step2，状态序列即为所求随机样本
特殊地，如果初始转移概率取为 $q(y \mid x) = q(x \mid y)$，则 $\alpha(y \mid x)=\min \left\{\displaystyle\frac{\pi(y) }{\pi(x) }, 1\right\}$；如果取 $q(x \mid y) = q(x)$，则 $\alpha(y \mid x)=\min \left\{\displaystyle\frac{\pi(y) q(x)}{\pi(x) q(y)}, 1\right\}$
示例：对一个 0-9 的离散分布进行 MH 采样，$q(x\mid y)=q(y\mid x)$ 为离散均匀分布
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22import random import numpy as np import pandas as pd import statsmodels.tsa.api as smt pi = np.arange(1, 11) pi = pi / np.sum(pi) x = 0 samples = [] for _ in range(10000): y = random.randint(0, 9) if random.random() < pi[y] / pi[x]: x = y samples.append(x) pd.Series(samples).value_counts(ascending=True, normalize=True).plot.bar(rot=0) # 观察样本间的序列相关性 smt.graphics.plot_acf(samples, lags=20) smt.graphics.plot_pacf(samples, lags=20)示例：用 MH 方法模拟 Beta(2, 4) 大于 0.8 的部分，$q(x\mid y) = q(y \mid x) \sim U[0.8, 1]$
1 2 3 4 5 6 7samples = [] a = random.random() * 0.2 + 0.8 for _ in range(10000): b = random.random() * 0.2 + 0.8 if random.random() < b * (1 - b) ** 3 / (a * (1 - a) ** 3): a = b samples.append(a)
Random Walk Metropolis-Hastings
令 $$ q(y \mid x) = x + \epsilon $$ 且 $\epsilon \sim N(0, \delta^2 I)$，这时候 $q(y \mid x) = q(x \mid y)$ .
示例：用 MH 方法模拟 $N\left(\begin{bmatrix}1 \\ 2\end{bmatrix}, \begin{bmatrix}1 & 0.5\\ 0.5 & 1\end{bmatrix}\right)$，$q(y \mid x) \sim x + N(0, 0.3I)$
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27import random import numpy as np import matplotlib.pyplot as plt from scipy import stats dist = stats.multivariate_normal( mean=np.array([1, 2]), cov=np.array([[1, 0.5], [0.5, 1]]) ) q = stats.multivariate_normal( mean=np.array([0, 0]), cov=np.array([[0.3, 0], [0, 0.3]]) ) x = np.array([0, 0]) samples = [] for _ in range(10000): y = x + q.rvs() if random.random() < dist.pdf(y) / dist.pdf(x): x = y samples.append(x) samples = np.array(samples) # samples = dist.rvs(10000) plt.scatter(x=samples[:, 0], y=samples[:, 1], marker='.')
Simulated Annealing
https://bookdown.org/rdpeng/advstatcomp/simulated-annealing.html