Generowanie Zmiennych Losowych

Metody, wybrane rozkłady i ocena jakości generatorów

Author

Dominika Grudzień

1.Rozkłady Ciągłe

1.Rozkład Normalny (CTG)

Metoda Superpozycji:

Zgodnie z Centralnym Twierdzeniem Granicznym (CTG), suma wielu losowych zdarzeń tworzy dzwon Gaussa. Jeśli zsumujemy dokładnie \(k=12\) liczb z płaskiego rozkładu jednostajnego \(U(0,1)\) i odejmiemy od wyniku 6, otrzymamy idealny standardowy rozkład normalny \(N(0,1)\).

Code 
n <- 10000

x_ctg <- replicate(n, sum(runif(12)) - 6)

hist(x_ctg, breaks=40, probability=TRUE, col="lightpink", 
     main="Metoda CTG (Suma 12 zmiennych U(0,1))")
curve(dnorm(x), add=TRUE, col="darkred", lwd=2)

2.Rozkład Normalny (Box-Muller)

Algorytm Dokładny:

Ponieważ CTG obcina ogony rozkładu, w praktyce stosuje się metodę biegunową Boxa-Muller’a. Losujemy współrzędne na płaszczyźnie, odrzucamy te poza kołem jednostkowym, a resztę transformujemy wzorem do płaszczyzny 3D.

Code 
#| layout-ncol: 2
#| fig-height: 3.5
box_muller_polar <- function(n) {
  x <- numeric(n); j <- 0
  for (i in 1:n) {
    repeat {
      v1 <- 2 * runif(1) - 1; v2 <- 2 * runif(1) - 1
      s <- v1^2 + v2^2
      if (s > 0 & s < 1) { 
        j <- j + 1; x[j] <- v1 * sqrt(-2 * log(s) / s); break 
      }
    }
  }
  return(x)
}

hist(box_muller_polar(2000), probability=TRUE, col="lightgreen", main="Box-Muller")
curve(dnorm(x), add=TRUE, col="red", lwd=2)

Code 
hist(rnorm(2000), probability=TRUE, col="lightgrey", main="Wbudowane rnorm()")
curve(dnorm(x), add=TRUE, col="red", lwd=2)

3.Rozkład Beta

Metoda Akceptacji-Odrzucenia:

Rozkład \(Beta(2,2)\) ma kształt góry uwięzionej między 0 a 1. Tworzone jest więc pudełko o wysokości 1.5 (maksimum gęstości), losujemy w nim punkty i zatrzymujemy tylko te, które wpadły pod krzywą funkcji \(f(x) = 6x(1-x)\).

Code 
rbeta_custom <- function(n) {
  res <- numeric(n); i <- 1
  while(i <= n) {
     x <- runif(1); y <- runif(1, 0, 1.5)      
     if (y < 6 * x * (1 - x)) { res[i] <- x; i <- i + 1 }
  }
  return(res)
}

hist(rbeta_custom(2000), probability=TRUE, col="pink", main="Akceptacja-Odrzucenie")
curve(dbeta(x, 2, 2), add=TRUE, col="purple", lwd=2)
hist(rbeta(2000, 2, 2), probability=TRUE, col="plum", main="Wbudowane rbeta()")
curve(dbeta(x, 2, 2), add=TRUE, col="purple", lwd=2)

4.Rozkład Cauchy’ego

Odwracanie Dystrybuanty:

Rozkład Cauchy’ego słynie z tzw. grubych ogonów (brak skończonej średniej). Najprostszą metodą jego wygenerowania jest klasyczne odwrócenie dystrybuanty. Podstawiamy losowy ułamek pod funkcję tangens’a: \(X = \tan(\pi(U - 0.5))\).

Code 
u <- runif(5000)
x_cauchy <- tan(pi * (u - 0.5))

hist(x_cauchy[x_cauchy > -10 & x_cauchy < 10], breaks=50, probability=TRUE, 
     col="salmon", main="Rozkład Cauchy'ego (Odwracanie Dystrybuanty)")
curve(dcauchy(x), add=TRUE, lwd=2, col="darkred")

2.Rozkłady Skokowe

1.Rozkład Dwumianowy:

Rozkład \(Binomial(n, p)\) mówi, ile sukcesów osiągniemy w \(n\) próbach. Losujemy \(n\) razy ułamek \(U(0,1)\). Każdy ułamek mniejszy niż nasze zadane \(p\) traktujemy jako wygraną i zliczamy je.

Code 
rbinom_custom <- function(symulacje, n, p) {
  wyniki <- numeric(symulacje)
  for(i in 1:symulacje) {
    wyniki[i] <- sum(runif(n) < p)
  }
  return(wyniki)
}

wlasne <- rbinom_custom(3000, n=10, p=0.5)
wbudowane <- rbinom(3000, size=10, prob=0.5)

barplot(table(wlasne)/3000, col="darkseagreen", main="Dwumianowy (n=10)")
barplot(table(wbudowane)/3000, col="darkolivegreen", main="Wbudowane rbinom()")

2. Rozkład Poissona

Opisuje rzadkie zjawiska (np. ilość awarii danej maszyny). Zgodnie z algorytmem Knutha, losujemy liczby jednostajne i następnie mnożymy je przez siebie tak długo, aż ich iloczyn spadnie poniżej progu \(e^{-\lambda}\). Ilość mnożeń to wynik z rozkładu Poissona.

Code 
rpois_custom <- function(n, lambda) {
  x <- numeric(n); L <- exp(-lambda)
  for(i in 1:n) {
    p <- 1; k <- 0
    while(p >= L) { k <- k + 1; p <- p * runif(1) }
    x[i] <- k - 1
  }
  return(x)
}

barplot(table(rpois_custom(3000, 4))/3000, col="khaki", main="Poisson (lambda=4)")
barplot(table(rpois(3000, 4))/3000, col="gold", main="Wbudowane rpois()")

3.Dowolny zadany rozkład

Koło fortuny:

Dla niestandardowej tabeli prawdopodobieństw (np. zdarzenie A: 20%, B: 50%, C: 30%) budowane są skumulowane przedziały. Losujemy ułamek \(U(0,1)\) i sprawdzane jest, w który przedział wpadł.

Code 
rdowolny <- function(n, wartosci, prawdopodobienstwa) {
  dystrybuanta <- cumsum(prawdopodobienstwa)
  u <- runif(n)
  sapply(u, function(x) wartosci[which.max(dystrybuanta >= x)])
}

w <- c("A (20%)", "B (50%)", "C (30%)"); p <- c(0.2, 0.5, 0.3)
wyniki <- rdowolny(5000, w, p)

barplot(table(wyniki)/5000, col="cyan", ylim=c(0, 0.6), main="Empiryczne vs Teoretyczne")
abline(h = p, col="red", lty=2, lwd=2)

3. Metoda ROU

Ratio of Uniforms (ROU):

Metoda ROU (iloraz zmiennych jednostajnych) to uniwersalna technika wprowadzona w 1977 r. Pozwala ona generować rozkłady ciągłe poprzez dzielenie przez siebie dwóch punktów \(V\) i \(U\). Dzięki temu, że \(U\) potrafi być bardzo bliskie zera, metoda ta świetnie radzi sobie z generowaniem ekstremalnych wartości (grubych ogonów). Poniżej implementacja ROU do wygenerowania rozkładu Cauchy’ego z obszaru półkola: \(U^2 + V^2 \le 1\).

Code 
rou_cauchy <- function(n) {
  x <- numeric(n); i <- 1
  while(i <= n) {
    u <- runif(1, 0, 1)      
    v <- runif(1, -1, 1)     
    if (u^2 + v^2 <= 1) { x[i] <- v / u; i <- i + 1 } 
  }
  return(x)
}

c_vals <- rou_cauchy(5000)
hist(c_vals[c_vals > -10 & c_vals < 10], breaks=50, probability=TRUE, 
     col="mediumpurple1", main="Rozkład Cauchy'ego zbudowany metodą ROU")

4. Ocena Jakości Generatora

1.Zły generator, który oszukuje:

Algorytmy komputerowe udają losowość. Budujemy więc historyczny, wadliwy algorytm (Liniowy Generator Kongruencyjny LCG), który dla niepoznaki nazywano RANDU. Gdy przeprowadzimy z nim porównanie z funkcjami R za pomocą histogramów (1D), stary algorytm wyda się idealny – słupki ułożą się równomiernie i płasko.

Code 
zly_generator <- function(n, seed=1) {
  x <- numeric(n); x[1] <- seed
  for(i in 2:n) x[i] <- (65539 * x[i-1]) %% (2^31)
  return(x / (2^31)) 
}

u_bad <- zly_generator(4000)
u_good <- runif(4000)

hist(u_bad, col="tomato", main="Zły Generator (który wydaje się dobry)", probability=TRUE)
hist(u_good, col="skyblue", main="Dobry Generator z R", probability=TRUE)

2.Widmo Generatora (Test Spektralny):

Skoro w jednym wymiarze błąd może nie być widoczny, musimy narysować punkty w 2D (kolejne losowania stają się współrzędnymi zarówno X jak i Y). Tutaj możemy zobaczyć, że wadliwy generator układa punkty w wyraźne “prążki” (zależności). Dobry generator wbudowany w R zachowuje się jak szum telewizyjny – tworzy chaos na całej płaszczyźnie, udowadniając swoją jakość.

Code 
plot(u_bad[-length(u_bad)], u_bad[-1], pch=20, cex=0.3, col="darkred", 
     main="Widmo słabego RANDU", xlab=expression(U[i]), ylab=expression(U[i+1]))
plot(u_good[-length(u_good)], u_good[-1], pch=20, cex=0.3, col="darkblue", 
     main="Widmo runif() z R", xlab=expression(U[i]), ylab=expression(U[i+1]))

Bibliografia:

1.“Non-Uniform Random Variate Generation”- Devroye, L. (1986)

https://link.springer.com/book/10.1007/978-1-4613-8643-8

  • Rozkład Normalny (Box-Muller i CTG): Rozdział IX, s. 380–381.

  • Rozkład Beta (Akceptacja-Odrzucenie): Rozdział II, s. 41–44, Rozdział IX, s. 428–430.

  • Odwracanie dystrybuanty (Rozkład Caucheg’o, Rozkłady Skokowe): Rozdział III, s. 85–87.

  • Metoda ROU: Rozdział VII, s. 194–197.

2.“The Art of Computer Programming, Volume 2: Seminumerical Algorithms” (3rd ed.)- Knuth, D. E. (1997)

https://seriouscomputerist.atariverse.com/media/pdf/book/Art%20of%20Computer%20Programming%20-%20Volume%202%20(Seminumerical%20Algorithms).pdf

  • Rozkład Poissona: Rozdział 3.4.1.C, s. 137.

  • Widmo generatora i wady RANDU (The Spectral Test): Rozdział 3.3.4, s. 93–104.

3.“A note on the generation of random normal deviates. The Annals of Mathematical Statistics” -Box, G. E. P., & Muller, M. E. (1958), 29(2), s. 610-611

https://projecteuclid.org/journals/annals-of-mathematical-statistics/volume-29/issue-2/A-Note-on-the-Generation-of-Random-Normal-Deviates/10.1214/aoms/1177706645.pdf

4. “Computer generation of random variables using the ratio of uniform deviates. ACM Transactions on Mathematical Software” -Kinderman, A. J. & Monahan, J. F. (1977), 3(3)

https://dl.acm.org/doi/epdf/10.1145/355744.355750

*ten link działa tylko przez sieć uczelnianą “Eduroam”

5.“Komputerowe generatory liczb losowych”- Wieczorowski & Zieliński (1997)

https://dokumen.pub/komputerowe-generatory-liczb-losowych.html

  • ocena jakości generatorów i test spektralny (widmo), s.45-55