numpy - 有效地过滤随机约束的随机样本

Question

我希望得到n个随机样本numpy，过滤以满足标准。我对目前的实施不满意;它对于大的N值（比如100,000）来说太慢了。如何更有效地过滤这些样品，以满足相关标准均匀随机样品小于f / g的条件？必须有更快的方法来实现此代码。

import numpy as np
from scipy.special import gamma
import matplotlib.pyplot as plt
def f(x): return 1. / gamma(3) * x * np.exp(-1 * x)
lambd = .2
c = 1 / lambd / gamma(3) * (2./(1-lambd)) ** 2 * np.exp(-1 * (1 - lambd) * (2. / (lambd - 1)))
def g(x): return c * lambd * np.exp(-1 * lambd * x)
x = np.linspace(0, 50, 1000)
samples = []
N = 100
while len(samples) < N:
    randou = np.random.uniform(0, 1)
    randoh = c * np.random.exponential(0.2)
    if randou <= f(randoh) / g(randoh): samples.append(randoh)
plt.hist(samples, 100, normed=True, label='Simulated PDF')
plt.plot(x, f(x), label='True PDF', lw=2)
plt.xlim(0, 10)
plt.show()

我还尝试一次性生成样本然后在while循环中过滤那些样本，但我不确定这个方法实际上有多快：

samples = np.random.uniform(0, 1, 100000)
hsamps = c * np.random.exponential(0.2, 100000)
N = 100
idx = np.array([True, False])
while len(idx[idx==True]) > 0:
    idx = samples > ( f(hsamps) / g(hsamps))
    samples[idx] = np.random.uniform(0, 1, len(idx[idx==True]))
    hsamps[idx] = c * np.random.exponential(0.2, len(idx[idx==True]))

Answer 1

要利用NumPy的速度，您需要使用大型数组而不是循环处理的单个标量。例如，你可以像这样生成N样本：

    randous = np.random.uniform(0, 1, size=N)
    randohs = c * np.random.exponential(0.2, size=N)

然后选择那些通过你的过滤器，如下所示：

    mask = randous <= f(randohs) / g(randohs)
    return randohs[mask]

唯一的问题是无法保证randohs[mask]具有所需数量的值（或任何值）。所以我们可能会重复这个，直到我们生成足够的样本：

while len(samples) < N:
    randohs = generate_samples()
    samples.extend(randohs)
samples = samples[:N]

尽管使用了while循环，但这仍然比一次生成一个样本快得多。

---

import numpy as np
from scipy.special import gamma
import matplotlib.pyplot as plt

def f(x): 
    return 1. / gamma(3) * x * np.exp(-1 * x)
def g(x): 
    return c * lambd * np.exp(-1 * lambd * x)

def generate_samples(N=10**5):
    randous = np.random.uniform(0, 1, size=N)
    randohs = c * np.random.exponential(0.2, size=N)
    mask = randous <= f(randohs) / g(randohs)
    return randohs[mask]

lambd = .2
c = (1 / lambd / gamma(3) * (2./(1-lambd)) ** 2 
     * np.exp(-1 * (1 - lambd) * (2. / (lambd - 1))))
x = np.linspace(0, 50, 1000)

samples = []
N = 10**5
while len(samples) < N:
    randohs = generate_samples()
    samples.extend(randohs)
samples = samples[:N]

plt.hist(samples, 100, density=True, label='Simulated PDF')
plt.plot(x, f(x), label='True PDF', lw=2)
plt.xlim(0, 10)
plt.show()