lundi 12 avril 2021

Rejection method on a landau pdf

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