Author

Wojciech Potrawa

Szacowanie liczby \(\pi\)

###Strzały do tarczy

Szacowanie wartości liczby \(\pi\) metodą Monte Carlo opiera się na geometrycznym prawdopodobieństwie trafienia losowego punktu w określony obszar. Proces ten wykorzystuje generowanie punktów z rozkładu jednostajnego.

Rozważmy kwadrat o boku \(2r\), w który wpisane jest koło o promieniu \(r\). Pole kwadratu wynosi \(P_{kwadratu} = (2r)^2 = 4r^2\), natomiast pole koła to \(P_{kola} = \pi r^2\). Stosunek pola koła do pola kwadratu wyraża się wzorem:

\[\frac{P_{kola}}{P_{kwadratu}} = \frac{\pi r^2}{4r^2} = \frac{\pi}{4}\]

Jeżeli wygenerujemy \(N\) niezależnych punktów o współrzędnych \((x, y)\) z rozkładu jednostajnego \(U(-r, r)\), to prawdopodobieństwo, że punkt znajdzie się wewnątrz koła (czyli spełni warunek \(x^2 + y^2 \le r^2\)), jest równe stosunkowi tych pól. Oznaczając liczbę punktów, które trafiły do wnętrza koła jako \(M\), otrzymujemy przybliżenie:

\[\pi \approx 4 \frac{M}{N}\]

Poniższy kod w języku R przeprowadza taką symulację dla \(r=1\) i \(N=10000\) iteracji, a następnie wizualizuje wynik:

Code 
N <- 10000

x <- runif(N, min = -1, max = 1)
y <- runif(N, min = -1, max = 1)

w_kole <- (x^2 + y^2) <= 1

M <- sum(w_kole)
pi_est <- 4 * M / N

cat("Oszacowana wartość pi dla N =", N, "wynosi:", pi_est, "\n")
Oszacowana wartość pi dla N = 10000 wynosi: 3.166
Code 
cat("Błąd bezwzględny:", abs(pi - pi_est), "\n")
Błąd bezwzględny: 0.02440735
Code 
plot(x, y, col = ifelse(w_kole, "steelblue", "red"), pch = 20, cex = 0.5,
     main = paste("Szacowanie liczby pi (Wynik:", pi_est, ")"),
     xlab = "x", ylab = "y", asp = 1)

theta <- seq(0, 2*pi, length.out = 100)
lines(cos(theta), sin(theta), col = "black", lwd = 2)

Wizualizacja szacowania liczby pi metodą Monte Carlo

Eksperyment Igły Buffona

Innym, klasycznym przykładem estymacji liczby \(\pi\) z wykorzystaniem symulacji jest tzw. problem igły Buffona. Zakłada on rzucanie igły o długości \(l\) na płaszczyznę pokrytą równoległymi liniami oddalonymi od siebie o odległość \(d\) (przy czym \(l \le d\)).

Zgodnie z klasycznymi dowodami, prawdopodobieństwo tego, że igła przetnie którąkolwiek z linii, wyraża się wzorem: \[P = \frac{2l}{\pi d}\] Przekształcając ten wzór i estymując prawdopodobieństwo na podstawie wyników symulacji komputerowej (stosunek rzutów przecinających linie do wszystkich rzutów), możemy oszacować wartość \(\pi\). Poniżej znajduje się implementacja tego eksperymentu:

Code 
set.seed(123)
N <- 100000   
l <- 1        
d <- 2        

alfa <- runif(N) * pi
x <- runif(N) * d

przeciecia <- sum((x + 0.5 * l * sin(alfa) > d) | (x - 0.5 * l * sin(alfa) < 0))

pi_buffon <- 2 * l * N / (d * przeciecia)

cat("Oszacowana wartość pi (Igła Buffona) dla N =", N, "wynosi:", pi_buffon, "\n")
Oszacowana wartość pi (Igła Buffona) dla N = 1e+05 wynosi: 3.132538

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

Współczesne symulacje numeryczne opierają się na programowych generatorach liczb pseudolosowych. Zgodnie z literaturą przedmiotu (T. Chwiej, Generatory liczb pseudolosowych, AGH), w przeciwieństwie do generatorów fizycznych (wykorzystujących np. zjawisko promieniotwórczości czy szumy elektroniczne), algorytmy komputerowe charakteryzują się powtarzalnością generowanych ciągów oraz dużą łatwością obsługi.

Jak w praktyce pracują takie algorytmy? W pierwszym etapie tworzony jest deterministyczny ciąg liczb naturalnych \(X_0, X_1, \dots, X_n\), które są ograniczone od góry przez określoną reprezentację maszynową \(m\) (np. \(m = 2^{32}\) dla 32-bitowej liczby całkowitej). Następnie, aby uzyskać pożądane wartości rzeczywiste mieszczące się w przedziale \((0, 1]\), stosuje się prostą transformację:

\[U = \frac{X}{m} \implies U_i \in (0, 1]\]

Warto w tym miejscu podkreślić fundamentalną różnicę między prawdziwą losowością a tą generowaną przez maszyny. Jak zauważają R. Wieczorkowski i R. Zieliński w swojej klasycznej monografii, ciągiem ściśle losowym jest taki ciąg, którego nie da się zapisać za pomocą algorytmu krótszego niż on sam. Komputery są urządzeniami deterministycznymi, dlatego generują tzw. liczby pseudolosowe.

Konsekwencją tego jest fakt, że źle zaprojektowany generator (np. o nieodpowiednio dobranych parametrach mnożnika) wyprodukuje ciąg silnie schematyczny. Tego typu “niepoprawne” generatory charakteryzują się bardzo krótkimi cyklami lub ukrytymi regularnościami – wygenerowane z nich liczby, po naniesieniu na wykres wielowymiarowy, zaczynają układać się w przewidywalne równoległe proste lub hiperpłaszczyzny, co całkowicie dyskwalifikuje je z poważnych symulacji Monte Carlo.

Metoda środków kwadratów

Dobrym przykładem obrazującym samą ideę generowania takich ciągów z pewnej wartości początkowej (ziarna) jest klasyczna metoda środków kwadratów. Polega ona na podnoszeniu liczby do potęgi drugiej, wyodrębnianiu jej środkowych cyfr, a następnie dzieleniu przez odpowiednią potęgę dziesiątki (pełniącą tu rolę naszego \(m\)), co sprowadza wynik do przedziału \((0, 1)\).

Code 
x <- 8234372336
y <- c()

for (i in 1:10) {
  x <- x^2
  # Zabezpieczenie przed notacją naukową w R
  x_str <- format(x, scientific = FALSE)
  x_str <- substr(x_str, 6, 15)
  x <- as.numeric(x_str) 
  y <- c(y, x / 10^10) 
}

print(y)
 [1] 0.88776788 0.18123108 0.70308935 0.46323960 0.09277493 0.88434369
 [7] 0.37675257 0.25024677 0.44729681 0.44405302

Liniowe generatory

Znacznie bardziej sformalizowaną i szeroko opisaną klasą algorytmów są generatory liniowe. Opierają się one na operacji dzielenia modulo (wyznaczania reszty z dzielenia), tworząc ciąg według rekurencyjnego schematu:

\[X_{n+1} = (a_1 X_n + a_2 X_{n-1} + \dots + a_k X_{n-k+1} + c) \pmod m\]

Gdzie parametry generatora (\(a_i, c, m\)) są z góry ustalonymi liczbami całkowitymi. Najprostszą i niezwykle popularną odmianą jest tzw. generator multiplikatywny, w którym stała \(c = 0\), a wzór redukuje się do postaci:

\[X_{i+1} = a X_{i-1} \pmod m\]

W tego typu algorytmach kluczowy jest dobór odpowiednich parametrów, ponieważ to od nich zależy tzw. okres generatora, czyli moment, w którym sekwencja liczb zacznie się powtarzać. Maksymalny możliwy okres dla generatora multiplikatywnego wynosi \(m-1\).

Poniżej przedstawiono przykład wygenerowania dużej próby takich liczb w środowisku R (przy użyciu domyślnych funkcji optymalizowanych):

Code 
set.seed(123)
u_vals <- runif(100000)

hist(u_vals, main="Histogram rozkładu jednostajnego", 
     xlab="Wartość z przedziału (0, 1)", ylab="Częstość", col="steelblue")

Warto również wspomnieć, że suma niezależnych zmiennych losowych o rozkładzie jednostajnym \(U(0, 1)\) tworzy rozkład, który zgodnie z Centralnym Twierdzeniem Granicznym (CTG) dąży do rozkładu normalnego. Należy przy tym pamiętać, że zbieżność ta następuje wraz ze wzrostem liczby sumowanych zmiennych losowych. W szczególnym przypadku, gdy sumujemy dokładnie dwie niezależne zmienne o rozkładzie \(U(0, 1)\), otrzymujemy precyzyjny rozkład trójkątny określony na przedziale \((0, 2)\):

Code 
set.seed(42)
hist(runif(100000) + runif(100000), 
     main = "Histogram sumy dwóch rozkładów jednostajnych",
     xlab = "Wartość z przedziału (0, 2)", 
     ylab = "Częstość", 
     col = "darkgreen")

Generatory wartości z zadanego rozkładu

Dysponując wartościami losowymi z rozkładu jednostajnego \(U(0,1)\), możemy wykorzystać odpowiednie metody transformacji, aby otrzymać wartości losowe z dowolnych innych rozkładów. Można tego dokonać na 4 różne metody:

Metoda odwracania dystrybuanty

Metoda ta opiera się na własnościach dystrybuanty \(F(x)\), która zawsze przyjmuje wartości z przedziału \([0, 1]\). Algorytm jest prosty: losujemy wartość \(U\) z rozkładu jednostajnego \(U(0, 1)\), a następnie przyrównujemy ją do dystrybuanty docelowego rozkładu i wyliczamy \(x\). Matematycznie sprowadza się to do użycia funkcji odwrotnej (tzw. funkcji kwantylowej): \[x = F^{-1}(U)\]

W środowisku R funkcja odwrotna do dystrybuanty jest często wbudowana jako np. qnorm. Poniżej przykład ręcznego odwrócenia dystrybuanty dla rozkładu wykładniczego, gdzie \(F(x) = 1 - e^{-x}\), co po przekształceniu daje wzór \(x = -\ln(1-U)\).

Code 
set.seed(123)
n <- 10000
U <- runif(n)
X <- -log(1 - U) 

hist(X, prob = TRUE, breaks = 40, col = "lightyellow", 
     main = "Rozkład wykładniczy (Odwracanie dystrybuanty)", xlab="x")

Metoda eliminacji (akceptacji-odrzucenia)

W momencie gdy znalezienie funkcji odwrotnej do dystrybuanty jest trudne lub niemożliwe, można skorzystać z metody eliminacji. Wymaga ona znajomości funkcji gęstości \(f(x)\) oraz jej maksimum na rozpatrywanym przedziale.

Doskonałym przykładem zastosowania tej metody jest generowanie wartości z rozkładu Beta o parametrach (2,2). Funkcja gęstości tego rozkładu na przedziale \([0, 1]\) ma postać \(f(x) = 6x(1-x)\), a jej maksymalna wartość (wierzchołek paraboli) wynosi \(1.5\) (czyli \(3/2\)).

Algorytm losuje punkt \((x, y)\) w prostokącie ograniczającym funkcję i akceptuje go tylko wtedy, gdy punkt “trafi” pod wykres gęstości. W przeciwnym razie losowanie powtarzane jest aż do skutku:

Code 
set.seed(123)
n <- 1000
wygenerowane_beta <- numeric(n)
licznik <- 1

while (licznik <= n) {
  x <- runif(1, 0, 1)
  y <- runif(1, 0, 3/2)      
  fx <- 6 * x * (1 - x)      
  
  if (y < fx) {
    wygenerowane_beta[licznik] <- x  
    licznik <- licznik + 1           
  }
}

hist(wygenerowane_beta, prob=TRUE, col="lightblue", 
     main="Metoda eliminacji (Rozkład Beta)", xlab="x", ylab="Gęstość")
curve(6*x*(1-x), add=TRUE, col="red", lwd=2)

Metoda odrzucania dla rozkładów warunkowych

Metoda ta jest specyficzną modyfikacją klasycznej eliminacji. Znajduje zastosowanie w sytuacji, gdy celem jest wygenerowanie zmiennej z tzw. rozkładu uciętego. Proces polega na losowaniu wartości z pełnego rozkładu bazowego, a następnie filtrowaniu wyników – odrzucane są wszystkie te, które nie spełniają zadanego logicznego warunku brzegowego.

Code 
set.seed(123)
warunkowe_x <- numeric(0)

while(length(warunkowe_x) < 2000) {
  kandydaci <- rnorm(500)
  zaakceptowani <- kandydaci[kandydaci > 0]
  warunkowe_x <- c(warunkowe_x, zaakceptowani)
}

warunkowe_x <- warunkowe_x[1:2000] 
hist(warunkowe_x, prob=TRUE, breaks=30, col="lavender",
     main="Warunkowy rozkład normalny (ucinka x > 0)", xlab="x", ylab="Gęstość")

Metoda superpozycji rozkładów

Metoda superpozycji jest stosowana, gdy docelowa gęstość prawdopodobieństwa stanowi kombinację liniową dwóch lub więcej innych gęstości. Sam proces symulacji jest dwuetapowy: najpierw z zadanym prawdopodobieństwem (wagą) losowany jest rozkład składowy, z którego w danej iteracji będzie pochodzić wynik, a w drugim etapie następuje właściwe wygenerowanie wartości z wybranego rozkładu.

Przykład praktyczny symulacji: Rozpatrzmy symulację czasu przybycia klientów do restauracji w ciągu dnia. Z empirycznych obserwacji wynika, że rozkład ten jest bimodalny i skupia się wokół dwóch szczytów natężenia ruchu: 1. Porannego (kawa wokół godziny 9:00, co stanowi 30% całkowitego ruchu). 2. Popołudniowego (obiad wokół godziny 14:00, co stanowi 70% całkowitego ruchu).

Code 
set.seed(123)
n_klientow <- 10000
czasy_przybycia <- numeric(n_klientow)

for(i in 1:n_klientow) {
  # Etap 1: Wybór rozkładu składowego (za pomocą zmiennej U(0,1))
  if(runif(1) < 0.30) {
    # 30% szans: Generowanie czasu dla szczytu porannego
    czasy_przybycia[i] <- rnorm(1, mean = 9, sd = 1.2)  
  } else {
    # 70% szans: Generowanie czasu dla szczytu popołudniowego
    czasy_przybycia[i] <- rnorm(1, mean = 14, sd = 2.0) 
  }
}

hist(czasy_przybycia, prob = TRUE, breaks = 40, col = "purple3",
     main = "Bimodalny czas przybycia klientów (Superpozycja)",
     xlab = "Godzina w ciągu dnia", ylab = "Częstość")
lines(density(czasy_przybycia), col = "darkred", lwd = 2)

Testy normalności

Zarówno przy badaniu rzeczywistych danych, jak i przy weryfikacji poprawności działania naszych generatorów, kluczowe jest sprawdzenie, czy otrzymana próba faktycznie pochodzi z rozkładu normalnego.

Do klasycznych metod wnioskowania statystycznego w tym zakresie zaliczamy testy normalności. Przyjmują one jako hipotezę zerową (\(H_0\)) założenie, że dane pochodzą z populacji o rozkładzie normalnym. Poniżej przedstawiono weryfikację rozkładu na przykładzie realnego zbioru danych mtcars, zmienna mpgjako spalanie w milach na galon z wykorzystaniem popularnych testów: Kołmogorowa-Smirnowa oraz Shapiro-Wilka.

Code 
dane_mpg <- mtcars$mpg

hist(dane_mpg, prob=TRUE, breaks=10, col="darkred", 
     main="Histogram dla zmiennej mpg ze zbioru mtcars", xlab="MPG")
lines(density(dane_mpg), col="darkgreen", lwd=2)

Code 
ks_wynik <- ks.test(dane_mpg, "pnorm", mean(dane_mpg), sd(dane_mpg))
print(ks_wynik)

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  dane_mpg
D = 0.1263, p-value = 0.687
alternative hypothesis: two-sided
Code 
shapiro_wynik <- shapiro.test(dane_mpg)
print(shapiro_wynik)

    Shapiro-Wilk normality test

data:  dane_mpg
W = 0.94756, p-value = 0.1229

Interpretacja: W obu przeprowadzonych testach wartość p-value jest większa od standardowego poziomu istotności \(\alpha = 0.05\) (dla testu Shapiro-Wilka wynosi 0.1229). Oznacza to, że nie mamy podstaw do odrzucenia hipotezy zerowej – możemy zatem przyjąć, że zmienna mpg w zbiorze mtcars ma rozkład zbliżony do normalnego.

Bibliografia

  1. Chwiej, T. (2022/2023). Generatory liczb pseudolosowych. Materiały do wykładu z Metod Numerycznych, Akademia Górniczo-Hutnicza w Krakowie (AGH). Dostępne online: https://galaxy.agh.edu.pl/~chwiej/mn/wyk/generatory_22_23.pdf
  2. Wieczorkowski, R., Zieliński, R. (1997). Komputerowe generatory liczb losowych. Warszawa: Wydawnictwa Naukowo-Techniczne (WNT), strony: 15-38.
  3. Bobrowski, D. (2004). Modele i metody matematyczne w symulacji komputerowej. Poznań: Wydawnictwo Naukowe UAM, strony: 55-82.