Generar números aleatorios siguiendo una distribución normal en C / C++

EDIT: Desde el 12 de agosto de 2011 tenemos C++11 que ofrece directamente std::normal_distribution, que es lo que haría hoy.

Aquí está la respuesta original:

Aquí hay algunas soluciones ordenadas por complejidad ascendente.

1. Añadir 12 números uniformes del 0 al 1 y restar 6. Esto coincidirá con la media y la desviación estándar de una variable normal. Un inconveniente obvio es que el rango está limitado a + / -6, a diferencia de una normal verdadera distribución.

2. Box-Muller transform - se listó anteriormente, y es relativamente fácil de implementar. Sin embargo, si necesita muestras muy precisas, tenga en cuenta que la transformación Box-Muller combinada con algunos generadores uniformes sufre de una anomalía llamada Efecto Neave.

   H. R. Neave, "On using the Box-Muller transformation with multiplicative congruential pseudorandom number generators," Applied Statistics, 22, 92-97, 1973

3. Para mejor precisión Sugiero dibujar uniformes y aplicar la distribución normal acumulada inversa para llegar a variaciones normalmente distribuidas. Puede encontrar un algoritmo muy bueno para la distribución normal acumulada inversa en

Https://web.archive.org/web/20151030215612/http://home.online.no/~pjacklam/notes/invnorm/

Esperanza que ayuda

Pedro

Un método rápido y fácil es simplemente sumar un número de números aleatorios distribuidos uniformemente y tomar su promedio. Vea el Teorema del Límite Central para una explicación completa de por qué esto funciona.

He creado un proyecto de código abierto en C++ para un benchmark de generación de números aleatorios normalmente distribuido.

Compara varios algoritmos, incluyendo

• Método del teorema del límite central
• Transformación de caja-Muller
• Método polar de Marsaglia
• Algoritmo de Zigurat
• Método de muestreo de transformada inversa.
• cpp11random usa C++11 std::normal_distribution con std::minstd_rand (en realidad es transformada de Box-Muller en clang).

Los resultados de versión de precisión simple (float) en iMac [email protected] , clang 6.1, 64-bit:

Para la corrección, el programa verifica la media, la desviación estándar, la asimetría y la curtosis de las muestras. Se encontró que el método CLT sumando 4, 8 o 16 números uniformes no tienen buena curtosis como los otros métodos.

El algoritmo Zigurat tiene mejor rendimiento que los otros. Sin embargo, no es adecuado para el paralelismo SIMD, ya que necesita búsqueda de tablas y ramas. Box-Muller con el conjunto de instrucciones SSE2/AVX es mucho más rápido (x1.79, x2.99) que la versión no SIMD del algoritmo zigurat.

Por lo tanto, sugeriré usar Box-Muller para arquitectura con conjuntos de instrucciones SIMD, y podría ser zigurat de lo contrario.

P.d. el benchmark utiliza un LCG PRNG más simple para generar números aleatorios distribuidos uniformes. Por lo tanto, puede no ser suficiente para algunas aplicaciones. Pero la comparación de rendimiento debe ser justa porque todas las implementaciones utilizan lo mismo PRNG, por lo que el benchmark prueba principalmente el rendimiento de la transformación.

Aquí hay un ejemplo de C++, basado en algunas de las referencias. Esto es rápido y sucio, es mejor no reinventar y usar la biblioteca boost.

#include "math.h" // for RAND, and rand
double sampleNormal() {
    double u = ((double) rand() / (RAND_MAX)) * 2 - 1;
    double v = ((double) rand() / (RAND_MAX)) * 2 - 1;
    double r = u * u + v * v;
    if (r == 0 || r > 1) return sampleNormal();
    double c = sqrt(-2 * log(r) / r);
    return u * c;
}

Puede usar una gráfica Q-Q para examinar los resultados y ver qué tan bien se aproxima a una distribución normal real (clasifique sus muestras 1..x, convertir los rangos en proporciones de la cuenta total de x es decir. cuántas muestras, obtener los valores z y trazarlos. Una línea recta hacia arriba es el resultado deseado).

Use std::tr1::normal_distribution.

El espacio de nombres std::tr1 no es parte de boost. Es el espacio de nombres que contiene las adiciones de la biblioteca del Informe Técnico de C++ 1 y está disponible en los compiladores de Microsoft y gcc actualizados, independientemente de boost.

Así es como se generan las muestras en un compilador moderno de C++.

#include <random>
...
std::mt19937 generator;
double mean = 0.0;
double stddev  = 1.0;
std::normal_distribution<double> normal(mean, stddev);
cerr << "Normal: " << normal(generator) << endl;

Puede usar el GSL. Algunos ejemplos completos se dan para demostrar cómo usarlo.

Echa un vistazo a: http://www.cplusplus.com/reference/random/normal_distribution / . Es la forma más sencilla de producir distribuciones normales.

Si está usando C++11, puede usar std::normal_distribution:

#include <random>

std::default_random_engine generator;
std::normal_distribution<double> distribution(/*mean=*/0.0, /*stddev=*/1.0);

double randomNumber = distribution(generator);

Hay muchas otras distribuciones que puede usar para transformar la salida del motor de números aleatorios.

He seguido la definición del PDF dada en http://www.mathworks.com/help/stats/normal-distribution.html y se le ocurrió esto:

const double DBL_EPS_COMP = 1 - DBL_EPSILON; // DBL_EPSILON is defined in <limits.h>.
inline double RandU() {
    return DBL_EPSILON + ((double) rand()/RAND_MAX);
}
inline double RandN2(double mu, double sigma) {
    return mu + (rand()%2 ? -1.0 : 1.0)*sigma*pow(-log(DBL_EPS_COMP*RandU()), 0.5);
}
inline double RandN() {
    return RandN2(0, 1.0);
}

Tal vez no sea el mejor enfoque, pero es bastante simple.

Implementación de Box-Muller:

#include <cstdlib>
#include <cmath>
#include <ctime>
#include <iostream>
using namespace std;
 // return a uniformly distributed random number
double RandomGenerator()
{
  return ( (double)(rand()) + 1. )/( (double)(RAND_MAX) + 1. );
}
 // return a normally distributed random number
double normalRandom()
{
  double y1=RandomGenerator();
  double y2=RandomGenerator();
  return cos(2*3.14*y2)*sqrt(-2.*log(y1));
}

int main(){
double sigma = 82.;
double Mi = 40.;
  for(int i=0;i<100;i++){
double x = normalRandom()*sigma+Mi;
    cout << " x = " << x << endl;
  }
  return 0;
}

1) Gráficamente forma intuitiva que puede generar números aleatorios gaussianos es mediante el uso de algo similar al método de Monte Carlo. Generaría un punto aleatorio en una caja alrededor de la curva gaussiana usando su generador de números pseudo-aleatorios en C. Puede calcular si ese punto está dentro o debajo de la distribución gaussiana usando la ecuación de la distribución. Si ese punto está dentro de la distribución gaussiana, entonces tienes tu número aleatorio gaussiano como el valor x del punto.

Este método no es perfecto porque técnicamente la curva gaussiana continúa hacia el infinito, y no se podría crear una caja que se acerca al infinito en la dimensión x. Pero la curva de Guassian se acerca a 0 en la dimensión y bastante rápido, así que no me preocuparía por eso. La restricción del tamaño de sus variables en C puede ser más un factor limitante para su precisión.

2) Otra forma sería utilizar el Teorema del Límite Central que establece que cuando las variables aleatorias independientes se añaden, forman una distribución normal. Manteniendo este teorema en mente, puede aproximar un número aleatorio gaussiano agregando una gran cantidad de variables aleatorias independientes.

Estos métodos no son los más prácticos, pero eso es de esperar cuando no desea usar una biblioteca preexistente. Tenga en cuenta que esta respuesta viene de alguien con poca o ninguna experiencia en cálculo o estadística.

Echa un vistazo a lo que encontré.

Esta biblioteca utiliza el algoritmo Zigurat.

El comp.lang.c FAQ list comparte tres formas diferentes de generar fácilmente números aleatorios con una distribución gaussiana.

Puedes echarle un vistazo: http://c-faq.com/lib/gaussian.html

Existen varios algoritmos para la distribución normal acumulada inversa. Los más populares en finanzas cuantitativas se prueban en http://chasethedevil.github.io/post/monte-carlo--inverse-cumulative-normal-distribution/

Además, muestra el inconveniente de los enfoques similares a Zigurat.

La computadora es un dispositivo determinista. No hay aleatoriedad en el cálculo. Por otra parte el dispositivo aritmético en CPU puede evaluar summ sobre un cierto sistema finito de números enteros (realizando la evaluación en campo finito) y sistema finito de números racionales reales. Y también realizó operaciones bitwise. Matemáticas tomar un acuerdo con más grandes conjuntos como [0.0, 1.0] con un número infinito de puntos.

Puede escuchar algún cable dentro de la computadora con algún controlador, pero ¿tendría distribuciones uniformes? Me no sé. Pero si se asume que su señal es el resultado de acumular valores gran cantidad de variables aleatorias independientes, entonces recibirá aproximadamente normal variable aleatoria distribuida (Se demostró en la Teoría de la Probabilidad)

Existen algoritmos llamados - pseudo generador aleatorio. Como me sentía el propósito de pseudo generador aleatorio es emular la aleatoriedad. Y el criterio de la bondad es: - la distribución empírica es convergente (en algún sentido, uniforme, L2) a teórico - los valores que recibe de random generator parecen ser independientes. Por supuesto que no es cierto desde el "punto de vista real", pero asumimos que es cierto.

Uno de los métodos populares - puede summ 12 i.r.v con distribuciones uniformes....Pero para ser honesto durante el Teorema de Límite Central de derivación con ayuda de la transformada de Fourier, Serie de Taylor, es necesario tener n - >+suposiciones inf veces. Así que por ejemplo teórico - Personalmente no entiendo cómo la gente realizar summ de 12 i.r. v. con distribución uniforme.

Tuve la teoría de la probilidad en la universidad. Y particularmente para mí es solo una pregunta matemática. En la universidad vi el siguiente modelo:

double generateUniform(double a, double b)
{
  return uniformGen.generateReal(a, b);
}

double generateRelei(double sigma)
{
  return sigma * sqrt(-2 * log(1.0 - uniformGen.generateReal(0.0, 1.0 -kEps)));
}
double generateNorm(double m, double sigma)
{
  double y2 = generateUniform(0.0, 2 * kPi);
  double y1 = generateRelei(1.0);
  double x1 = y1 * cos(y2);
  return sigma*x1 + m;
}

De tal manera que todo era solo un ejemplo, supongo que existen otras formas de implementarlo.

Prueba de que es correcto se puede encontrar en este libro "Moscow, BMSTU, 2004: XVI Probability Theory, Example 6.12, p. 246-247" of Krishchenko Alexander Petrovich ISBN 5-7038-2485-0

Desafortunadamente no conozco la existencia de la traducción de este libro al inglés.