如何根据自定义概率密度函数(Python)生成随机数？

Question

我有一个包含随机变量X和它们出现的分数的列表，所以如果我绘制这些，我得到一个概率密度函数。我想知道如何利用这个概率密度函数来产生一些随机数？

我使用scipy.interpolate.CubicSpline来获得这个数据的Python函数。如何使用这个函数生成随机数？

Answer 1

为了重新表达你的问题，你已经想出了一个pdf (一个包含随机变量X和它们发生的分数的列表)，并想知道如何从有这个pdf的分布中抽取随机样本。根据你想要的正式程度，有两种方法(我知道)可以做到这一点。

TLDR:对于简单的情况，可以使用NumPy实现，因为它干净、简单和快速。如果您想要一个更正式的版本，因为您使用的是更大的统计框架，那么也许SciPy版本更适合。

SciPy

如果您希望它适合SciPy分发框架，那么您可以使用rv_discrete类并对其进行扩展。在您的例子中，这看起来应该是：

from scipy.stats import rv_discrete

# these are your variables X
vals = [1, 2, 3]
# these are the fractions they occur
probs = [0.2, 0.5, 0.3]

# define discrete distribution
distrib = rv_discrete(values=(range(len(vals)), probs)) 

# sample 10 values from this distribution
distrib.rvs(size=10)
array([1, 0, 1, 2, 1, 1, 0, 1, 1, 1])

# distrib outputs indices in vals, not actual vals
[vals[x] for x in distrib.rvs(size=10)]
[3, 2, 3, 2, 2, 2, 1, 1, 2, 2]

并进行快速的速度测试，以获得良好的效果：

%timeit [vals[x] for x in distrib.rvs(size=10000)]
2.34 ms ± 195 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

NumPy

正如注释中提到的，您只需直接使用NumPy函数即可。

import numpy as np
np.random.choice(vals, size=10, p=probs)
array([2, 2, 1, 2, 2, 2, 2, 3, 1, 2])

虽然它不是SciPy发行框架的一部分，但它简单而干净，如下所示，速度更快：

%timeit np.random.choice(vals, size=10000, p=probs)
639 µs ± 204 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

NumPy比SciPy的速度提高的部分原因是它们产生的伪随机数的不同，而伪随机数是采样过程的基础。NumPy已将它们的默认PRNG进程升级为随机数生成器的PCG家族，而SciPy仍在使用默森-图斯特尔。他们宣布了这个这里，如果你对它的工作方式感到好奇的话，我已经写了一个简单的解释这里。专家之间还有更多更详细的这里。

通过将NumPy PRNG传递给SciPy，我们可以看到这种速度提高的影响：

# default SciPy
distrib = rv_discrete(values=(range(len(vals)), probs))
%timeit [vals[x] for x in distrib.rvs(size=1000000)]
358 ms ± 204 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

# SciPy but we pass NumPy's new PCG PRNG
np_seed = np.random.default_rng(123)
distrib = rv_discrete(values=(range(len(vals)), probs), seed=np_seed)
%timeit [vals[x] for x in distrib.rvs(size=1000000)]
221 ms ± 15.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

因此，使用NumPy的PRNG是关于1.5x更快。