9

我试图将泊松连续误差条放在我用 matplotlib 制作的直方图上,但我似乎找不到一个 numpy 函数,假设泊松数据会给我一个 95% 的置信区间。理想情况下,该解决方案不依赖于 scipy,但任何事情都会奏效。有这样的功能吗?我发现了很多关于引导的信息,但在我的情况下这似乎有点过分了。

4

3 回答 3

11

我最终根据在 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

这将返回连续(而不是离散)边界,这在我的领域中更为标准。

于 2013-02-12T12:25:45.487 回答
8

使用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])
于 2013-02-11T15:00:31.573 回答
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 即可得到这些。

请查阅论文以了解这些公式对于您需要的置信区间的准确性,但对于大多数实际应用来说,这些公式应该足够准确。

于 2017-04-05T15:34:17.700 回答