如何正确采样截断的分布？
Ale*_*ska 4 python random numpy probability mcmc
我正在尝试学习如何对截短的分布进行采样。首先，我决定尝试一个简单的示例，在这里找到示例
我不太了解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)
编辑：
我根据此方案通过对解析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.
这通常是大都会-黑斯廷斯（Hastings）的问题。参见下图：
我也检查了 shift=150
最重要的是，改变移位大小肯定会改善收敛。痛苦是为什么，因为高斯不受限制。
小智 5
您说您想学习采样截断分布的基本思想，但是您的来源是有关Metropolis-Hastings算法的博客文章 ？您实际上是否需要这种“从直接采样困难的概率分布中获取随机样本序列的方法”？以此为起点，就像通过阅读莎士比亚来学习英语。
截断法线
对于截断的正态，您只需要基本的拒绝采样：生成原始分布的样本，拒绝超出范围的样本。正如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并仅保留两个结果均为非负数的行（这意味着我们有一个> = 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]])