如何生成具有预定义概率分布的随机数？

Question

我想在 python 中实现一个函数（使用 numpy），它将数学函数（例如下面的 p(x) = e^(-x)）作为输入并生成随机数，是根据该数学函数的概率分布分布的。我需要绘制它们，以便我们可以看到分布。

我实际上需要一个随机数生成器函数来准确地输入以下 2 个数学函数，但如果它可以采用其他函数，为什么不呢：

1) p(x) = e^(-x)
2) g(x) = (1/sqrt(2*pi)) * e^(-(x^2)/2)

有谁知道这在 python 中是如何实现的？

Answer 1

NumPy 提供 a wide range of probability distributions。

第一个函数是带有参数 1 的 exponential distribution。

np.random.exponential(1)

第二个是 normal distribution，均值为 0，方差为 1。

np.random.normal(0, 1)

请注意，在这两种情况下，参数都是可选的，因为这些是这些分布的默认值。

作为旁注，您还可以在 random 模块中找到这些分布，分别为 random.expovariate 和 random.gauss。

更一般的分布

虽然 NumPy 可能会满足您的所有需求，但请记住，您始终可以计算逆 cumulative distribution function of your distribution and input values from a uniform distribution。

inverse_cdf(np.random.uniform())

例如，如果 NumPy 没有提供 指数分布，您可以这样做。

def exponential():
    return -np.log(-np.random.uniform())

如果遇到CDF不好计算的分布，那么考虑.

Answer 2

对于您需要的简单分布，或者如果您有一个易于反转的封闭形式 CDF，您可以在 NumPy 中找到大量采样器，正如 Olivier 的回答中正确指出的那样。

对于任意分布，您可以使用马尔可夫链蒙特卡洛抽样方法。

这些算法的最简单且可能更容易理解的变体是 Metropolis 采样。

基本思路是这样的：

从一个随机点开始x并采取随机步骤xnew = x + delta
在起点 p(x) 和新起点 p(xnew)
如果新点更有可能 p(xnew)/p(x) >= 1 接受移动
如果新点的可能性较小，则根据新点的可能性¹随机决定是接受还是拒绝
从此点开始新步骤并重复循环

可以显示，例如Sokal²，即用该方法采样的点服从接受概率分布。

Python 中蒙特卡洛方法的广泛实现可以在 PyMC3 包中找到。

示例实现

这是一个玩具示例，只是为了向您展示基本思想，并不意味着以任何方式作为参考实现。任何认真的工作请参考成熟的包。

def uniform_proposal(x, delta=2.0):
    return np.random.uniform(x - delta, x + delta)

def metropolis_sampler(p, nsamples, proposal=uniform_proposal):
    x = 1 # start somewhere

    for i in range(nsamples):
        trial = proposal(x) # random neighbour from the proposal distribution
        acceptance = p(trial)/p(x)

        # accept the move conditionally
        if np.random.uniform() < acceptance:
            x = trial

        yield x

让我们看看它是否适用于一些简单的发行版

高斯混合

def gaussian(x, mu, sigma):
    return 1./sigma/np.sqrt(2*np.pi)*np.exp(-((x-mu)**2)/2./sigma/sigma)

p = lambda x: gaussian(x, 1, 0.3) + gaussian(x, -1, 0.1) + gaussian(x, 3, 0.2)
samples = list(metropolis_sampler(p, 100000))

柯西

def cauchy(x, mu, gamma):
    return 1./(np.pi*gamma*(1.+((x-mu)/gamma)**2))

p = lambda x: cauchy(x, -2, 0.5)
samples = list(metropolis_sampler(p, 100000))

任意函数

您实际上不必从适当的概率分布中抽样。您可能只需要强制执行一个有限的域，在该域中对随机步骤进行采样³

p = lambda x: np.sqrt(x)
samples = list(metropolis_sampler(p, 100000, domain=(0, 10)))

p = lambda x: (np.sin(x)/x)**2
samples = list(metropolis_sampler(p, 100000, domain=(-4*np.pi, 4*np.pi)))

结论

关于提案分布、收敛性、相关性、效率、应用程序、贝叶斯形式主义、其他 MCMC 采样器等，还有太多话要说。我不认为这是合适的地方，而且有很多比我在这里写的在线可用的更好的东西 material。

这里的想法是支持概率较高的探索，但仍然关注低概率区域，因为它们可能会导致其他峰值。基础是提议分布的选择，即你如何选择新的点来探索。步长太小可能会将您限制在分布的有限区域，步长太大可能会导致探索效率非常低下。
物理方向。如今，贝叶斯形式主义 (Metropolis-Hastings) 是首选，但恕我直言，它对初学者来说有点难掌握。网上有很多教程，请参阅例如this one 来自杜克大学。
未显示的实现不会增加太多混乱，但很简单，您只需在域边缘包装试验步骤或使所需函数在域外变为零。

如何生成具有预定义概率分布的随机数？

How to generate random numbers with predefined probability distribution?

python

numpy

probability

probability-density

probability-distribution

更一般的分布

示例实现

高斯混合

柯西

任意函数

结论