python - 从 N 个元素的拒绝采样中获得 N 个元素输出

Question

我正在对给定概率密度函数 pdf 的 N 个元素应用拒绝采样。当对 N 个元素应用此方法时，您可能会返回一个值数组，该数组与您正在评估的 N 个数字相比，其元素数量较少，这是由于应用拒绝方法而不循环条件 condition 的值而导致的 是 False。

为了协调这一点，我可以尝试循环不满足条件的值，直到它们为True。但是，我不确定如何循环条件，直到数组中的元素数量与给定我定义的函数的值数量 N 具有相同的长度。

import numpy

N = 1000 # number of elements
x = np.linspace(0, 200, N) 
pdf = pdf(x) # some pdf


# Rejection Method #1
# -------------------
fx = np.random.random_sample(size=N) * x.max() # uniform random samples scaled out
u = np.random.random_sample(size=N) # uniform random sample
condition = np.where(u <= pdf/pdf.max())[0] # Run first rejection criterion that returns bool values
x_arr = fx[condition] 

# Here, len(x_arr) < N, so I want to fix it until len(x_arr) == N
while len(x_arr) < N:
    ...
    if len(x_arr) == N:
        break

此后，我无法形成迭代方法，直到 len(x_arr) = N。

Answer 1

这是使用 bool 和高级索引的一种方法。它保留一个被拒绝的值的索引列表，并重新绘制这些值，直到列表为空。

采样和接受/拒绝函数示例:

def sample(N):
    return np.random.uniform(-3, 3, (N,))

def accept(v):
    return np.random.rand(v.size) < stats.norm().pdf(v)

主循环:

def draw(N, f_sample, f_accept):
    out = f_sample(N)
    mask = f_accept(out)
    reject, = np.where(~mask)
    while reject.size > 0:
        fill = f_sample(reject.size)
        mask = f_accept(fill)
        out[reject[mask]] = fill[mask]
        reject = reject[~mask]
    return out

健全性检查:

>>> counts, bins = np.histogram(draw(100000, sample, accept))
>>> counts / stats.norm().pdf((bins[:-1] + bins[1:]) / 2)
array([65075.50020815, 65317.17811578, 60973.84255365, 59440.53739031,
       58969.62310004, 59267.33983256, 60565.1928325 , 61108.60840388,
       64303.2863583 , 68293.86441234])

看起来大致平坦，所以还可以。