使用scipy的rv_continuous方法创建自定义连续分布

Question

我正在尝试计算 E[f(x)]，这个值是我从数据中生成或估计的概率密度函数（pdf）。

文档中提到：

子类化

可以通过创建一个新的类，继承自 rv_continuous 类，并重新定义至少一个 _pdf 或 _cdf 方法（这些方法需要标准化到位置 0 和尺度 1），这样可以接收干净的参数（在 a 和 b 之间），并通过参数检查方法。

如果你的随机变量的正参数检查不正确，那么你还需要重新定义 _argcheck 方法。

所以我创建了一个子类，并定义了 _pdf，但每当我尝试调用：

print my_continuous_rv.expect(lambda x: x)

scipy 就会给我报错：

AttributeError: 'your_continuous_rv' object has no attribute 'a'

这很合理，因为我猜它在试图确定积分的下限，错误信息中也打印了：

lb = loc + self.a * scale

我尝试定义属性 self.a 和 self.b（我认为这是随机变量定义的范围/区间）：

self.a = float("-inf")
self.b = float("inf")

然而，当我这样做时，它又报错说：

if N > self.numargs:
AttributeError: 'your_continuous_rv' object has no attribute 'numargs'

我不太确定 numargs 应该是什么，但在查看了 scipy 的代码后，我发现有这么一行代码：

if not hasattr(self, 'numargs'):
    # allows more general subclassing with *args
    self.numargs = len(shapes)

我猜这就是我的函数应该接受的随机变量的形状。

目前我只在做一个非常简单的随机变量，它的可能值是一个浮点数。所以我决定把 numargs 硬编码为 1。但这又导致了 scipy 的更多报错。

所以，归根结底，我觉得文档中没有清楚说明我在子类化时需要做什么，因为我按照他们的说法重写了 _pdf，但之后它又让我定义 self.a，我硬编码了这个值，然后又让我定义 numargs，此时我觉得我并不清楚他们希望我如何子类化 rv_continuous。有没有人知道？我可以从我想要拟合的数据中生成 pdf，然后从 pdf 中获取期望值和其他东西，我还需要在 rv_continuous 中初始化什么才能让它正常工作？

Answer 1

我的猜测是，你定义了一个构造函数，使用了 __init__，但没有调用基类 rv_continuous 的构造函数。

如果你参考文档中的例子，并添加一个像这样的构造函数

from scipy.stats import rv_continuous
class gaussian_gen(rv_continuous):
    "Gaussian distribution"
    def __init__(self, **kwargs):
        pass
        #rv_continuous.__init__(self, **kwargs)
    def _pdf(self, x):
        return np.exp(-x**2 / 2.) / np.sqrt(2.0 * np.pi)

gaussian = gaussian_gen(name='gaussian')
print(gaussian.pdf(3))

那么调用 pdf 就会失败。如果你把调用 rv_continous.__init__ 的那一行注释掉，它就能正常工作，因为那样基类的构造函数会为这个实例添加所有必要的字段。

Answer 2

SciPy的 rv_histogram 方法可以让你输入一些数据，然后它会帮你计算出概率密度函数（pdf）、累积分布函数（cdf）以及生成随机数的方法。

Answer 3

由于历史原因，scipy中的分布是以实例的形式存在的，所以你需要有一个你自己子类的实例。举个例子：

>>> class MyRV(stats.rv_continuous):
...    def _pdf(self, x, k):
...      return k * np.exp(-k*x)
>>> my_rv = MyRV(name='exp', a=0.)     # instantiation

注意需要指定支持的范围：默认值是 a=-inf 和 b=inf。

>>> my_rv.a, my_rv.b
(0.0, inf)
>>> my_rv.numargs        # gets figured out automagically
1

一旦你指定了，比如说 _pdf，你就得到了一个可以使用的分布实例：

>>> my_rv.cdf(4, k=3)
0.99999385578764677
>>> my_rv.rvs(k=3, size=4)
array([ 0.37696127,  1.10192779,  0.02632473,  0.25516446])
>>> my_rv.expect(lambda x: 1, args=(2,))    # k=2 here
0.9999999999999999