使用numpy计算泊松分布的置信区间

Question

我正在使用matplotlib制作直方图，并尝试在上面放置泊松连续误差条，但是我似乎找不到一个numpy函数，可以给我一个95％的置信区间，假设是泊松数据。理想情况下，解决方案不依赖于scipy，但任何东西都可以。是否存在这样的功能？我发现很多关于引导的内容，但在我的情况下似乎有点过度。

- Shep

Answer 1

我最终根据维基百科上的一些特性编写了自己的函数。

def poisson_interval(k, alpha=0.05): 
    """
    uses chisquared info to get the poisson interval. Uses scipy.stats 
    (imports in function). 
    """
    from scipy.stats import chi2
    a = alpha
    low, high = (chi2.ppf(a/2, 2*k) / 2, chi2.ppf(1-a/2, 2*k + 2) / 2)
    if k == 0: 
        low = 0.0
    return low, high

这将返回连续（而不是离散）的边界，这在我的领域中更为常见。

- Shep

Answer 2

使用scipy.stats.poisson和interval方法：

>>> scipy.stats.poisson.interval(0.95, [10, 20, 30])
(array([  4.,  12.,  20.]), array([ 17.,  29.,  41.]))

虽然对于非整数值计算泊松分布只有有限的意义，但是可以按照以下方式计算OP请求的精确置信区间：

>>> data = np.array([10, 20, 30])
>>> scipy.stats.poisson.interval(0.95, data)
(array([  4.,  12.,  20.]), array([ 17.,  29.,  41.]))
>>> np.array(scipy.stats.chi2.interval(.95, 2 * data)) / 2 - 1
array([[  3.7953887 ,  11.21651959,  19.24087402],
       [ 16.08480345,  28.67085357,  40.64883744]])

---

同样可以使用ppf方法：

>>> data = np.array([10, 20, 30])
>>> scipy.stats.poisson.ppf([0.025, 0.975], data[:, None])
array([[  4.,  17.],
       [ 12.,  29.],
       [ 20.,  41.]])

但是由于分布是离散的，返回值将为整数，置信区间将不会完全跨越95％：

>>> scipy.stats.poisson.ppf([0.025, 0.975], 10)
array([  4.,  17.])
>>> scipy.stats.poisson.cdf([4, 17], 10)
array([ 0.02925269,  0.98572239])

- Jaime

---

你知道如何获取精确的返回值吗？ - Shep

@Shep，我在我的回答中添加了一个基于卡方检验的版本，但是使用了“interval”。 - Jaime

@Jaime 精确公式是不正确的。例如，对于较小的 k 值，它返回负值。它与 Shep 的答案不等价，其中 chi2.ppf 是使用 $2k$ 和 $2k+2$ 计算的。 - Anton Mellit

Answer 3

这个问题在天文学（我的领域！）中经常出现，而这篇论文是这些置信区间的首选参考：Gehrels 1980 对于带有Poisson统计的任意置信区间，这篇论文中提供了大量数学内容，但对于双侧95%置信区间（相当于高斯置信区间的2-sigma，或在本文的上下文中为S = 2），当测量N个事件时，一些简单的解析公式可用于计算上限和下限的置信度限制。

upper = N + 2. * np.sqrt(N + 1) + 4. / 3.
lower = N * (1. - 1. / (9. * N) - 2. / (3. * np.sqrt(N))) ** 3.

我已经将它们以Python格式为您准备好。您只需要使用numpy或其他喜欢的平方根模块即可。请记住，这些将为您提供事件的上限和下限 - 而不是+/-值。您只需从两者中减去N即可获得那些值。

请参考论文以获取所需置信区间的公式精度，但对于大多数实际应用而言，这些应该已经足够准确了。

- Steven Ehlert

---

感谢 @firelynx 的编辑。这样更易读了。由于我更多地从事科学而非软件工程，所以我经常忘记遵循 PEP8。 - Steven Ehlert