Scipy rv_continuous 错误地从分布中生成样本

Question

from scipy import stats
import numpy as np 

class your_distribution(stats.rv_continuous):
    def _pdf(self, x):
        p0 = 10.9949
        p1 = 0.394447
        p2 = 12818.4
        p3 = 2.38898

        return ((p1*p3)/(p3*p0+p2*p1))*((p0*np.exp(-1.0*p1*x))+(p2*np.exp(-1.0*p3*x)))

distribution = your_distribution(a=0.15, b=10.1)
sample = distribution.rvs(size=50000)

上面的代码从 0.15 到 10.1 范围内的标准化 pdf 生成 50000 个样本。但是，在上限b=10.1 处生成了不成比例的大量样本。这没有意义，正如绘制 pdf 时所看到的那样。

我该如何解决这个问题？

Answer 1

PDF 已针对整个分布范围正确标准化。然而，设置a 和b 只是简单地剪切PDF 而不进行任何重新规范化。使用(a=0.15, b=10.1)，PDF 不再集成为 1，并且由于 scipy 实现的一个怪癖，剩余的密度显然被添加到了范围的末尾。这会导致上界的样本数量过多。

我们可以通过绘制 a=0 和 a=0.15 的累积密度函数 (CDF) 来可视化正在发生的事情：

x = np.linspace(0, 15, 1000)

distribution = your_distribution(a=0.0, b=10.1)
plt.plot(x, distribution.cdf(x), label='a=0')

distribution = your_distribution(a=0.15, b=10.1)
plt.plot(x, distribution.cdf(x), label='a=0.15')

plt.legend()

为了消除 CDF 中的跳跃和上限范围内的虚假样本，我们需要针对 a..b 范围重新规范化 PDF。我懒得去分析算出正确的因数，所以让 scipy 做艰苦的工作：

from scipy import stats
from scipy.integrate import quad
import numpy as np

# I pulled the definition of the PDF out of the class so we can use it to
# compute the scale factor.
def pdf(x):
    p0 = 10.9949
    p1 = 0.394447
    p2 = 12818.4
    p3 = 2.38898

    return ((p1*p3)/(p3*p0+p2*p1))*((p0*np.exp(-1.0*p1*x))+(p2*np.exp(-1.0*p3*x)))    

class your_distribution(stats.rv_continuous):        
    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)

        # integrate area of the PDF in range a..b
        self.scale, _ = quad(pdf, self.a, self.b)

    def _pdf(self, x):
        return pdf(x) / self.scale  # scale PDF so that it integrates to 1 in range a..b 

distribution = your_distribution(a=0.15, b=10.1)
sample = distribution.rvs(size=1000)

如果你碰巧知道积分的解析解，你可以使用它而不是调用quad。