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

Respuestas:

92

Hay muchos métodos para generar números distribuidos en Gauss a partir de un RNG regular .

La transformada de Box-Muller se usa comúnmente. Produce correctamente valores con una distribución normal. Las matemáticas son fáciles. Generas dos números aleatorios (uniformes) y, al aplicarles una fórmula, obtienes dos números aleatorios distribuidos normalmente. Devuelve uno y guarda el otro para la próxima solicitud de un número aleatorio.

S.Lott
fuente
10
Sin embargo, si necesita velocidad, entonces el método polar es más rápido. Y el algoritmo Ziggurat aún más (aunque mucho más complejo de escribir).
Joey
2
Encontré una implementación del Zigurat aquí people.sc.fsu.edu/~jburkardt/c_src/ziggurat/ziggurat.html Es bastante completo.
dwbrito
24
Tenga en cuenta que C ++ 11 agrega std::normal_distributionque hace exactamente lo que pide sin profundizar en los detalles matemáticos.
3
No se garantiza que std :: normal_distribution sea coherente en todas las plataformas. Estoy haciendo las pruebas ahora y MSVC proporciona un conjunto de valores diferente de, por ejemplo, Clang. Los motores de C ++ 11 parecen generar las mismas secuencias (dada la misma semilla), pero las distribuciones de C ++ 11 parecen implementarse utilizando diferentes algoritmos en diferentes plataformas.
Arno Duvenhage
47

C ++ 11

Ofertas de C ++ 11 std::normal_distribution, que es el camino que seguiría hoy.

C o anterior C ++

Aquí hay algunas soluciones en orden ascendente de complejidad:

  1. Sume 12 números aleatorios uniformes de 0 a 1 y reste 6. Esto coincidirá con la desviación estándar y media de una variable normal. Un inconveniente obvio es que el rango está limitado a ± 6, a diferencia de una verdadera distribución normal.

  2. La transformación de Box-Muller. Esto se enumera anteriormente y es relativamente sencillo de implementar. Sin embargo, si necesita muestras muy precisas, tenga en cuenta que la transformada Box-Muller combinada con algunos generadores uniformes sufre una anomalía llamada Neave Effect 1 .

  3. Para obtener la mejor precisión, sugiero dibujar uniformes y aplicar la distribución normal acumulativa inversa para llegar a variantes distribuidas normalmente. Aquí hay un muy buen algoritmo para distribuciones normales acumulativas inversas.

1. HR Neave, "Sobre el uso de la transformación de Box-Muller con generadores de números pseudoaleatorios congruentes multiplicativos", Estadística aplicada, 22, 92-97, 1973

Peter G.
fuente
por casualidad, ¿tendría otro enlace al pdf sobre el efecto Neave? o la referencia del artículo de revista original? gracias
pyCthon
2
@stonybrooknick Se agrega la referencia original. Comentario interesante: mientras buscaba en Google "box muller neave" para encontrar la referencia, ¡esta misma pregunta de stackoverflow apareció en la primera página de resultados!
Peter G.
sí, no es muy conocido fuera de ciertas comunidades pequeñas y grupos de interés
pyCthon
@Peter G. ¿Por qué alguien rechazaría tu respuesta? - posiblemente la misma persona también hizo mi comentario a continuación, lo cual estoy bien, pero pensé que tu respuesta fue muy buena. Sería bueno que SO hiciera votos en contra para forzar un comentario real. Sospecho que la mayoría de los votos en contra de temas antiguos son simplemente frívolos y triviales.
Pete855217
"Suma 12 números uniformes del 0-1 y resta 6." - ¿La distribución de esta variable tendrá distribución normal? ¿Puede proporcionar un vínculo con la derivación, porque durante la derivación el teorema del límite central, n -> + inf es una suposición muy necesaria?
bruziuz
31

Un método rápido y sencillo consiste simplemente en sumar un número de números aleatorios distribuidos uniformemente y calcular su promedio. Consulte el Teorema del límite central para obtener una explicación completa de por qué funciona.

Paul R
fuente
+1 Enfoque muy interesante. ¿Está verificado para dar realmente subconjuntos distribuidos normalmente para grupos más pequeños?
Morlock
4
@Morlock Cuanto mayor sea el número de muestras que promedie, más se acercará a una distribución gaussiana. Si su aplicación tiene requisitos estrictos para la precisión de la distribución, es mejor que utilice algo más riguroso, como Box-Muller, pero para muchas aplicaciones, por ejemplo, generar ruido blanco para aplicaciones de audio, puede salirse con la suya con un número bastante pequeño. de muestras promediadas (por ejemplo, 16).
Paul R
2
Además, ¿cómo se parametriza esto para obtener una cierta cantidad de varianza, digamos que desea una media de 10 con una desviación estándar de 1?
Morlock
1
@Ben: ¿podrías señalarme un algoritmo eficiente para esto? Solo he usado la técnica de promediado para generar aproximadamente ruido gaussiano para el procesamiento de audio e imágenes con restricciones en tiempo real; si hay una forma de lograr esto en menos ciclos de reloj, entonces podría ser muy útil.
Paul R
1
@Petter: probablemente tengas razón en el caso general, para valores de punto flotante. Sin embargo, todavía hay áreas de aplicación como el audio, donde desea ruido gaussiano entero rápido (o punto fijo) y la precisión no es demasiado importante, donde el método de promedio simple es más eficiente y útil (especialmente para aplicaciones integradas, donde puede que ni siquiera haya ser soporte de hardware de punto flotante).
Paul R
24

Creé un proyecto de código abierto C ++ para un punto de referencia de generación de números aleatorios distribuidos normalmente .

Compara varios algoritmos, incluidos

  • Método del teorema del límite central
  • Transformada de Box-Muller
  • Método polar de Marsaglia
  • Algoritmo Zigurat
  • Método de muestreo por transformación inversa.
  • cpp11randomusa C ++ 11 std::normal_distributioncon std::minstd_rand(en realidad es la transformación de Box-Muller en clang).

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

normaldistf

Para su 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 al sumar 4, 8 o 16 números uniformes no tiene una buena curtosis como los otros métodos.

El algoritmo Ziggurat tiene un mejor rendimiento que los demás. 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 sin SIMD del algoritmo ziggurat.

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


PD: el punto de referencia utiliza un LCG PRNG más simple para generar números aleatorios distribuidos uniformemente. Por tanto, puede que no sea suficiente para algunas aplicaciones. Pero la comparación de rendimiento debería ser justa porque todas las implementaciones utilizan el mismo PRNG, por lo que el punto de referencia prueba principalmente el rendimiento de la transformación.

Milo Yip
fuente
2
"Pero la comparación de rendimiento debería ser justa porque todas las implementaciones usan el mismo PRNG". Excepto que BM usa un RN de entrada por salida, mientras que CLT usa muchos más, etc ... así que el tiempo para generar un # aleatorio uniforme es importante.
greggo
14

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 de impulso.

#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 utilizar un gráfico QQ para examinar los resultados y ver qué tan bien se aproxima a una distribución normal real (clasifique sus muestras 1..x, convierta los rangos en proporciones del recuento total de x, es decir, cuántas muestras, obtenga los valores z y trazarlos. Una línea recta hacia arriba es el resultado deseado).

Pete855217
fuente
1
¿Qué es sampleNormalManual ()?
resolvingPuzzles
@solvingPuzzles: lo siento, se corrigió el código. Es una llamada recursiva.
Pete855217
1
Esto está destinado a fallar en algún evento raro (¿mostrar la aplicación a tu jefe te suena?). Esto debe implementarse mediante un bucle, no mediante recursividad. El método parece desconocido. ¿Cuál es la fuente / cómo se llama?
los cerdos
Box-Muller transcrito de una implementación java. Como dije, es rápido y sucio, siéntete libre de arreglarlo.
Pete855217
1
FWIW, muchos compiladores podrán convertir esa llamada recursiva en particular en un 'salto a la parte superior de la función'. La pregunta es si desea contar con ella :-) Además, la probabilidad de que se necesiten> 10 iteraciones es de 1 en 4,8 millones. p (> 20) es el cuadrado de que, etc.
Greggo
12

Utilice std::tr1::normal_distribution.

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

JoeG
fuente
25
No pidió estándar, pidió "no impulso".
JoeG
12

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

#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;
Petter
fuente
el generatorrealmente debería ser sembrado.
Walter
Siempre está sembrado. Hay una semilla predeterminada.
Petter
4

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 utilizar para transformar la salida del motor de números aleatorios.

Dibujó Noakes
fuente
Ben ya lo ha mencionado ( stackoverflow.com/a/11977979/635608 )
Mat
3

Seguí la definición del PDF dada en http://www.mathworks.com/help/stats/normal-distribution.html y se me 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);
}

Quizás no sea el mejor enfoque, pero es bastante simple.

MJVC
fuente
-1 No funciona para, por ejemplo, RANDN2 (0.0, d + 1.0). Las macros son conocidas por esto.
Petter
La macro fallará si rand()de RANDUdevuelve un cero, puesto que no está definido Ln (0).
interDist
¿De verdad has probado este código? Parece que ha creado una función que genera números distribuidos por Rayleigh . Compare con la transformada Box-Muller , donde se multiplican con cos(2*pi*rand/RAND_MAX), mientras que usted multiplica con (rand()%2 ? -1.0 : 1.0).
Hola
1

La lista de preguntas frecuentes de comp.lang.c comparte tres formas diferentes de generar fácilmente números aleatorios con una distribución gaussiana.

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

Delgan
fuente
1

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;
}
Administrador de sistemas
fuente
1

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

En mi opinión, no hay mucho incentivo para usar algo más que el algoritmo AS241 de Wichura : es una máquina de precisión, confiable y rápida. Los cuellos de botella rara vez se encuentran en la generación de números aleatorios gaussianos.

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

La respuesta principal aquí aboga por Box-Müller, debe tener en cuenta que tiene deficiencias conocidas. Cito https://www.sciencedirect.com/science/article/pii/S0895717710005935 :

En la literatura, Box-Muller a veces se considera ligeramente inferior, principalmente por dos razones. Primero, si se aplica el método de Box-Muller a números de un generador congruencial lineal incorrecto, los números transformados proporcionan una cobertura extremadamente pobre del espacio. Se pueden encontrar tramas de números transformados con colas en espiral en muchos libros, sobre todo en el libro clásico de Ripley, quien probablemente fue el primero en hacer esta observación "

jherek
fuente
0

1) La forma gráficamente intuitiva de generar números aleatorios gaussianos es utilizando algo similar al método Monte Carlo. Generaría un punto aleatorio en un cuadro alrededor de la curva gaussiana usando su generador de números pseudoaleatorios 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 avanza hacia el infinito y no se puede crear una caja que se acerque 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 usar el Teorema del límite central que establece que cuando se agregan variables aleatorias independientes, forman una distribución normal. Teniendo 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 es de esperar cuando no desea utilizar una biblioteca preexistente. Tenga en cuenta que esta respuesta proviene de alguien con poca o ninguna experiencia en cálculo o estadística.

dan dan
fuente
0

Método de Monte Carlo La forma más intuitiva de hacer esto sería utilizar un método de Monte Carlo . Tome un rango adecuado -X, + X. Los valores más grandes de X darán como resultado una distribución normal más precisa, pero tarda más en converger. a. Elija un número aleatorio z entre -X y X. b. Mantener con una probabilidad de N(z, mean, variance)donde N es la distribución gaussiana. De lo contrario, deje caer y vuelva al paso (a).

Jagat
fuente
-1

Eche un vistazo a lo que encontré.

Esta biblioteca utiliza el algoritmo Ziggurat.

dwbrito
fuente
-3

La computadora es un dispositivo determinista. No hay aleatoriedad en el cálculo. Además, el dispositivo aritmético en la CPU puede evaluar la suma sobre un conjunto finito de números enteros (realizando la evaluación en un campo finito) y un conjunto finito de números racionales reales. Y también realizó operaciones bit a bit. Las matemáticas toman un trato con conjuntos más grandes 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? No lo sé. Pero si se supone que su señal es el resultado de acumular valores de una gran cantidad de variables aleatorias independientes, recibirá una variable aleatoria distribuida aproximadamente normal (se demostró en la teoría de la probabilidad)

Existen algoritmos llamados - generador pseudoaleatorio. Como sentí, el propósito del generador pseudoaleatorio es emular la aleatoriedad. Y el criterio de bondad es: - la distribución empírica es convergente (en cierto sentido - puntual, uniforme, L2) a teórica - los valores que recibe del generador aleatorio 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 más populares: puede sumar 12 irv con distribuciones uniformes ... Pero para ser honesto durante la derivación Teorema del límite central con la ayuda de la Transformada de Fourier, Serie de Taylor, es necesario tener n -> + infosupuestos un par de veces. Entonces, por ejemplo, teóricamente, personalmente no entiendo cómo la gente realiza una suma de 12 irv con distribución uniforme.

Tenía teoría de la capacidad en la universidad. Y en particular para mí, es solo una cuestión de matemáticas. 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;
}

Así como todo fue solo un ejemplo, supongo que existen otras formas de implementarlo.

La prueba de que es correcta se puede encontrar en este libro "Moscú, BMSTU, 2004: XVI Teoría de la probabilidad, ejemplo 6.12, p.246-247" de Krishchenko Alexander Petrovich ISBN 5-7038-2485-0

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

bruziuz
fuente
Tengo varios votos negativos. Déjame saber qué hay de malo aquí.
bruziuz
La pregunta es cómo generar números pseudoaleatorios en la computadora (lo sé, el lenguaje está suelto aquí), no es una cuestión de existencia matemática.
user2820579
Sí tienes razón. Y la respuesta es cómo generar un número pseudoaleatorio con distribución normal basado en un generador que tiene distribución uniforme. Se ha proporcionado el código fuente, puede reescribirlo en cualquier idioma.
bruziuz
Claro, creo que el tipo está buscando, por ejemplo, "Recetas numéricas en C / C ++". Por cierto, sólo para complementar nuestra discusión, los autores de este último libro dan referencias interesantes para un par de generadores pseudoaleatorios que cumplen con los estándares de ser generadores "decentes".
user2820579
1
Hice una copia de seguridad aquí: sites.google.com/site/burlachenkok/download
bruziuz