몬테카를로 시뮬레이션 in R

몬테카를로 시뮬레이션(Monte Carlo Simulation)은 확률을 계산할 때 사용하는 반복 계산법이다. 주사위를 던졌을 때 1이 나올 확률은 1/6이라고 이미 알고 있지만, 만약 모른다고 했을 때 어떻게 계산할 수 있을까? 주사위를 10,000번 던져서 1이 나온 비율(1이 나온 경우/전체 경우)을 계산하면 근사치를 구해낼 수 있다. 

 

예제 1) 동전을 20번 던져서 앞 면이 4번 이상 나올 확률 구하기

 

1. 동전의 앞면을 1, 뒷면을 0이라고 설정
2. 20번을 무작위로 추출해서 1을 모두 더함(=앞 면이 나온 횟수, 변수 S에 저장)
3. 위 과정을 10만번 반복
4. "S가 4이상인 경우/전체 경우"를 계산

 

runs <- 100000 # 총 횟수

# 20번 시행하여 앞면이 4번 이상 나온 경우 계산
single.run <- function(){ 
  sum(sample(c(0,1),20,replace=T)) > 3
}

# 위 함수를 10만 번 반복 시행 
# 앞면이 4번 이상 나온 경우/전체 경우
sum(replicate(runs,single.run()))/runs

 

예제 2) 원주율 구하기

그림 1. 반지름, 원의 면적, 사각형의 면적

반지름이 r인 원의 면적은 pi * r2

한 변의 길이가 2r인 사각형의 면적은 4r2(2*2r)

원의 면적/사각형의 면적 = (pi * r2)/4r2

pi=4*원의 면적/사각형의 면적

 

즉, 원의 면적과 사각형의 면적의 비율을 구할 수 있으면 pi를 구할 수 있게 된다. 그림 1의 사각형을 가로(x) 10만 칸, 세로(y) 10만 칸이 있는 바둑판이라고 가정하자. 원 내부에 있는 칸의 수를 더하면 원의 면적이라고 볼 수 있다. x2+y2=r2의 공식을 이용하여 x, y를 10만 번 추출하여 공식에 대입한다. 결과가, r2보다 작은 경우에는 원의 내부에 해당되고 그렇지 많으면 원의 외부에 해당한다.

 

#반지름 1
r <- 1
runs <- 100000

x <- runif(runs,min=-r,max=r)
y <- runif(runs,min=-r,max=r)

#원의 내부
in.circle <- x^2 + y^2 <= r^2

#(원의 면적/전체 면적)/4
(sum(in.circle)/runs)*4 

그림 2. 몬테카를로 시뮬레이션 결과

그림 2는 원 내부의 칸을 노란색으로 외부의 칸은 녹색으로 표현한 것이다. R 코드의 실행결과는 다음과 같다. 우리가 알고 있는 원주율은 3.141592... 라 시뮬레이션 값은 근사치라는 것을 알 수 있다.

> (sum(in.circle)/runs)*4
[1] 3.14056

 

 

참고자료

www.countbayesie.com/blog/2015/3/3/6-amazing-trick-with-monte-carlo-simulations