Gerar números aleatórios seguindo uma distribuição normal em C / C ++

Respostas:

92

Existem muitos métodos para gerar números com distribuição gaussiana a partir de um RNG regular .

A transformação Box-Muller é comumente usada. Ele produz valores corretamente com uma distribuição normal. A matemática é fácil. Você gera dois números aleatórios (uniformes) e, ao aplicar uma fórmula a eles, obtém dois números aleatórios normalmente distribuídos. Retorne um e salve o outro para a próxima solicitação de um número aleatório.

S.Lott
fonte
10
Se você precisa de velocidade, o método polar é mais rápido. E o algoritmo Zigurate ainda mais (embora seja muito mais complexo de escrever).
Joey
2
encontrei uma implementação do Zigurate aqui people.sc.fsu.edu/~jburkardt/c_src/ziggurat/ziggurat.html Está bastante completo.
dwbrito
24
Observe, C ++ 11 adiciona o std::normal_distributionque faz exatamente o que você pede sem se aprofundar em detalhes matemáticos.
3
std :: normal_distribution não tem garantia de consistência em todas as plataformas. Estou fazendo os testes agora e o MSVC fornece um conjunto diferente de valores de, por exemplo, Clang. Os motores C ++ 11 parecem gerar as mesmas sequências (dada a mesma semente), mas as distribuições C ++ 11 parecem ser implementadas usando algoritmos diferentes em plataformas diferentes.
Arno Duvenhage de
47

C ++ 11

C ++ 11 oferece std::normal_distribution, que é o caminho que eu faria hoje.

C ou C ++ mais antigo

Aqui estão algumas soluções em ordem de complexidade crescente:

  1. Adicione 12 números aleatórios uniformes de 0 a 1 e subtraia 6. Isso corresponderá à média e ao desvio padrão de uma variável normal. Uma desvantagem óbvia é que o intervalo é limitado a ± 6 - ao contrário de uma distribuição normal verdadeira.

  2. A transformação Box-Muller. Isso está listado acima e é relativamente simples de implementar. Se você precisar de amostras muito precisas, no entanto, esteja ciente de que a transformada Box-Muller combinada com alguns geradores uniformes sofre de uma anomalia chamada Efeito Neave 1 .

  3. Para melhor precisão, sugiro desenhar uniformes e aplicar a distribuição normal cumulativa inversa para chegar às variáveis ​​normalmente distribuídas. Aqui está um algoritmo muito bom para distribuições normais cumulativas inversas.

1. HR Neave, "Sobre o uso da transformação de Box-Muller com geradores de números pseudo-aleatórios congruenciais multiplicativos", Applied Statistics, 22, 92-97, 1973

Peter G.
fonte
por acaso você teria outro link para o pdf sobre o efeito Neave? ou a referência do artigo de periódico original? obrigado
pyCthon
2
@stonybrooknick A referência original é adicionada. Observação legal: enquanto pesquisava "box muller neave" no Google para encontrar a referência, esta questão stackoverflow apareceu na primeira página de resultados!
Peter G.
sim, não é muito conhecido fora de certas pequenas comunidades e grupos de interesse
pyCthon
@Peter G. Por que alguém votaria negativamente em sua resposta? - possivelmente a mesma pessoa fez meu comentário abaixo também, o que não me incomoda, mas achei sua resposta muito boa. Seria bom se os votos negativos do SO fizessem um comentário real ... Suspeito que a maioria dos votos negativos de tópicos antigos são apenas frívolos e triviais.
Pete855217
"Adicione 12 números uniformes de 0-1 e subtraia 6." - a distribuição desta variável terá distribuição normal? Você pode fornecer um link com a derivação, porque durante a derivação o teorema do limite central, n -> + inf é uma suposição muito necessária.
bruziuz
31

Um método rápido e fácil é somar um número de números aleatórios distribuídos uniformemente e obter sua média. Veja o Teorema do Limite Central para uma explicação completa de por que isso funciona.

Paul R
fonte
1 Abordagem muito interessante. É verificado que realmente dá subconjuntos normalmente distribuídos para grupos menores?
Morlock
4
@Morlock Quanto maior o número de amostras que você calcula, mais perto você chega de uma distribuição gaussiana. Se seu aplicativo tem requisitos estritos para a precisão da distribuição, então você pode ser melhor usar algo mais rigoroso, como Box-Muller, mas para muitos aplicativos, por exemplo, gerar ruído branco para aplicativos de áudio, você pode se safar com um número bastante pequeno de amostras médias (por exemplo, 16).
Paul R
2
Além disso, como você parametriza isso para obter uma certa variação, digamos que você queira uma média de 10 com um desvio padrão de 1?
Morlock
1
@Ben: você poderia me apontar um algoritmo eficiente para isso? Eu só usei a técnica de média para gerar aproximadamente ruído gaussiano para processamento de áudio e imagem com restrições de tempo real - se houver uma maneira de conseguir isso em menos ciclos de clock, isso poderia ser muito útil.
Paul R
1
@Petter: você provavelmente está certo no caso geral, para valores de ponto flutuante. Ainda existem áreas de aplicação como áudio, onde você quer ruído gaussiano de número inteiro rápido (ou ponto fixo) e a precisão não é muito importante, onde o método de média simples é mais eficiente e útil (especialmente para aplicativos incorporados, onde pode nem mesmo ser suporte de ponto flutuante de hardware).
Paul R
24

Eu criei um projeto de código aberto C ++ para benchmark de geração de número aleatório normalmente distribuído .

Ele compara vários algoritmos, incluindo

  • Método do teorema do limite central
  • Transformada Box-Muller
  • Método polar Marsaglia
  • Algoritmo Zigurate
  • Método de amostragem por transformação inversa.
  • cpp11randomusa C ++ 11 std::normal_distributioncom std::minstd_rand(na verdade é uma transformação Box-Muller em clang).

Os resultados da floatversão de precisão única ( ) no iMac [email protected], clang 6.1, 64 bits:

normaldistf

Para correção, o programa verifica a média, desvio padrão, assimetria e curtose das amostras. Verificou-se que o método CLT pela soma de 4, 8 ou 16 números uniformes não apresenta boa curtose como os demais métodos.

O algoritmo Zigurate tem melhor desempenho que os outros. No entanto, ele não é adequado para paralelismo SIMD, pois precisa de consulta de tabela e ramificações. Box-Muller com conjunto de instruções SSE2 / AVX é muito mais rápido (x1,79, x2,99) do que a versão não SIMD do algoritmo zigurate.

Portanto, vou sugerir o uso de Box-Muller para arquitetura com conjuntos de instruções SIMD, e pode zigurar caso contrário.


PS o benchmark usa um LCG PRNG mais simples para gerar números aleatórios distribuídos uniformes. Portanto, pode não ser suficiente para alguns aplicativos. Mas a comparação de desempenho deve ser justa porque todas as implementações usam o mesmo PRNG, portanto, o benchmark testa principalmente o desempenho da transformação.

Milo Yip
fonte
2
"Mas a comparação de desempenho deve ser justa porque todas as implementações usam o mesmo PRNG" .. Exceto que o BM usa um RN de entrada por saída, enquanto o CLT usa muitos mais, etc ... então o tempo para gerar um # aleatório uniforme é importante.
greggo
14

Aqui está um exemplo de C ++, baseado em algumas das referências. Isso é rápido e sujo, é melhor você não reinventar e usar a 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;
}

Você pode usar um gráfico QQ para examinar os resultados e ver o quão bem ele se aproxima de uma distribuição normal real (classifique suas amostras 1..x, transforme as classificações em proporções de contagem total de x, ou seja, quantas amostras, obtenha os valores z e plote-os. Uma linha reta para cima é o resultado desejado).

Pete855217
fonte
1
O que é sampleNormalManual ()?
resolvendo enigmas de
@solvingPuzzles - desculpe, corrigiu o código. É uma chamada recursiva.
Pete855217
1
Isso está fadado a falhar em algum evento raro (mostrar o aplicativo para seu chefe soa um sino?). Isso deve ser implementado usando um loop, não usando recursão. O método parece desconhecido. Qual é a fonte / como é chamada?
o suíno de
Box-Muller transcrito de uma implementação java. Como eu disse, é rápido e sujo, fique à vontade para consertar.
Pete855217
1
FWIW, muitos compiladores serão capazes de transformar essa chamada recursiva em um 'salto para o início da função'. A questão é se você quer contar com isso :-) Além disso, a probabilidade de que sejam> 10 iterações é de 1 em 4,8 milhões. p (> 20) é o quadrado disso, etc.
greggo
12

Use std::tr1::normal_distribution.

O namespace std :: tr1 não faz parte do boost. É o namespace que contém as adições de biblioteca do C ++ Technical Report 1 e está disponível em compiladores Microsoft e gcc atualizados, independentemente do boost.

JoeG
fonte
25
Ele não pediu padrão, ele pediu 'não impulso'.
JoeG
12

É assim que você gera os exemplos em um compilador 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
fonte
o generatordeve realmente ser semeado.
Walter
É sempre semeado. Existe uma semente padrão.
Petter
4

Se estiver usando C ++ 11, você pode 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);

Existem muitas outras distribuições que você pode usar para transformar a saída do mecanismo de números aleatórios.

Drew Noakes
fonte
Isso já foi mencionado por Ben ( stackoverflow.com/a/11977979/635608 )
Mat
3

Segui a definição do PDF fornecida em http://www.mathworks.com/help/stats/normal-distribution.html e cheguei a isto:

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

Talvez não seja a melhor abordagem, mas é bastante simples.

MJVC
fonte
-1 Não funciona, por exemplo, RANDN2 (0,0, d + 1,0). As macros são notórias por isso.
Petter
A macro falhará se rand()de RANDUretornar um zero, pois Ln (0) é indefinido.
interDist
Você realmente tentou este código? Parece que você criou uma função que gera números que são distribuídos por Rayleigh . Compare com a transformação Box-Muller , onde eles se multiplicam por cos(2*pi*rand/RAND_MAX), enquanto você multiplica por (rand()%2 ? -1.0 : 1.0).
HelloGoodbye
1

A lista de perguntas frequentes do comp.lang.c compartilha três maneiras diferentes de gerar facilmente números aleatórios com uma distribuição Gaussiana.

Você pode dar uma olhada: http://c-faq.com/lib/gaussian.html

Delgan
fonte
1

Implementação do 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;
}
Sysadmin
fonte
1

Existem vários algoritmos para a distribuição normal cumulativa inversa. Os mais populares em finanças quantitativas são testados em http://chasethedevil.github.io/post/monte-carlo--inverse-cumulative-normal-distribution/

Na minha opinião, não há muito incentivo para usar algo diferente do algoritmo AS241 de Wichura : é uma máquina de precisão, confiável e rápido. Os gargalos raramente estão na geração de números aleatórios gaussianos.

Além disso, mostra a desvantagem das abordagens do tipo Zigurate.

A principal resposta aqui defende a Box-Müller, você deve estar ciente de que ela tem deficiências conhecidas. Cito https://www.sciencedirect.com/science/article/pii/S0895717710005935 :

na literatura, Box-Muller é às vezes considerado ligeiramente inferior, principalmente por duas razões. Primeiro, se alguém aplicar o método Box-Muller a números de um gerador de congruência linear ruim, os números transformados fornecem uma cobertura extremamente pobre do espaço. Lotes de números transformados com caudas em espiral podem ser encontrados em muitos livros, mais notavelmente no livro clássico de Ripley, que foi provavelmente o primeiro a fazer esta observação "

Jherek
fonte
0

1) A maneira graficamente intuitiva de gerar números aleatórios gaussianos é usando algo semelhante ao método de Monte Carlo. Você geraria um ponto aleatório em uma caixa ao redor da curva gaussiana usando seu gerador de números pseudo-aleatórios em C. Você pode calcular se esse ponto está dentro ou abaixo da distribuição gaussiana usando a equação da distribuição. Se esse ponto estiver dentro da distribuição gaussiana, você obteve seu número aleatório gaussiano como o valor x do ponto.

Este método não é perfeito porque tecnicamente a curva gaussiana segue em direção ao infinito, e você não poderia criar uma caixa que se aproxima do infinito na dimensão x. Mas a curva Guassiana se aproxima de 0 na dimensão y muito rápido, então eu não me preocuparia com isso. A restrição do tamanho de suas variáveis ​​em C pode ser mais um fator limitante para sua precisão.

2) Outra forma seria usar o Teorema do Limite Central, que afirma que quando variáveis ​​aleatórias independentes são adicionadas, elas formam uma distribuição normal. Mantendo esse teorema em mente, você pode aproximar um número aleatório gaussiano adicionando uma grande quantidade de variáveis ​​aleatórias independentes.

Esses métodos não são os mais práticos, mas isso é esperado quando você não deseja usar uma biblioteca pré-existente. Tenha em mente que esta resposta vem de alguém com pouca ou nenhuma experiência em cálculo ou estatística.

dan dan
fonte
0

Método de Monte Carlo A maneira mais intuitiva de fazer isso seria usar um método de monte carlo. Pegue um intervalo adequado -X, + X. Valores maiores de X resultarão em uma distribuição normal mais precisa, mas leva mais tempo para convergir. uma. Escolha um número aleatório z entre -X a X. b. Mantenha a probabilidade de N(z, mean, variance)onde N é a distribuição gaussiana. Solte caso contrário e volte para a etapa (a).

Jagat
fonte
-1

Dê uma olhada no que encontrei.

Esta biblioteca usa o algoritmo Zigurate.

dwbrito
fonte
-3

O computador é um dispositivo determinístico. Não há aleatoriedade no cálculo. Além disso, o dispositivo aritmético na CPU pode avaliar a soma sobre algum conjunto finito de números inteiros (realizando avaliação em campo finito) e conjunto finito de números racionais reais. E também executou operações bit a bit. A matemática trata de conjuntos maiores como [0,0, 1,0] com número infinito de pontos.

Você pode ouvir algum fio dentro do computador com algum controlador, mas teria distribuição uniforme? Eu não sei. Mas se for assumido que seu sinal é o resultado de acumular valores enormes quantidades de variáveis ​​aleatórias independentes, então você receberá variáveis ​​aleatórias distribuídas aproximadamente normais (Isso foi provado na Teoria da Probabilidade)

Existem algoritmos chamados - gerador pseudo aleatório. Como eu senti, o objetivo do gerador pseudo aleatório é emular a aleatoriedade. E o critério de goodnes é: - a distribuição empírica é convergente (em algum sentido - pontual, uniforme, L2) para teórico - os valores que você recebe do gerador aleatório parecem ser independentes. Claro que não é verdade do 'ponto de vista real', mas presumimos que seja verdade.

Um dos métodos populares - você pode somar 12 irv com distribuições uniformes ... Mas para ser honesto durante a derivação Teorema do Limite Central com ajuda da Transformada de Fourier, Série de Taylor, é necessário ter n -> + inf suposições algumas vezes. Então, por exemplo, teoricamente - Pessoalmente, não entendo como as pessoas realizam a soma de 12 IRV com distribuição uniforme.

Tive teoria da probabilidade na universidade. E, particularmente para mim, é apenas uma questão de matemática. Na universidade, vi o seguinte 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;
}

Assim como fazer foi apenas um exemplo, acho que existem outras formas de implementá-lo.

Prova de que está correto pode ser encontrada neste livro "Moscow, BMSTU, 2004: XVI Probability Theory, Example 6.12, p.246-247" de Krishchenko Alexander Petrovich ISBN 5-7038-2485-0

Infelizmente não sei da existência de tradução deste livro para o inglês.

Bruziuz
fonte
Eu tenho vários votos negativos. Deixe-me saber o que é ruim aqui?
bruziuz
A questão é como gerar números pseudoaleatórios no computador (eu sei, a linguagem é solta aqui), não é uma questão de existência matemática.
user2820579
Sim, você está certo. E a resposta é como gerar um número pseudoaleatório com distribuição normal baseado em gerador que possui distribuição uniforme. O código-fonte foi fornecido, você pode reescrevê-lo em qualquer idioma.
bruziuz
Claro, acho que o cara está procurando, por exemplo, "Receitas numéricas em C / C ++". A propósito, apenas para complementar nossa discussão, os autores deste último livro fornecem referências interessantes para alguns geradores pseudoaleatórios que atendem aos padrões de serem geradores "decentes".
user2820579
1
Fiz backup aqui: sites.google.com/site/burlachenkok/download
bruziuz