python - tabla - sumatoria en qgis

Muestreo multinomial a partir de pequeños vectores de probabilidad de registro en números/caracteres

Question

¿Existe una función en numpy / scipy que le permite muestrear multinomial de un vector de probabilidades de registro pequeño, sin perder precisión? ejemplo:

# sample element randomly from these log probabilities
l = [-900, -1680]

el método ingenuo falla debido al desbordamiento:

import scipy
import numpy as np
# this makes a all zeroes
a = np.exp(l) / scipy.misc.logsumexp(l)
r = np.random.multinomial(1, a)

este es un intento

def s(l):
    m = np.max(l)
    norm = m + np.log(np.sum(np.exp(l - m)))
    p = np.exp(l - norm)
    return np.where(np.random.multinomial(1, p) == 1)[0][0]

¿Es este el método mejor / más rápido y puede np.exp() en el último paso?

Answer 1

En primer lugar, creo que el problema que está encontrando es porque está normalizando sus probabilidades de manera incorrecta. Esta línea es incorrecta:

a = np.exp(l) / scipy.misc.logsumexp(l)

Estás dividiendo una probabilidad por una probabilidad de registro, lo que no tiene sentido. En su lugar, probablemente quieras

a = np.exp(l - scipy.misc.logsumexp(l))

Si lo hace, encontrará a = [1, 0] y su muestreador multinomial funciona como se espera hasta la precisión de punto flotante en la segunda probabilidad.

Una solución para pequeñas N: Histogramas

Dicho esto, si aún necesita más precisión y el rendimiento no es tan preocupante, una forma en que podría progresar es implementando un muestreador multinomial desde cero, y luego modificándolo para que funcione con una mayor precisión.

La función multinomial de NumPy se implementa en Cython , y esencialmente realiza un bucle sobre varias muestras binomiales y las combina en una muestra multinomial. y puedes llamarlo así:

np.random.multinomial(10, [0.1, 0.2, 0.7])
# [0, 1, 9]

(Tenga en cuenta que los valores de salida precisos aquí y abajo son aleatorios, y cambiarán de una llamada a otra).

Otra forma de implementar un muestreador multinomial es generar N valores aleatorios uniformes y luego calcular el histograma con intervalos definidos por las probabilidades acumuladas:

def multinomial(N, p):
    rand = np.random.uniform(size=N)
    p_cuml = np.cumsum(np.hstack([[0], p]))
    p_cuml /= p_cuml[-1]
    return np.histogram(rand, bins=p_cuml)[0]

multinomial(10, [0.1, 0.2, 0.7])
# [1, 1, 8]

Con este método en mente, podemos pensar en hacer las cosas con mayor precisión manteniendo todo en el espacio de registro. El truco principal es darse cuenta de que el registro de desviaciones aleatorias uniformes es equivalente al negativo de desviaciones aleatorias exponenciales, por lo que puede hacer todo lo anterior sin dejar espacio de registro:

def multinomial_log(N, logp):
    log_rand = -np.random.exponential(size=N)
    logp_cuml = np.logaddexp.accumulate(np.hstack([[-np.inf], logp]))
    logp_cuml -= logp_cuml[-1]
    return np.histogram(log_rand, bins=logp_cuml)[0]

multinomial_log(10, np.log([0.1, 0.2, 0.7]))
# [1, 2, 7]

Los dibujos multinomiales resultantes mantendrán la precisión incluso para valores muy pequeños en la matriz p . Desafortunadamente, estas soluciones basadas en histogramas serán mucho más lentas que la función numpy.multinomial nativa, por lo que si el rendimiento es un problema, es posible que necesite otro enfoque. Una opción sería adaptar el código Cython vinculado arriba para trabajar en el espacio de registro, usando trucos matemáticos similares a los que usé aquí.

Una solución para N grande: aproximación de Poisson

El problema con la solución anterior es que a medida que N crece, se vuelve muy lento. Estaba pensando en esto y me di cuenta de que hay una manera más eficiente de avanzar, a pesar de np.random.multinomial falla para probabilidades más pequeñas que 1E-16 o algo así.

Aquí hay un ejemplo de ese error: en una máquina de 64 bits, esto siempre dará cero para la primera entrada debido a la forma en que se implementa el código, cuando en realidad debería dar algo cerca de 10:

np.random.multinomial(1E18, [1E-17, 1])
# array([                  0, 1000000000000000000])

Si profundiza en la fuente, puede rastrear este problema hasta la función binomial sobre la que se construye la función multinomial. El código cython internamente hace algo como esto:

def multinomial_basic(N, p, size=None):
    results = np.array([np.random.binomial(N, pi, size) for pi in p])
    results[-1] = int(N) - results[:-1].sum(0)
    return np.rollaxis(results, 0, results.ndim)

multinomial_basic(1E18, [1E-17, 1])
# array([                  0, 1000000000000000000])

El problema es que la función binomial ahoga en valores muy pequeños de p ; esto se debe a que el algoritmo calcula el valor (1 - p) , por lo que el valor de p está limitado por la precisión del punto flotante.

Entonces, ¿qué podemos hacer? Bueno, resulta que para valores pequeños de p, la distribución de Poisson es una aproximación extremadamente buena de la distribución binomial, y la implementación no tiene estos problemas. Así que podemos construir una función multinomial robusta basada en un muestreador binomial robusto que cambia a un muestreador Poisson en la pequeña p:

def binomial_robust(N, p, size=None):
    if p < 1E-7:
        return np.random.poisson(N * p, size)
    else:
        return np.random.binomial(N, p, size)

def multinomial_robust(N, p, size=None):
    results = np.array([binomial_robust(N, pi, size) for pi in p])
    results[-1] = int(N) - results[:-1].sum(0)
    return np.rollaxis(results, 0, results.ndim)

multinomial_robust(1E18, [1E-17, 1])
array([                 12, 999999999999999988])

La primera entrada es distinta de cero y cerca de 10 como se esperaba. Tenga en cuenta que no podemos usar N mayor que 1E18 , ya que desbordará el entero largo. Pero podemos confirmar que nuestro enfoque funciona para probabilidades más pequeñas utilizando el parámetro de size y promediando los resultados:

p = [1E-23, 1E-22, 1E-21, 1E-20, 1]
size = int(1E6)
multinomial_robust(1E18, p, size).mean(0)
# array([  1.70000000e-05,   9.00000000e-05,   9.76000000e-04,
#          1.00620000e-02,   1.00000000e+18])

Vemos que incluso para estas muy pequeñas probabilidades, los valores multinomiales están apareciendo en la proporción correcta. El resultado es una aproximación muy robusta y muy rápida a la distribución multinomial para p pequeñas.