C#: algoritmo numerico per generare numeri dalla distribuzione binomiale

Question

Devo generare numeri casuali dalla distribuzione binomiale (n, p).

Una variabile casuale binomiale (n, p) è la somma di n variabili uniformi che prendono 1 con probabilità p. In pseudo codice, x=0; for(i=0; i<n; ++i) x+=(rand()<p?1:0); genererà un binomio (n, p).

Ho bisogno di generare questo per n piccoli e molto grandi, per esempio n = 10^6 e p = 0.02. C'è un algoritmo numerico veloce per generarlo?

EDIT -

In questo momento questo è quello che ho come approssimazione (con funzioni per la precisa distribuzione di Poisson e Normale) -

public long Binomial(long n, double p) { 
     // As of now it is an approximation 
     if (n < 1000) { 
      long result = 0; 
      for (int i=0; i<n; ++i) 
       if (random.NextDouble() < p) result++; 
      return result; 
     } 
     if (n * p < 10) return Poisson(n * p); 
     else if (n * (1 - p) < 10) return n - Poisson(n * p); 
     else { 
      long v = (long)(0.5 + nextNormal(n * p, Math.Sqrt(n * p * (1 - p)))); 
      if (v < 0) v = 0; 
      else if (v > n) v = n; 
      return v; 
     } 
    } 

Answer 1

Se si è disposti a pagare, dare un'occhiata a NMath di Centerspace.

In caso contrario, il codice C utilizzato dal programma Stats R è here e dovrebbe essere semplice da portare a C#.

MODIFICA: Vi sono dettagli (codice inc.) Sulla creazione di un metodo per questo su p178 di Practical Numerical Methods with C# di Jack Xu.

UN ALTRO MODIFICO: A free C# library che fa quello che vuoi.

Answer 2

Non c'è modo ovvio per farlo in modo efficiente Per la piccola n, potresti anche solo noi la formula per calcolare il PDF inverso. Per n maggiore, probabilmente stai meglio usando uno degli approximations to other distributions che sono più facili da calcolare.

Answer 3

Un'altra opzione sarebbe campionare da Normal o Poisson come si fa e quindi aggiungere un passaggio Metropolis-Hastings per accettare o rifiutare il campione. Se accetti di aver finito, se rifiuti, devi ricampionare completamente. La mia ipotesi è che, poiché l'approssimazione è così vicina, si otterrà quasi sempre un passo di accettazione, di tanto in tanto si potrebbe rifiutare.

Anche Luc Devroye's book ha alcuni grandi algoritmi per il campionamento binomiale.

PS Se si finisce con un buon algoritmo; ti dispiacerebbe condividerlo a Math.Net Numerics?