蒙特卡罗法抽样方法
1 前言
蒙特卡罗法有 直接抽样法、接收-拒绝抽样法、重要性抽样法。
2 直接抽样法
假如 $p(x)$ 是要采样的分布。若可以得到 $x$ 的概率密度函数 $p(x)$, 对 $p(x)$ 求积分得到累积分布函数 $F(x)$,显然 $x$ 和 $y$ 是一 一对应的。
先随机抽取 $y$,然后通过求 $F(x)$ 的反函数 $F^{-1}(y)$ 得到随机变量 $x$ 的值:
$\boldsymbol{x}=\boldsymbol{F}^{-1}(\boldsymbol{y})$
随机变量 $y $ 在区间 $[0,1]$ 上均匀分布 $\Rightarrow$ 利用 $[0,1]$ 区 间上均匀分布随机数产生器抽取。
同理可以抽样出 $\mathrm{N}$ 个样本点 $x^{(1)},x^{(2)}, \ldots, x^{(N)}$
该方法的缺点是现实中很难得到概率密度函数,因此累积分布函数是求不出来的,所以很难通过这种方法去抽样得到样本点。这种方法只能适用于一 些非常有限的分布。
3 接收-拒绝抽样法
其中,$p(x)$ 是概率密度函数,$q(x)$ 为建议分布函数,$c$ 为 $q(x)$ 前的参数。
由于概率密度累积和为 $1$ ,故建议分布 $q(x)$ 不可能处处大于 $p(x)$ ,因此在 $q(x)$ 前面乘上一个参数 $c$,使其位于 $p(x)$ 之上,使其满足 $\forall x_{i}, c q\left(x_{i}\right) \geq p\left(x_{i}\right)$ 如图所示。
引入建议分布 $q(x)$ 的原因是在现实中 $p(x)$ 往往比较复杂,我们无法直接从 $p(x)$ 中采样。而 $q(x)$ 是相对比较简单的分布我们可以直接从里面采样, 再看看抽样的样本点是否能够代表 $p(x)$ 。
具体步骤如下:
输入:抽样的目标概率分布的概率密度函数 $p(x)$;
输出:概率分布的随机样本 $x_{1} \cdot x_{2}, \ldots, x_{n} $;
参数:样本数 $n$;(1) 选择合适的概率密度函数作为建议分布 $q(x)$,使其对任 $x$ 满足 $c q(x) \geq p(x)$ ,其中 $c>0$;
(2) 在建议分布 $q(x)$ 随机抽样得到样本 $x^{*}$ ,再按照均匀分布在 $(0,1)$ 范围内抽样得到 $u$;
(3) 如果 $u \leq \frac{p\left(x^{*}\right)}{c q\left(x^{*}\right)}$ ,则将 $x^{*}$ 作为抽样结果;否则,回到步骤 (2);
(4) 直至得到 $\mathrm{n}$ 个随机样本,结束。
解释一下这里 $\mathrm{u}$ 的用处。当 $p(x)$ 和 $c q(x)$ 接近时,也就是图中拒绝域比较窄时, $\frac{p\left(x^{*}\right)}{c q\left(x^{*}\right)}$ 比例接近 $1$,这时 $\mathrm{u}$ 在 $(0,1)$ 中取值满足小于 $\frac{p\left(x^{*}\right)}{c q\left(x^{*}\right)}$ 的概率比较大,接受率高,抽样效率比较高; 而当 $\frac{p\left(x^{*}\right)}{c q\left(x^{*}\right)}$ 比较小时,就会导致拒绝比例比较高,抽样效率较低。这样做能达到在 $p(x)$ 和 $c q(x)$ 接近的地方多抽样,在差异大的地方少抽样的效果。
接受-拒绝法的优点是容易实现,缺点是效率可能不高。如果 $p(x)$ 的涵盖体积占 $c q(x)$ 的涵盖体积的比例很低,就会导致拒绝的比例很高,抽样效率很低。注意,一般是在高维空间进行抽样,即使 $p(x)$ 与 $c q(x)$ 很接近,两者涵盖体积的差异也可能很大(与我们在三维空间 的直观不同)。
例子:
假设要抽样的分布 $p(x)$ 是:
$p(x)=\frac{1}{1.2113}\left(0.3 \exp \left(-(z-0.3)^{2}\right)+0.7 \exp \left(-(z-2)^{2} / 0.3\right)\right.$
在这个例子中,选择均值为 $1.4$,方差为 $1.2$ 的正态分布作为我们的建议分布,取常数 $C =2.5$。
代码片段:
from scipy.stats import norm
import seaborn
import numpy as np
seaborn.set()
import matplotlib.pyplot as plt
# 目标采样分布的概率密度函数
def p(x):
return (0.3 * np.exp(-(x - 0.3) ** 2) + 0.7 * np.exp(-(x - 2.) ** 2 / 0.3)) / 1.2113
# 建议分布G
norm_rv = norm(loc=1.4, scale=1.2)
# C值
C = 2.5
x = np.arange(-4., 6., 0.01)
plt.plot(x, p(x), color='r', lw=5, label='p(x)')
plt.plot(x, C*norm_rv.pdf(x), color='b', lw=5, label='C*g(x)')
plt.legend()
plt.show()
输出结果:
接下来具体使用接受-拒绝算法
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import uniform, norm
import seaborn
seaborn.set()
# 目标采样分布的概率密度函数
def p(x):
return (0.3 * np.exp(-(x - 0.3) ** 2) + 0.7 * np.exp(-(x - 2.) ** 2 / 0.3)) / 1.2113
# 建议分布
norm_rv = norm(loc=1.4, scale=1.2)
# C值
C = 2.5
uniform_rv = uniform(loc=0, scale=1)
sample = []
for i in range(100000):
Y = norm_rv.rvs(1)[0]
U = uniform_rv.rvs(1)[0]
if p(Y) >= U * C * norm_rv.pdf(Y):
sample.append(Y)
x = np.arange(-3., 5., 0.01)
plt.gca().axes.set_xlim(-3, 5)
plt.plot(x, p(x), color='r')
plt.hist(sample, color='b', bins=150, normed=True, edgecolor='k')
plt.show()
输出结果:
我们通过接受-拒绝采样方法,进行了10万次的数据采样,最后把采样结果绘制成。通过运行结果我们发现,柱状图的形状几乎和原始目标分布的概率密度函数曲线完全拟合。由此通过实践证明了这个采样的方法确实是有效的,他高度近似的拟合了原始样本的分布特性。
3 重要性抽样法
上面两种采样方法或多或少都有一些问题。
上面描述了两种从另一个分布获取指定分布的采样样本的算法,对于1. 在实际工作中,一般来说我们需要采样的分布都及其复杂, 太可能求解出它的反函数,但 $p(x)$ 的值也许还是可以计算的。对于 2 .找到一个合适的 $c q(x)$ 往往很困难,接受概率有可能会很低。 那我们回过头来看我们采样的目的: 其实是想求得 $E[f(x)]$, $x \sim p$ , 也就是
$E[f(x)]=\int_{x} f(x) p(x) d x$
如果符合 $\mathrm{p}(\mathrm{x})$ 分布的样本不太好生成,我们可以引入另一个分布 $q(x)$ ,可以很方便地生成样本。使得
$\begin{aligned}\int_{x} f(x) p(x) d x &=\int_{x} f(x) \frac{p(x)}{q(x)} q(x) d x \\&=\int_{x} g(x) q(x) d x\end{aligned}$
其中,$g(x)=f(x) \frac{p(x)}{q(x)}=f(x) w(x) $
我们将问题转化为了求 $g(x)$ 在 $q(x)$ 分布下的期望! ! !
我们称其中的 $w(x)=\frac{p(x)}{q(x)}$ 叫做Importance Weight.
重要性抽样法解决的问题:
首先当然是我们本来没办法直接从 $p$ 进行采样,这个是我们看到的,重要性采样将之转化为了从 $q$ 分布进行采样;同时重要性采样有时候还可以改进原来的采样,比如说:
可以看到如果我们直接从 $p$ 进行采样,而实际上这些样本对应的 $f(x)$ 都很小,采样数量有限的情况下很有可能都无法获得 $f(x)$ 值较大的样 本,这样评估出来的期望偏差会较大;
而如果我们找到一个 $q$ 分布,使得它能在 $f(x) * p(x)$ 较大的地方采集到样本,则能更好地逼近 $[E f(x)]$ ,因为有Importance Weight控制其比重,所以也不会导致结果出现过大偏差。
所以选择一个好的q也能帮助你sample出来的效率更高,要使得 $f(x) p(x)$ 较大的地方能被sample出来。
https://www.bilibili.com/read/cv4939412
https://blog.csdn.net/lukas_ten/article/details/114458591
因上求缘,果上努力~~~~ 作者:图神经网络,转载请注明原文链接:https://www.cnblogs.com/BlairGrowing/p/16395842.html