Zagadnienia z Wykładu

Antoni Dałek

1. Szacowanie liczby pi

1.1. Igła Buffona

Problem igły Buffona polega na rzucaniu igły o długości \(l\) na płaszczyznę z narysowanymi równoległymi liniami w odstępie \(a\) (przy założeniu \(l \le a\)). Prawdopodobieństwo \(P\), że igła przetnie linię, wyraża się wzorem:\[P = \frac{2l}{\pi a}\]Z tego wynika, że liczbę \(\pi\) możemy oszacować jako:\[\pi \approx \frac{2lN}{na}\]gdzie:
\(N\) – całkowita liczba rzutów,
\(n\) – liczba rzutów, w których igła przecięła linię.

1.2 Igła Buffona w R

Poniższy kod przedstawia symulację rzutów igłą Buffona.

Code

set.seed(67)
N = 10000 #Liczba rzutów igłą         
l = 0.8 #Długość igły           
a = 1.0 #Odległość między liniami           
x = runif(N, min = 0, max = a / 2) #Losowanie odległości środka igły od najbliższej linii
ang = runif(N, min = 0, max = pi / 2) #Losowanie kąta nachylenia igły
p = x <= (l / 2) * sin(ang) #Sprawdzenie, czy igła przecina linię
kump = cumsum(p) #Wyliczenie skumulowanej sumy przecięć
nr_rzutu = 1:N
pi_szac = (2 * l * nr_rzutu) / (a * kump) #Szacowanie pi
pi_kon = pi_szac[N] #Końcowy wynik szacowania

cat("Oszacowana wartość pi:", pi_kon, "\n")
## Oszacowana wartość pi: 3.13664
cat("Prawdziwa wartość pi:", pi, "\n")
## Prawdziwa wartość pi: 3.141593

Szacowane pi w zależności od liczby wykonanych rzutów:

1.3 Strzelanie do tarczy

Strzelanie do tarczy to metoda estymacji liczby \(\pi\) oparta na polach figur. Przebiega ona następująco:

Losowanie punktów: Generowanie losowych współrzędnych \((x, y)\) na polu o postaci kwadratu (z rozkładu jednostajnego \(U[-1, 1]\)).

Warunek trafienia: Sprawdzenie, czy wylosowany punkt wpadł do koła wpisanego w ten kwadrat, czyli czy spełnia warunek odległości od środka: \(x^2 + y^2 \le 1\).

Estymacja \(\pi\): Obliczenie stosunku liczby punktów znajdujących się w kole do liczby wszystkich punktów. Ponieważ stosunek pola koła do pola kwadratu wynosi \(\pi/4\), wartość \(\pi\) szacujemy mnożąc otrzymany stosunek przez 4.

1.4 Strzelanie do tarczy w R

Poniższy kod przedstawia estymację liczby \(\pi\) za pomocą strzelania do tarczy.

Code

set.seed(67)
N = 5000  #Liczba punktów     
x = runif(N, min = -1, max = 1) #Losowanie współrzędnych x i y z rozkładu jednostajnego U[-1, 1]
y = runif(N, min = -1, max = 1)
traf = (x^2 + y^2) <= 1 #Sprawdzenie, czy punkt trafił w koło
pi_szac = 4 * sum(traf) / N #Szacowanie liczby pi
cat("Oszacowana wartość pi:", pi_szac, "\n")
## Oszacowana wartość pi: 3.1528
cat("Prawdziwa wartość pi:", pi, "\n")
## Prawdziwa wartość pi: 3.141593

Wizualizacja wylosowanych punktów:

2. Generatory wartości

2.1 Przekształcanie rozkładu \(U[0, 1]\) na dowolny rozkład \(U[a, b]\)

Z wykorzystaniem rozkładu jednostajnego \(U[0, 1]\) można wygenerować liczby z dowolnego przedziału \([a, b]\). Aby to zrobić, stosujemy prostą transformację liniową: \[Y = (b-a)U + a\] Prosty przykład przekształcenia 10 liczb z rozkładu \(U[0, 1]\) na takie z przedziału [10, 14] z wykorzystaniem R.

set.seed(67)
x = runif(10) #Generujemy 10 wartości z rozkładu U[0, 1]
y = x * 4 + 10 #Transformacja Y = (b-a)U + a
cat("Pierwotne wartości:", x, "\n")
## Pierwotne wartości: 0.8886542 0.1627505 0.4504991 0.9213891 0.1691578 0.635835 0.4684489 0.5786171 0.1264339 0.7146106
cat("Przekształcone wartości:", y, "\n") #Wartości są z przedziału [10, 4]
## Przekształcone wartości: 13.55462 10.651 11.802 13.68556 10.67663 12.54334 11.8738 12.31447 10.50574 12.85844

2.2 Metoda odwracania dystrybuanty

Rozkład prawdopodobieństwa można wygenerować, jeśli znamy funkcję odwrotną do jego dystrybuanty (\(F^{-1}\)).

Kroki algorytmu:

Losujemy wartość \(u\) z rozkładu jednostajnego \(U(0, 1)\).
Wyznaczamy wartość zmiennej \(z\) poprzez funkcję odwrotną:\[z = F^{-1}(u)\]

2.3 Odwracanie dystrybuanty, rozkład normalny

Komenda ‘qnorm’ w programie R to funkcja \(F^{-1}\) dla rozkładu normalnego. Poniższy kod przedstawia wygenerowanie 10000 wartości z rozkładu jednostajnego na przedziale \((0, 1)\), a następnie przekształcenie je w zmienne o rozkładzie normalnym \(N(0, 1)\).

Code

set.seed(70)
N = 10000
u = runif(N) #Generowanie wartości z U(0,1)
z = qnorm(u) #Odwrócenie dystrybuanty dla rozkładu N(0,1)

Wartości przedstawione jako histogram:

2.4 Odwracanie dystrybuanty, rozkład wykładniczy

Poniższy kod przedstawia wygenerowanie 10000 wartości z rozkładu jednostajnego na przedziale \((0, 1)\), a następnie przekształcenie je w zmienne o rozkładzie wykładniczym

Code

set.seed(123)
N = 10000
lambda = 2  

#Generowanie wartości z U(0, 1)
u = runif(N)

#Odwrócenie dystrybuanty dla rozkładu wykładniczego
z <- -log(1 - u) / lambda

Wartości przedstawione jako histogram:

2.5 Metoda odrzucania dla rozkładów warunkowych

Jest to generowanie zmiennych losowych z rozkładu opartego na pewnym warunku.

Zasada działania:

Generowane są wartości z rozkładu bezwarunkowego (bazowego).
Odrzucane są wszystkie te wartości, które nie spełniają narzuconego warunku.

2.6 Metoda odrzucania dla rozkładów warunkowych w R

Poniższy skrypt generuje zmienne z rozkładu \(N(0, 1)\) pod warunkiem \(X \in (-2, 2)\).

Code

set.seed(110)
N = 1000
baza = rnorm(N) #Generowanie wartości z rozkładu bezwarunkowego N(0, 1)
warunek = (baza < 2 & baza > -2) #Warunek wartości z przedziału (-2, 2)
x = baza[warunek] #Zachowanie wartości spełniających warunek
cat("Liczba prób wyjściowych:", N, "\n")
## Liczba prób wyjściowych: 1000
cat("Liczba zaakceptowanych prób:", length(x), "\n")
## Liczba zaakceptowanych prób: 965

Poniższy histogram przedstawia gęstość wygenerowanego rozkładu warunkowego. Rozkład jest ostro ucięty, a wszystkie wygenerowane wartości mieszczą się ściśle w zadanym przedziale \((-2, 2)\).

2.7 Metoda eliminacji

Celem jest wygenerowanie zmiennej losowej o gęstości \(f(x)\). Dysponujemy rozkładem pomocniczym o gęstości \(g(x)\), z którego potrafimy generować wartości, oraz stałą \(c \ge 1\), taką że dla każdego \(x\) spełniony jest warunek:\[f(x) \le c \cdot g(x)\] Kroki algorytmu:

Generujemy wartość \(Y\) z rozkładu o gęstości \(g(x)\).
Generujemy wartość \(U\) z rozkładu jednostajnego \(U(0, 1)\).
Sprawdzamy warunek akceptacji:Jeśli \(U \le \frac{f(Y)}{c \cdot g(Y)}\), to akceptujemy \(Y\) i przyjmujemy \(X = Y\). W przeciwnym razie odrzucamy \(Y\) i powtarzamy kroki 1-3.

2.8 Metoda eliminacji w R

W poniższym kodzie generujemy wartości z rozkładu \(Beta(2, 2)\), którego gęstość to \(f(x) = 6x(1-x)\) na przedziale \([0,1]\). Jako rozkład pomocniczy \(g(x)\) przyjęto rozkład jednostajny \(U(0, 1)\), dla którego \(g(x) = 1\). Stała \(c = 1.5\)

Code

set.seed(67)
N = 5000
c = 1.5

f = function(x) { 6 * x * (1 - x) } #Definicja gęstości docelowej f(x)

#Generowanie propozycji Y z rozkładu pomocniczego g(x) = U(0,1)
Y = runif(N)

#Generowanie zmiennej decyzyjnej U z rozkładu U(0,1)
U = runif(N)

a = (U * c) <= f(Y) #Weryfikacja warunku akceptacji

Xa = Y[a] #Wyodrębnienie zaakceptowanych wartości

Wizualizacja algorytmu (wartości zaakceptowane są oznaczone na zielono i tworzą rozkład \(Beta(2, 2)\)):

2.9 Metoda superpozycji rozkładów

Dla \(k\) rozkładów składowych gęstość wyraża się wzorem: \[f(x) = \sum_{i=1}^{k} p_i f_i(x)\] gdzie \(p_i > 0\) to wagi (prawdopodobieństwa wyboru danego rozkładu), a \(\sum_{i=1}^{k} p_i = 1\).

Algorytm generacji:

Generujemy dyskretną zmienną losową \(I \in \{1, 2, \dots, k\}\), przyjmującą wartość \(i\) z prawdopodobieństwem \(p_i\). Zmienna ta pełni rolę wskaźnika wyboru rozkładu składowego.
Generujemy wartość zmiennej losowej \(X\) z wybranego w pierwszym kroku rozkładu o gęstości \(f_I(x)\).

2.10 Metoda superpozycji rozkładów w R

Jako przykład pokazuję generowanie zmiennych z dwumodalnego rozkładu będącego mieszaniną dwóch rozkładów normalnych: \(0.3 \cdot N(-3, 1) + 0.7 \cdot N(3, 1.5)\).

Code

set.seed(67)
N = 10000

p = c(0.3, 0.7) #Wagi składowe (prawdopodobieństwa p_i)

#Losowanie wektora wskaźników 
wsk = sample(1:2, size = N, replace = TRUE, prob = p)

#Zliczenie, ile zmiennych wygenerować z poszczególnych rozkładów
n1 = sum(wsk == 1)
n2 = sum(wsk == 2)

#Generowanie z rozkładów składowych
x1 = rnorm(n1, mean = -3, sd = 1)
x2 = rnorm(n2, mean = 3, sd = 1.5)

#Złożenie wyników w jeden końcowy wektor 
Xfin = c(x1, x2)

Histogram obrazuje dwumodalny kształt wygenerowanych danych.

3. Bibliografia

Haugh, Martin. “Generating Random Variables and Stochastic Processes.” Columbia University, www.columbia.edu/~mh2078/MonteCarlo/MCS_Generate_RVars.pdf

Hwang, Chi-Ok, et al. “Buffon’s Needle Algorithm to Estimate π.” SCIRP, Scientific Research Publishing, 3 Mar. 2017, www.scirp.org/journal/paperinformation?paperid=74541

Rolski, Tomasz. “Symulacje Stochastyczne I Teoria Monte Carlo.” Uniwersytet Wrocławski, Instytut Matematyczny, www.math.uni.wroc.pl/~rolski/Zajecia/sym17.pdf
(54 - 57)

Dobra, Adrian. “General Sampling Methods.” University of Washington, Department of Statistics, sites.stat.washington.edu/mmp/courses/stat534/spring19/Handouts/534sampling.pdf