importar - scipy.stats python

Intervalo de confianza de Poisson con numpy

Question

Estoy tratando de poner barras de error continuas de Poisson en un histograma que estoy haciendo con matplotlib, pero parece que no puedo encontrar una función numpy que me dé un intervalo de confianza del 95% asumiendo datos poissonianos. Idealmente, la solución no depende de nada, pero cualquier cosa funcionará. ¿Existe tal función? He encontrado mucho sobre bootstrapping pero esto parece un poco excesivo en mi caso.

Answer 1

Terminé escribiendo mi propia función basada en algunas propiedades que encontré en 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

Esto devuelve límites continuos (en lugar de discretos), que es más estándar en mi campo.

Answer 2

Usando scipy.stats.poisson , y el método de interval :

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

Aunque solo tiene sentido limitado calcular la distribución de Poisson para valores no enteros, los intervalos de confianza exactos solicitados por el OP se pueden calcular de la siguiente manera:

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

También es posible usar el método ppf :

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

Pero debido a que la distribución es discreta, los valores de retorno serán enteros, y el intervalo de confianza no abarcará exactamente el 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 3

Este problema surge mucho en astronomía (¡mi campo!) Y este documento es la referencia a seguir para estos intervalos de confianza: Gehrels 1980

Tiene un montón de matemática para un intervalo de confianza arbitrario con las estadísticas de Poisson, pero para un intervalo de confianza del 95% bilateral (correspondiente a un intervalo de confianza gaussiano de 2 sigma, o S = 2 en el contexto de este documento) algunos las fórmulas analíticas simples para los límites de confianza superior e inferior para cuando se miden N eventos son

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

donde ya los puse en formato Python. Todo lo que necesitas es numpy o tu otro módulo de raíz cuadrada favorito. Tenga en cuenta que estos le darán los límites superior e inferior para los eventos, no los valores +/-. Simplemente resta N de ambos para obtenerlos.

Consulte en el documento la precisión de estas fórmulas para el intervalo de confianza que necesita, pero estas deberían ser más que suficientemente precisas para la mayoría de las aplicaciones prácticas.