CODEWITHSUNDEEP
pythonrandomprobabilityprobability-distribution
user11847849
user11847849

Reputation:

Write a random number generator that, based on uniformly distributed numbers between 0 and 1, samples from a Lévy-distribution?

I'm completely new to Python. Could someone show me how can I write a random number generator which samples from the Levy Distribution? I've written the function for the distribution, but I'm confused about how to proceed further! The random numbers generated by this distribution I want to use them to simulate a 2D random walk.

I'm aware that from scipy.stats I can use the Levy class, but I want to write the sampler myself.

import numpy as np
import matplotlib.pyplot as plt

# Levy distribution
"""
    f(x) = 1/(2*pi*x^3)^(1/2) exp(-1/2x)
"""
def levy(x):
    return 1 / np.sqrt(2*np.pi*x**3) * np.exp(-1/(2*x))

N = 50
foo = levy(N)

Upvotes: 2

Views: 883

Answers (2)

Severin Pappadeux
Severin Pappadeux

Reputation: 20080

@pjs code looks ok to me, but there is a discrepancy between his code and what SciPy thinks about Levy - basically, sampling is different from PDF.

Code, Python 3.8 Windows 10 x64

import numpy as np
from scipy.stats import levy
from scipy.stats import norm
import matplotlib.pyplot as plt

rng = np.random.default_rng(312345)

# Arguments
#   u: a uniform[0,1) random number
#   c: scale parameter for Levy distribution (defaults to 1)
#   mu: location parameter (offset) for Levy (defaults to 0)
def my_levy(u, c = 1.0, mu = 0.0):
    return mu + c / (2.0 * (norm.ppf(1.0 - u))**2)

fig, ax = plt.subplots()

rnge=(0, 20.0)

x = np.linspace(rnge[0], rnge[1], 1001)

N = 200000
q = np.empty(N)

for k in range(0, N):
    u = rng.random()
    q[k] = my_levy(u)

nrm = levy.cdf(rnge[1])
ax.plot(x, levy.pdf(x)/nrm, 'r-', lw=5, alpha=0.6, label='levy pdf')
ax.hist(q, bins=100, range=rnge, density=True, alpha=0.2)
plt.show()

produce graph

enter image description here

UPDATE

Well, I tried to use home-made PDF, same output, same problem

# replace levy.pdf(x) with PDF(x)
def PDF(x):
    return np.where(x <= 0.0, 0.0, 1.0 / np.sqrt(2*np.pi*x**3) * np.exp(-1./(2.*x)))

UPDATE II

After applying @pjs corrected sampling routine, sampling and PDF are aligned perfectly. New graph

enter image description here

Upvotes: 2

pjs
pjs

Reputation: 19855

Here's a straightforward implementation of the generating algorithm for the Levy distribution found on Wikipedia:

import random
from scipy.stats import norm

# Arguments
#   u: a uniform[0,1) random number
#   c: scale parameter for Levy distribution (defaults to 1)
#   mu: location parameter (offset) for Levy (defaults to 0)
def my_levy(u, c = 1.0, mu = 0.0):
    return mu + c / (2 * norm.ppf(1.0 - u)**2)

# Generate a handful of samples
for _ in range(10):
    print(my_levy(random.random()))

I don't normally use Python, so please suggest improvements.

ADDENDUM

Kudos to Severin Pappadeux for the work in his response. I had already noted that a simpler answer would be to take the inverse of a squared Gaussian, but Advaita had asked for an explicit function of U ~ Uniform(0,1) so I didn't pursue that. It turns out that I should have. The Wikipedia cite mentions that, but without the scale factor of 2 in the denominator. When I take the 2 out of the implementation of Wikipedia's generating algorithm, i.e. change the implemention to

def my_levy(u, c = 1.0, mu = 0.0):
    return mu + c / (norm.ppf(1.0 - u)**2)

the resulting histogram aligns beautifully with the normalized plot of the pdf. (Note - I've now also edited the incorrect Wikipedia entry to correct the formula.)

Upvotes: 1

Related Questions