python - 优化生成变量的拒绝方法

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

My first question: Is my probability plot made properly?

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

And the second, what is in the title. How to optimise that method? Is it possible to reduce the time and number of rejections?

当然，您掌握了 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)

你应该得到类似的东西