python - 编写一个随机数生成器,它基于 0 和 1 之间的均匀分布数,从 Lévy 分布中采样?
问题描述
我对 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 对 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 完美对齐。新图表
推荐阅读
- linux - 使用“getopts”和“case”错误处理输入
- android - 如何将 ValueAnimation 放入线程中?
- ios - 如何使用 CPTLegen 作为按钮切换图形可见性
- bash - 如何在shell中提取由分隔符分隔的位置值
- java - 如何清理 Mina FTP 服务器中的会话
- octave - 无法在 Octave 中加载除 MATLAB 文件以外的文件格式
- r - 删除一些值后矩阵的 R 最大距离
- ios - 用于购买无期限云存储的应用内购买类型是什么
- javascript - 主机更改后 jQuery 无法正常运行。旧的服务器代码运行良好
- javascript - 在模态上使用 jqGrid 有什么条件吗?