NMAD: Generowanie liczb losowych z zadanego rozkładu

Autor

Julianna Kukuła

1. Wstęp

Celem tego sprawozdania jest prezentacja metod losowania liczb pseudolosowych z wybranych rozkładów prawdopodobieństwa oraz weryfikacja ich jakości.

Kod

# Wczytanie bibliotek
options(repos = c(CRAN = "[https://cloud.r-project.org/](https://cloud.r-project.org/)"))
if(!require(tidyverse)) install.packages("tidyverse")
library(tidyverse)

2. Rozkłady ciągłe

2.1. Rozkład normalny

Fundamentem wszelkich symulacji stochastycznych jest generator liczb pseudolosowych z rozkładu jednostajnego $U(0,1)$.

Klasyczny algorytm Boxa-Mullera (1958) opiera się na transformacji dwóch niezależnych zmiennych z rozkładu $U(0,1)$ do układu współrzędnych biegunowych. Po odpowiednim przekształceniu z wykorzystaniem funkcji trygonometrycznych, otrzymujemy dwie niezależne zmienne o standardowym rozkładzie normalnym $N(0,1)$.

Ze względu na wysoki koszt obliczeniowy funkcji trygonometrycznych w architekturze procesorów, zaimplementowano metodę biegunową (zaproponowaną przez G. Marsaglię w 1964 r.). Algorytm ten losuje punkty z kwadratu $[-1, 1] \times [-1, 1]$, po czym drogą eliminacji odrzuca te, które znajdują się poza kołem jednostkowym. Pozwala to zastąpić sinus i cosinus znacznie szybszymi operacjami algebraicznymi przy zachowaniu identycznego rozkładu wyjściowego.

Aby wygenerować wartości ze standardowego rozkładu normalnego $N(0,1)$ , możemy posłużyć się różnymi podejściami. Poniżej porównamy trzy z nich:

1. Wbudowaną funkcję języka R (`rnorm`).

2. Klasyczny algorytm Boxa-Mullera (wykorzystujący funkcje trygonometryczne).

3. Algorytm biegunowy Marsaglii (ulepszonie Boxa-Mullera).

Kod

n_losowan = 5000

# 1. Wbudowana funkcja R
wyniki_rnorm = rnorm(n_losowan)

# 2. Klasyczny algorytm Boxa-Mullera
box_muller = function(n) {
  u1 = runif(n)
  u2 = runif(n)
  z0 = sqrt(-2 * log(u1)) * cos(2 * pi * u2)
  return(z0) 
}
wyniki_bm = box_muller(n_losowan)

# 3. Algorytm biegunowy Marsaglii
marsaglia = function(n) {
  z = numeric(n)
  j = 0
  
  while(j < n) {
    v1 = runif(1, -1, 1)
    v2 = runif(1, -1, 1)
    s = v1^2 + v2^2
    
    if (s < 1 & s > 0) {
      j = j + 1
      mnoznik = sqrt(-2 * log(s) / s)
      z[j] = v1 * mnoznik
    }
  }
  return(z)
}
wyniki_marsaglia = marsaglia(n_losowan)

# Wizualizacja i porównanie metod
df_norm_porownanie = bind_rows(
  data.frame(Metoda = "1. Wbudowane rnorm()", Wartosc = wyniki_rnorm),
  data.frame(Metoda = "2. Box-Muller (klasyczny)", Wartosc = wyniki_bm),
  data.frame(Metoda = "3. Algorytm Marsaglii", Wartosc = wyniki_marsaglia)
)

ggplot(df_norm_porownanie, aes(x = Wartosc)) +
  geom_histogram(aes(y = after_stat(density)), bins = 40, fill = "pink", color = "white", alpha = 0.8) +
  stat_function(fun = dnorm, color = "red", linewidth = 1) +
  facet_wrap(~Metoda, ncol = 3) +
  labs(
    title = "Porównanie metod generowania rozkładu N(0,1)", 
    subtitle = "Histogramy wygenerowanych prób na tle teoretycznej gęstości (czerwona linia)", 
    x = "Wartość", 
    y = "Gęstość"
  ) +
  theme_minimal() +
  theme(strip.text = element_text(face = "bold", size = 11))

2.2. Rozkład Beta

Metoda akceptacji-odrzucenia, sformułowana przez J. von Neumanna (1951), to fundamentalne narzędzie symulacyjne pozwalające na generowanie zmiennych z rozkładów o skomplikowanej funkcji gęstości $f(x)$, dla których analityczne odwrócenie dystrybuanty jest trudne lub niemożliwe.

Generowanie wartości z rozkładu Beta polega na wylosowaniu punktu $(X, Y)$ ze zdefiniowanego obszaru prostokątnego, który ogranicza z góry naszą funkcję gęstości (z wykorzystaniem tzw. stałej majorantyzującej). Jeżeli wylosowana rzędna $Y$ spełnia warunek $Y \le f(X)$, odcięta $X$ jest akceptowana jako realizacja poszukiwanej zmiennej losowej; w przeciwnym razie wektor jest odrzucany, a proces iterowany.

Kod

# Metoda akceptacji-odrzucenia dla Beta(2,2) 
rbeta_odrzucenie = function(n) {
  wyniki = numeric(n)
  j = 0
  while(j < n) {
    x = runif(1, 0, 1)
    y = runif(1, 0, 1.5) 
    fx = 6 * x * (1 - x)
    if (y < fx) {
      j = j + 1
      wyniki[j] = x
    }
  }
  return(wyniki)
}

wyniki_beta = rbeta_odrzucenie(5000)

#Wizualizacja
df_beta = data.frame(Wartosc = wyniki_beta)

ggplot(df_beta, aes(x = Wartosc)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "orchid", color = "white", alpha=0.8) +
  stat_function(fun = function(x) dbeta(x, 2, 2), color = "red", linewidth = 1) +
  labs(title = "Rozkład Beta(2,2)", subtitle = "Metoda akceptacji-odrzucenia", x = "Wartość", y = "Gęstość") +
  theme_minimal()

2.3. Rozkład Cauchy’ego

Rozkład Cauchy’ego charakteryzuje się tzw. grubymi ogonami, co skutkuje brakiem skończonej nadzei matematycznej oraz wariancji.

Zjawisko to symulujemy na podstawie znanej własności statystycznej: iloraz dwóch niezależnych zmiennych losowych o standardowym rozkładzie normalnym

$$X, Y \sim N(0,1) \implies \frac{X}{Y} \sim \text{Cauchy}(0,1)$$

Kod

# Generowanie rozkładu Cauchy'ego przez iloraz normalnych 
generuj_cauchy = function(n) {
  x = rnorm(n)
  y = rnorm(n)
  return(x / y)
}

wyniki_cauchy = generuj_cauchy(5000)

df_cauchy = data.frame(Wartosc = wyniki_cauchy)

ggplot(df_cauchy, aes(x = Wartosc)) +
  geom_histogram(aes(y = after_stat(density)), breaks = seq(-10, 10, by = 0.5), fill = "pink", color = "white", alpha=0.8) +
  stat_function(fun = dcauchy, color = "red", linewidth = 1, xlim =c(-10,10), n = 500) +
  coord_cartesian(xlim = c(-10, 10)) +
  labs(title = "Rozkład Cauchy'ego", subtitle = "Iloraz dwóch rozkładów N(0,1)", x = "Wartość", y = "Gęstość") +
  theme_minimal()

3. Rozkłady skokowe

3.1. Rozkład Dwumianowy

Rozkład dwumianowy $B(n, p)$ reprezentuje liczbę sukcesów w $n$ próbach Bernoulliego. W implementacji przeprowadzamy symulację poprzez wygenerowanie wektora n zmiennych z rozkładu $U(0,1)$. Każda wartość mniejsza niż zadane prawdopodobieństwo p jest traktowana jako sukces. Suma tych zdarzeń w jednej iteracji stanowi pojedynczą realizację zmiennej z rozkładu dwumianowego.

Kod

# Generowanie rozkładu dwumianowego (n=10, p=0.5) 
generuj_dwumianowy = function(N_prob, n, p) {
  wyniki = numeric(N_prob)
  for(i in 1:N_prob) {
    u = runif(n)
    sukcesy = sum(u < p)
    wyniki[i] = sukcesy
  }
  return(wyniki)
}

wyniki_binom = generuj_dwumianowy(5000, 10, 0.5)

3.2. Rozkład Poissona

Rozkład Poissona modeluje liczbę zdarzeń w ustalonym przedziale czasu. Polega na mnożeniu niezależnych zmiennych jednostajcnych $U_{i}$ i sprawdzaniu kiedy spadnie poniżej $e^-\lambda$

Kod

generuj_poissona = function(N_prob, lambda) {
  wyniki = numeric(N_prob)
  L = exp(-lambda)
  for(i in 1:N_prob) {
    k = 0
    p = 1
    repeat {
      k = k + 1
      p = p * runif(1)
      if (p <= L) break
    }
    wyniki[i] = k - 1
  }
  return(wyniki)
}

wyniki_pois = generuj_poissona(5000, 4)

Porównanie rozkładów skokowych

Kod

df_skokowe = bind_rows(
  data.frame(Rozklad = "Dwumianowy B(10, 0.5)", Wartosc = wyniki_binom),
  data.frame(Rozklad = "Poissona (lambda = 4)", Wartosc = wyniki_pois)
)

ggplot(df_skokowe, aes(x = Wartosc, fill = Rozklad)) +
  geom_bar(position = position_dodge(width = 0.8), color = "black", alpha = 0.8) +
  scale_fill_manual(values = c("pink", "purple")) +
  
  scale_x_continuous(breaks = seq(0, max(df_skokowe$Wartosc), by = 1)) +
  
  labs(
    title = "Generowanie zmiennych skokowych", 
    subtitle = "Zestawienie prób z rozkładu dwumianowego i Poissona",
    x = "Wylosowana wartość (k)", 
    y = "Liczebność"
  ) +
  theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.title = element_blank(),   
    panel.grid.minor = element_blank()
  )

3.3. Dowolny zadany rozkład dyskretny (Metoda odwracania dystrybuanty)

Generowanie zmiennych z arbitralnie zadanego rozkładu dyskretnego opiera się na adaptacji metody odwracania dystrybuanty. Posiadając zbiór wartości $X = \{x_1, x_2, \dots, x_k\}$ oraz odpowiadający im wektor prawdopodobieństw $P = \{p_1, p_2, \dots, p_k\}$, wyznaczamy dystrybuantę (skumulowane sumy prawdopodobieństw). Następnie losujemy zmienną $U \sim U(0,1)$ i przypisujemy jej wartość $x_i$, dla której $U$ trafia w wyznaczony przedział.

Zbudujmy własny, abstrakcyjny rozkład. Załóżmy, że nasza zmienna przyjmuje wartości: 10, 20, 30 i 40 z prawdopodobieństwami odpowiednio: 0.1, 0.4, 0.3 oraz 0.2.

Kod

generuj_zadany_skokowy = function(n, wartosci, prawdopodobienstwa) { 
  if(abs(sum(prawdopodobienstwa) - 1) > 1e-6) { stop("Prawdopodobieństwa muszą sumować się do 1!") }

dystrybuanta = cumsum(prawdopodobienstwa) 
wyniki = numeric(n)

for(i in 1:n) { 
  u = runif(1) 
  indeks = which(u <= dystrybuanta)[1] 
  wyniki[i] = wartosci[indeks] }

return(wyniki) }
# Parametry naszego zadanego rozkładu
moje_wartosci = c(10, 20, 30, 40)
moje_prawd = c(0.1, 0.4, 0.3, 0.2)

wyniki_zadany = generuj_zadany_skokowy(5000, moje_wartosci, moje_prawd)

#Wizualizacja
# Przygotowanie danych do wykresu
df_zadany = data.frame(Wartosc = factor(wyniki_zadany, levels = moje_wartosci))

df_teoria = data.frame(
  Wartosc = factor(moje_wartosci),
  Prawdopodobienstwo_teoretyczne = moje_prawd * 5000 # Przeskalowane do liczby losowań
)

ggplot(df_zadany, aes(x = Wartosc)) +
  geom_bar(fill = "pink", color = "black", alpha = 0.8, width = 0.6) +
  geom_point(data = df_teoria, aes(x = Wartosc, y = Prawdopodobienstwo_teoretyczne), 
             color = "red", size = 4, shape = 18) +
  labs(
    title = "Dowolnie zadany rozkład skokowy",
    subtitle = "Słupki: empiryczne losowania | Czerwone romby: teoretyczne oczekiwania",
    x = "Zadana wartość (x)",
    y = "Liczebność (liczba wystąpień)"
  ) +
  theme_minimal() +
  theme(
    panel.grid.minor = element_blank()
  )

4. Metoda ROU

Metoda Ratio of Uniforms (zaproponowana w 1977 r. przez A. Kindermana i J. Monahana) stanowi zaawansowaną alternatywę dla metody eliminacji. Opiera się na matematycznym twierdzeniu orzekającym, że obszar pod krzywą gęstości można przekształcić w przestrzeń dwuwymiarową.

Algorytm polega na jednostajnym próbkowaniu punktów $(U, V)$ z wielokąta krzywoliniowego (tzw. obszaru akceptacji) zdefiniowanego nierównością $0 \le U \le \sqrt{f(V/U)}$. Jeśli wylosowany punkt leży wewnątrz tego obszaru, to iloraz $X = V/U$ ma rozkład o funkcji gęstości ściśle proporcjonalnej do $f(x)$. Istotną zaletą analityczną tej metody jest brak konieczności znajomości stałej normalizującej (całki z funkcji gęstości).

Kod

generuj_ROU_normalny = function(n) {
  wyniki = numeric(n)
  v_max = sqrt(2/exp(1))
  v_min = -v_max
  j = 0
  
  while(j < n) {
    u = runif(1, 0, 1)
    v = runif(1, v_min, v_max)
    x = v / u
    
    # Warunek akceptacji
    if (u^2 <= exp(-x^2 / 2)) {
      j = j + 1
      wyniki[j] = x
    }
  }
  return(wyniki)
}

wyniki_rou = generuj_ROU_normalny(5000)

df_rou = data.frame(Wartosc = wyniki_rou)

ggplot(df_rou, aes(x = Wartosc)) +
  geom_histogram(aes(y = after_stat(density)), bins = 40, fill = "orchid", color = "white", alpha = 0.8) +
  stat_function(fun = dnorm, color = "black", linewidth = 1, linetype="dashed") +
  labs(title = "Metoda ROU", subtitle = "Generowanie rozkładu N(0,1)", x = "Wartość", y = "Gęstość") +
  theme_minimal()

5. Ocena jakości generatora

5.1. Widmo generatora

Warto zbadać zachowanie punktów parami w postaci wykresu rozrzutu $(U_i, U_{i+1})$, zwanego widmem generatora. Dobry generator wypełnia przestrzeń losowo i bez wyraźnych wzorów.

Kod

# 1. Definiujemy celowo BARDZO SŁABY generator (Liniowy Generator Kongruencyjny)
zly_generator = function(n, seed = 123) {
  u = numeric(n)
  m = 256
  a = 137
  c = 0
  x = seed
  for(i in 1:n) {
    x = (a * x + c) %% m
    u[i] = x / m
  }
  return(u)
}

n_prob = 1500
dane_zle = zly_generator(n_prob) 
dane_r = runif(n_prob)           # Dobry generator (Wbudowany w R)

df_widmo = bind_rows(
  data.frame(X = dane_zle[1:(n_prob-1)], Y = dane_zle[2:n_prob], Generator = "Słaby generator (LCG)"),
  data.frame(X = dane_r[1:(n_prob-1)], Y = dane_r[2:n_prob], Generator = "Dobry generator (runif)")
)

#Wizualizacja
ggplot(df_widmo, aes(x = X, y = Y)) +
  geom_point(alpha = 0.5, size = 1, color = "purple") +
  facet_wrap(~Generator) +
  labs(
    title = "Widmo generatora (Wykres rozrzutu)", 
    subtitle = "U[i] względem U[i+1] - wizualny test niezależności",
    x = "Wartość U[i]", 
    y = "Wartość U[i+1]"
  ) +
  theme_minimal() +
  theme(
    strip.text = element_text(face = "bold", size = 12),
    panel.border = element_rect(color = "gray", fill = NA, linewidth = 0.5) 
  )

Wykres prawy (Słaby LCG): Punkty nie są losowe, tylko układają się w siatkę/prążki. Oznacza to, że znając poprzednią liczbę, możemy w miarę dokładnie przewidzieć następną. To dyskwalifikuje taki generator.
Wykres lewy (runif): Punkty tworzą równomierną “chmurę” (szum), bez pustych przestrzeni i bez wzorów. Wskazuje to na brak autokorelacji i dobrą jakość generatora.

6. Bibliografia

Takashi Shinzato “Box Muller Method” https://d1wqtxts1xzle7.cloudfront.net/37993877/boxmuller-libre.pdf?1435166544=&response-content-disposition=inline%3B+filename%3DBox_Muller_Method.pdf&Expires=1775907394&Signature=YPJOy6ryDjl8qEPyffxonGc7dptVwYGzPRLFllcmeqxhAWkt17r4vffAAi5UnLHW
Marsaglia George “Random Number Generators” https://web.archive.org/web/20180720073146id_/https://digitalcommons.wayne.edu/cgi/viewcontent.cgi?article=1725&context=jmasm
https://mst.mimuw.edu.pl/lecture.php?lecture=sst&part=Ch3
https://dl.acm.org/doi/pdf/10.1145/76738.76801
https://www.ceeol.com/search/article-detail?id=885076https://www.ceeol.com/search/article-detail?id=885076 strony 6,8
https://bibliotekanauki.pl/articles/212494.pdf - rozdział 5.1

--- title: "NMAD: Generowanie liczb losowych z zadanego rozkładu" author: "Julianna Kukuła" format: html: theme: cosmo toc: true toc-title: "Nawigacja" code-fold: show code-tools: true smooth-scroll: true number-sections: false lang: pl editor: visual --- ## 1. Wstęp Celem tego sprawozdania jest prezentacja metod losowania liczb pseudolosowych z wybranych rozkładów prawdopodobieństwa oraz weryfikacja ich jakości. ```{r} #| label: setup #| include: false set.seed(47) ``` ```{r} #| message: false #| warning: false # Wczytanie bibliotek options(repos = c(CRAN = "[https://cloud.r-project.org/](https://cloud.r-project.org/)")) if(!require(tidyverse)) install.packages("tidyverse") library(tidyverse) ``` ## 2. Rozkłady ciągłe ### 2.1. Rozkład normalny Fundamentem wszelkich symulacji stochastycznych jest generator liczb pseudolosowych z rozkładu jednostajnego $U(0,1)$. Klasyczny algorytm Boxa-Mullera (1958) opiera się na transformacji dwóch niezależnych zmiennych z rozkładu $U(0,1)$ do układu współrzędnych biegunowych. Po odpowiednim przekształceniu z wykorzystaniem funkcji trygonometrycznych, otrzymujemy dwie niezależne zmienne o standardowym rozkładzie normalnym $N(0,1)$. Ze względu na wysoki koszt obliczeniowy funkcji trygonometrycznych w architekturze procesorów, zaimplementowano metodę biegunową (zaproponowaną przez G. Marsaglię w 1964 r.). Algorytm ten losuje punkty z kwadratu $[-1, 1] \times [-1, 1]$, po czym drogą eliminacji odrzuca te, które znajdują się poza kołem jednostkowym. Pozwala to zastąpić sinus i cosinus znacznie szybszymi operacjami algebraicznymi przy zachowaniu identycznego rozkładu wyjściowego. Aby wygenerować wartości ze standardowego rozkładu normalnego $N(0,1)$ , możemy posłużyć się różnymi podejściami. Poniżej porównamy trzy z nich: 1\. Wbudowaną funkcję języka R (\`rnorm\`). 2\. Klasyczny algorytm Boxa-Mullera (wykorzystujący funkcje trygonometryczne). 3\. Algorytm biegunowy Marsaglii (ulepszonie Boxa-Mullera). ```{r} #| code-fold: show n_losowan = 5000 # 1. Wbudowana funkcja R wyniki_rnorm = rnorm(n_losowan) # 2. Klasyczny algorytm Boxa-Mullera box_muller = function(n) { u1 = runif(n) u2 = runif(n) z0 = sqrt(-2 * log(u1)) * cos(2 * pi * u2) return(z0) } wyniki_bm = box_muller(n_losowan) # 3. Algorytm biegunowy Marsaglii marsaglia = function(n) { z = numeric(n) j = 0 while(j < n) { v1 = runif(1, -1, 1) v2 = runif(1, -1, 1) s = v1^2 + v2^2 if (s < 1 & s > 0) { j = j + 1 mnoznik = sqrt(-2 * log(s) / s) z[j] = v1 * mnoznik } } return(z) } wyniki_marsaglia = marsaglia(n_losowan) # Wizualizacja i porównanie metod df_norm_porownanie = bind_rows( data.frame(Metoda = "1. Wbudowane rnorm()", Wartosc = wyniki_rnorm), data.frame(Metoda = "2. Box-Muller (klasyczny)", Wartosc = wyniki_bm), data.frame(Metoda = "3. Algorytm Marsaglii", Wartosc = wyniki_marsaglia) ) ggplot(df_norm_porownanie, aes(x = Wartosc)) + geom_histogram(aes(y = after_stat(density)), bins = 40, fill = "pink", color = "white", alpha = 0.8) + stat_function(fun = dnorm, color = "red", linewidth = 1) + facet_wrap(~Metoda, ncol = 3) + labs( title = "Porównanie metod generowania rozkładu N(0,1)", subtitle = "Histogramy wygenerowanych prób na tle teoretycznej gęstości (czerwona linia)", x = "Wartość", y = "Gęstość" ) + theme_minimal() + theme(strip.text = element_text(face = "bold", size = 11)) ``` ### 2.2. Rozkład Beta Metoda akceptacji-odrzucenia, sformułowana przez J. von Neumanna (1951), to fundamentalne narzędzie symulacyjne pozwalające na generowanie zmiennych z rozkładów o skomplikowanej funkcji gęstości $f(x)$, dla których analityczne odwrócenie dystrybuanty jest trudne lub niemożliwe. Generowanie wartości z rozkładu Beta polega na wylosowaniu punktu $(X, Y)$ ze zdefiniowanego obszaru prostokątnego, który ogranicza z góry naszą funkcję gęstości (z wykorzystaniem tzw. stałej majorantyzującej). Jeżeli wylosowana rzędna $Y$ spełnia warunek $Y \le f(X)$, odcięta $X$ jest akceptowana jako realizacja poszukiwanej zmiennej losowej; w przeciwnym razie wektor jest odrzucany, a proces iterowany. ```{r} #| code-fold: show # Metoda akceptacji-odrzucenia dla Beta(2,2) rbeta_odrzucenie = function(n) { wyniki = numeric(n) j = 0 while(j < n) { x = runif(1, 0, 1) y = runif(1, 0, 1.5) fx = 6 * x * (1 - x) if (y < fx) { j = j + 1 wyniki[j] = x } } return(wyniki) } wyniki_beta = rbeta_odrzucenie(5000) #Wizualizacja df_beta = data.frame(Wartosc = wyniki_beta) ggplot(df_beta, aes(x = Wartosc)) + geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "orchid", color = "white", alpha=0.8) + stat_function(fun = function(x) dbeta(x, 2, 2), color = "red", linewidth = 1) + labs(title = "Rozkład Beta(2,2)", subtitle = "Metoda akceptacji-odrzucenia", x = "Wartość", y = "Gęstość") + theme_minimal() ``` ### 2.3. Rozkład Cauchy'ego Rozkład Cauchy'ego charakteryzuje się tzw. grubymi ogonami, co skutkuje brakiem skończonej nadzei matematycznej oraz wariancji. Zjawisko to symulujemy na podstawie znanej własności statystycznej: iloraz dwóch niezależnych zmiennych losowych o standardowym rozkładzie normalnym \$\$X, Y \\sim N(0,1) \\implies \\frac{X}{Y} \\sim \\text{Cauchy}(0,1)\$\$ ```{r} #| code-fold: show # Generowanie rozkładu Cauchy'ego przez iloraz normalnych generuj_cauchy = function(n) { x = rnorm(n) y = rnorm(n) return(x / y) } wyniki_cauchy = generuj_cauchy(5000) df_cauchy = data.frame(Wartosc = wyniki_cauchy) ggplot(df_cauchy, aes(x = Wartosc)) + geom_histogram(aes(y = after_stat(density)), breaks = seq(-10, 10, by = 0.5), fill = "pink", color = "white", alpha=0.8) + stat_function(fun = dcauchy, color = "red", linewidth = 1, xlim =c(-10,10), n = 500) + coord_cartesian(xlim = c(-10, 10)) + labs(title = "Rozkład Cauchy'ego", subtitle = "Iloraz dwóch rozkładów N(0,1)", x = "Wartość", y = "Gęstość") + theme_minimal() ``` ## 3. Rozkłady skokowe ### 3.1. Rozkład Dwumianowy Rozkład dwumianowy $B(n, p)$ reprezentuje liczbę sukcesów w $n$ próbach Bernoulliego. W implementacji przeprowadzamy symulację poprzez wygenerowanie wektora n zmiennych z rozkładu $U(0,1)$. Każda wartość mniejsza niż zadane prawdopodobieństwo p jest traktowana jako sukces. Suma tych zdarzeń w jednej iteracji stanowi pojedynczą realizację zmiennej z rozkładu dwumianowego. ```{r} #| code-fold: show # Generowanie rozkładu dwumianowego (n=10, p=0.5) generuj_dwumianowy = function(N_prob, n, p) { wyniki = numeric(N_prob) for(i in 1:N_prob) { u = runif(n) sukcesy = sum(u < p) wyniki[i] = sukcesy } return(wyniki) } wyniki_binom = generuj_dwumianowy(5000, 10, 0.5) ``` ### 3.2. Rozkład Poissona Rozkład Poissona modeluje liczbę zdarzeń w ustalonym przedziale czasu. Polega na mnożeniu niezależnych zmiennych jednostajcnych $U_{i}$ i sprawdzaniu kiedy spadnie poniżej $e^-\lambda$ ```{r} #| code-fold: show generuj_poissona = function(N_prob, lambda) { wyniki = numeric(N_prob) L = exp(-lambda) for(i in 1:N_prob) { k = 0 p = 1 repeat { k = k + 1 p = p * runif(1) if (p <= L) break } wyniki[i] = k - 1 } return(wyniki) } wyniki_pois = generuj_poissona(5000, 4) ``` Porównanie rozkładów skokowych ```{r} #| code-fold: show df_skokowe = bind_rows( data.frame(Rozklad = "Dwumianowy B(10, 0.5)", Wartosc = wyniki_binom), data.frame(Rozklad = "Poissona (lambda = 4)", Wartosc = wyniki_pois) ) ggplot(df_skokowe, aes(x = Wartosc, fill = Rozklad)) + geom_bar(position = position_dodge(width = 0.8), color = "black", alpha = 0.8) + scale_fill_manual(values = c("pink", "purple")) + scale_x_continuous(breaks = seq(0, max(df_skokowe$Wartosc), by = 1)) + labs( title = "Generowanie zmiennych skokowych", subtitle = "Zestawienie prób z rozkładu dwumianowego i Poissona", x = "Wylosowana wartość (k)", y = "Liczebność" ) + theme_minimal() + theme( legend.position = "bottom", legend.title = element_blank(), panel.grid.minor = element_blank() ) ``` ### 3.3. Dowolny zadany rozkład dyskretny (Metoda odwracania dystrybuanty) Generowanie zmiennych z arbitralnie zadanego rozkładu dyskretnego opiera się na adaptacji metody odwracania dystrybuanty. Posiadając zbiór wartości $X = \{x_1, x_2, \dots, x_k\}$ oraz odpowiadający im wektor prawdopodobieństw $P = \{p_1, p_2, \dots, p_k\}$, wyznaczamy dystrybuantę (skumulowane sumy prawdopodobieństw). Następnie losujemy zmienną $U \sim U(0,1)$ i przypisujemy jej wartość $x_i$, dla której $U$ trafia w wyznaczony przedział. Zbudujmy własny, abstrakcyjny rozkład. Załóżmy, że nasza zmienna przyjmuje wartości: 10, 20, 30 i 40 z prawdopodobieństwami odpowiednio: 0.1, 0.4, 0.3 oraz 0.2. ```{r} #| code-fold: show generuj_zadany_skokowy = function(n, wartosci, prawdopodobienstwa) { if(abs(sum(prawdopodobienstwa) - 1) > 1e-6) { stop("Prawdopodobieństwa muszą sumować się do 1!") } dystrybuanta = cumsum(prawdopodobienstwa) wyniki = numeric(n) for(i in 1:n) { u = runif(1) indeks = which(u <= dystrybuanta)[1] wyniki[i] = wartosci[indeks] } return(wyniki) } # Parametry naszego zadanego rozkładu moje_wartosci = c(10, 20, 30, 40) moje_prawd = c(0.1, 0.4, 0.3, 0.2) wyniki_zadany = generuj_zadany_skokowy(5000, moje_wartosci, moje_prawd) #Wizualizacja # Przygotowanie danych do wykresu df_zadany = data.frame(Wartosc = factor(wyniki_zadany, levels = moje_wartosci)) df_teoria = data.frame( Wartosc = factor(moje_wartosci), Prawdopodobienstwo_teoretyczne = moje_prawd * 5000 # Przeskalowane do liczby losowań ) ggplot(df_zadany, aes(x = Wartosc)) + geom_bar(fill = "pink", color = "black", alpha = 0.8, width = 0.6) + geom_point(data = df_teoria, aes(x = Wartosc, y = Prawdopodobienstwo_teoretyczne), color = "red", size = 4, shape = 18) + labs( title = "Dowolnie zadany rozkład skokowy", subtitle = "Słupki: empiryczne losowania | Czerwone romby: teoretyczne oczekiwania", x = "Zadana wartość (x)", y = "Liczebność (liczba wystąpień)" ) + theme_minimal() + theme( panel.grid.minor = element_blank() ) ``` ## 4. Metoda ROU Metoda Ratio of Uniforms (zaproponowana w 1977 r. przez A. Kindermana i J. Monahana) stanowi zaawansowaną alternatywę dla metody eliminacji. Opiera się na matematycznym twierdzeniu orzekającym, że obszar pod krzywą gęstości można przekształcić w przestrzeń dwuwymiarową. Algorytm polega na jednostajnym próbkowaniu punktów $(U, V)$ z wielokąta krzywoliniowego (tzw. obszaru akceptacji) zdefiniowanego nierównością $0 \le U \le \sqrt{f(V/U)}$. Jeśli wylosowany punkt leży wewnątrz tego obszaru, to iloraz $X = V/U$ ma rozkład o funkcji gęstości ściśle proporcjonalnej do $f(x)$. Istotną zaletą analityczną tej metody jest brak konieczności znajomości stałej normalizującej (całki z funkcji gęstości). ```{r} #| code-fold: show generuj_ROU_normalny = function(n) { wyniki = numeric(n) v_max = sqrt(2/exp(1)) v_min = -v_max j = 0 while(j < n) { u = runif(1, 0, 1) v = runif(1, v_min, v_max) x = v / u # Warunek akceptacji if (u^2 <= exp(-x^2 / 2)) { j = j + 1 wyniki[j] = x } } return(wyniki) } wyniki_rou = generuj_ROU_normalny(5000) df_rou = data.frame(Wartosc = wyniki_rou) ggplot(df_rou, aes(x = Wartosc)) + geom_histogram(aes(y = after_stat(density)), bins = 40, fill = "orchid", color = "white", alpha = 0.8) + stat_function(fun = dnorm, color = "black", linewidth = 1, linetype="dashed") + labs(title = "Metoda ROU", subtitle = "Generowanie rozkładu N(0,1)", x = "Wartość", y = "Gęstość") + theme_minimal() ``` ## 5. Ocena jakości generatora ### 5.1. Widmo generatora Warto zbadać zachowanie punktów parami w postaci wykresu rozrzutu $(U_i, U_{i+1})$, zwanego widmem generatora. Dobry generator wypełnia przestrzeń losowo i bez wyraźnych wzorów. ```{r} #| code-fold: show # 1. Definiujemy celowo BARDZO SŁABY generator (Liniowy Generator Kongruencyjny) zly_generator = function(n, seed = 123) { u = numeric(n) m = 256 a = 137 c = 0 x = seed for(i in 1:n) { x = (a * x + c) %% m u[i] = x / m } return(u) } n_prob = 1500 dane_zle = zly_generator(n_prob) dane_r = runif(n_prob) # Dobry generator (Wbudowany w R) df_widmo = bind_rows( data.frame(X = dane_zle[1:(n_prob-1)], Y = dane_zle[2:n_prob], Generator = "Słaby generator (LCG)"), data.frame(X = dane_r[1:(n_prob-1)], Y = dane_r[2:n_prob], Generator = "Dobry generator (runif)") ) #Wizualizacja ggplot(df_widmo, aes(x = X, y = Y)) + geom_point(alpha = 0.5, size = 1, color = "purple") + facet_wrap(~Generator) + labs( title = "Widmo generatora (Wykres rozrzutu)", subtitle = "U[i] względem U[i+1] - wizualny test niezależności", x = "Wartość U[i]", y = "Wartość U[i+1]" ) + theme_minimal() + theme( strip.text = element_text(face = "bold", size = 12), panel.border = element_rect(color = "gray", fill = NA, linewidth = 0.5) ) ``` - **Wykres prawy (Słaby LCG):** Punkty nie są losowe, tylko układają się w siatkę/prążki. Oznacza to, że znając poprzednią liczbę, możemy w miarę dokładnie przewidzieć następną. To dyskwalifikuje taki generator. - **Wykres lewy (runif):** Punkty tworzą równomierną "chmurę" (szum), bez pustych przestrzeni i bez wzorów. Wskazuje to na brak autokorelacji i dobrą jakość generatora. ## 6. Bibliografia 1. Takashi Shinzato "Box Muller Method" <https://d1wqtxts1xzle7.cloudfront.net/37993877/boxmuller-libre.pdf?1435166544=&response-content-disposition=inline%3B+filename%3DBox_Muller_Method.pdf&Expires=1775907394&Signature=YPJOy6ryDjl8qEPyffxonGc7dptVwYGzPRLFllcmeqxhAWkt17r4vffAAi5UnLHW> 2. Marsaglia George "Random Number Generators" <https://web.archive.org/web/20180720073146id_/https://digitalcommons.wayne.edu/cgi/viewcontent.cgi?article=1725&context=jmasm> 3. <https://mst.mimuw.edu.pl/lecture.php?lecture=sst&part=Ch3> 4. <https://dl.acm.org/doi/pdf/10.1145/76738.76801> 5. [https://www.ceeol.com/search/article-detail?id=885076https://www.ceeol.com/search/article-detail?id=885076](https://www.ceeol.com/search/article-detail?id=885076) strony 6,8 6. <https://bibliotekanauki.pl/articles/212494.pdf> - rozdział 5.1