优化生成变量的拒绝方法

Question

我对生成连续随机变量的拒绝方法的优化有疑问。我有一个密度：f(x) = 3/2 (1-x^2)。这是我的代码：

import random
import matplotlib.pyplot as plt
import numpy  as np
import time
import scipy.stats as ss

a=0   # xmin
b=1   # xmax

m=3/2 # ymax
variables = [] #list for variables

def f(x):
    return 3/2 * (1 - x**2)  #probability density function

reject = 0   # number of rejections
start = time.time()
while len(variables) < 100000:  #I want to generate 100 000 variables
    u1 = random.uniform(a,b)
    u2 = random.uniform(0,m)

    if u2 <= f(u1):
        variables.append(u1)
    else:
        reject +=1
end = time.time()

print("Time: ", end-start)
print("Rejection: ", reject)
x = np.linspace(a,b,1000)
plt.hist(variables,50, density=1)
plt.plot(x, f(x))
plt.show()

ss.probplot(variables, plot=plt)
plt.show()

我的第一个问题：我的概率图是否正确？ 第二，标题中的内容。如何优化该方法？我想得到一些建议来优化代码。现在该代码大约需要 0.5 秒，并且大约有 50 000 次拒绝。是否可以减少拒绝的时间和次数？如果需要，我可以使用不同的生成变量的方法进行优化。

Answer 1

我的第一个问题：我的概率图是否正确？

没有。它与默认的正态分布相比较。您必须将您的函数f(x) 打包到派生自 stats.rv_continuous 的类中，使其成为 _pdf 方法，并将其传递给 probplot

第二个，标题是什么。如何优化该方法？是否可以减少拒绝的时间和次数？

当然，您拥有 NumPy 矢量能力的强大功能。永远不要编写显式循环 - vectoriz、vectorize 和 vectorize！

看看下面修改过的代码，不是一个循环，一切都是通过 NumPy 向量完成的。 100000 个样本（Xeon、Win10 x64、Anaconda Python 3.7）在我的计算机上的时间从 0.19 下降到 0.003。

import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
import time

a = 0.  # xmin
b = 1.  # xmax

m = 3.0/2.0 # ymax

def f(x):
    return 1.5 * (1.0 - x*x)  # probability density function

start  = time.time()

N = 100000
u1 = np.random.uniform(a, b, N)
u2 = np.random.uniform(0.0, m, N)

negs = np.empty(N)
negs.fill(-1)
variables = np.where(u2 <= f(u1), u1, negs) # accepted samples are positive or 0, rejected are -1

end = time.time()

accept = np.extract(variables>=0.0, variables)
reject = N - len(accept)

print("Time: ", end-start)
print("Rejection: ", reject)

x = np.linspace(a, b, 1000)
plt.hist(accept, 50, density=True)
plt.plot(x, f(x))
plt.show()

ss.probplot(accept, plot=plt) # against normal distribution
plt.show()

关于减少拒绝的数量，您可以使用 0 拒绝进行反向采样，它是三次方程，因此可以轻松使用

更新

这是用于probplot的代码：

class my_pdf(ss.rv_continuous):
    def _pdf(self, x):
        return 1.5 * (1.0 - x*x)

ss.probplot(accept, dist=my_pdf(a=a, b=b, name='my_pdf'), plot=plt)

你应该得到类似的东西

Answer 2

关于您的第一个问题，scipy.stats.probplot 将您的样本与正态分布的分位数进行比较。如果您希望它与您的f(x) 分布的分位数进行比较，请查看probplot 的dist 参数。

就加快此采样过程而言，避免循环通常是要走的路。将 start = ... 和 end = ... 之间的代码替换为以下代码，对我来说速度提高了 20 倍以上。

n_before_accept_reject = 150000
u1 = np.random.uniform(a, b, size=n_before_accept_reject)
u2 = np.random.uniform(0, m, size=n_before_accept_reject)
variables = u1[u2 <= f(u1)]
reject = n_before_accept_reject - len(variables)

请注意，每次运行它时，这将为您提供大约 100000 个接受的样本。您可以稍微提高n_before_accept_reject 的值，以有效保证variables 将始终具有>100000 个可接受的值，然后在必要时将变量的大小限制为恰好返回100000。

Answer 3

其他人谈到了概率图，我将讨论拒绝算法的效率。

接受/拒绝方案基于 m(x)，一个“主要函数”。一个majorizing 函数应该有两个性质： 1) m(x)≥ f(x) ∀ x； 2) m(x)，当缩放为分布时，应该很容易从中生成值。 你使用了常数函数 m = 3/2，它满足了这两个要求，但没有非常紧密地绑定 f(x)。从零积分到一，面积为3/2。您的 f(x) 是一个有效的密度函数，其面积为 1。因此，∫f(x)) / ∫m(x)) = 1 / (3/2) = 2/3。换句话说，您从专业化函数生成的值中有 2/3 被接受，而您拒绝了 1/3 的尝试。

您需要一个 m(x)，它为 f(x) 提供更严格的界限。我选择了一条在 x = 1/2 处与 f(x) 相切的线。通过一点微积分得到斜率，我推导出了m(x) = 15/8 - 3x/2。

这个选择的 m(x) 的面积是 9/8，所以只有 1/9 的值会被拒绝。更多的微积分产生了基于这个 m(x) 的 x 的逆变换生成器是 x = (5 - sqrt(25 - 24U)) / 4，其中 U 是一个统一的 (0,1) 随机变量。

这是一个基于您的原始版本的实现。我将拒绝方案包装在一个函数中，并使用列表理解而不是附加到列表来创建值。如果你运行它，你会看到，它产生的拒绝比你原来的版本少得多。

import random
import matplotlib.pyplot as plt
import numpy  as np
import time
import math
import scipy.stats as ss

a = 0   # xmin
b = 1   # xmax

reject = 0   # number of rejections

def f(x):
    return 3.0 / 2.0 * (1.0 - x**2)  #probability density function

def m(x):
    return 1.875 - 1.5 * x

def generate_x():
    global reject
    while True:
        x = (5.0 - math.sqrt(25.0 - random.uniform(0.0, 24.0))) / 4.0
        u = random.uniform(0, m(x))
        if u <= f(x):
            return x 
        reject += 1    

start = time.time()
variables = [generate_x() for _ in range(100000)]
end = time.time()

print("Time: ", end-start)
print("Rejection: ", reject)
x = np.linspace(a,b,1000)
plt.hist(variables,50, density=1)
plt.plot(x, f(x))
plt.show()