Implementazione Box-Mueller generatore di numeri casuali in C #

Question

questa domanda: Random generatore di numeri che gravita numeri per un dato numero in gamma? ho fatto qualche ricerca da quando ho incontrato un tale generatore di numeri casuali. Tutto quello che ricordo è stato il nome di "Mueller", quindi credo che ho trovato, qui:

• Box-Mueller trasformare

riesco a trovare numerose implementazioni di esso in altre lingue, ma io non riesco a implementare correttamente in C #.

In questa pagina, per esempio, Il Metodo Box-Muller per i numeri casuali Generazione gaussiana dice che il codice dovrebbe essere così (questo non è 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;
}

Ora, ecco la mia attuazione di quanto sopra in C #. Si noti che la trasformazione produce 2 numeri, quindi il trucco con la "fase" di cui sopra. Ho semplicemente scartare il secondo valore e restituire il primo.

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;
}

La mia domanda è la seguente scenario specifico, in cui il mio codice non restituisce un valore nella gamma di 0-1, e non riesco a capire come il codice originale possibile.

• U = 0,5, v = 0,1
• S diventa 0.5*0.5 + 0.1*0.1 = 0.26
• fac diventa ~ 3.22
• il valore di ritorno è quindi ~ ~ 0.5 * 3.22 o 1.6

Non è all'interno 0 .. 1.

Che cosa sto facendo di sbagliato / non capire?

Se io modificare il mio codice in modo che invece di moltiplicare fac con u, ho moltiplicare per S, ottengo un valore che varia da 0 a 1, ma ha la distribuzione sbagliata (sembra avere una distribuzione massima intorno 0,7-0,8 e poi si assottiglia in entrambe le direzioni).

Answer 1

Il tuo codice va bene. Il tuo errore è pensare che deve restituire i valori esclusivamente all'interno [0, 1]. (Standard) distribuzione normale è una distribuzione di peso diverso da zero sull'intera retta reale. Vale a dire, i valori al di fuori di [0, 1] sono possibili. Infatti, i valori all'interno [-1, 0] sono altrettanto probabile come valori all'interno [0, 1], e inoltre, il complemento di [0, 1] ha circa 66% del peso della distribuzione normale. Pertanto, il 66% del tempo ci aspettiamo un valore al di fuori di [0, 1].

Inoltre, credo che questo non è il Box-Mueller trasformare, ma è in realtà il metodo polare Marsaglia.

Answer 2

Non sono un matematico o statistico, ma se penso a questo non mi aspetto una distribuzione gaussiana per tornare i numeri in un intervallo esatto. Data la sua attuazione la media è 0 e la deviazione standard è 1 così che ci si aspetterebbe valori distribuiti sulla curva a campana con 0 al centro e quindi riducendo i numeri da 0 discostano su entrambi i lati. Quindi la sequenza sarebbe sicuramente coprire entrambe +/- numeri.

Quindi dal momento che è statistica, perché sarebbe difficile limitata a -1..1 solo perché lo std.dev è 1? Non ci può essere statisticamente certo gioco su entrambi i lati ed ancora soddisfare il requisito statistico.

Answer 3

Il variata casuale uniforme è davvero all'interno 0..1, ma la variata casuale gaussiana (che è ciò che l'algoritmo Box-Muller genera) può essere in qualsiasi punto della retta reale. Vedere wiki / NormalDistribution per i dettagli.

Answer 4

Credo che la funzione restituisce le coordinate polari. Quindi è necessario entrambi i valori per ottenere risultati corretti.

Inoltre, distribuzione gaussiana non è tra 0 .. 1. Si può facilmente finire come 1000, ma la probabilità di tale evento è estremamente bassa.

Answer 5

Questo è un metodo Monte Carlo così non si può bloccare il risultato, ma che cosa si può fare è ignorare campioni.

// 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;
    }
}