Convertendo uma distribuição uniforme em uma distribuição normal

106

Como posso converter uma distribuição uniforme (como a maioria dos geradores de números aleatórios produzem, por exemplo, entre 0,0 e 1,0) em uma distribuição normal? E se eu quiser uma média e um desvio padrão de minha escolha?

Terhorst
fonte
3
Você tem uma especificação de linguagem ou esta é apenas uma questão geral de algoritmo?
Bill the Lizard
3
Questão geral do algoritmo. Eu não me importo com qual idioma. Mas eu preferiria que a resposta não dependesse da funcionalidade específica que apenas aquela linguagem oferece.
Terhorst,

Respostas:

47

O algoritmo Zigurate é bastante eficiente para isso, embora a transformação Box-Muller seja mais fácil de implementar do zero (e não muito lenta).

Tyler
fonte
7
Os avisos usuais sobre geradores congruentes lineares se aplicam a ambos os métodos, então use um gerador subjacente decente. Felicidades.
dmckee --- ex-moderador gatinho
3
Como Mersenee Twister, ou você tem outras sugestões?
Gregg Lind
47

Existem muitos métodos:

  • fazer não usar Box Muller. Especialmente se você desenhar muitos números gaussianos. Box Muller produz um resultado que fica entre -6 e 6 (assumindo precisão dupla. As coisas pioram com flutuações.). E é realmente menos eficiente do que outros métodos disponíveis.
  • O Zigurate é bom, mas precisa de uma consulta de tabela (e alguns ajustes específicos da plataforma devido a problemas de tamanho do cache)
  • Proporção de uniformes é o meu favorito, apenas algumas adições / multiplicações e um log 1/50 do tempo (por exemplo, olhe lá ).
  • Inverter o CDF é eficiente (e esquecido, por quê?), Você tem implementações rápidas disponíveis se você pesquisar no google. É obrigatório para números quase aleatórios.
Alexandre C.
fonte
2
Tem certeza sobre a fixação [-6,6]? Este é um ponto bastante significativo se verdadeiro (e digno de uma nota na página da wikipedia).
redcalx
1
@locster: isso é o que um professor meu me disse (ele estudou esses geradores, e confio em sua palavra). Posso encontrar uma referência para você.
Alexandre C.
7
@locster: esta propriedade indesejável também é compartilhada pelo método CDF inverso. Consulte cimat.mx/~src/prope08/randomgauss.pdf . Isso pode ser aliviado usando um RNG uniforme que tem probabilidade diferente de zero de produzir um número de ponto flutuante muito próximo de zero. A maioria dos RNG não, pois eles geram um inteiro (normalmente de 64 bits) que é mapeado para [0,1]. Isso torna esses métodos inadequados para amostragem de caudas de variáveis ​​gaussianas (pense em precificar opções de ataque alto / baixo em finanças computacionais).
Alexandre C.
6
@AlexandreC. Só para ficar claro em dois pontos, usando números de 64 bits, as caudas vão para 8,57 ou 9,41 (o valor mais baixo corresponde à conversão para [0,1) antes de tirar o log). Mesmo se fixado em [-6, 6], as chances de estar fora desse intervalo são cerca de 1,98e-9, bom o suficiente para a maioria das pessoas, mesmo na ciência. Para os números de 8,57 e 9,41, torna-se 1,04e-17 e 4,97e-21. Esses números são tão pequenos que a diferença entre uma amostragem Box Muller e uma amostragem verdadeira gaussiana em termos de dito limite é quase puramente acadêmica. Se precisar de algo melhor, basta somar quatro deles e dividir por 2.
CrazyCasta
6
Acho que a sugestão de não usar a transformação Box Muller é enganosa para uma grande porcentagem de usuários. É ótimo saber sobre a limitação, mas, como CrazyCasta aponta, para a maioria dos aplicativos que não dependem muito de outliers, você provavelmente não precisa se preocupar com isso. Por exemplo, se você já dependeu da amostragem de um normal usando numpy, você dependeu da transformação Box Muller (forma de coordenada polar) github.com/numpy/numpy/blob/… .
Andreas Grivas
30

Alterar a distribuição de qualquer função para outra envolve o uso do inverso da função desejada.

Em outras palavras, se você almeja uma função de probabilidade específica p (x), você obtém a distribuição integrando-a -> d (x) = integral (p (x)) e usa seu inverso: Inv (d (x)) . Agora use a função de probabilidade aleatória (que tem distribuição uniforme) e lance o valor do resultado por meio da função Inv (d (x)). Você deve obter valores aleatórios lançados com distribuição de acordo com a função escolhida.

Esta é a abordagem matemática genérica - ao usá-la agora você pode escolher qualquer probabilidade ou função de distribuição que você tenha, desde que tenha uma aproximação inversa ou boa.

Espero que tenha ajudado e obrigado pelo pequeno comentário sobre o uso da distribuição e não a probabilidade em si.

Adi
fonte
4
+1 Este é um método esquecido para gerar variáveis ​​gaussianas que funciona muito bem. O CDF inverso pode ser calculado com eficiência com o método de Newton neste caso (a derivada é e ^ {- t ^ 2}), uma aproximação inicial é fácil de obter como uma fração racional, então você precisa de 3-4 avaliações de erf e exp. É obrigatório se você usar números quase aleatórios, um caso em que você deve usar exatamente um número uniforme para obter um gaussiano.
Alexandre C.
9
Observe que você precisa inverter a função de distribuição cumulativa, não a função de distribuição de probabilidade. Alexandre sugere isso, mas achei que mencionar mais explicitamente não faria mal - já que a resposta parece sugerir o PDF
ltjax
Você pode usar o PDF se estiver preparado para selecionar aleatoriamente uma direção em relação à média; eu entendo isso certo?
Mark McKenna
2
Isso é chamado de amostragem de transformação inversa
dashesy
1
Aqui está uma questão relacionada em SE com uma resposta mais generalizada com uma boa explicação.
rápido
23

Aqui está uma implementação de javascript usando a forma polar da transformação Box-Muller.

/*
 * Returns member of set with a given mean and standard deviation
 * mean: mean
 * standard deviation: std_dev 
 */
function createMemberInNormalDistribution(mean,std_dev){
    return mean + (gaussRandom()*std_dev);
}

/*
 * Returns random number in normal distribution centering on 0.
 * ~95% of numbers returned should fall between -2 and 2
 * ie within two standard deviations
 */
function gaussRandom() {
    var u = 2*Math.random()-1;
    var v = 2*Math.random()-1;
    var r = u*u + v*v;
    /*if outside interval [0,1] start over*/
    if(r == 0 || r >= 1) return gaussRandom();

    var c = Math.sqrt(-2*Math.log(r)/r);
    return u*c;

    /* todo: optimize this algorithm by caching (v*c) 
     * and returning next time gaussRandom() is called.
     * left out for simplicity */
}
user5084
fonte
5

Use o teorema do limite central, entrada da Wikipédia, entrada mathworld a seu favor.

Gere n dos números uniformemente distribuídos, some-os, subtraia n * 0,5 e você terá a saída de uma distribuição aproximadamente normal com média igual a 0 e variância igual a (1/12) * (1/sqrt(N))(ver wikipedia sobre distribuições uniformes para esse último)

n = 10 dá a você algo meio decente rápido. Se você quiser algo mais da metade decente, vá para a solução Tylers (conforme observado no entrada Wikipédia em distribuições normais )

jilles de wit
fonte
1
Isso não fornecerá um normal particularmente próximo (as "caudas" ou pontos finais não estarão próximos da distribuição normal real). Box-Muller é melhor, como outros sugeriram.
Peter K.
1
Box Muller também tem caudas erradas (retorna um número entre -6 e 6 em precisão dupla)
Alexandre C.
n = 12 (soma 12 números aleatórios no intervalo de 0 a 1 e subtrai 6) resulta em desd padrão = 1 e média = 0. Isso pode então ser usado para gerar qualquer distribuição normal. Simplesmente multiplique o resultado pelo stddev desejado e adicione a média.
JerryM
3

Eu usaria o Box-Muller. Duas coisas sobre isso:

  1. Você acaba com dois valores por iteração
    Normalmente, você armazena em cache um valor e retorna o outro. Na próxima chamada para uma amostra, você retorna o valor em cache.
  2. Box-Muller dá uma pontuação Z
    Você deve então dimensionar a pontuação Z pelo desvio padrão e adicionar a média para obter o valor total na distribuição normal.
Hughdbrown
fonte
Como você escala o Z-score?
Terhorst,
3
scaled = mean + stdDev * zScore // dá-lhe normal (mean, stdDev ^ 2)
yoyoyoyosef
2

Onde R1, R2 são números uniformes aleatórios:

DISTRIBUIÇÃO NORMAL, com SD de 1: sqrt (-2 * log (R1)) * cos (2 * pi * R2)

Isso é exato ... não há necessidade de fazer todos aqueles loops lentos!

Erik Aronesty
fonte
Antes que alguém me corrigisse ... aqui está a aproximação que eu fiz: (1,5- (R1 + R2 + R3)) * 1,88. Eu gosto disso também.
Erik Aronesty,
2

Parece incrível que eu pudesse adicionar algo a isso depois de oito anos, mas para o caso de Java, gostaria de apontar aos leitores o método Random.nextGaussian () , que gera uma distribuição gaussiana com média 0,0 e desvio padrão 1,0 para você.

Uma simples adição e / ou multiplicação mudará a média e o desvio padrão de acordo com suas necessidades.

Pepijn Schmitz
fonte
1

O módulo de biblioteca Python padrão aleatório tem o que você deseja:

normalvariate (mu, sigma)
Distribuição normal. mu é a média e sigma é o desvio padrão.

Para o algoritmo em si, dê uma olhada na função em random.py na biblioteca Python.

A entrada manual está aqui

Brent.Longborough
fonte
2
Infelizmente, a biblioteca de python usa Kinderman, AJ e Monahan, JF, "Computer generation of random variables using the ratio of uniform desvio", ACM Trans Math Software, 3, (1977), pp257-260. Isso usa duas variáveis ​​aleatórias uniformes para gerar o valor normal, ao invés de uma única, então não é óbvio como usá-lo como o mapeamento que o OP queria.
Ian,
1

Esta é minha implementação JavaScript do Algoritmo P ( Método Polar para desvios normais ) da Seção 3.4.1 do livro de Donald Knuth, The Art of Computer Programming :

function normal_random(mean,stddev)
{
    var V1
    var V2
    var S
    do{
        var U1 = Math.random() // return uniform distributed in [0,1[
        var U2 = Math.random()
        V1 = 2*U1-1
        V2 = 2*U2-1
        S = V1*V1+V2*V2
    }while(S >= 1)
    if(S===0) return 0
    return mean+stddev*(V1*Math.sqrt(-2*Math.log(S)/S))
}
Alessandro Jacopson
fonte
0

Acho que você deve tentar isso no EXCEL: =norminv(rand();0;1) . Isso produzirá os números aleatórios que devem ser normalmente distribuídos com a média zero e a variância unitária. "0" pode ser fornecido com qualquer valor, de modo que os números tenham a média desejada e, mudando "1", você obterá a variância igual ao quadrado de sua entrada.

Por exemplo: =norminv(rand();50;3)resultará em números normalmente distribuídos com MEAN = 50 VARIANCE = 9.

Hipopótamo
fonte
0

P Como posso converter uma distribuição uniforme (como a maioria dos geradores de números aleatórios produzem, por exemplo, entre 0,0 e 1,0) em uma distribuição normal?

  1. Para implementação de software, conheço alguns nomes de geradores aleatórios que fornecem uma sequência aleatória pseudo-uniforme em [0,1] (Mersenne Twister, Linear Congruate Generator). Vamos chamá-lo de U (x)

  2. Existe uma área matemática que se denomina teoria da probabilidade. Primeira coisa: se você quiser modelar RV com distribuição integral F, então você pode tentar avaliar apenas F ^ -1 (U (x)). Na teoria pr. Foi provado que tal va terá distribuição integral F.

  3. A etapa 2 pode ser aplicada para gerar rv ~ F sem o uso de quaisquer métodos de contagem quando F ^ -1 pode ser derivado analiticamente sem problemas. (por exemplo, exp.distribution)

  4. Para modelar a distribuição normal, você pode calcular y1 * cos (y2), onde y1 ~ é uniforme em [0,2pi]. e y2 é a distribuição relei.

P: E se eu quiser uma média e um desvio padrão de minha escolha?

Você pode calcular sigma * N (0,1) + m.

Pode ser mostrado que tal mudança e escala levam a N (m, sigma)

Bruziuz
fonte
0

Esta é uma implementação Matlab usando a forma polar da transformação Box-Muller :

Função randn_box_muller.m:

function [values] = randn_box_muller(n, mean, std_dev)
    if nargin == 1
       mean = 0;
       std_dev = 1;
    end

    r = gaussRandomN(n);
    values = r.*std_dev - mean;
end

function [values] = gaussRandomN(n)
    [u, v, r] = gaussRandomNValid(n);

    c = sqrt(-2*log(r)./r);
    values = u.*c;
end

function [u, v, r] = gaussRandomNValid(n)
    r = zeros(n, 1);
    u = zeros(n, 1);
    v = zeros(n, 1);

    filter = r==0 | r>=1;

    % if outside interval [0,1] start over
    while n ~= 0
        u(filter) = 2*rand(n, 1)-1;
        v(filter) = 2*rand(n, 1)-1;
        r(filter) = u(filter).*u(filter) + v(filter).*v(filter);

        filter = r==0 | r>=1;
        n = size(r(filter),1);
    end
end

E invocar histfit(randn_box_muller(10000000),100);este é o resultado: Box-Muller Matlab Histfit

Obviamente, é realmente ineficiente em comparação com o randn integrado do Matlab .

madx
fonte
0

Tenho o seguinte código que talvez possa ajudar:

set.seed(123)
n <- 1000
u <- runif(n) #creates U
x <- -log(u)
y <- runif(n, max=u*sqrt((2*exp(1))/pi)) #create Y
z <- ifelse (y < dnorm(x)/2, -x, NA)
z <- ifelse ((y > dnorm(x)/2) & (y < dnorm(x)), x, z)
z <- z[!is.na(z)]
mentes brilhantes pensam igual
fonte
0

Também é mais fácil usar a função implementada rnorm (), pois é mais rápida do que escrever um gerador de números aleatórios para a distribuição normal. Veja o código a seguir como prova

n <- length(z)
t0 <- Sys.time()
z <- rnorm(n)
t1 <- Sys.time()
t1-t0
peterweethetbeter
fonte
-2
function distRandom(){
  do{
    x=random(DISTRIBUTION_DOMAIN);
  }while(random(DISTRIBUTION_RANGE)>=distributionFunction(x));
  return x;
}

fonte
Não há garantia de retorno, não é? ;-)
Peter K.
5
Os números aleatórios são muito importantes para serem deixados ao acaso.
Drew Noakes em
Não responde à pergunta - a distribuição normal tem um domínio infinito.
Matt,