samedi 20 avril 2019

How to optimise the rejection method for generating variables?

I have a problem with optimisation the rejection method of generating continuous random variables. I've got a density: f(x) = 3/2 (1-x^2). Here's my code:

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()

My first question: Is my probability plot made properly? And the second, what is in the title. How to optimise that method? I would like to get some advice to optimise the code. Now that code take about 0.5 second and there are about 50 000 rejections. Is it possible to reduce the time and number of rejections? If it's needed I can optimise using different method of generating variables. Thanks in advance!




Aucun commentaire:

Enregistrer un commentaire