如何正确地截取截断的分布？

Question

我正在尝试学习如何对截断的分布进行采样。首先，我决定尝试一个我在这里找到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，根据@Crazy Ivan和@Leandro Caniglia给出的答案，通过对分析pdf进行规范来改变代码，以便删除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.

这通常是大都市 - 黑斯廷斯的问题。见下图：

我还检查了shift=150

底线是改变移位大小肯定会改善收敛。痛苦是为什么，因为高斯是无界的。

Answer 1

你说你想学习采样截断分布的基本思想，但你的来源是关于Metropolis–Hastings algorithm的博客文章？您是否真的需要这种“从直接采样难以获得的概率分布中获取随机样本序列的方法”？以此为出发点就像阅读莎士比亚一样学习英语。

Truncated normal

对于截断的正常，基本拒绝采样就是您所需要的：为原始分布生成样本，拒绝那些超出边界的样本。正如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()

Truncated multivariate normal

现在我们要保持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]])

Answer 2

要从整个密度函数pdf_t计算截断密度函数pdf，请执行以下操作：

1. 让[a, b]成为截断区间; （x轴）
2. 让A := cdf(a)和B := cdf(b); （cdf =非截断累积分布函数）
3. 然后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”部分。