Блог

Я совершенно новичок в Python. Может кто-нибудь показать мне, как я могу написать генератор случайных чисел, который выбирает из распределения Леви? Я написал функцию для распределения, но я не знаю, как действовать дальше! Случайные числа, сгенерированные этим распределением, я хочу использовать их для имитации случайного блуждания 2D.

Я знаю, что из scipy.stats я могу использовать класс Levy, но я хочу сам написать сэмплер.

 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)
  

код @pjs выглядит нормально для меня, но есть несоответствие между его кодом и тем, что SciPy думает о Леви — в основном, выборка отличается от PDF.

Код, 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()
  

создайте график

Обновить

Ну, я попытался использовать самодельный PDF, тот же результат, та же проблема

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

ОБНОВЛЕНИЕ II

После применения исправленной процедуры выборки @pjs выборка и PDF идеально выровнены. Новый график

Вот простая реализация алгоритма генерации для распределения Леви, найденного в Википедии:

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

Обычно я не использую Python, поэтому, пожалуйста, предложите улучшения.

ДОБАВЛЕНИЕ

Благодарность Северину Паппаде за работу в его ответе. Я уже отмечал, что более простым ответом было бы взять обратную квадрату гауссова, но Адвайта попросил явную функцию U ~ Uniform(0,1), поэтому я не стал этого делать. Оказывается, я должен был. В Википедии упоминается об этом, но без масштабного коэффициента 2 в знаменателе. Когда я беру 2 из реализации алгоритма генерации Википедии, т. Е. Изменяю реализацию на

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

результирующая гистограмма прекрасно совпадает с нормализованным графиком pdf. (Примечание — теперь я также отредактировал неправильную запись в Википедии, чтобы исправить формулу.)