Apr-20-2019, 01:42 PM
Hello, 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!