1

采样（三）：重要性采样与接受拒绝采样

 2 years ago
source link: https://allenwind.github.io/blog/10466/
Go to the source link to view the article. You can view the picture content, updated content and better typesetting reading experience. If the link is broken, please click the button below to view the snapshot at that time.

本文介绍重要性采样与接受拒绝采样，并给出使用Laplace分布作为参考分布，采样正态分布的例子。

重要性采样（importance sampling）

在概率分布p(x)p(x)下计算f(x)f(x)均值，

∫Ωp(x)f(x)dx=Ex∼p(x)[f(x)]≈1nn∑i=1f(xi)∫Ωp(x)f(x)dx=Ex∼p(x)[f(x)]≈1n∑i=1nf(xi)

其中xi∼p(x)xi∼p(x)，然而在很多情况下p(x)p(x)并不是简单的分布，无法像均匀分布、正态分布一样直接采样。这种情况下重要性采样引入一个方便采样的分布q(x)q(x)，把均值计算改写为，

∫Ωp(x)f(x)dx=Ex∼p(x)[f(x)]=∫Ωq(x)p(x)q(x)f(x)dx=Ex∼q(x)[p(x)q(x)f(x)]≈1nn∑i=1p(xi)q(xi)f(xi)∫Ωp(x)f(x)dx=Ex∼p(x)[f(x)]=∫Ωq(x)p(x)q(x)f(x)dx=Ex∼q(x)[p(x)q(x)f(x)]≈1n∑i=1np(xi)q(xi)f(xi)

其中xi∼q(x)xi∼q(x)比较容易​，而p(xi)q(xi)f(xi)p(xi)q(xi)f(xi)​的计算没有任何难度。注意到q(xi)≠0q(xi)≠0​。p(xi)q(xi)p(xi)q(xi)​可以看着是样本的权重，这正是“重要性”来源。

p(x)p(x)为不容易采样的分布，计算f(x)f(x)的均值有，

Ex∼p(x)[f(x)]=∑xp(x)f(x)=∑xq(x)p(x)q(x)f(x)=Ex∼q(x)[p(x)q(x)f(x)]=Ex∼q(x)[w(x)⋅f(x)]Ex∼p(x)[f(x)]=∑xp(x)f(x)=∑xq(x)p(x)q(x)f(x)=Ex∼q(x)[p(x)q(x)f(x)]=Ex∼q(x)[w(x)⋅f(x)]

对每个样本加权，

w(x)=p(x)q(x)w(x)=p(x)q(x)

可以验证，权重符合归一性，

Ex∼q(x)[w(x)]=∑xq(x)p(x)q(x)=∑xp(x)=1Ex∼q(x)[w(x)]=∑xq(x)p(x)q(x)=∑xp(x)=1

接受-拒绝采样（accept-reject sampling）

假设有概率分布p(x)p(x)，称为目标分布，需要从中采样往往比较难，于是引入较简单的易于采样的分布q(x)q(x)，称为参考分布或建议分布，并定义如下函数，

α(x)=p(x)q(x)∈[0,1]α(x)=p(x)q(x)∈[0,1]

现在采样样本xi∼q(x)xi∼q(x)​​，采样随机数εi∼U[0,1]εi∼U[0,1]​​，如果εi≤α(xi)εi≤α(xi)​​，则接受样本xixi​​否则拒绝该样本，继续以上流程采样新样本，如此下去获得的样本集x1,x2,…,xnx1,x2,…,xn​​​，服从目标分布p(x)p(x)​。

q(x)q(x)为建议分布，p(x)p(x)为目标分布，满足如下关系，

M⋅q(x)≥p(x)M⋅q(x)≥p(x)

即M×q(x)M×q(x)​​要包络p(x)p(x)​​。

1=∫p(x)dx≈n∑i=1α(xi)q(xi)1=∫p(x)dx≈∑i=1nα(xi)q(xi)

因此，归一化后可获得p(x)p(x)的采样集，

α(xi)q(xi)n∑i=1α(xi)q(xi)α(xi)q(xi)∑i=1nα(xi)q(xi)

那么接受概率是多少呢？直觉上这明显取决于M⋅q(x)M⋅q(x)​ 包络 p(x)p(x)​​ 的质量，考虑到接受样本需要满足，

εi≤p(xi)M×q(xi),xi∼q(x),εi∼U[0,1]εi≤p(xi)M×q(xi),xi∼q(x),εi∼U[0,1]

注意下式是从q(x)q(x)​​中采样，接受概率为，

P(接受)=P(U<p(x)M⋅q(x))=∫q(x)p(x)M⋅q(x)dx=1MP(接受)=P(U<p(x)M⋅q(x))=∫q(x)p(x)M⋅q(x)dx=1M

因此，高质量的采样，需要 MM 尽可能小，否则从q(x)q(x)​采样出来的样本大部分都会被拒绝掉，提高了程序运行时的空转时间。

利用贝叶斯公式可以知道最后的采样分布（接受样本条件下得到的分布）为，

q(x|接受)=q(x)×P(接受|X=x)P(接受)=p(x)q(x|接受)=q(x)×P(接受|X=x)P(接受)=p(x)

接下来我们给出一个例子，使用Laplace分布作为参考分布，采样正态分布。

拉普拉斯分布采样正态分布

以拉普拉斯分布为参考分布，采样目标分布为正态分布p(x)p(x)

f(x;μ,σ)=1√2πσe−(x−μ)22σ2∣∣μ=0f(x;μ,σ)=12πσe−(x−μ)22σ2|μ=0

拉普拉斯分布的概率密度q(x)q(x)，

f(x∣μ,b)=12bexp(−|x−μ|b)∣∣μ=0f(x∣μ,b)=12bexp⁡(−|x−μ|b)|μ=0

均值和方差分别为μ,2b2μ,2b2。

计算M⋅q(x)M⋅q(x)​ 包络 p(x)p(x)​​ 中MM的最小值，

1√2πσexp(−x22σ2)≤M×12bexp(−|x|b)12πσexp⁡(−x22σ2)≤M×12bexp⁡(−|x|b)

取特殊值，解得包络值的下界，

M≥2b√2π⋅exp(12b2)M≥2b2π⋅exp⁡(12b2)

取等号获得最优包络的MM​值。

均匀分布采样正态分布

类似地，以均匀分布为参考分布，

1√2πσexp(−x22σ2)≤M×1b−a12πσexp⁡(−x22σ2)≤M×1b−a

这里正态分布参数σ=1σ=1​，认为x∈(−5σ,5σ)x∈(−5σ,5σ)足够（正态分布的kσkσ​法则），解得

M≥10√2πM≥102π

Python实现如下，

import numpy as np

def accept_reject_sampling(pf, qf, q_gen, M):
    """
    pf:目标采样概率密度函数
    qf:参考分布概率密度函数
    q_gen:参考分布采样生成器
    M:使得pf(x)<=qf(x)*M，qf(x)*M称为pf(x)的包络函数
    """
    for x in q_gen():
        u = np.random.uniform()
        if u < pf(x) / (M * qf(x)):
            yield x

以Laplace分布作为参考分布，采样正态分布，

以均匀分布作为参考分布，采样正态分布，

以上相关分析的代码实现已经开源。更详细的实验代码实现更新到：sampling-from-distribution 。

以上介绍重要性采样与接受拒绝采样。对于接受拒绝采样，给出了两个例子：

  • 均匀分布作为参考分布，采样正态分布
  • 拉普拉斯分布作为参考分布，采样正态分布

转载请包括本文地址：https://allenwind.github.io/blog/10466
更多文章请参考：https://allenwind.github.io/blog/archives/


report

Recommend

About Joyk


Aggregate valuable and interesting links.
Joyk means Joy of geeK