如何创建线性分数分布作为自定义离散概率分布？

Question

我定义了以下自定义概率分布：

import scipy.stats as st

# parameters
a = 3 / 16
b = 1

class linear_fractional(st.rv_discrete):
    def _pdf(self, n):
        if (n == 0):
            return (a + b - 1) / (a + b)
        else:
            return (a * b ** (n - 1)) / (a + b) ** (n + 1)

LF = linear_fractional()
LF.rvs()

当我让我的脚本运行时，我收到一条冗长的错误消息：

Traceback (most recent call last):
File "C:/Users/thoma/PycharmProjects/Host_Parasite_Coevolution/Asymptotics.py", line 17, in <module> LF.rvs()
File "C:\Users\thoma\AppData\Local\Programs\Python\Python37-32\lib\site-packages\scipy\stats\_distn_infrastructure.py", line 2969, in rvs
    return super(rv_discrete, self).rvs(*args, **kwargs)

...

RecursionError: maximum recursion depth exceeded while calling a Python object

如果我改为使用LF.mean()，我会得到

Fatal Python error: Cannot recover from stack overflow.

有谁知道这是为什么以及如何解决这个问题？我必须为我的概率分布定义一个上限吗？

Answer 1

按照the docs 和this post 给出的示例，该方法需要进行一些修改。重要的是，因为它是一个离散分布，所以应该使用_pmf 而不是_pdf。此外，_pmf 将被 n 的 numpy 样式数组调用，n == 0 的工作方式完全不同。

当n == 0 时(a * b ** (n - 1)) / (a + b) ** (n + 1) 等于(a + b - 1) / (a + b)，我们可以对所有n 使用第一个表达式。但是，当b 是整数且n = -1 时，numpy 会产生错误。将b 与1.0 相乘会将其更改为浮点数，numpy 不会给出此类错误。如果多次使用相同的参数a 和b，可能会生成冻结分布。

这是一个示例，它创建生成样本的直方图，并将其与pmf 进行比较。

import scipy.stats as st
import numpy as np
from matplotlib import pyplot as plt

class linear_fractional(st.rv_discrete):
    def _pmf(self, n, a, b):
        return (a * (1.0 * b) ** (n - 1)) / (a + b) ** (n + 1)

# parameters
a = 3 / 16
b = 1

LF = linear_fractional()

N = 10000
plt.hist(LF.rvs(a, b, size=N), bins=np.arange(-0.5, 50), ec='w', label='histogram of samples')
plt.plot(LF.pmf(np.arange(50), a, b) * N, 'ro', label='probability mass function (scaled)')
plt.legend(title=f'$a={a}; b={b}$')
plt.autoscale(enable=True, axis='x', tight=True)
plt.show()

LF.mean(a, b) 输出5.33333333333286

散点图是说明分布样本的另一种方法：

plt.scatter(np.random.uniform(0, 1, N), LF.rvs(a, b, size=N), marker=',', alpha=0.2, lw=0, s=1, color='crimson')

PS：当b=1 时，这个分布的公式等于geometric distribution 和p = a/(a+1) 减去1。这要快得多，因为它完全是在numpy 内部计算的。

samples = np.random.geometric(a/(a+1), size=1000) - 1