Landau 上的拒绝方法 pdf

Question

我一直在尝试使用拒绝接受算法为landau PDF 创建一个随机数生成器。问题是编译似乎需要很长时间，并且生成的值似乎不遵循 PDF。 这是我尝试使用的代码：

import matplotlib.pyplot as plt
import numpy  as np
import time
from sympy import sqrt,exp,pi

a=0   # xmin
b=150 # xmax

h=0.25 # ymax
variables = [] #list for variables
def f(x):
    return 1/(sqrt(2*pi))*exp(-1/2*((x-25)+exp(-(x-25))))  #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,h)

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

print("Time: ", end-start)
print("Rejection: ", reject)
xx = np.linspace(a,b,150)

plt.hist(variables,50, density=1)
plt.show()

欢迎提出任何意见

Answer 1

为什么同情？？？拒绝方法总是有性能问题，如果你把它和 sympy 结合起来......这不是一个好主意。我建议您将上面脚本的第一个块替换为：

import matplotlib.pyplot as plt
import numpy  as np
import time
import random
from scipy import sqrt,exp,pi

然后至少它会运行（一些您可以轻松修复的警告）并产生正确的输出：

Time:  61.50497007369995
Rejection:  3671570

为什么慢？

因为您丢弃概率为1-f(u1)/0.25 的随机值，该概率仅在您的分布峰值u1 ~ 25 附近很小。它实际上在数字上变成了1.0，远远早于你到达u1=150，因为f2(150)=12.8675100380087668e-28！这意味着您会丢弃 100% 的随机数尝试以获取更大的 u1 值，这是非常低效的。

加快代码速度的一个好方法是缩小u1 的范围。只需更改 b=50 即可：

Time:  20.477684020996094
Rejection:  1145527

--> 快三倍，结果相同。

就这些了吗？

不，如果你真的关心速度，你不应该总是在 0 到 0.25 的区间上生成 u2。您可以使用任何上限，只要它保持小于实际的朗道分布即可。它必须是一个“信封”。由于您降低了拒绝概率，这可能会大大增加计算时间。对于朗道犬来说，这并不简单，因为它的尾巴很长。但是您肯定可以至少在u1 的特定范围内找到一些东西。

我刚刚测试的一个相当愚蠢的方法是将您的问题分成四个区间 [a,b) = [0,27), [27,30), [30,35), [35,50) 有上限在 0.25、0.11、0.05、0.01 的 u2 上，获得几乎 2 倍加速的另一个因子：

Time:  12.69986867904663
Rejection:  661561

当然，同样，统计结果相同：

通过一些工作，您可以获得比这更好的。编辑后的示例原理验证代码是（这可以做得更好/更好）：

import matplotlib.pyplot as plt
import numpy  as np
import time
import random
from scipy import sqrt,exp,pi

a=0   # xmin
b=50 # xmax

#h=0.25 # ymax
variables = [] #list for variables
envelope = []
def f(x):
    return 1/(sqrt(2*pi))*exp(-1/2*((x-25)+exp(-(x-25))))  #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)
    h = 0.01
    if (u1<35): h=0.05
    if (u1<30): h=0.11
    if (u1<27): h=0.25
    
    u2 = random.uniform(0,h)

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

print("Time: ", end-start)
print("Rejection: ", reject)
xx = np.linspace(a,b,150)

plt.hist(variables, 50, weights=envelope,  density=1)
plt.show()

打败它

看完下面的cmets我想分享一下我充分利用numpy向量化的终极表现：

import matplotlib.pyplot as plt
import numpy  as np
import time

a=0   # xmin
b=50 # xmax

def f(x):
    return 1/(np.sqrt(2*np.pi))*np.exp(-1/2*((x-25)+np.exp(-(x-25))))

start = time.time()

n = 100000
u1 = np.random.rand(n)*(b-a)+a
h = np.full(n, 0.01)
h[u1<35] = 0.05
h[u1<30] = 0.11
h[u1<27] = 0.25
u2 = np.random.rand(n) * h

accept = u2 <= f(u1)
reject = n - np.count_nonzero(accept)

end = time.time()

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

plt.hist(u1[accept], 50, weights=h[accept],  density=1)
plt.show()

这个输出：

Time:  0.015883445739746094
Rejection:  86907

与我的初始版本相比，这是 3892 倍的性能提升，我相信与 sympy 相比，它甚至不止于此。了解我们编码的内容、原因和方式非常重要。

其实我很感兴趣：有没有人可以进一步改进这个最新版本的代码？