Attuazione Box-Mueller generatore di numeri casuali in C#

Da questa domanda: generatore di numeri Casuali che gravita numeri per un dato numero di serie? ho fatto qualche ricerca dal momento che ho incontrato ad un generatore di numero casuale prima. Tutto quello che ricordo era il nome di “Mueller”, quindi credo che ho trovato qui:

Posso trovare numerose implementazioni in altre lingue, ma non riesco ad implementare correttamente in C#.

Questa pagina, per esempio, Il Box-Muller Metodo per la Generazione di Numeri Casuali Gaussiane dice che il codice dovrebbe essere simile a questo (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 il mio 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 ignorare il secondo valore e il ritorno 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, dove il mio codice non restituisce un valore nell’intervallo 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 è così ~0.5 * 3.22 o ~1.6

Che non entro 0 .. 1.

Che cosa sto facendo di sbagliato/non capire?

Se posso modificare il mio codice in modo che invece di moltiplicare fac con u, ho moltiplicare per S, ho un valore che va da 0 a 1, ma ha il torto di distribuzione (sembra di avere un massimo di distribuzione di circa 0,7-0,8 e poi si assottiglia in entrambe le direzioni).

  • Nota che ho controllato un paio di esempi di codice di cui sopra, di solito in C o Java, e sembrano tutti più o meno lo stesso.
  • Sei sicuro che il codice C genera esattamente quello che vuoi?



5 Replies
  1. 9

    Il tuo codice va bene. Il tuo errore è pensare che si deve restituire valori esclusivamente all’interno [0, 1]. (Standard) distribuzione normale è una distribuzione con peso diverso da zero su tutta la retta reale. Che è, valori al di fuori di [0, 1] sono possibili. Infatti, i valori all’interno di [-1, 0] sono altrettanto probabile che, come valori all’interno di [0, 1], e inoltre, il complemento di [0, 1] ha circa il 66% del peso della distribuzione normale. Pertanto, il 66% del tempo ci si aspetta un valore al di fuori del [0, 1].

    Inoltre, penso che questo non è il Box-Mueller di trasformazione, ma in realtà è la Marsaglia polar metodo.

  2. 2

    Io non sono un matematico o statistico, ma se penso che su questo non mi aspetterei una distribuzione Gaussiana per tornare numeri in un intervallo esatto. Data di implementazione la media è 0 e la deviazione standard è 1 quindi mi aspetto che i valori distribuiti su una curva a campana con 0 al centro e quindi riducendo i numeri deviare dal 0 su entrambi i lati. Così la sequenza sarebbe sicuramente coprire sia +/- numeri.

    Quindi dato che è di tipo statistico, perché sarebbe difficile limitato a -1..1 solo perché la std.dev è 1? Non ci può essere statisticamente un po ‘ di gioco su entrambi i lati e ancora soddisfare le statistiche requisito.

  3. 2

    Uniforme casuale variante è infatti all’interno di 0..1, ma la gaussiana casuale variante (che è quello che di Box-Muller algoritmo genera) può essere qualsiasi punto della retta reale. Vedere wiki/NormalDistribution per i dettagli.

  4. 1

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

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

  5. 0

    Questo è un metodo monte carlo, così non è possibile bloccare il risultato, ma quello che si può fare è ignorare i 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;
        }
    }

Lascia un commento