Implementando el generador de números aleatorios Box-Mueller en C#

Question

De this question: Random number generator which gravitates numbers to any given number in range? Investigué un poco ya que me encontré con un generador de números aleatorios. Todo lo que recuerdo es el nombre "Mueller", así que supongo que lo encontré, aquí:

• Box-Mueller transform

puedo encontrar numerosas implementaciones de la misma en otros idiomas, pero me parece que no puede implementarlo correctamente en C#.

Esta página, por ejemplo, The Box-Muller Method for Generating Gaussian Random Numbers dice que el código debería tener este aspecto (esto no es C#):

#include <stdio.h> 
#include <stdlib.h> 
#include <math.h> 
#include <time.h> 

double 
gaussian(void) 
{ 
    static double v, fac; 
    static int phase = 0; 
    double S, Z, U1, U2, u; 

    if (phase) 
     Z = v * fac; 
    else 
    { 
     do 
     { 
     U1 = (double)rand()/RAND_MAX; 
     U2 = (double)rand()/RAND_MAX; 

     u = 2. * U1 - 1.; 
     v = 2. * U2 - 1.; 
     S = u * u + v * v; 
     } while(S >= 1); 

     fac = sqrt (-2. * log(S)/S); 
     Z = u * fac; 
    } 

    phase = 1 - phase; 

    return Z; 
} 

Ahora, aquí está mi aplicación de lo anterior en C#. Tenga en cuenta que la transformación produce 2 números, de ahí el truco con la "fase" anterior. Simplemente descarto el segundo valor y devuelvo el primero.

public static double NextGaussianDouble(this Random r) 
{ 
    double u, v, S; 

    do 
    { 
     u = 2.0 * r.NextDouble() - 1.0; 
     v = 2.0 * r.NextDouble() - 1.0; 
     S = u * u + v * v; 
    } 
    while (S >= 1.0); 

    double fac = Math.Sqrt(-2.0 * Math.Log(S)/S); 
    return u * fac; 
} 

Mi pregunta es la siguiente escenario específico, donde mi código no devuelve un valor en el rango de 0-1, y no puedo entender cómo el código original sea posible.

• u = 0,5, v = 0,1
• S se convierte en 0.5*0.5 + 0.1*0.1 = 0.26
• FAC se convierte en ~ 3.22
• el valor de retorno es, pues, ~ 0.5 * 3.22 o ~ 1.6

Eso no está dentro de 0 .. 1.

¿Qué estoy haciendo mal/no estoy entendiendo?

Si modifico mi código para que en vez de multiplicar fac con u, multiplico por S, puedo obtener un valor que oscila entre 0 y 1, pero tiene la mala distribución (parece tener una distribución máxima alrededor de 0,7- 0.8 y luego se estrecha en ambas direcciones.)

Answer 1

Su código es correcto. Su error es pensar que debe devolver valores exclusivamente dentro de [0, 1]. La distribución normal (estándar) es una distribución con un peso distinto de cero en toda la línea real. Es decir, son posibles valores fuera de [0, 1]. De hecho, los valores dentro de [-1, 0] son tan probables como los valores dentro de [0, 1], y además, el complemento de [0, 1] tiene aproximadamente el 66% del peso de la distribución normal. Por lo tanto, el 66% de las veces esperamos un valor fuera de [0, 1].

Además, creo que esta no es la transformación de Box-Mueller, pero en realidad es el método polar de Marsaglia.

Answer 2

Creo que la función devuelve coordenadas polares. Entonces necesita ambos valores para obtener resultados correctos.

Además, la distribución gaussiana no está entre 0 .. 1. Puede terminar fácilmente como 1000, pero la probabilidad de que ocurra es extremadamente baja.

Answer 3

No soy matemático o estadístico, pero si lo pienso no esperaría que una distribución gaussiana devuelva números en un rango exacto. Dada su implementación, la media es 0 y la desviación estándar es 1, por lo que esperaría que los valores se distribuyeran en la curva de la campana con 0 en el centro y luego reduciendo a medida que los números se desvían de 0 en cada lado. Entonces la secuencia definitivamente cubriría ambos números +/-.

Luego, dado que es estadístico, ¿por qué sería difícil limitarlo a -1..1 solo porque std.dev es 1? No puede haber algo de juego estadísticamente en ninguno de los lados y aún así cumplir con los requisitos estadísticos.

Answer 4

La variante aleatoria uniforme está efectivamente dentro de 0..1, pero la variable aleatoria gaussiana (que es lo que genera el algoritmo Box-Muller) puede estar en cualquier lugar de la línea real.Ver wiki/NormalDistribution para más detalles.

Answer 5

Este es un método de monte carlo, por lo que no puede fijar el resultado, pero lo que puede hacer es ignorar las muestras.

// return random value in the range [0,1]. 
double gaussian_random() 
{ 
    double sigma = 1.0/8.0; // or whatever works. 
    while (1) { 
     double z = gaussian() * sigma + 0.5; 
     if (z >= 0.0 && z <= 1.0) 
      return z; 
    } 
}