Log-computaciones en Python

Estoy buscando para calcular algo como:

fórmula

Donde f(i) es una función que devuelve un número real en [-1,1] para cualquier i en {1,2,...,5000} .

Obviamente, el resultado de la sum está en algún lugar en [-1,1] , pero cuando parece que no puedo calcularlo en Python usando encoding directa, como 0.5 5000 convierte en 0 y comb(5000,2000) convierte en inf , que da como resultado que la sum calculada se convierta en NaN .

La solución requerida es usar el registro en ambos lados.

Eso es usar la identidad a × b = 2 log(a) + log(b) , si pudiera calcular log(a) y log(b) podría calcular la sum, incluso si a es grande y b es casi 0 .

Así que supongo que lo que pregunto es si hay una forma fácil de computar

 log2(scipy.misc.comb(5000,2000)) 

Así podría calcular mi sum simplemente por

 sum([2**(log2comb(5000,i)-5000) * f(i) for i in range(1,5000) ]) 

La solución de @abarnert, mientras trabajaba para la cifra 5000, soluciona el problema aumentando la precisión con la que se calcula el peine. Esto funciona para este ejemplo, pero no se escala, ya que la memoria requerida boostía significativamente si en lugar de 5000 tuviéramos 1e7 por ejemplo.

Actualmente, estoy usando una solución que es fea, pero que mantiene bajo el consumo de memoria:

 log2(comb(5000,2000)) = sum([log2 (x) for x in 1:5000])-sum([log2 (x) for x in 1:2000])-sum([log2 (x) for x in 1:3000]) 

¿Hay una manera de hacerlo en una expresión legible?

La sum

fórmula

es la expectativa de f con respecto a una distribución binomial con n = 5000 p = 0.5 .

Puede calcular esto con scipy.stats.binom.expect :

 import scipy.stats as stats def f(i): return i n, p = 5000, 0.5 print(stats.binom.expect(f, (n, p), lb=0, ub=n)) # 2499.99999997 

También tenga en cuenta que a medida que n va hacia el infinito, con p fijo, la distribución binomial se aproxima a la distribución normal con media np y varianza np*(1-p) . Por lo tanto, para grandes n puede en cambio calcular:

 import math print(stats.norm.expect(f, loc=n*p, scale=math.sqrt((n*p*(1-p))), lb=0, ub=n)) # 2500.0 

Por defecto, comb le da un float64 , que se desborda y le da inf .

Pero si pasas exact=True , te da un int Python de tamaño variable, que no puede desbordarse (a menos que te vuelvas tan grande que te quedes sin memoria).

Y, aunque no puede usar np.log2 en un int , puede usar math.log2 de Python.

Asi que:

 math.log2(scipy.misc.comb(5000, 2000, exact=True)) 

Como alternativa, tu pariente de que n elige k se define como n!k / k! , ¿derecho? Puede reducir eso a ∏(i=1...k)((n+1-i)/i) , que es fácil de calcular.

O, si desea evitar el desbordamiento, puede hacerlo alternando * (ni) y / (ki) .

Lo que, por supuesto, también puede reducir al agregar y restar registros. Creo que el bucle en Python y la computación de logaritmos 4000 va a ser más lento que el bucle en C y la computación de 4000 multiplicaciones, pero siempre podemos vectorizarlo, y luego, podría ser más rápido. Vamos a escribirlo y probarlo:

 In [1327]: n, k = 5000, 2000 In [1328]: %timeit math.log2(scipy.misc.comb(5000, 2000, exact=True)) 100 loops, best of 3: 1.6 ms per loop In [1329]: %timeit np.log2(np.arange(n-k+1, n+1)).sum() - np.log2(np.arange(1, k+1)).sum() 10000 loops, best of 3: 91.1 µs per loop 

Por supuesto, si estás más preocupado por la memoria en lugar del tiempo … bueno, esto obviamente lo empeora. Tenemos 2000 flotadores de 8 bytes en lugar de un entero de 608 bytes a la vez. Y si sube a 100000, 20000, obtiene 20000 flotantes de 8 bytes en lugar de un entero de 9K. Y a 1000000, 200000, son 200000 flotantes de 8 bytes frente a un entero de 720K.

No estoy seguro de por qué de cualquier manera es un problema para ti. Especialmente teniendo en cuenta que está utilizando una lista comp en lugar de un genexpr, y por lo tanto, está creando una lista innecesaria de 5000, 100000 o 1000000 flotadores de Python: 24 MB no es un problema, pero ¿720 K lo es? Pero si lo es, obviamente podemos hacer lo mismo de forma iterativa, a costa de cierta velocidad:

 r = sum(math.log2(ni) - math.log2(ki) for i in range(nk)) 

Esto no es mucho más lento que la solución scipy , y nunca usa más que un pequeño número constante de bytes (un puñado de flotadores de Python). (A menos que estés en Python 2, en cuyo caso … solo usa xrange lugar de range y vuelve a ser constante).


Como nota al margen, ¿por qué utiliza una comprensión de lista en lugar de una matriz NumPy con operaciones vectorizadas (para velocidad y también un poco de compacidad) o una expresión de generador en lugar de una comprensión de lista (para ningún uso de memoria, en absoluto) costo a la velocidad)?

EDIT: @unutbu ha respondido la pregunta real , pero dejaré esto aquí en caso de que log2comb(n, k) sea ​​útil para cualquiera.


comb(n, k) es n! / ((nk)! k!), y n! se puede calcular utilizando la función gamma(n+1) . Scipy proporciona la función scipy.special.gamma . Scipy también proporciona gammaln , que es el registro (log natural, es decir) de la función Gamma.

Por lo tanto, log(comb(n, k)) se puede calcular como gammaln(n+1) - gammaln(n-k+1) - gammaln(k+1)

Por ejemplo, log (comb (100, 8)) (después de ejecutar from scipy.special import gammaln ):

 In [26]: log(comb(100, 8)) Out[26]: 25.949484949043022 In [27]: gammaln(101) - gammaln(93) - gammaln(9) Out[27]: 25.949484949042962 

y registro (peine (5000, 2000)):

 In [28]: log(comb(5000, 2000)) # Overflow! Out[28]: inf In [29]: gammaln(5001) - gammaln(3001) - gammaln(2001) Out[29]: 3360.5943053174142 

(Por supuesto, para obtener el logaritmo en base 2, simplemente divídalo por log(2) .

Por conveniencia, puede definir:

 from math import log from scipy.special import gammaln def log2comb(n, k): return (gammaln(n+1) - gammaln(n-k+1) - gammaln(k+1)) / log(2)