如何正确采样截断分布?
How to properly sample truncated distributions?
我正在尝试学习如何对截断分布进行采样。首先,我决定尝试在这里找到的一个简单示例 example
我不太理解 CDF 的除法,因此我决定稍微调整一下算法。被采样的是值 x>0
的指数分布 这是一个示例 python 代码:
# Sample exponential distribution for the case x>0
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
def pdf(x):
return x*np.exp(-x)
xvec=np.zeros(1000000)
x=1.
for i in range(1000000):
a=x+np.random.normal()
xs=x
if a > 0. :
xs=a
A=pdf(xs)/pdf(x)
if np.random.uniform()<A :
x=xs
xvec[i]=x
x=np.linspace(0,15,1000)
plt.plot(x,pdf(x))
plt.hist([x for x in xvec if x != 0],bins=150,normed=True)
plt.show()
Ant 的输出是:
上面的代码似乎只在使用条件 if a > 0. :
时工作正常,即正 x
,选择另一个条件(例如 if a > 0.5 :
)会产生错误的结果。
由于我的最终目标是在截断的间隔上对 2D-Gaussian - pdf 进行采样,因此我尝试使用指数分布扩展简单示例(请参见下面的代码)。不幸的是,由于简单的情况不起作用,我假设下面给出的代码会产生错误的结果。
我假设所有这些都可以使用 python 的高级工具来完成。但是,由于我的主要想法是了解背后的原理,我将非常感谢您帮助理解我的错误。
谢谢您的帮助。
编辑:
# code updated according to the answer of CrazyIvan
from scipy.stats import multivariate_normal
RANGE=100000
a=2.06072E-02
b=1.10011E+00
a_range=[0.001,0.5]
b_range=[0.01, 2.5]
cov=[[3.1313994E-05, 1.8013737E-03],[ 1.8013737E-03, 1.0421529E-01]]
x=a
y=b
j=0
for i in range(RANGE):
a_t,b_t=np.random.multivariate_normal([a,b],cov)
# accept if within bounds - all that is neded to truncate
if a_range[0]<a_t and a_t<a_range[1] and b_range[0]<b_t and b_t<b_range[1]:
print(dx,dy)
编辑:
我根据 this scheme 规范解析 pdf 并根据 @Crazy Ivan 和 @Leandro Caniglia 给出的答案更改了代码,针对删除了 pdf 底部的情况。这是除以 (1-CDF(0.5)),因为我的接受条件是 x>0.5
。这似乎再次显示出一些差异。谜团再次盛行..
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
def pdf(x):
return x*np.exp(-x)
# included the corresponding cdf
def cdf(x):
return 1. -np.exp(-x)-x*np.exp(-x)
xvec=np.zeros(1000000)
x=1.
for i in range(1000000):
a=x+np.random.normal()
xs=x
if a > 0.5 :
xs=a
A=pdf(xs)/pdf(x)
if np.random.uniform()<A :
x=xs
xvec[i]=x
x=np.linspace(0,15,1000)
# new part norm the analytic pdf to fix the area
plt.plot(x,pdf(x)/(1.-cdf(0.5)))
plt.hist([x for x in xvec if x != 0],bins=200,normed=True)
plt.savefig("test_exp.png")
plt.show()
似乎可以通过选择更大的移位大小来解决这个问题
shift=15.
a=x+np.random.normal()*shift.
这通常是 Metropolis - Hastings 的一期。见下图:
我也查了shift=150
最重要的是,改变移位大小肯定会提高收敛性。痛苦就是为什么,因为高斯是无界的。
要从整个密度函数 pdf
计算截断密度函数 pdf_t
,请执行以下操作:
- 设
[a, b]
为截断间隔; (x
轴)
- 让
A := cdf(a)
和B := cdf(b)
; (cdf
= 非截断累积分布函数)
- 然后
pdf_t(x) := pdf(x) / (B - A)
如果 x
在 [a, b]
和 0
其他地方。
在 a = -infinity
(resp. b = +infinity
)的情况下,取 A := 0
(resp. B := 1
)。
至于你看到的“玄机”
请注意,您的蓝色曲线错误。它不是截断分布的 pdf
,它只是未截断分布的 pdf
,按正确的量缩放(除以 1-cdf(0.5)
)。实际截断的 pdf
曲线从 x = 0.5
上的一条垂直线开始,一直上升直到到达您当前的蓝色曲线。换句话说,您只缩放了曲线但忘了截断它,在本例中是向左截断。这样的截断对应于上述算法中步骤 3 的“0
其他地方”部分。
你说你想学习对截断分布进行采样的基本思想,但你的来源是一个博客 post 关于
Metropolis–Hastings algorithm?你真的需要这个 "method for obtaining a sequence of random samples from a probability distribution for which direct sampling is difficult" 吗?以此为出发点,就好比通过阅读莎士比亚来学习英语。
正常截断
对于截断正态,基本拒绝抽样 是您所需要的:为原始分布生成样本,拒绝边界外的样本。正如 Leandro Caniglia 指出的那样,您不应期望截断分布具有相同的 PDF,除非间隔较短——这显然是不可能的,因为 PDF 图形下方的面积始终为 1。如果您从侧面切掉东西,则必须多在中间; PDF 被重新缩放。
当你需要100000个时,一个一个地收集样本是非常低效的。我会一次抓取100000个正常样本,只接受那些适合的;然后重复直到我有足够的。 amin 和 amax 之间截断正态采样的示例:
import numpy as np
n_samples = 100000
amin, amax = -1, 2
samples = np.zeros((0,)) # empty for now
while samples.shape[0] < n_samples:
s = np.random.normal(0, 1, size=(n_samples,))
accepted = s[(s >= amin) & (s <= amax)]
samples = np.concatenate((samples, accepted), axis=0)
samples = samples[:n_samples] # we probably got more than needed, so discard extra ones
这是与 PDF 曲线的比较,重新缩放 除以 cdf(amax) - cdf(amin)
,如上所述。
from scipy.stats import norm
_ = plt.hist(samples, bins=50, density=True)
t = np.linspace(-2, 3, 500)
plt.plot(t, norm.pdf(t)/(norm.cdf(amax) - norm.cdf(amin)), 'r')
plt.show()
截断多元正态分布
现在我们要保留amin 和amax 之间的第一个坐标,bmin 和bmax 之间的第二个坐标。同样的故事,除了将有一个 2 列数组并且以相对偷偷摸摸的方式完成与边界的比较:
(np.min(s - [amin, bmin], axis=1) >= 0) & (np.max(s - [amax, bmax], axis=1) <= 0)
这意味着:从每行中减去 amin、bmin 并仅保留两个结果均为非负的行(意味着我们有 a >= amin 和 b >= bmin)。也对 amax、bmax 做类似的事情。只接受同时满足这两个条件的行。
n_samples = 10
amin, amax = -1, 2
bmin, bmax = 0.2, 2.4
mean = [0.3, 0.5]
cov = [[2, 1.1], [1.1, 2]]
samples = np.zeros((0, 2)) # 2 columns now
while samples.shape[0] < n_samples:
s = np.random.multivariate_normal(mean, cov, size=(n_samples,))
accepted = s[(np.min(s - [amin, bmin], axis=1) >= 0) & (np.max(s - [amax, bmax], axis=1) <= 0)]
samples = np.concatenate((samples, accepted), axis=0)
samples = samples[:n_samples, :]
不打算绘制,但这里有一些值:自然地,在范围内。
array([[ 0.43150033, 1.55775629],
[ 0.62339265, 1.63506963],
[-0.6723598 , 1.58053835],
[-0.53347361, 0.53513105],
[ 1.70524439, 2.08226558],
[ 0.37474842, 0.2512812 ],
[-0.40986396, 0.58783193],
[ 0.65967087, 0.59755193],
[ 0.33383214, 2.37651975],
[ 1.7513789 , 1.24469918]])
我正在尝试学习如何对截断分布进行采样。首先,我决定尝试在这里找到的一个简单示例 example
我不太理解 CDF 的除法,因此我决定稍微调整一下算法。被采样的是值 x>0
的指数分布 这是一个示例 python 代码:
# Sample exponential distribution for the case x>0
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
def pdf(x):
return x*np.exp(-x)
xvec=np.zeros(1000000)
x=1.
for i in range(1000000):
a=x+np.random.normal()
xs=x
if a > 0. :
xs=a
A=pdf(xs)/pdf(x)
if np.random.uniform()<A :
x=xs
xvec[i]=x
x=np.linspace(0,15,1000)
plt.plot(x,pdf(x))
plt.hist([x for x in xvec if x != 0],bins=150,normed=True)
plt.show()
Ant 的输出是:
上面的代码似乎只在使用条件 if a > 0. :
时工作正常,即正 x
,选择另一个条件(例如 if a > 0.5 :
)会产生错误的结果。
由于我的最终目标是在截断的间隔上对 2D-Gaussian - pdf 进行采样,因此我尝试使用指数分布扩展简单示例(请参见下面的代码)。不幸的是,由于简单的情况不起作用,我假设下面给出的代码会产生错误的结果。
我假设所有这些都可以使用 python 的高级工具来完成。但是,由于我的主要想法是了解背后的原理,我将非常感谢您帮助理解我的错误。 谢谢您的帮助。
编辑:
# code updated according to the answer of CrazyIvan
from scipy.stats import multivariate_normal
RANGE=100000
a=2.06072E-02
b=1.10011E+00
a_range=[0.001,0.5]
b_range=[0.01, 2.5]
cov=[[3.1313994E-05, 1.8013737E-03],[ 1.8013737E-03, 1.0421529E-01]]
x=a
y=b
j=0
for i in range(RANGE):
a_t,b_t=np.random.multivariate_normal([a,b],cov)
# accept if within bounds - all that is neded to truncate
if a_range[0]<a_t and a_t<a_range[1] and b_range[0]<b_t and b_t<b_range[1]:
print(dx,dy)
编辑:
我根据 this scheme 规范解析 pdf 并根据 @Crazy Ivan 和 @Leandro Caniglia 给出的答案更改了代码,针对删除了 pdf 底部的情况。这是除以 (1-CDF(0.5)),因为我的接受条件是 x>0.5
。这似乎再次显示出一些差异。谜团再次盛行..
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
def pdf(x):
return x*np.exp(-x)
# included the corresponding cdf
def cdf(x):
return 1. -np.exp(-x)-x*np.exp(-x)
xvec=np.zeros(1000000)
x=1.
for i in range(1000000):
a=x+np.random.normal()
xs=x
if a > 0.5 :
xs=a
A=pdf(xs)/pdf(x)
if np.random.uniform()<A :
x=xs
xvec[i]=x
x=np.linspace(0,15,1000)
# new part norm the analytic pdf to fix the area
plt.plot(x,pdf(x)/(1.-cdf(0.5)))
plt.hist([x for x in xvec if x != 0],bins=200,normed=True)
plt.savefig("test_exp.png")
plt.show()
似乎可以通过选择更大的移位大小来解决这个问题
shift=15.
a=x+np.random.normal()*shift.
这通常是 Metropolis - Hastings 的一期。见下图:
我也查了shift=150
最重要的是,改变移位大小肯定会提高收敛性。痛苦就是为什么,因为高斯是无界的。
要从整个密度函数 pdf
计算截断密度函数 pdf_t
,请执行以下操作:
- 设
[a, b]
为截断间隔; (x
轴) - 让
A := cdf(a)
和B := cdf(b)
; (cdf
= 非截断累积分布函数) - 然后
pdf_t(x) := pdf(x) / (B - A)
如果x
在[a, b]
和0
其他地方。
在 a = -infinity
(resp. b = +infinity
)的情况下,取 A := 0
(resp. B := 1
)。
至于你看到的“玄机”
请注意,您的蓝色曲线错误。它不是截断分布的 pdf
,它只是未截断分布的 pdf
,按正确的量缩放(除以 1-cdf(0.5)
)。实际截断的 pdf
曲线从 x = 0.5
上的一条垂直线开始,一直上升直到到达您当前的蓝色曲线。换句话说,您只缩放了曲线但忘了截断它,在本例中是向左截断。这样的截断对应于上述算法中步骤 3 的“0
其他地方”部分。
你说你想学习对截断分布进行采样的基本思想,但你的来源是一个博客 post 关于 Metropolis–Hastings algorithm?你真的需要这个 "method for obtaining a sequence of random samples from a probability distribution for which direct sampling is difficult" 吗?以此为出发点,就好比通过阅读莎士比亚来学习英语。
正常截断
对于截断正态,基本拒绝抽样 是您所需要的:为原始分布生成样本,拒绝边界外的样本。正如 Leandro Caniglia 指出的那样,您不应期望截断分布具有相同的 PDF,除非间隔较短——这显然是不可能的,因为 PDF 图形下方的面积始终为 1。如果您从侧面切掉东西,则必须多在中间; PDF 被重新缩放。
当你需要100000个时,一个一个地收集样本是非常低效的。我会一次抓取100000个正常样本,只接受那些适合的;然后重复直到我有足够的。 amin 和 amax 之间截断正态采样的示例:
import numpy as np
n_samples = 100000
amin, amax = -1, 2
samples = np.zeros((0,)) # empty for now
while samples.shape[0] < n_samples:
s = np.random.normal(0, 1, size=(n_samples,))
accepted = s[(s >= amin) & (s <= amax)]
samples = np.concatenate((samples, accepted), axis=0)
samples = samples[:n_samples] # we probably got more than needed, so discard extra ones
这是与 PDF 曲线的比较,重新缩放 除以 cdf(amax) - cdf(amin)
,如上所述。
from scipy.stats import norm
_ = plt.hist(samples, bins=50, density=True)
t = np.linspace(-2, 3, 500)
plt.plot(t, norm.pdf(t)/(norm.cdf(amax) - norm.cdf(amin)), 'r')
plt.show()
截断多元正态分布
现在我们要保留amin 和amax 之间的第一个坐标,bmin 和bmax 之间的第二个坐标。同样的故事,除了将有一个 2 列数组并且以相对偷偷摸摸的方式完成与边界的比较:
(np.min(s - [amin, bmin], axis=1) >= 0) & (np.max(s - [amax, bmax], axis=1) <= 0)
这意味着:从每行中减去 amin、bmin 并仅保留两个结果均为非负的行(意味着我们有 a >= amin 和 b >= bmin)。也对 amax、bmax 做类似的事情。只接受同时满足这两个条件的行。
n_samples = 10
amin, amax = -1, 2
bmin, bmax = 0.2, 2.4
mean = [0.3, 0.5]
cov = [[2, 1.1], [1.1, 2]]
samples = np.zeros((0, 2)) # 2 columns now
while samples.shape[0] < n_samples:
s = np.random.multivariate_normal(mean, cov, size=(n_samples,))
accepted = s[(np.min(s - [amin, bmin], axis=1) >= 0) & (np.max(s - [amax, bmax], axis=1) <= 0)]
samples = np.concatenate((samples, accepted), axis=0)
samples = samples[:n_samples, :]
不打算绘制,但这里有一些值:自然地,在范围内。
array([[ 0.43150033, 1.55775629],
[ 0.62339265, 1.63506963],
[-0.6723598 , 1.58053835],
[-0.53347361, 0.53513105],
[ 1.70524439, 2.08226558],
[ 0.37474842, 0.2512812 ],
[-0.40986396, 0.58783193],
[ 0.65967087, 0.59755193],
[ 0.33383214, 2.37651975],
[ 1.7513789 , 1.24469918]])