使用numpy计算泊松置信区间

Question

我正在用matplotlib制作一个直方图，想要在上面加上泊松分布的连续误差条，但我找不到一个numpy的函数可以给我95%的置信区间，前提是数据是泊松分布的。理想情况下，我希望这个解决方案不依赖于scipy，但其他的也可以。有没有这样的函数呢？我看到很多关于自助法（bootstrapping）的资料，但在我的情况下，这似乎有点过于复杂了。

Answer 1

这个问题在天文学中经常出现（这是我的专业！），而这篇论文是关于置信区间的权威参考资料：Gehrels 1980

里面有很多数学公式，主要是针对任意置信区间和泊松统计的。不过，如果你只需要一个双侧的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，就能得到±值。

请查看论文以确认这些公式在你需要的置信区间下的准确性，但对于大多数实际应用来说，这些公式应该足够准确了。

Answer 2

我最后根据我在维基百科上找到的一些特性写了自己的函数。

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

这个函数返回的是连续的范围，而不是离散的，这在我的领域里更常见。

Answer 3

使用 scipy.stats.poisson 这个库，还有 interval 方法：

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

虽然对于非整数值计算泊松分布的意义有限，但可以按照以下方式计算出提问者所要求的精确置信区间：

>>> 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])