从自定义连续概率密度函数生成随机数

Question

正如标题所述，我试图从自定义连续概率密度函数生成随机数，即：

0.001257 *x^4 * e^(-0.285714 *x)

为此，我使用（在Python3上）scipy.stats.rv_continuous然后rvs()来生成它们

from decimal import Decimal
from scipy import stats
import numpy as np

class my_distribution(stats.rv_continuous):
    def _pdf(self, x):
        return (Decimal(0.001257) *Decimal(x)**(4)*Decimal(np.exp(-0.285714 *x)))

distribution = my_distribution()
distribution.rvs()

注意，我使用Decimal来摆脱OverflowError: (34, 'Result too large')

尽管如此，我还是得到了一个错误RuntimeError: Failed to converge after 100 iterations

那里发生了什么事？实现我需要做的事情的正确方法是什么

Answer 1

我已经找到了你问题的原因

rvs默认情况下使用数值积分，这是一个缓慢的过程，在某些情况下可能会失败。您的PDF可能是其中一种情况，左侧增长时没有边界

因此，您应该按如下方式指定发行版的支持（以下示例显示支持在区间[-4,4]）：

distribution = my_distribution(a = -4, b = 4)

在这个时间间隔内，PDF将从上方进行绑定，从而允许积分（从而生成随机数）正常工作。请注意，默认情况下，rv_continuous假定整个实数行都支持该分布

但是，这只适用于您在此处提供的特定PDF，不一定适用于任意PDF

---

通常，当你只给你的rv_continuous子类一个PDF时，子类的rvs、mean等会非常慢，因为每次需要生成随机数或计算统计数据时，该方法都需要集成PDF。例如，随机数生成需要使用数值积分来积分PDF，并且此过程可能无法收敛，这取决于PDF

在将来处理任意分布的情况下，尤其是在速度很快的情况下，您将需要添加到使用自己的采样器的_rvs方法中。这方面的一个例子是我自己的^{}，当只给定PDF和采样域时，它通过数值反转生成随机数。另一个例子是answer to a related question中给出的更简单的拒绝采样器

另见我的章节“Random Numbers from an Arbitrary Distribution”