聰明不如鈍筆
총명불여둔필
몬테카를로 카지노 앞에서 열린 포뮬러원(F1) 모나코 그랑프리 경주 장면. F1 홈페이지
세상을 살다 보면 '시뮬레이션'이라는 걸 해보고 싶을 때가 있습니다.
확률을 기반으로 시뮬레이션을 진행할 때는 '몬테카를로 방법'이 거의 글로벌 스탠더드라고 해도 과언이 아닙니다.
이 블로그와 아주 좋은 친구 사이인 위피디디아는 몬테카를로 방법을 이렇게 설명합니다.
몬테카를로 방법(Monte Carlo method)은 난수를 이용하여 함수의 값을 확률적으로 계산하는 알고리즘을 부르는 용어이다.
수학이나 물리학 등에 자주 사용되며, 계산하려는 값이 닫힌 형식으로 표현되지 않거나 복잡한 경우에 근사적으로 계산할 때 사용된다.
스타니스와프 울람이 모나코의 유명한 도박의 도시 몬테카를로의 이름을 본따 명명하였다.
1930년 엔리코 페르미가 중성자의 특성을 연구하기 위해 이 방법을 사용한 것으로 유명하다.
무슨 소리인지 이해 가시나요?
이렇게 설명만 보고 무슨 뜻인지 잘 모를 때는 직접 해보는 게 오히려 도움이 될 수 있습니다.
이번 포스트에서는 몬테카를로 방법으로 원주율($\pi$) 값을 구하는 방법에 대해 알아보겠습니다.
일단 개념을 설명하는 차원에서 반지름 길이(r)가 1인 원과 이 원에 외접하는 정사각형이 있다고 가정해 보겠습니다.
그러면 이 사각형 넓이에서 원 넓이가 차지하는 비율은 얼마나 될까요?
일단 사각형 넓이를 구하는 공식은 '가로 길이 $\times$ 세로 길이'입니다.
이 정사각형은 가로와 세로가 모두 $2r=2$이니까 넓이는 $4$입니다.
원 면적을 구하는 공식은 '$\pi r^2$'이고 여기서 $r=1$이니까 이 원 넓이는 $\pi$가 됩니다.
그러면 이 사각형 넓이에서 원이 차지하는 비율은 $\frac{\pi}{4}$가 됩니다.
따라서 이 비율을 알아낸 다음 그 비율에 4를 곱하면 $\pi$ 값을 얻을 수 있습니다.
그러면 몬테카를로 시뮬레이션으로는 이 비율을 어떻게 알아낼 수 있을까요?
몬테카를로 시뮬레이션 제일 중요한 개념은 난수(亂數) 그러니까 (어떤 범위 안에서) 무작위로 뽑아낸 값입니다.
난수 두 개를 발생시킨 다음 이 숫자 두 개를 좌표로 점을 만들어 찍습니다. 이런 작업을 여러 번 반복합니다.
점을 찍을 때는 이 점이 원 안에 있는지 바깥에 있는지 판정합니다.
그러고 나서 원 안에 있는 점 개수와 원 바깥에 있는 점 개수를 따로 따로 셉니다.
그러면 이 점 숫자로 비율을 계산할 수 있습니다.
이 비율에 4를 곱하면 그 값이 바로 몬테카를로 시뮬레이션을 통해 추정한 $\pi$ 값입니다.
기본 개념을 확인했으니 이제 실제로 R로 코딩을 해보겠습니다.
제일 먼저 할 일은 언제나 그런 것처럼 tidyverse 패키지 (설치하고) 불러오기.
#install.packages('tidyverse')
library('tidyverse')## -- Attaching packages ------------------------------------------------------------------------------------------------------- tidyverse 1.3.0 --
## √ ggplot2 3.3.2 √ purrr 0.3.4 ## √ tibble 3.0.3 √ dplyr 1.0.2 ## √ tidyr 1.1.2 √ stringr 1.4.0 ## √ readr 1.3.1 √ forcats 0.5.0
## -- Conflicts ---------------------------------------------------------------------------------------------------------- tidyverse_conflicts() -- ## x dplyr::filter() masks stats::filter() ## x dplyr::lag() masks stats::lag()
이제 난수를 발생할 차례입니다.
R는 각 확률 분포에 난수를 만들 수 있는 함수를 마련해두고 있습니다.
우리는 균등하게 퍼져 있는 숫자를 추출할 거니까 runif() 함수를 쓰면 됩니다.
위에 있는 표에서 확인할 수 있는 것처럼 runif() 함수는 기본적으로 최솟값 0, 최댓값 1 사이에서 난수를 발생합니다.
이 범위를 -1과 1 사이로 바꿔서 숫자를 하나 만들어보겠습니다.
runif(n=1, min=-1, max=1)## [1] 0.7058747
※ 문자 그대로 난수이기 때문에 앞으로 나오는 숫자는 여러분 결과와 다를 수 있습니다.
이런 숫자를 10개 만들고 싶을 때는 'n=1'만 'n=10'으로 바꾸면 됩니다.
runif(n=10, min=-1, max=1)## [1] 0.71997636 -0.44805362 0.06111398 0.50521892 -0.13111627 -0.24473839 ## [7] 0.06546944 0.24709610 -0.36177674 -0.42482012
계속해서 이렇게 뽑은 난수 10개가 각각 x, y열에 들어 있는 데이터 프레임을 만들어 봅니다.
이 때는 tibble() 함수를 활용하시면 됩니다.
tibble(
x=runif(n=10, min=-1, max=1),
y=runif(n=10, min=-1, max=1)
)## # A tibble: 10 x 2 ## x y ## <dbl> <dbl> ## 1 -0.269 0.423 ## 2 0.0545 0.655 ## 3 0.518 -0.475 ## 4 0.210 0.615 ## 5 -0.319 -0.981 ## 6 0.783 -0.377 ## 7 -0.566 0.654 ## 8 0.269 -0.694 ## 9 0.829 -0.786 ## 10 -0.992 0.0944
이제 두 숫자를 좌표로 만들어 점을 찍을 겁니다.
위에서 점을 찍을 때는 이 점이 원 안에 있는지 바깥에 있는지 구분을 해줘야 한다고 했습니다.
수학에서 원은 $x^2 + y^2 = r^2$라는 공식으로 표시하기도 합니다.
그래서 어떤 좌표가 $x^2 + y^2 \le r^2$ 조건을 만족하면 그 지점은 반지름 길이 안에 있는 점이 됩니다.
다시 숫자 10개를 발생한 다음 그 점이 원 안에 있는지 바깥에 있는지 구분해 보도록 하겠습니다.
이렇게 조건을 걸 때는 ifelse() 함수를 쓰면 됩니다.
ifelse() 함수는 마이크로소프트(MS) 엑셀 if() 함수처럼 ifelse(조건, 참일 때, 거짓일 때) 형태로 씁니다.
여기서 우리가 확인하고 싶은 조건은 'x^2+y^2<=1'이 될 겁니다.
tibble(
x=runif(n=10, min=-1, max=1),
y=runif(n=10, min=-1, max=1),
in_circle=ifelse(x^2+y^2<=1, 1, 0)
)## # A tibble: 10 x 3 ## x y in_circle ## <dbl> <dbl> <dbl> ## 1 0.893 -0.285 1 ## 2 0.701 -0.417 1 ## 3 -0.167 -0.738 1 ## 4 -0.759 0.113 1 ## 5 0.0941 -0.813 1 ## 6 0.821 0.243 1 ## 7 -0.856 -0.258 1 ## 8 -0.645 -0.955 0 ## 9 0.534 0.258 1 ## 10 0.815 -0.946 0
10개를 만들 수 있다는 건 10만 개도 가능하다는 뜻이겠죠?
tibble(
x=runif(n=100000, min=-1, max=1),
y=runif(n=100000, min=-1, max=1),
in_circle=ifelse(x^2+y^2<=1, 1, 0)
) ## # A tibble: 100,000 x 3 ## x y in_circle ## <dbl> <dbl> <dbl> ## 1 -0.280 -0.0716 1 ## 2 0.938 -0.641 0 ## 3 0.324 0.531 1 ## 4 0.00239 -0.410 1 ## 5 0.630 -0.0997 1 ## 6 -0.702 -0.0535 1 ## 7 -0.456 -0.793 1 ## 8 -0.976 -0.261 0 ## 9 0.0308 -0.456 1 ## 10 0.861 -0.181 1 ## # ... with 99,990 more rows
이제 ggplot2로 실제 점을 찍어 보도록 하겠습니다.
tibble(
x=runif(n=100000, min=-1, max=1),
y=runif(n=100000, min=-1, max=1),
in_circle=ifelse(x^2+y^2<=1, 1, 0)
) %>%
ggplot(aes(x=x, y=y, color=as_factor(in_circle))) +
geom_point() +
scale_color_manual(values=c('#cdcdcd', '#79c3d6')) +
coord_equal() +
guides(color=FALSE)위에서 우리가 확인했던 그림과 아주 닮은 그림이 나왔습니다.
혹시 이 그림을 그리는 코드가 이해가 잘 가지 않으시는 분은 '최대한 친절하게 쓴 R로 그래프 그리기(feat. ggplot2)' 포스트가 도움이 될 수 있습니다.
이제 점이 몇 개씩인지 세어봐야겠죠?
원 안에 있는 점은 in_circle 열이 1입니다.
in_circle이 1인 열만 골라내서 몇 개인지 n() 함수로 셉니다.
이런 점이 몇 개인지 비율을 알아봐야 하니 이 숫자를 전체 점 숫자(10만 개)로 나눕니다.
이 숫자에 4를 곱하면 $\pi$가 나옵니다.
format() 함수를 써서 소수점 이하를 조금 더 열어보는 작업까지 진행하겠습니다.
tibble(
x=runif(n=100000, min=-1, max=1),
y=runif(n=100000, min=-1, max=1),
in_circle=ifelse(x^2+y^2<=1, 1, 0)
) %>%
filter(in_circle==1) %>%
summarise(pi=format(4*n()/100000, digits=5))## # A tibble: 1 x 1 ## pi ## <chr> ## 1 3.1437
실제로 $\pi$는 3.14159265359…니까 이 값이 아주 정확하다고 할 수는 없습니다.
그래도 우리가 초등학교 때 배운 것처럼 $\pi = 3.14$라고 표현하는 데도 큰 무리가 없습니다.
이 정도면 꽤 쓸 만한 결과라고 할 수 있는 셈입니다.
이번에도 filter(), summarise() 함수가 이해가 가지 않으시는 분이 계시다면 '최대한 친절하게 쓴 R로 데이터 뽑아내기(feat. dplyr)' 포스트가 도움이 될 수 있습니다.
우리는 이번 시뮬레이션에서 점을 10만 개 찍었습니다.
점 개수가 늘어날 때마다 이 결과가 얼마나 실제로 $\pi$에 가까워지는지 아래 코드로 확인할 수 있습니다.
tibble(
x=runif(n=100000, min=-1, max=1),
y=runif(n=100000, min=-1, max=1),
in_circle=ifelse(x^2+y^2<=1, 1, 0)
) %>%
mutate(prop=cumsum(in_circle)*4/row_number()) %>%
ggplot(aes(x=c(1:100000), y=prop)) +
geom_line() +
geom_hline(yintercept=pi, lwd=1, color='#79c3d6', linetype='dotted')$\pi$ 값을 알아볼 때는 약 1만2500개 정도 점을 찍으면 비교적 안정적인 값을 얻는다는 사실을 확인할 수 있습니다.
그럼 여러분 모두 Happy Tidyversing!
댓글,
