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

Question

我想在 python 中实现一个函数（使用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 的最佳答案。