Symulacje to coś, do czego komputery nadają się najlepiej. Zamiast liczyć na kartce wykorzystujemy maszynę i dostajemy wynik. Dzisiaj policzymy liczbę Pi i przeanalizujemy pewien teleturniej telewizyjny.
Liczba Pi i metoda Monte Carlo
Policzymy wartość liczby Pi… upuszczając ołówek na kartkę. Możecie zrobić to samodzielnie, albo wykorzystując metodę Monte Carlo.
Metoda Monte Carlo to numeryczny sposób rozwiązywania problemów. Wykonujemy bardzo dużo obliczeń (opartych o losowe symulacje) i na podstawie ich wyniku coś określamy.
Metoda została opracowana przez Stanisława Ulama, który należał do tak zwanej lwowskiej szkoły matematycznej i był między innymi współtwórcą amerykańskiej bomby termojądrowej (co akurat powinno wzbudzić szacun a nie jakiś postrach). Ulam był też jednym z pierwszych naukowców, którzy wykorzystywali w swych pracach komputer. Możecie (zachęcam!) obejrzeć film o lwowskiej szkole matamatycznej pochodzący z kanału Historia Bez Cenzury – mam nadzieję, że wzbudzi Wasz podziw (i odkrycie ludzkich cech u matematyków) dla polskich przedwojennych matematyków. Każdy kto miał na studiach matematykę (wszystkie uczelnie techniczne) zapewne słyszał o przestrzeni Banacha. Tak, ma to związek z lwowską szkołą, a Banach zresztą pochowany jest na Cmentarzu Łyczakowskim, podobnie jak wielu innych znanych Polaków (ze znanych ze szkoły nazwisk: Maria Konopnicka i Gabriela Zapolska).
Wróćmy jednak do naszego eksperymentu. Polega on na losowym upuszczeniu (na przykład) ołówka, tak aby zaznaczył punkt na kartce. Wcześniej na tej kartce rysujemy koło i opisany na nim kwadrat. Ołówkiem rzucamy z góry tak, aby trafić w kwadrat. Liczymy liczbę rzutów i liczbę trafień w koło.
Zamiast robić to ręcznie wykorzystamy komputer (i oczywiście R). Zamiast rzucania wylosujemy punkt:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 |
kolo <- data_frame(x = numeric(), y = numeric(), k = logical()) # promień koła r <- 100 # rzucamy 5000 razy ołówkiem for(i in 1:5000) { # losowe współrzędne punktu x <- (-1)^sample(c(1, 2), 1) * sample(0:r, 1) y <- (-1)^sample(c(1, 2), 1) * sample(0:r, 1) # w jakim miejscu upadł ołówek? kolo[i, "x"] <- x kolo[i, "y"] <- y # czy wylosowany punkt wypadł w kole? kolo[i, "k"] <- ifelse(x^2 + y^2 <= r^2, TRUE, FALSE) } |
Gdzie upadł ołówek?
|
1 2 3 4 5 6 7 |
library(tidyverse) kolo %>% ggplot() + geom_point(aes(x, y, color = k)) + coord_equal() + labs(color = "Czy punkt\njest w kole?") |
Widać, że mniej więcej wszystkie punkty rozłożone są równomiernie. Liczba Pi to liczba punktów w kole w stosunku do liczby wszystkich punktów – daje nam to 3.0896, czyli całkiem niezły wynik. Dla przypomnienia Pi to 3.1415926 (tyle liczb pamiętam).
A teraz to samo ćwiczenie powtórzymy miliony razy. (Prawie) za każdym razem z różną liczbą upuszczanych ołówków – tak, aby sprawdzić dokładność metody i czy z większą liczbą losowań (ołówków) dostaniemy lepszy wynik?
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 |
df <- data.frame(i = numeric(), p = numeric(), pi_mc = numeric()) r <- 1000 k <- 0 max_p <- 200 max_j <- 16 for(j in 1:max_j) { # określamy liczbę losowań - kolejne potęgi dwójki max_i <- 2^j # prosty "progress bar" :) cat(paste("\rj:", j, "/", max_j, "\tmax_i: ", max_i, "\n")) # dla każdego z max_p losowych kółek for(p in 1:max_p) { # losujemy max_i punktów i sprawdzamy czy wypadają w kole kolo <- 0 for(i in 1:max_i) { # losowe współrzędne koła x <- (-1)^sample(1:max_i, 1) * sample(0:r, 1) y <- (-1)^sample(1:max_i, 1) * sample(0:r, 1) # czy wylosowany punkt wypadł w kole? if(x^2 + y^2 <= r^2) kolo <- kolo + 1 } # obliczamy Pi Pi <- 4 * kolo / max_i # 4, bo kwadrat ma pole 4r^2: (2r)^2, gdzie r jest od -r do +r # zapisujemy wyniki w tabeli k <- k + 1 df[k, "i"] <- max_i df[k, "pi_mc"] <- Pi df[k, "p"] <- p } } |
Obliczenia trwają długo – mamy po 2, 4, 8, …, 216 losowanych punktów dla jednego koła, a kół będzie po 200 dla każdej z liczby losowanych punktów. Łącznie daje to sum(200*2^(1:16)) losowanych punktów w kole (czyli jakieś 26.2 milionów losowań).
obliczenia są bardzo długie i poprawnie byłoby je zrównoleglić (przeczytaj tekst), a jeszcze lepiej napisać w C/C++ (z pomocą pakietu Rcpp – polecam tekst oraz wpis). Tutaj tego nie użyłem ze względu na czytelność kodu, a także charakter wpisu – chodzi mi o wyjaśnienie algorytmu i komplikowanie całości przez wprowadzenie dodatkowych rozwiązań do niczego dobrego nie prowadzi.
Zobaczmy co nam wyszło – jak wygląda wartość Pi w zależności od liczby losowań?
|
1 2 3 4 5 6 |
df %>% ggplot() + geom_boxplot(aes(i, pi_mc, group=i), fill = "lightgreen", color = "darkred") + geom_hline(yintercept = pi, color="red") + geom_jitter(aes(i, pi_mc), width = 0.05, alpha=0.1, color = "blue") + scale_x_log10() |
Widzimy, że im więcej powtórzeń tym bardziej obliczona liczba Pi bliższa jest wartości prawdziwej, jednocześnie z coraz mniejszym rozrzutem. W którymś momencie dokładność jest wystarczająca.
Do czego można to wykorzystać? Całkowanie! W metodzie MC możemy zrealizować to poprzez losowanie dla każdego X jakiejś wartości (wielokrotnie) i sprawdzanie czy jest pod czy nad krzywą, której całkę chcemy policzyć. To co jest pomiędzy osią OX (y = 0) a krzywą jest wartością naszej całki. Im więcej losowań dla danego X i im mniejsze odstępy pomiędzy kolejnymi Xami tym lepszy (dokładniejszy) wynik. Mechanika jest identyczna jak w przypadku kół i liczby Pi.
Będąc jeszcze przy liczbie Pi: wiecie, że rozkład wystąpienia poszczególnych cyfr jest mniej więcej równy? To znaczy, że biorąc pod uwagę na przykład tysiąc miejsc po przecinku otrzymamy zbiór liczb z 100 jedynek, 100 dwójek, 100 trójek i tak dalej? Dokładniej możecie to zobaczyć na poniższej animacji (via Reddit)
a jeśli samodzielnie chcecie sobie to sprawdzić to na stronie PiDay.org znajdziecie liczbę Pi w rozwinięciu do miliona miejsc po przecinku. Proszę policzyć rozkład poszczególnych cyfr. W pamięci. No dobra, może być na palcach ;)
Więcej o liczbie Pi chociażby na Wikipedii. Ciekawostką jest angielskie zdanie How I want a drink, alcoholic of course, after the heavy lectures involving quantum mechanics!. Dlaczego jest ciekawe? Ten kod da odpowiedź:
|
1 2 3 4 5 |
"How I want a drink, alcoholic of course, after the heavy lectures involving quantum mechanics!" %>% gsub("[,!]", "", .) %>% strsplit(" ") %>% unlist() %>% nchar() |
Teleturniej “The Wall”
W piątkowy wieczór trafiłem w TVP na teleturniej “The Wall” (na postawie amerykańskiej wersji). Pokrótce polega on na odpowiadaniu na pytania, a nagrodą jest kwota określona na podstawie tego, gdzie wyląduje puszczona po ścianie kula. Plansza (owa ściana) jest mniej więcej prostokątna (to chyba raczej trapez, nie przyglądałem się aż tak dokładnie), uczestnik wybiera kolumnę z której wypada kula. Na każdym poziomie kula odbija się od wystających ze ściany kołków, co powoduje jej odbicie w lewo lub w prawo, czasem kula spada bez zmiany kierunku.
Czy możemy określić z którego miejsca należy puszczać kulę, aby mieć największe szanse na najwyższe nagrody?
Przygotujmy mechanizm symulujący drogę przez planszę dla pojedynczej kulki. Zaczniemy od wybranej kolumny, a przy przejściu do kolejnego rzędu wylosujemy kierunek zmiany (w lewo, w prawo albo bezpośrednio w dół). W zależności od wylosowanego kierunku zmieniamy numer kolumny w której będzie kulka (o jeden więcej, mniej albo bez zmiany) w kolejnym rzędzie. Założenie dodatkowe: przy dojściu do bocznej ściany kulka odbija się od niej.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 |
rm(list = ls()) # rozmiar planszy wysokosc_planszy <- 20 szerokosc_planszy <- 10 # startowa kolumna kolumna <- 5 # definiujemy pustą planszę plansza <- matrix(0, nrow = wysokosc_planszy, ncol = szerokosc_planszy) # przejście przez planszę od góry w dół for(rzad in 1:nrow(plansza)) { # w odpowiednim rzędzie i kolumnie zapisujemy położenie kulki # w tym przypadku jest to numer kolejnego kroku równy rzędowi plansza[rzad, kolumna] <- rzad # w lewo, w prawo czy w dół? losujemy kierunek switch(sample(c("L", "P", "S"), size = 1), "L" = kolumna <- kolumna + 1, "P" = kolumna <- kolumna - 1, "S" = kolumna <- kolumna) # odbicie od prawej krawędzi planszy if(kolumna > szerokosc_planszy) kolumna <- ncol(plansza) # odbicie od lewej krawędzi planszy if(kolumna == 0) kolumna <- 1 } |
Co z tego wyszło? Zobaczmy na obrazku (wpierw robiąc trochę gimnastyki z macierzą tak, aby łatwo ją narysować z użyciem ggplot):
|
1 2 3 4 5 6 7 8 9 10 11 12 |
plansza %>% as.data.frame() %>% mutate(Rzad = row_number()) %>% gather(key = Kolumna, value = Krok, -Rzad) %>% mutate(Kolumna = as.numeric(gsub("V", "", Kolumna))) %>% mutate(Kolumna = factor(Kolumna, levels = 1:szerokosc_planszy), Rzad = factor(Rzad, levels = wysokosc_planszy:1)) %>% ggplot() + geom_tile(aes(Kolumna, Rzad, fill = ifelse(Krok == 0, "white", "black")), color = "black", show.legend = FALSE) + geom_text(aes(Kolumna, Rzad, label = ifelse(Krok == 0, "", Krok))) + scale_fill_manual(values = c("white" = "white", "black" = "orange")) |
Na pomarańczowo widzimy drogę jaką przebyła kula, a liczby w pomarańczowych polach to kolejny krok.
Wiecie, że to wszystko to łancuch Markova (w najprostszym wariancie zwanym spacerem losowym)? To znaczy, że stan przyszły (w naszym przypadku numer kolumny w kolejnym rzędzie) zależy tylko od stanu obecnego (numer kolumny w aktualnym rzędzie). Kulka nie pamięta swojej historii (i nie wie, że jeśli ostatnio było w prawo to teraz na przykład powinno być w lewo). Istnieje oczywiście pakiet do tego typu zagadnień – mcmc, ale szczerze mówiąc nie zagłębiałem się w niego i nie wiem jak to samo zrobić za jego pomocą. Na pewno krócej i szybciej.
Teraz zrobimy dokładnie to samo co z rzucaniem ołówkiem – powtórzymy to samo ćwiczenie tysiące razy. Dodatkowo zaczynając od każdej z kolumn.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 |
df <- data_frame() iteracje <- 10000 powtorki <- 1000 # powtorka wszystkich przejść for(t in 1:powtorki) { # amatorski progress bar :) cat(paste0("\rt = ", t, "/", powtorki)) # losowy generator losowy :) set.seed(floor(runif(1, min = 1, max = 5000))) # dla każdej kolumny startowej... for(kolumna_startowa in 1:szerokosc_planszy) { # ... definiujemy czystą planszę sumaryczną... plansza_full <- matrix(0, nrow = wysokosc_planszy, ncol = szerokosc_planszy) # ...robimy 'iteracje' powtorzeń przejścia przez planszę for(i in 1:iteracje) { # czyścimy planszę plansza <- matrix(0, nrow = wysokosc_planszy, ncol = szerokosc_planszy) # zaczynamy od kolumny startowej kolumna <- kolumna_startowa # przejście przez planszę od góry w dół for(rzad in 1:nrow(plansza)) { # jedynką oznaczamy pozycję kulki plansza[rzad, kolumna] <- 1 # w lewo, w prawo czy w dół? switch(sample(c("L", "P", "S"), size = 1), "L" = kolumna <- kolumna + 1, "P" = kolumna <- kolumna - 1, "S" = kolumna <- kolumna) # odbicie od prawej strony if(kolumna > szerokosc_planszy) kolumna <- ncol(plansza) # odbicie od lewej strony if(kolumna == 0) kolumna <- 1 } # heatmapa wszystkich przejść = suma wszystkich plansz plansza_full <- plansza_full + plansza } # największe prawdopodobieństwo zakończenia # w której kolumnie ostatniego wiersza było najwięcej kulek na planszy sumarycznej? df <- bind_rows(df, data_frame(start = kolumna_startowa, finisz = which.max(plansza_full[wysokosc_planszy,]))) } } |
Znowu trwało to bardzo długo, znowu nalepiej byłoby to zrobić korzystając ze zrównoleglenia obliczeń i opierając je na kodzie w C/C++. Zysk na wydajności byłby pewnie z 20-30 krotny.
Ale zobaczmy gdzie wpadała kula (oś pionowa) w zależności od tego skąd została wyrzucona (oś pozioma):
|
1 2 3 |
df %>% ggplot() + geom_jitter(aes(start, finisz), alpha = 0.1) |
Widać coś, ale przeliczmy liczbę punktów na prawdopodobieństwo i nieco uporządkujmy sam wykres:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 |
df %>% mutate(start = as.factor(start), finisz = as.factor(finisz)) %>% count(start, finisz) %>% ungroup() %>% group_by(start) %>% mutate(p = 100*n/sum(n)) %>% mutate(maks = n == max(n)) %>% ungroup() %>% ggplot() + geom_tile(aes(start, finisz, fill = n), color = "gray80", show.legend = FALSE) + geom_text(aes(start, finisz, label = sprintf("%.1f%%", p), size = maks), show.legend = FALSE) + scale_fill_gradient(low = "darkgreen", high = "yellow") + scale_size_manual(values = c("TRUE" = 5, "FALSE" = 3)) + labs(x = "Kolumna początkowa", y = "Kolumna końcowa") |
Prawdopodobieństwo liczone jest w pione powyższego wykresu. To znaczy, że upuszczając kulę z kolumny 4 w 50.1% przypadków wylądowała ona w kolumnie 1, w 37.6% w kolumnie 2, a w 12.3% – w kolumnie 3. Widać, że największy rozrzut jest dla kolumn środkowych. Prawie pewni możemy być, że puszczając kulkę z kolumn 2, 8 i 10 wylądujemy w kolumnach skrajnych (odpowiednio 1 i 10).
Przygotujmy jeszcze tabelę kolumna startowa – kolumna końcowa, która powinna graczom ułatwić wybór kolumn:
|
1 2 3 4 5 6 7 8 9 |
df %>% count(start, finisz) %>% ungroup() %>% group_by(start) %>% mutate(p = 100*n/sum(n)) %>% mutate(maks = n == max(n)) %>% ungroup() %>% filter(maks) %>% select(start, finisz) |
| Kolumna startowa | Kolumna końcowa |
|---|---|
| 1 | 1 |
| 2 | 1 |
| 3 | 1 |
| 4 | 1 |
| 5 | 4 |
| 6 | 7 |
| 7 | 10 |
| 8 | 10 |
| 9 | 10 |
| 10 | 10 |
Można to uprościć do reguł:
- 1 wypadnie, kiedy zaczniemy od kolumn 1 do 4
- dla początkowej 5 wypadnie 4
- dla 6 wypadnie 7
- dla pozostałych (8-10) wypadnie 10
Teraz tylko wystarczy zobaczyć gdzie jest największa kasa (interesują nas w sumie tylko cztery kolumny: 1, 4, 7 i 10) i odpowiednio wypuścić kulkę.
Oczywiście to gra losowa i przy jednej próbie nie musi wcale tak być. Dodatkowo ważny aspekt naszej symulacji: uprościliśmy tablicę do kwadratowej (w teleturnieju jest trapezem), uprościliśmy rozmieszczenie kołków odbijających kulę oraz uprościliśmy liczbę kolumn startowych do takiej samej dla każdej rundy (zdaje się, że w teleturnieju z każdą rundą zmienia się liczba kolumn z której może wypaść kula).
Szukając materiałów o teleturnieju znalazłem podobną analizę na blogu, który może zainteresować fanów NBA (pod kątem analitycznym). Tam znajdziecie lepiej rozpoznaną planszę i więcej liczb (dlaczego autor nie zrobił z tego heatmapy? porównywanie liczb z czterema miejscami po przecinku to koszmar; mógł chociaż wyróżnić największe wartości w rzędzie w ostatniej tabeli…). Idea pozostaje jednak taka sama.
W piątek pewnie obejrzę teleturniej jeszcze raz (powtórka jest bodaj w niedziele) – synowi starszemu strasznie te kule się podobają, a Was zachęcam do eksperymentów. Policzcie jaka jest oczekiwana wartość wygranej przy poprawnych odpowiedziach na wszystkie pytania.
Ciekawy temat, ale mam kilka uwag.
Pierwsza sprawa – Losowanie współrzędnych.
x <- (-1)^sample(c(1, 2), 1) * sample(0:r, 1)
Ewidentnie pierwszy sample służy do losowania znaku. Nie można prościej?
sample(c(-1,1),1) ?
Jeszcze śmieszniej jest w drugim przykładzie:
x <- (-1)^sample(1:max_i, 1) * sample(0:r, 1)
Losować liczbę z przedziału 1..65536 i podnosić potem -1 do tak wysokiej potęgi tylko po to, żeby wylosować znak? To nie może być wydajne. :)
Z tym znakiem jest też inny problem. Zero jest losowane dwa razy częściej niż inne wartości. Widać to jak r będzie małe, np.2. Wtedy losowanie jest ze zbioru (-1)*(0,1,2) + (+1)*(0,1,2). Widać, że mamy zbiór (0,0,1,-1,2,-2). Owszem, przy r=1000 różnica jest już mała, ale przecież zależy nam na precyzji Pi.
Można pozbyć się obu powyższych problemów losując liczby tylko z pierwszej ćwiartki. Pomijamy wtedy mnożenie przez 4 we wzorze Pi <- 4 * kolo / max_i
Druga sprawa – dokładność
"Widzimy, że im więcej powtórzeń tym bardziej obliczona liczba Pi bliższa jest wartości prawdziwej" To nie jest cała prawda. Dokładność obliczenia Pi zależy od przyjętego R. Wraz ze zwiększaniem liczby losowań dla r=1000 zbliżamy się do Pi=3,141549. Dlaczego? Bo w kwadracie o boku 2000 wewnątrz koła jest 3141549 punktów o współrzędnych całkowitych. Do zwiększenia dokładności trzeba zwiększyć r.
Trzecia sprawa – zasadność MC dla Pi
MC jest całkiem ciekawą metodą, ale akurat dla liczenia Pi jest chyba najgorszą możliwą. :) Po milionach operacji uzyskaliśmy wartość Pi z dokładnością do 4 cyfr po przecinku. Czyli mniej niż w ułamku 335/113. :)
Ciekawe było by obliczyć/sprawdzić ile powtórzeń trzeba by zastosować dla r dajmy na to 10^9, żeby uzyskać pi z sensowną dokładnością. To mogło by posłużyć do oszacowania ile razy losować przy użyciu MC do rzeczywistych zastosowań.
Warto też podkreślać, że MC jest bardzo kosztowną metodą i przed użyciem warto sprawdzić, czy nie ma innego sposobu. Do tego trzeba się upewnić czy rozkład badanego zjawiska jest równomierny.
Brawo, takie komentarze lubię! :)
I wszystko się zgadza. Tutaj (jak w większości moich wpisów) chodziło o pokazanie metody.
Co do losowania znaku – pisząc kod zafiksowałem się na wzorach matematycznych całkowicie zapominając o wydajności. Tak więc znak n-tej liczby określamy jako (-1)^n, co przy dużym n robi się problemem.
Losowanie dla całego koła zamiast ćwiartki – chodzi tylko o obrazek (w końcu rzucamy ołówkiem w koło, więc powinno być je widać).
„Liczba Pi to liczba punktów w kole w stosunku do liczby wszystkich punktów” ale wszystkich punktów jest więcej więc dlaczego wychodzi większa od jeden?
Stąd mnożenie przez 4:
# obliczamy PiPi < - 4 * kolo / max_i # 4, bo kwadrat ma pole 4r^2: (2r)^2, gdzie r jest od -r do +r
Ciekawy wpis, a uwagi z komentarzy trafne. Mi przydało się w szkole :P