python - 编写一个随机数生成器，基于 0 到 1 之间的均匀分布数，从 Lévy 分布中采样？

Question

我对 Python 完全陌生。有人可以告诉我如何编写一个从 Levy Distribution 中采样的随机数生成器吗？ ？我已经编写了分发函数，但我对如何进一步进行感到困惑! 我想用这个分布生成的随机数来模拟 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)

Answer 1

@pjs 代码对我来说看起来不错，但他的代码与 SciPy 对 Levy 的看法之间存在差异 - 基本上，采样与 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)))

更新二

应用@pjs修正的采样例程后，采样和PDF完美对齐。新图表