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

Question

我想在python(using numpy)中实现一个函数,它接受一个数学函数(p(x) = e^(-x)例如下面的例子)作为输入并生成随机数,这些随机数根据数学函数的概率分布进行分配.我需要绘制它们,所以我们可以看到分布.

实际上我需要一个随机数生成器函数来完成以下两个数学函数作为输入,但如果它可以采取其他函数,为什么不:

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

有没有人知道如何在python中这样做？

Answer 1

对于像您需要的那些简单分布，或者如果您可以轻松地以封闭形式的 CDF 求反，您可以在 NumPy 中找到大量采样器，正如 Olivier 的回答中正确指出的那样。

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

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

基本思路是这样的：

• 从一个随机点开始x并采取随机步骤xnew = x + delta
• 评估起点p(x)和新起点的期望概率分布p(xnew)
• 如果新点更有可能p(xnew)/p(x) >= 1接受移动
• 如果新点是不太可能的随机决定是否接受或拒绝取决于如何可能¹新点
• 从这一点开始新的一步并重复循环

可以证明，例如，参见Sokal ²，用这种方法采样的点遵循接受概率分布。

可以在PyMC3包中找到在 Python 中广泛实现的 Montecarlo 方法。

示例实现

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

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 采样器等，还有太多要说的。我不认为这是合适的地方，有很多比什么更好的材料我可以在这里写在线可用。

1. 这里的想法是支持概率较高的探索，但仍然关注低概率区域，因为它们可能会导致其他峰值。基本的是提案分布的选择，即您如何选择新的探索点。步长太小可能会将您限制在分布的有限区域，太大可能会导致探索效率非常低。

2. 物理导向。贝叶斯形式主义（Metropolis-Hastings）现在是首选，但恕我直言，初学者有点难以掌握。网上有很多教程，例如杜克大学的这个教程。

3. 未显示的实现不会增加太多混淆，但很简单，您只需在域边缘包装试验步骤或使所需的函数在域外归零。

Answer 2

NumPy 提供了广泛的概率分布。

第一个函数是参数为 1的指数分布。

np.random.exponential(1)

第二个是一个正态分布均值为0，方差为1。

np.random.normal(0, 1)

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

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

更一般的分布

虽然 NumPy 可能会满足您的所有需求，但请记住，您始终可以计算分布的逆累积分布函数和来自均匀分布的输入值。

inverse_cdf(np.random.uniform())

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

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

如果您遇到不容易计算 CDF 的分布，请考虑filippo 的出色答案。