编写一个随机数生成器，它基于 0 和 1 之间的均匀分布的数字，从 Lévy 分布中抽取样本？

Question

我是 Python 的新手。有人可以告诉我如何编写一个从 Levy Distribution 中采样的随机数生成器吗？我已经为分发编写了函数，但我对如何进一步进行感到困惑！这个分布生成的随机数我想用它们来模拟2D随机游走。

我知道从 scipy.stats 开始我可以使用 Levy class，但我想自己编写采样器。

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

这是在维基百科上找到的 generating algorithm for the Levy distribution 的直接实现：

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，所以请提出改进建议。

附录

感谢 Severin Pappadeux 在他的回复中所做的工作。我已经注意到，一个更简单的答案是取平方高斯的倒数，但 Advaita 要求 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 的归一化图完美对齐。（注意 - 我现在还编辑了不正确的维基百科条目以更正公式。）

Answer 2

@pjs 代码在我看来还不错，但他的代码与 SciPy 对 Levy 的看法存在差异 - 基本上，抽样与 PDF 不同。

代码，Python3.8Windows10x64

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 完美对齐。新图表

编写一个随机数生成器，它基于 0 和 1 之间的均匀分布的数字，从 Lévy 分布中抽取样本？

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

python

random

probability

probability-distribution