Symulowanie losowań z rozkładu jednolitego przy użyciu losowań z rozkładu normalnego

Niedawno kupiłem zasób wywiadu danych, w którym jedno z pytań prawdopodobieństwa było następujące:

Biorąc pod uwagę losowania z rozkładu normalnego o znanych parametrach, jak można symulować losowania z rozkładu jednolitego?

Mój pierwotny proces myślowy polegał na tym, że dla dyskretnej zmiennej losowej możemy rozbić rozkład normalny na K unikalnych podsekcji, gdzie każda podsekcja ma równe pole pod krzywą normalną. Następnie możemy określić, które z wartości K przyjmuje zmienna, rozpoznając, w którym obszarze krzywej normalnej spada zmienna.

Ale działałoby to tylko w przypadku dyskretnych zmiennych losowych. Przeprowadziłem badania nad tym, jak moglibyśmy zrobić to samo dla ciągłych zmiennych losowych, ale niestety mogłem znaleźć tylko techniki takie jak próbkowanie z transformacją odwrotną, które wykorzystywałyby jako dane wejściowe jednolitą zmienną losową i mogłyby generować zmienne losowe z innego rozkładu. Myślałem, że być może moglibyśmy wykonać ten proces odwrotnie, aby uzyskać jednolite zmienne losowe?

Zastanawiałem się również nad wykorzystaniem normalnych zmiennych losowych jako danych wejściowych do liniowego generatora kongruencjalnego, ale nie jestem pewien, czy to zadziała.

Masz jakieś przemyślenia na temat tego, jak mogę podejść do tego pytania?

W duchu stosowania prostych obliczeń algebraicznych niezwiązanych z obliczeniami rozkładu normalnego skłaniam się ku następującemu. Są one uporządkowane tak, jak o nich myślałem (i dlatego potrzebowałem, aby uzyskać coraz więcej kreatywności), ale zachowałem to, co najlepsze - i najbardziej zaskakujące - na trwałe.

1. Odwróć technikę Boxa-Muellera : z każdej pary normalnych można zbudować dwa niezależne mundury jako atan2 ( Y , X ) (w przedziale [ - π , π ] ) i exp ( - ( X 2 + Y 2 ) / 2 ) (w przedziale [ 0 , 1 ] ).$(X,Y)$$\text{atan2}(Y,X)$$[-\pi, \pi]$$\exp(-(X^2+Y^2)/2)$$[0,1]$

2. Podjąć normalne w grupach po dwa, a suma ich kwadratów, aby uzyskać sekwencję zmiennych towarzyszących Y 1 , Y 2 , ... , Y i , ... . Wyrażenia uzyskane z par$\chi^2_2$$Y_1, Y_2, \ldots, Y_i, \ldots$

   $$X_i = \frac{Y_{2i}}{Y_{2i-1}+Y_{2i}}$$

   będzie miał rozkład , który jest jednolity.$\text{Beta}(1,1)$

   Że wymaga to tylko podstawowej, prostej arytmetyki, powinno być jasne.

3. Ponieważ dokładny rozkład współczynnika korelacji Pearsona dla czteroparowej próbki ze standardowej dwuwymiarowej Rozkład normalny jest równomiernie rozłożony na , możemy po prostu wziąć normalne w grupach po cztery pary (to znaczy, osiem wartości w każdego zestawu) i zwraca współczynnik korelacji tych par. (Obejmuje to prostą arytmetykę plus dwie operacje pierwiastkowe).$[-1,1]$

4. Od czasów starożytnych wiadomo, że cylindryczny rzut kuli (powierzchnia w trzech przestrzeniach) ma równą powierzchnię . Oznacza to, że w rzucie równomiernego rozkładu na kulę zarówno współrzędna pozioma (odpowiadająca długości geograficznej), jak i współrzędna pionowa (odpowiadająca szerokości geograficznej) będzie miała równomierne rozkłady. Ponieważ trójwartościowy standardowy rozkład normalny jest sferycznie symetryczny, jego rzut na kulę jest jednorodny. Uzyskanie długości geograficznej jest zasadniczo tym samym obliczeniem co kąt w metodzie Boxa-Muellera ( qv ), ale rzutowana szerokość geograficzna jest nowa. Rzut na sferę normalizuje jedynie potrójną współrzędną iw tym punkcie z jest rzutowaną szerokością geograficzną. Zatem weź zmienne Normalne w grupach trzyosobowych, X 3 i - 2 , X 3 i - 1 , X 3 i , i oblicz$(x,y,z)$$z$$X_{3i-2}, X_{3i-1}, X_{3i}$

   $$\frac{X_{3i}}{\sqrt{X_{3i-2}^2 + X_{3i-1}^2 + X_{3i}^2}}$$

   dla .$i=1, 2, 3, \ldots$

5. Ponieważ większość systemów komputerowych reprezentacji liczb w formacie binarnym , jednolity numer generacji zazwyczaj zaczyna się od produkcji równomiernie rozłożonych liczb całkowitych pomiędzy i 2 32 - 1 (lub wysokiej mocy 2 związane z długością słowa komputer) i przeskalowanie je w razie potrzeby. Takie liczby całkowite są reprezentowane wewnętrznie jako ciągi 32 cyfr binarnych. Możemy uzyskać niezależne losowe bity poprzez porównanie zmiennej Normal z jej medianą. Wystarczy zatem podzielić zmienne normalne na grupy wielkości równe pożądanej liczbie bitów, porównać każdą ze średnią i złożyć wynikowe sekwencje wyników prawda / fałsz w liczbę binarną. Pisanie $0$$2^{32}-1$$2$$32$$k$dla liczby bitów i dla znaku (to znaczy H ( x ) = 1 gdy x > 0 i H ( x ) = 0 w przeciwnym razie) możemy wyrazić wynikową znormalizowaną jednolitą wartość w [ 0 , 1 ) za pomocą wzoru$H$$H(x)=1$$x\gt 0$$H(x)=0$$[0, 1)$

   $$\sum_{j=0}^{k-1} H(X_{ki - j})2^{-j-1}.$$

   Zmienne można wyciągnąć z dowolnego ciągłego rozkładu, którego mediana wynosi 0 (na przykład standardowa Normalna); są przetwarzane w grupach k, przy czym każda grupa wytwarza jedną taką pseudo-jednolitą wartość.$X_n$$0$$k$

6. Próbkowanie przy odrzuceniu jest standardowym, elastycznym i potężnym sposobem rysowania zmiennych losowych z dowolnych rozkładów. Załóżmy, że rozkład docelowy ma format PDF . Wartość Y jest rysowana zgodnie z innym rozkładem w formacie PDF g . W etapie odrzucania rysowana jest jednolita wartość U leżąca między 0 a g ( Y ) niezależnie od Y i porównywana z f ( Y ) : jeśli jest mniejsza, $f$$Y$$g$$U$$0$$g(Y)$$Y$$f(Y)$$Y$zostaje zachowany, ale w przeciwnym razie proces się powtarza. To podejście wydaje się jednak okrągłe: w jaki sposób generujemy zmienną jednolitą z procesem, który wymaga na początku zmiennego jednolitego?

   Odpowiedź jest taka, że ​​tak naprawdę nie potrzebujemy jednolitej wariacji, aby wykonać krok odrzucenia. Zamiast tego (zakładając ), możemy rzucić uczciwą monetą, aby uzyskać losowo 0 lub 1 . Będzie to interpretowane jako pierwszy bit w binarnej reprezentacji zmiennej zmiennej U w przedziale [ 0 , 1 ) . Gdy wynik jest 0 , to znaczy 0 ≤ U < 1 / 2 ; w przeciwnym razie, 1 / 2 ≤ U < 1 . $g(Y)\ne 0$$0$$1$$U$$[0,1)$$0$$0 \le U \lt 1/2$$1/2\le U \lt 1$Połowę czasu, to wystarczy, aby zdecydować, etap odrzucenia: jeśli a moneta jest 0 , Y należy przyjąć; jeśli f ( T ) / g ( Y ) < 1 / 2 a moneta jest 1 , Y należy odrzucić; w przeciwnym razie, musimy ponownie odwrócić monetę w celu uzyskania następny bit U . Ponieważ - bez względu na wartość f ( $f(Y)/g(Y) \ge 1/2$$0$$Y$$f(Y)/g(Y) \lt 1/2$$1$$Y$$U$ posiada - nie jest 1 / 2 szansa zatrzymania po każdej klapki, oczekiwana liczba koziołki jest tylko 1 / 2 ( 1 ) + 1 / 4 ( 2 ) + 1 / 8 ( 3 ) + ⋯ + 2 - n ( n ) + ⋯ = 2 .$f(Y)/g(Y)$$1/2$$1/2(1)+1/4(2)+1/8(3)+\cdots+2^{-n}(n)+\cdots=2$

   Pobieranie próbek odrzucenia może być opłacalne (i wydajne), pod warunkiem, że spodziewana liczba odrzuceń jest niewielka. Możemy to osiągnąć, dopasowując największy możliwy prostokąt (reprezentujący jednolity rozkład) pod normalnym plikiem PDF.

   

   Korzystanie Rachunek zoptymalizować obszar prostokąta, można zauważyć, że jego punkty końcowe powinny leżeć w , gdzie jego wysokość jest równa exp ( - 1 / 2 ) / $\pm 1$, co czyni jego powierzchnię nieco większą niż0,48. Dzięki zastosowaniu tego standardowego gęstość zwykle Jakg, odrzucając wszystkie wartości spoza przedziału[-1,1],automatycznie lub w inny sposób, stosując procedurę odrzucenia, otrzymamy w jednolity zmiennych towarzyszących[-1,1]skutecznie:$\exp(-1/2)/\sqrt{2\pi}\approx 0.241971$$0.48$$g$$[-1,1]$$[-1,1]$

   • W ułamku czasu, Zmienna normalna leży poza [ - 1 , 1 ] i jest natychmiast odrzucana. ( Φ jest standardowym normalnym CDF.)$2\Phi(-1) \approx 0.317$$[-1,1]$$\Phi$

   • W pozostałym ułamku czasu należy zastosować procedurę odrzucenia binarnego, wymagającą średnio dwóch dodatkowych zmiennych normalnych.

   • Ogólna procedura wymaga średnio kroki.$1/(2\exp(-1/2)/\sqrt{2\pi}) \approx 2.07$

   Oczekiwana liczba zmiennych normalnych potrzebnych do uzyskania każdego jednolitego wyniku działa

   $$\sqrt{2 e \pi}\left(1-2\Phi(-1)\right) \approx 2.82137.$$

   Chociaż jest to dość wydajne, należy pamiętać, że (1) obliczenie normalnego PDF wymaga obliczenia wykładniczego i (2) wartość musi zostać obliczona raz na zawsze. To wciąż trochę mniej obliczeń niż metoda Boxa-Muellera ( qv ).$\Phi(-1)$

7. Te statystyki zamówień o równomiernym rozkładzie mają wykładnicze luk. Ponieważ suma kwadratów dwóch normalnych (średniej zerowej) jest wykładnicza, możemy wygenerować realizację niezależnych mundurów przez zsumowanie kwadratów par takich normalnych, obliczenie skumulowanej ich sumy, przeskalowanie wyników w celu uzyskania przedziału [ 0 , 1 ] i upuszczenie ostatniego (który zawsze będzie równy 1 ). Jest to przyjemne podejście, ponieważ wymaga tylko podniesienia do kwadratu, zsumowania i (na końcu) pojedynczego podziału.$n$$[0,1]$$1$

   Wartości będą automatycznie rosły w porządku rosnącym. Jeśli takie sortowanie jest pożądane, metoda ta jest lepsza pod względem obliczeniowym od wszystkich innych, o ile pozwala uniknąć kosztu O ( n log ( n ) ) . Jeśli jednak potrzebna jest sekwencja niezależnych mundurów, losowe posortowanie tych n wartości załatwi sprawę. Ponieważ (jak widać w metodzie Boxa-Muellera, qv ) stosunki każdej pary normalnych są niezależne od sumy kwadratów każdej pary, mamy już środki do uzyskania tej losowej permutacji: uporządkuj skumulowane sumy według odpowiednich stosunków . (Jeśli $n$$O(n\log(n))$$n$$n$jest bardzo duży, proces ten można przeprowadzić w mniejszych grupach przy niewielkiej utracie wydajności, ponieważ każda grupa potrzebuje tylko 2 ( k + 1 ) normałów, aby utworzyć k jednolitych wartości. W przypadku ustalonego k asymptotyczny koszt obliczeniowy wynosi zatem O ( n log ( k ) ) = O ( n ) , wymagając 2 n ( 1 + 1 / k ) Zmiennych normalnych w celu wygenerowania n jednolitych wartości.)$k$$2(k+1)$$k$$k$$O(n\log(k))$$O(n)$$2n(1+1/k)$$n$

8. W doskonałym przybliżeniu każda zmienna Normal z dużym odchyleniem standardowym wygląda jednolicie w zakresach o wiele mniejszych wartości. Po zwinięciu tego rozkładu do zakresu (biorąc tylko ułamkowe części wartości), uzyskujemy w ten sposób rozkład, który jest jednolity dla wszystkich praktycznych celów. Jest to niezwykle wydajne, wymagające jednej z najprostszych operacji arytmetycznych ze wszystkich: wystarczy zaokrąglić każdą wartość normalną do najbliższej liczby całkowitej i zatrzymać nadmiar. Prostota tego podejścia staje się przekonująca, gdy przeanalizujemy praktyczne wdrożenie:$[0,1]$R

   rnorm(n, sd=10) %% 1
   

   niezawodnie wytwarza njednolite wartości w przedziale kosztem tylko zmiennych normalnych i prawie żadnych obliczeń.$[0,1]$n

   (Nawet jeśli odchylenie standardowe wynosi , plik PDF z tym przybliżeniem różni się od jednolitego pliku PDF, jak pokazano na poniższym rysunku, o mniej niż jedną część na 10 8 ! Aby to niezawodnie wykryć wymagałoby próbki 10 16 wartości - to już nie jest w stanie wykonać żadnego standardowego testu losowości. Przy większym odchyleniu standardowym niejednorodność jest tak mała, że ​​nie można jej nawet obliczyć. Na przykład przy SD 10 jak pokazano w kodzie, maksymalne odchylenie od jednolitości PDF ma tylko 10 - 857 ).$1$$10^8$$10^{16}$$10$$10^{-857}$

   

W każdym przypadku zmienne normalne „ze znanymi parametrami” można łatwo zaktualizować i przeskalować do przyjętych wyżej norm normalnych. Następnie powstałe równomiernie rozłożone wartości można zaktualizować i przeskalować w celu uwzględnienia dowolnego pożądanego przedziału. Wymagają one jedynie podstawowych operacji arytmetycznych.

Łatwość tych konstrukcji dowodzi następujący Rkod, który używa tylko jednej lub dwóch linii dla większości z nich. Ich poprawności, świadczy otrzymanego prawie jednolite histogramy na podstawie niezależne wartości w każdym przypadku (wymaga około 12 sekund dla wszystkich siedmiu symulacji). Dla odniesienia - w przypadku, gdy martwisz się ilością zmian pojawiających się na którejkolwiek z tych wykresów - na końcu dołączony jest histogram wartości jednolitych symulowanych za pomocą generatora jednolitych liczb losowych.$100,000$R

$\chi^2$$1000$$3\%$R

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")

You can use a trick very similar to what you mention. Let's say that $X \sim N(\mu, \sigma^2)$ is a normal random variable with known parameters. Then we know its distribution function, $\Phi_{\mu,\sigma^2}$, and $\Phi_{\mu,\sigma^2}(X)$ will be uniformly distributed on $(0,1)$. To prove this, note that for $d \in (0,1)$ we see that

$P(\Phi_{\mu,\sigma^2}(X) \leq d) = P(X \leq \Phi_{\mu,\sigma^2}^{-1}(d)) = d$.

The above probability is clearly zero for non-positive $d$ and $1$ for $d \geq 1$. This is enough to show that $\Phi_{\mu,\sigma^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 $\sigma$-algebra on $\mathbb{R}$. Thus, you can just tranform the normally distributed data by the distribution function and you'll get uniformly distributed data.