NMAD

Author

Aneta Drabik

Published

March 26, 2026

Zaawansowane Symulacje i Generatory w R

1. Szacowanie liczby \(\pi\) (Metoda Monte Carlo)

Metoda ta wykorzystuje geometryczną definicję prawdopodobieństwa. Losujemy punkty w kwadracie i sprawdzamy, jaki ich odsetek mieści się w ćwiartce koła. Jest to klasyczny przykład wykorzystania naśladowania losu do rozwiązywania problemów deterministycznych.

Code 
library(ggplot2)
Warning: pakiet 'ggplot2' został zbudowany w wersji R 4.5.2
Code 
set.seed(2026)
n <- 5000
df_pi <- data.frame(x = runif(n), y = runif(n))
df_pi$inside <- ifelse(df_pi$x^2 + df_pi$y^2 <= 1, "Wewnątrz", "Zewnątrz")

pi_est <- 4 * sum(df_pi$inside == "Wewnątrz") / n

ggplot(df_pi, aes(x, y, color = inside)) +
  geom_point(alpha = 0.5, size = 0.8) +
  annotate("path", x = cos(seq(0, pi/2, length.out = 1000)), 
           y = sin(seq(0, pi/2, length.out = 1000)), color = "black", size = 1) +
  labs(title = paste("oszacowane Pi =", pi_est), subtitle = paste("liczba prób:", n)) +
  coord_fixed() + theme_minimal()
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

Wizualizacja pokazuje rozsyp punktów trafiających w tarczę. Stosunek punktów niebieskich (wewnątrz koła) do wszystkich wylosowanych pozwala nam wyznaczyć stałą \(\pi\) z dużą dokładnością.

2. Generatory wartości z rozkładu U[0, 1]

Najczęściej stosowaną metodą komputerowego naśladowania losu jest Liniowy Generator Kongruencyjny (LCG). Autorzy podręcznika wskazują, że wybór odpowiednich parametrów (mnożnika i modułu) jest kluczowy dla jakości uzyskanego nieładu.

Code 
lcg_gen <- function(n, seed = 123) {
  a <- 16807; m <- 2^31 - 1; c <- 0
  res <- numeric(n); curr <- seed
  for(i in 1:n) {
    curr <- (a * curr + c) %% m
    res[i] <- curr / m
  }
  return(res)
}

u_data <- data.frame(val = lcg_gen(2000))

ggplot(u_data, aes(x = val)) +
  geom_histogram(bins = 20, fill = "blue", color = "yellow") +
  labs(title = "histogram Generatora", x = "wrtość", y = "częstość") +
  theme_light()

Histogram potwierdza, że algorytm LCG generuje liczby rozłożone równomiernie na odcinku (0, 1). Choć ciąg jest czysto deterministyczny, spełnia on statystyczne testy jakości.

3. Metoda odwracania dystrybuanty

Jeśli dysponujemy rozkładem jednostajnym, możemy go przekształcić w dowolny inny rozkład klasyczny. Metoda ta polega na wyliczeniu wartości funkcji odwrotnej do dystrybuanty \(F^{-1}(u)\).

Code 
n <- 2000
lambda <- 1.5
u <- runif(n)
x_exp <- -log(1 - u) / lambda
df_exp <- data.frame(x = x_exp)

ggplot(df_exp, aes(x)) +
  geom_histogram(aes(y = ..density..), bins = 25, fill = "seagreen", alpha = 0.6) +
  stat_function(fun = dexp, args = list(rate = lambda), color = "purple", size = 1) +
  labs(title = "rozkład wykładniczy: Metoda Odwracania", x = "x", y = "gęstość") +
  theme_minimal()
Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
ℹ Please use `after_stat(density)` instead.

Histogram wygenerowanych danych pokrywa się z fioletową linią gęstości teoretycznej, co dowodzi poprawności matematycznego przekształcenia rozkładu U(0,1).

4. Metoda eliminacji

Stosowana, gdy notacja ułatwiająca implementację bezpośrednią dystrybuanty jest zbyt trudna. Losujemy punkty pod prostą ograniczającą i „wycinamy” te, które nie mieszczą się pod funkcją gęstości celu.

Code 
f <- function(x) dbeta(x, 2, 5)
M <- 2.6
n_try <- 3210
df_rej <- data.frame(x = runif(n_try), y = runif(n_try, 0, M))
df_rej$accept <- df_rej$y <= f(df_rej$x)

ggplot(df_rej, aes(x, y, color = accept)) +
  geom_point(alpha = 0.4) +
  stat_function(fun = f, color = "black", size = 1) +
  scale_color_manual(values = c("red", "purple")) +
  labs(title = "Metoda eliminacji: Rozkład Beta(2, 5)") +
  theme_bw()

Punkty fioletowe to te, które zostały zaakceptowane. Widzimy, że tworzą one kształt gęstości Beta pod czarną krzywą, podczas gdy czerwone punkty zostały odrzucone.

5. Metoda odrzucania dla rozkładów warunkowych

Pozwala generować wartości z „uciętych” fragmentów populacji. Jest to technika szczególnie przydatna przy modelowaniu zjawisk ograniczonych, np. fizycznie niemożliwych wartości ujemnych.

Code 
n_war <- 1234
res_war <- c()
while(length(res_war) < n_war) {
  val <- rnorm(1, 0, 1)
  if(val > 1) res_war <- c(res_war, val)
}

ggplot(data.frame(x = res_war), aes(x)) +
  geom_density(fill = "pink", alpha = 0.4) +
  labs(title = "Ucięty rozkład normalny (X > 1)") +
  theme_minimal()

Gęstość zaczyna się gwałtownie od wartości 1. Wszystkie losowania, które nie spełniały warunku, zostały odrzucone, pozostawiając jedynie „ogon” rozkładu normalnego.

6. Metoda superpozycji rozkładów

Metoda polega na składaniu skomplikowanej gęstości z prostszych elementów (mieszanek). Najpierw wybieramy komponent (np. rzutem monetą), a potem losujemy z niego wartość.

Code 
n_mix <- 2000
waga <- 0.4
u_mix <- runif(n_mix)
mix_data <- ifelse(u_mix < waga, rnorm(n_mix, -2, 1), rnorm(n_mix, 2, 1.5))

ggplot(data.frame(x = mix_data), aes(x)) +
  geom_histogram(bins = 30, fill = "gold", color = "brown", alpha = 0.7) +
  geom_density(size = 1) +
  labs(title = "Superpozycja dwóch rozkładów normalnych") +
  theme_minimal()

Wykres pokazuje dwa wyraźne „garby” (mody). Jest to efekt nałożenia na siebie dwóch różnych procesów losowych w jedną populację.

7. Spacer Pijaka (Błądzenie Losowe)

W kontekście komputerowego naśladowania losu, o którym piszą Wieczorkowski i Zieliński, spacer pijaka jest fascynującym testem dla generatorów. Choć algorytm jest deterministyczny (tzw. jawna grzeszność według von Neumanna), to przy dużej liczbie kroków generuje strukturę, która statystycznie zachowuje się jak procesy naturalne, np. ruchy cząsteczek w cieczy.

Sam termin został po raz pierwszy użyty przez Karla Pearsona w 1905 roku na łamach czasopisma „Nature”. Pearson szukał rozwiązania matematycznego dla problemu „pijaka”, który po każdym kroku zapomina, skąd przyszedł i wybiera kolejny kierunek całkowicie losowo. Nazwa ta przyjęła się w nauce, ponieważ idealnie obrazuje proces, w którym każdy kolejny ruch jest niezależny od poprzedniego, co prowadzi do chaotycznej i nieprzewidywalnej ścieżki.

Code 
n_krokow <- 10000

# Symulacja błądzenia: losujemy kąt (0-360 stopni) i długość kroku
katy <- runif(n_krokow, 0, 2*pi)
dx <- cos(katy)
dy <- sin(katy)

# Tworzymy ścieżkę jako sumę skumulowaną przesunięć
df_pijak <- data.frame(
  x = cumsum(dx),
  y = cumsum(dy),
  czas = 1:n_krokow
)

ggplot(df_pijak, aes(x, y, color = czas)) +
  geom_path(alpha = 0.7) +
  scale_color_viridis_c(option = "magma") +
  labs(title = "Spacer Pijaka: Symulacja błądzenia na płaszczyźnie",
       subtitle = "Każdy krok ma losowy kierunek, ale stałą długość",
       x = "Pozycja X", y = "Pozycja Y") +
  theme_minimal()

BIBLIOGRAFIA

  1. Ryszard Magiera, MODELE I METODY STATYSTYKI MATEMATYCZNEJ Część II WNIOSKOWANIE STATYSTYCZNE, Wydanie III, rozszerzone (rodział 5, strony 221-231)

  2. Robert, C. P., Casella, G. (2010). Introducing Monte Carlo Methods with R. Springer.

  3. Szekli, R. (2021). Rachunek Prawdopodobieństwa: Plan wykładu i skrypt. Wrocław: Instytut Matematyczny Uniwersytetu Wrocławskiego (strony 34-59)