Estadísticas SciPy BINOM vuelve cdf nan

votos
2

Si he entendido bien, el cdfde una scipy.statsdistribución discreta debe devolver la suma de las probabilidades de los valores hasta el parámetro dado.

Por lo tanto, scipy.stats.binom(7000000000, 0.5).cdf(6999999999)debe devolver algo casi exactamente 1, ya que en 7 mil millones de ensayos, con una posibilidad de 50/50, la probabilidad de obtener un éxito en 7 mil millones menos 1 de ellos o menos es más o menos cierto. En su lugar, consigo np.nan. De hecho, para cualquier valor proporcionado a .cdfSALVO 7 mil millones en sí mismo (o más), le regreso np.nan.

¿Que esta pasando aqui? ¿Hay algún límite a los números que scipy.statslas distribuciones pueden manejar que no está en los documentos?

Publicado el 07/11/2018 a las 22:48
fuente por usuario
En otros idiomas...                            


1 respuestas

votos
2

TL; DR

La falta de precisión de coma flotante durante los cálculos internos. Aunque scipy es una biblioteca de Python, su núcleo está escrito en C y C utiliza tipos numéricos.


Dejame mostrarte un ejemplo:

import scipy.stats

for i in range (13):
    trials = 10 ** i
    print(f"i: {i}\tprobability: {scipy.stats.binom(trials, 0.5).cdf(trials - 1)}")

Y la salida es:

i: 0    probability: 0.5
i: 1    probability: 0.9990234375
i: 2    probability: 0.9999999999999999
i: 3    probability: 0.9999999999999999
i: 4    probability: 0.9999999999999999
i: 5    probability: 0.9999999999999999
i: 6    probability: 0.9999999999999999
i: 7    probability: 0.9999999999999999
i: 8    probability: 0.9999999999999999
i: 9    probability: 0.9999999999999999
i: 10   probability: nan
i: 11   probability: nan
i: 12   probability: nan

La razón radica en la fórmula para la distribución binomial CDF (no puedo insertar una imagen asi que aquí hay un enlace a wiki: https://en.wikipedia.org/wiki/Binomial_distribution

Dentro de las fuentes de SciPy veríamos un refence a esta aplicación: http://www.netlib.org/cephes/doubldoc.html#bdtr

Profundamente dentro de ella implica la división por trials( incbet.c, line 375: ai = 1.0 / a;aquí se llama a, pero nwm). Y si su trialses demasiado grande, el resultado de esta división es tan pequeña que cuando añadimos este pequeño número a otro, no tan poco número, no cambia realmente porque nos falta precisión de coma flotante aquí (sólo 64 bits de hasta ahora) . Entonces, después de algunos más aritmética, tratamos de obtener el logaritmo de un número, pero es igual a cero, ya que no tiene cambios cuando debe. Y log(0)no está definido, lo que equivale a np.nan.

Respondida el 08/11/2018 a las 00:02
fuente por usuario

Cookies help us deliver our services. By using our services, you agree to our use of cookies. Learn more