I have been trying to create a random number generator for a landau PDF using the rejection acceptance algorithm. The problem is it seems to take a very long time to compile and the generated values don't seem to follow the PDF. This is the code I tried to use:
import matplotlib.pyplot as plt
import numpy as np
import time
import math
import scipy.stats as ss
a = 0 # xmin
b = 150 # xmax
reject = 0 # number of rejections
def f(x):
return 1/(sqrt(2*pi))*exp(-1/2*((x-25)+exp(-(x-25)))) #probability density function
def m(x):
return 134168.643260437*sqrt(2)*(0.5*exp(25 - x) - 0.5)*exp(-0.5*x - 0.5*exp(25 - x))/sqrt(pi)
def generate_x():
global reject
while True:
x = random.uniform(0.0, 150)
u = random.uniform(0, m(x))
if m(x) <= f(x):
return x
reject += 1
start = time.time()
variables = [generate_x() for _ in range(1000)]
end = time.time()
print("Time: ", end-start)
print("Rejection: ", reject)
plt.hist(variables,50, density=1)
plt.plot(x, f(x))
plt.show()
Any thoughts would be welcome
Aucun commentaire:
Enregistrer un commentaire