具有numpy的泊松置信区间

Question

我想将泊松连续误差条放在我用matplotlib制作的直方图上，但我似乎无法找到一个numpy函数，它会给出95％的置信区间，假设泊松数据。理想情况下，解决方案不依赖于scipy，但任何东西都可以工作。这样的功能是否存在？我已经发现了很多关于bootstrapping的内容，但在我看来这似乎有点过分。

Answer 1

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

Answer 2

我最后写基于我自己的函数some properties I found on Wikipedia。

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

此问题出现在天文学很多（我的领域！），这纸是去到这些置信区间参考：Gehrels 1980

它有很多数学在它的任意置信区间与泊松统计量，但对于双侧95％置信区间（对应于2西格玛高斯置信区间，或本文中上下文中的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格式中为你alr伊迪。所有你需要的是numpy或你最喜欢的平方根模块。请记住，这些会给你事件的上限和下限 - 而不是+/-值。你只需从这两个中减去N就可以得到这些。

有关这些公式的准确性，请查阅您需要的置信区间，但这些公式对于大多数实际应用而言应该足够精确。