variables de generación de Poisson en C++

Question

I implementan esta función para generar una variable aleatoria de Poisson

typedef long unsigned int luint; 
luint poisson(luint lambda) { 
    double L = exp(-double(lambda)); 
    luint k = 0; 
    double p = 1; 
    do { 
     k++; 
     p *= mrand.rand(); 
    } while(p > L); 
    return (k-1); 
} 

donde mrand es el generador de números aleatorios Mersenne Twister. Encuentro que, al aumentar lambda, la distribución esperada va a ser incorrecta, con una media que se satura alrededor de 750. ¿Se debe a aproximaciones numéricas o cometí algún error?

Answer 1

De another question I asked earlier, parece que también podría aproximada poisson(750) como poisson(375) + poisson(375).

Answer 2

exp (-750) es un número muy pequeño, muy cerca del doble más pequeño posible, por lo que su problema es numérico. En cualquier caso, tu complejidad será lineal en lambda, por lo que el algoritmo no es muy eficiente para lambda alta. A menos que tenga una buena razón para codificar esto usted mismo, probablemente sea lógico usar una implementación de biblioteca existente, ya que estos algoritmos numéricos tienden a ser delicado precisamente para los problemas de precisión con los que se encuentra.

Answer 3

Dado que solo usa L en la expresión (p>L), básicamente está probando (log(p) > -lambda). Esa no es una transformación muy útil. Claro, ya no necesita exp (-750), pero simplemente desbordará p.

Ahora, p es sólo Π (mrand.rand()), y el registro (p) es el registro (Π (mrand.rand())) es Σ (log (mrand.rand()). Eso le da la necesaria transformación:

double logp = 0; 
do { 
    k++; 
    logp += log(mrand.rand()); 
} while(logp > -lambda); 

double tiene sólo 11 bits de exponente, pero unas 52 bits de la mantisa tanto, esto es un aumento masivo de la estabilidad numérica el precio pagado es que se necesita un log en cada iteración, en lugar de.. un solo exp por adelantado

Answer 4

En situaciones como estas, no necesita invocar el generador de números aleatorios más de una vez. u necesidad es una tabla de probabilidades acumuladas:

double c[k] = // the probability that X <= k (k = 0,...) 

luego generar un número aleatorio 0 <= r < 1, y tomar el primer número entero tal que Xc[X] > r. Puede encontrar este X con una búsqueda binaria.

Para generar esta tabla, necesitamos las probabilidades individuales

p[k] = lambda^k/(k! e^lambda) // // the probability that X = k 

Si lambda es grande, esto se convierte en muy impreciso, ya que ha encontrado. Pero podemos usar un truco aquí: comenzar en (o cerca) el valor más grande, con k = floor[lambda], y pretender por el momento que p[k] es igual a 1. Luego calcular p[i] para i > k usando la relación de recurrencia

p[i+1] = (p[i]*lambda)/(i+1) 

y para i < k usando

p[i-1] = (p[i]*i)/lambda 

Esto asegura que las mayores probabilidades tienen la mayor precisión posible.

Ahora acaba de calcular utilizando c[i]c[i+1] = c[i] + p[i+1], hasta el punto de que c[i+1] es lo mismo que c[i]. Entonces puede normalizar la matriz dividiendo por este valor límite c[i]; o puede dejar la matriz tal como está, y usar un número aleatorio 0 <= r < c[i].

Ver: http://en.wikipedia.org/wiki/Inverse_transform_sampling

Answer 5

Si va a la ruta "biblioteca existente", su compilador ya puede admitir el paquete C++ 11 std :: random. Así es como se lo utiliza:

#include <random> 
#include <ctime> 
#include <iostream> 

std::mt19937 mrand(std::time(0)); // seed however you want 

typedef long unsigned int luint; 

luint poisson(luint lambda) 
{ 
    std::poisson_distribution<luint> d(lambda); 
    return d(mrand); 
} 

int main() 
{ 
    std::cout << poisson(750) << '\n'; 
    std::poisson_distribution<luint> d(750); 
    std::cout << d(mrand) << '\n'; 
    std::cout << d(mrand) << '\n'; 
} 

Lo he utilizado dos formas anteriores:

1. he tratado de imitar la interfaz existente.

2. Si crea un std :: poisson_distribution con una media, es más eficiente usar esa distribución una y otra vez por la misma media (como se hizo en main()).

Aquí está un ejemplo de salida para mí:

751 
730 
779