Simulação de desenhos de uma distribuição uniforme usando desenhos de uma distribuição normal

15

Recentemente, comprei um recurso de entrevista em ciência de dados no qual uma das perguntas de probabilidade era a seguinte:

Dados os desenhos de uma distribuição normal com parâmetros conhecidos, como você pode simular os desenhos de uma distribuição uniforme?

Meu processo de pensamento original era que, para uma variável aleatória discreta, poderíamos dividir a distribuição normal em K subseções únicas, nas quais cada subseção tem uma área igual abaixo da curva normal. Então poderíamos determinar em qual dos valores K a variável leva reconhecendo em qual área da curva normal a variável acaba caindo.

Mas isso funcionaria apenas para variáveis ​​aleatórias discretas. Eu fiz alguma pesquisa sobre como podemos fazer o mesmo para variáveis aleatórias contínuas, mas infelizmente eu só poderia encontrar técnicas como transformada inversa amostragem que usar como entrada uma variável aleatória uniforme, e poderia saída de variáveis aleatórias de alguma outra distribuição. Eu estava pensando que talvez pudéssemos fazer esse processo ao contrário para obter variáveis ​​aleatórias uniformes?

Também pensei em possivelmente usar as variáveis ​​aleatórias normais como entradas em um gerador congruencial linear, mas não tenho certeza se isso funcionaria.

Alguma idéia de como eu poderia abordar essa questão?

Wellington
fonte

Respostas:

30

No espírito de usar cálculos algébricos simples que não estão relacionados ao cálculo da distribuição Normal , eu me inclinaria para o seguinte. Eles são ordenados como eu pensava neles (e, portanto, precisava ser cada vez mais criativo), mas guardei o melhor - e o mais surpreendente - para durar.

  1. Inverta a técnica de Box-Mueller : de cada par de normais (X,Y) , dois uniformes independentes podem ser construídos como atan2(Y,X) (no intervalo [π,π] ) e exp((X2+Y2)/2) (no intervalo [0,1] ).

  2. Pegue as normais em grupos de dois e some seus quadrados para obter uma sequência de χ22 variáveis Y1,Y2,,Yi, . As expressões obtidas dos pares

    Xi=Y2iY2i1+Y2i

    terá uma distribuição Beta(1,1) uniforme.

    Que isso exija apenas aritmética básica e simples deve ficar claro.

  3. Como a distribuição exata do coeficiente de correlação de Pearson de uma amostra de quatro pares de uma distribuição normal bivariada padrão é distribuída uniformemente em [1,1] , podemos simplesmente tomar as normais em grupos de quatro pares (ou seja, oito valores em cada conjunto) e retorne o coeficiente de correlação desses pares. (Isso envolve operações aritméticas simples mais operações com duas raízes quadradas.)

  4. Sabe-se desde os tempos antigos que uma projeção cilíndrica da esfera (uma superfície em três espaços) é de área igual . Isso implica que, na projeção de uma distribuição uniforme na esfera, tanto a coordenada horizontal (correspondente à longitude) quanto a vertical (correspondente à latitude) terão distribuições uniformes. Como a distribuição normal padrão trivariada é esfericamente simétrica, sua projeção na esfera é uniforme. A obtenção da longitude é essencialmente o mesmo cálculo do ângulo no método Box-Mueller ( qv ), mas a latitude projetada é nova. A projeção na esfera apenas normaliza um triplo de coordenadas (x,y,z) e nesse pontoz é a latitude projetada. Portanto, pegue as variáveis ​​Normal em grupos de três,X3i2,X3i1,X3i e calcule

    X3iX3i22+X3i12+X3i2

    para i=1,2,3, .

  5. Uma vez que a maioria dos sistemas de computação representar números em binário , geração de número uniforme geralmente começa por produzir uniformemente distribuídas números inteiros entre 0 e 2321 (ou algum alto poder de 2 relacionada com o comprimento da palavra de computador) e rescaling-los, conforme necessário. Tais números inteiros são representados internamente como cadeias de caracteres de 32 dígitos binários. Podemos obter bits aleatórios independentes comparando uma variável Normal com sua mediana. Assim, basta dividir as variáveis ​​normais em grupos de tamanho igual ao número desejado de bits, comparar cada uma com sua média e montar as seqüências resultantes de resultados verdadeiro / falso em um número binário. Escrevendo kpara o número de bits e H para o sinal (ou seja, H(x)=1 quando x>0 e H(x)=0 caso contrário), podemos expressar o valor uniforme normalizado resultante em [0,1) com a fórmula

    j=0k1H(Xkij)2j1.

    As variáveis podem ser obtidas de qualquer distribuição contínua cuja mediana é 0 (como uma Normal normal); eles são processados ​​em grupos de k com cada grupo produzindo um desses valores pseudo-uniformes.Xn0k

  6. A amostragem por rejeição é uma maneira padrão, flexível e poderosa de extrair variáveis ​​aleatórias de distribuições arbitrárias. Suponha que a distribuição de destino tenha PDF . Um valor Y é desenhado de acordo com outra distribuição no PDF g . Na etapa de rejeição, um valor uniforme U entre 0 e g ( Y ) é desenhado independentemente de Y e comparado a f ( Y ) : se for menor, YfYgU0g(Y)Yf(Y)Yé retido, mas, caso contrário, o processo é repetido. Essa abordagem parece circular: como podemos gerar uma variável uniforme com um processo que precisa de uma variável uniforme para começar?

    A resposta é que, na verdade, não precisamos de uma variável uniforme para executar a etapa de rejeição. Em vez disso (assumindo ), podemos jogar uma moeda justa para obter um 0 ou 1 aleatoriamente. Isso será interpretado como o primeiro bit na representação binária de uma variável uniforme U no intervalo [ 0 , 1 ) . Quando o resultado é 0 , o que significa 0 L < 1 / 2 ; de outro modo, um / 2 L < 1 . g(Y)001U[0,1)00U<1/21/2U<1Metade do tempo, este é suficiente para decidir o passo rejeição: se , mas a moeda é 0 , Y deve ser aceite; Se f ( Y ) / g ( Y ) < 1 / 2 , mas a moeda é 1 , Y deve ser rejeitado; caso contrário, precisamos jogar a moeda novamente, a fim de obter o próximo bit de U . Porque - não importa qual valor f ( Yf(Y)/g(Y)1/20Yf(Y)/g(Y)<1/21YU contém - há um 1 / 2 possibilidade de parar após cada aleta, o número esperado de flips só é 1 / 2 ( 1 ) + 1 / 4 ( 2 ) + 1 / 8 ( 3 ) + + 2 - n ( n ) + = 2 .f(Y)/g(Y)1/21/2(1)+1/4(2)+1/8(3)++2n(n)+=2

    A amostragem de rejeição pode valer a pena (e ser eficiente), desde que o número esperado de rejeições seja pequeno. Podemos fazer isso ajustando o maior retângulo possível (representando uma distribuição uniforme) abaixo de um PDF normal.

    Normal and Uniform PDFs

    Usando cálculo para otimizar a área do retângulo, você vai achar que seus endpoints deve ficar em , onde a sua altura é igual a exp ( - 1 / 2 ) / ±1, tornando sua área um pouco maior que0,48. Usando essa densidade normal padrão comogerejeitando todos os valores fora do intervalo[-1,1]automaticamente e aplicando o procedimento de rejeição, obteremos variáveis ​​uniformes em[-1,1] de formaeficiente:exp(1/2)/2π0.2419710.48g[1,1][1,1]

    • Em uma fração do tempo, a variável Normal está além de [ - 1 , 1 ] e é imediatamente rejeitada. ( Φ é o CDF normal padrão.)2Φ(1)0.317[1,1]Φ

    • Na fração restante do tempo, o procedimento de rejeição binária deve ser seguido, exigindo mais duas variáveis ​​normais em média.

    • O processo global requer uma média de etapas.1/(2exp(1/2)/2π)2.07

    O número esperado de variáveis ​​normais necessárias para produzir cada resultado uniforme funciona para

    2eπ(12Φ(1))2.82137.

    Embora isso seja bastante eficiente, observe que (1) o cálculo do PDF normal exige o cálculo de um exponencial e (2) o valor deve ser pré-computado de uma vez por todas. Ainda é um pouco menos de cálculo do que o método Box-Mueller ( qv ).Φ(1)

  7. As estatísticas de pedidos de uma distribuição uniforme têm lacunas exponenciais. Como a soma dos quadrados de duas Normais (com média zero) é exponencial, podemos gerar uma realização de uniformes independentes somando os quadrados dos pares dessas Normais, calculando a soma acumulada destes, redimensionando os resultados para cair no intervalo [ 0 , 1 ] e soltando o último (que sempre será igual a 1 ). Essa é uma abordagem agradável porque requer apenas quadratura, soma e (no final) uma única divisão.n[0,1]1

    The n values will automatically be in ascending order. If such a sorting is desired, this method is computationally superior to all the others insofar as it avoids the O(nlog(n)) cost of a sort. If a sequence of independent uniforms is needed, however, then sorting these n values randomly will do the trick. Since (as seen in the Box-Mueller method, q.v.) the ratios of each pair of Normals are independent of the sum of squares of each pair, we already have the means to obtain that random permutation: order the cumulative sums by the corresponding ratios. (If n is very large, this process could be carried out in smaller groups of k with little loss of efficiency, since each group needs only 2(k+1) Normals to create k uniform values. For fixed k, the asymptotic computational cost is therefore O(nlog(k)) = O(n), needing 2n(1+1/k) Normal variates to generate n uniform values.)

  8. To a superb approximation, any Normal variate with a large standard deviation looks uniform over ranges of much smaller values. Upon rolling this distribution into the range [0,1] (by taking only the fractional parts of the values), we thereby obtain a distribution that is uniform for all practical purposes. This is extremely efficient, requiring one of the simplest arithmetic operations of all: simply round each Normal variate down to the nearest integer and retain the excess. The simplicity of this approach becomes compelling when we examine a practical R implementation:

    rnorm(n, sd=10) %% 1
    

    reliably produces n uniform values in the range [0,1] at the cost of just n Normal variates and almost no computation.

    (Even when the standard deviation is 1, the PDF of this approximation varies from a uniform PDF, as shown in the following figure, by less than one part in 108! To detect it reliably would require a sample of 1016 values--that's already beyond the capability of any standard test of randomness. With a larger standard deviation the non-uniformity is so small it cannot even be calculated. For instance, with an SD of 10 as shown in the code, the maximum deviation from a uniform PDF is only 10857.)

    Approximate PDF


In every case Normal variables "with known parameters" can easily be recentered and rescaled into the Standard Normals assumed above. Afterwards, the resulting uniformly distributed values can be recentered and rescaled to cover any desired interval. These require only basic arithmetic operations.

The ease of these constructions is evidenced by the following R code, which uses only one or two lines for most of them. Their correctness is witnessed by the resulting near-uniform histograms based on 100,000 independent values in each case (requiring around 12 seconds for all seven simulations). For reference--in case you are worried about the amount of variation appearing in any of these plots--a histogram of uniform values simulated with R's uniform random number generator is included at the end.

Histograms

All these simulations were tested for uniformity using a χ2 test based on 1000 bins; none could be considered significantly non-uniform (the lowest p-value was 3%--for the results generated by R's actual uniform number generator!).

set.seed(17)
n <- 1e5
y <- matrix(rnorm(floor(n/2)*2), nrow=2)
x <- c(atan2(y[2,], y[1,])/(2*pi) + 1/2, exp(-(y[1,]^2+y[2,]^2)/2))
hist(x, main="Box-Mueller")

y <- apply(array(rnorm(4*n), c(2,2,n)), c(3,2), function(z) sum(z^2))
x <- y[,2] / (y[,1]+y[,2])
hist(x, main="Beta")

x <- apply(array(rnorm(8*n), c(4,2,n)), 3, function(y) cor(y[,1], y[,2]))
hist(x, main="Correlation")

n.bits <- 32; x <-  (2^-(1:n.bits)) %*% matrix(rnorm(n*n.bits) > 0, n.bits)
hist(x, main="Binary")

y <- matrix(rnorm(n*3), 3)
x <- y[1, ] / sqrt(apply(y, 2, function(x) sum(x^2)))
hist(x, main="Equal area")

accept <- function(p) { # Using random normals, return TRUE with chance `p`
  p.bit <- x <- 0
  while(p.bit == x) {
    p.bit <- p >= 1/2
    x <- rnorm(1) >= 0
    p <- (2*p) %% 1
  }
  return(x == 0)
}
y <- rnorm(ceiling(n * sqrt(exp(1)*pi/2))) # This aims to produce `n` uniforms
y <- y[abs(y) < 1]
x <- y[sapply(y, function(x) accept(exp((x^2-1)/2)))]
hist(x, main="Rejection")

y <- matrix(rnorm(2*(n+1))^2, 2)
x <- cumsum(y)[seq(2, 2*(n+1), 2)]
x <- x[-(n+1)] / x[n+1]
x <- x[order(y[2,-(n+1)]/y[1,-(n+1)])] 
hist(x, main="Ordered")

x <- rnorm(n) %% 1 # (Use SD of 5 or greater in practice)
hist(x, main="Modular")

x <- runif(n)      # Reference distribution
hist(x, main="Uniform")
whuber
fonte
2
(+1) If I were asking this question in an interview, I'd modify it to ask about the case where the parameters are fixed, but unknown, which strikes me as more interesting. The Pearson correlation approach (#3) goes through unchanged, but is perhaps slightly esoteric. The Beta approach (#2) requires only slight modification by considering the squares of differences of disjoint pairs. Similarly, Z=(X1X2)/(X3X4) is standard Cauchy (regardless of the mean and variance of X), which has a nice cdf.
cardinal
1
More generally, the principle is to find a pivotal quantity from the sample with a computationally amenable cdf. This ties in nicely with constructing confidence intervals and hypothesis tests, with the twist that we might seek to optimize the number of elements used rather than the latter cases which focus more on maximizing the information from a fixed sample size.
cardinal
@Cardinal Thank you for the interesting comments, as well as the ninth method (Cauchy). Even finding a pivotal quantity is unnecessary when only a good approximation is sought. For instance, (8) works perfectly well if you reserve a small number of initial results to establish a rough scale.
whuber
8

You can use a trick very similar to what you mention. Let's say that XN(μ,σ2) is a normal random variable with known parameters. Then we know its distribution function, Φμ,σ2, and Φμ,σ2(X) will be uniformly distributed on (0,1). To prove this, note that for d(0,1) we see that

P(Φμ,σ2(X)d)=P(XΦμ,σ21(d))=d.

The above probability is clearly zero for non-positive d and 1 for d1. This is enough to show that Φμ,σ2(X) has a uniform distribution on (0,1) as we have shown that the corresponding measures are equal for a generator of the Borel σ-algebra on R. Thus, you can just tranform the normally distributed data by the distribution function and you'll get uniformly distributed data.

swmo
fonte
4
It's the inverse of inverse transform sampling!
Roger Fan
Could you possibly elaborate on the second sentence of your second paragraph? I am having trouble understanding the following: "This is enough to show that Φμ,σ2(X) has a uniform distribution on (0,1) as we have shown that the corresponding measures are equal for a generator of the Borel σ-algebra on ℝ."
wellington
To show that some real random variable, X, has a uniform distribution, we should show that its corresponding measure, X(P) equals that of the uniform distribution for all measurable sets of the real line. However, it's actually enough to consider some generator of the σ-algebra, due to a uniqueness of measures-theorem. If they are equal on sets of the generator, they'll be equal for all measurable sets. This is just a measure-theoretic attachment to the answer.
swmo