生成具有给定概率密度函数的随机数答案

【问题标题】：Generating random numbers with a given probability density function生成具有给定概率密度函数的随机数
【发布时间】：2020-10-22 09:15:25
【问题描述】：

我想指定一个分布的probability density function，然后在 Python 中从该分布中提取 N 个随机数。我该怎么做呢？

【问题讨论】：

你可以看看scipy.stats库。有很多有用的知名发行版可供您抽样。
如果您的分布来自直方图或非参数，那么我可以发布一个配方，您如何从这些分布中采样。

标签： python numpy scipy probability

【解决方案1】：

一般来说，您希望获得逆累积概率密度函数。一旦你有了它，那么沿着分布生成随机数就很简单了：

import random

def sample(n):
    return [ icdf(random.random()) for _ in range(n) ]

或者，如果你使用 NumPy：

import numpy as np

def sample(n):
    return icdf(np.random.random(n))

在这两种情况下，icdf 都是逆累积分布函数，它接受 0 到 1 之间的值并从分布中输出相应的值。

为了说明icdf 的性质，我们将以值 10 和 12 之间的简单均匀分布为例：

概率分布函数在 10 到 12 之间为 0.5，其他地方为零
累积分布函数为 0 低于 10（没有样本低于 10），1 高于 12（没有样本高于 12）并且在值之间线性增加（PDF 的积分）
逆累积分布函数只定义在0到1之间。0时为10，12时为1，值之间呈线性变化

当然，困难的部分是获得逆累积密度函数。这真的取决于你的分布，有时你可能有一个分析函数，有时你可能想求助于插值。数值方法可能很有用，因为数值积分可用于创建 CDF，而插值可用于反转它。

【讨论】：

【解决方案2】：

这是我的函数，用于检索根据给定概率密度函数分布的单个随机数。我使用了类似蒙特卡洛的方法。当然调用这个函数n次可以生成n个随机数。

    """
    Draws a random number from given probability density function.

    Parameters
    ----------
        pdf       -- the function pointer to a probability density function of form P = pdf(x)
        interval  -- the resulting random number is restricted to this interval
        pdfmax    -- the maximum of the probability density function
        integers  -- boolean, indicating if the result is desired as integer
        max_iterations -- maximum number of 'tries' to find a combination of random numbers (rand_x, rand_y) located below the function value calc_y = pdf(rand_x).

    returns a single random number according the pdf distribution.
    """
    def draw_random_number_from_pdf(pdf, interval, pdfmax = 1, integers = False, max_iterations = 10000):
        for i in range(max_iterations):
            if integers == True:
                rand_x = np.random.randint(interval[0], interval[1])
            else:
                rand_x = (interval[1] - interval[0]) * np.random.random(1) + interval[0] #(b - a) * random_sample() + a

            rand_y = pdfmax * np.random.random(1) 
            calc_y = pdf(rand_x)

            if(rand_y <= calc_y ):
                return rand_x

        raise Exception("Could not find a matching random number within pdf in " + max_iterations + " iterations.")

如果您不必检索大量随机变量，我认为此解决方案的性能优于其他解决方案。另一个好处是您只需要 PDF 并避免计算 CDF、逆 CDF 或权重。

【讨论】：