몬테카를로(Monte-Carlo) 시뮬레이션에 대한 이론적인 설명

1. 목표

 

직사각형 안에 어떤 도형을 그려놓자.

 

그림1. 직사각형안에 그린 도형

 

빨간색 영역의 넓이는 얼마인지 알고 싶다.

 

 

2. 기본적인 원리

 

만약, 위와 같은 직사각형에서 임의의 난수를 하나 뽑는다고 하자.

 

그 난수가 빨간색 영역인 HIT에 들어갈 확률은 얼마인가?

 

직사각형의 넓이는 $c(b-a)$이고

 

빨간색 영역의 넓이를 $S$라고 하면, 기하학적 확률의 원리에 의해 $$p= \frac{(난수가 \; 목표로 \; 하는 \; 빨간색 \; 영역의 \; 넓이)}{(난수가 \; 있을 \; 수 \; 있는 \; 전체 \; 영역의 \; 넓이)} = \frac{S}{c(b-a)}$$

 

그러나 $S$를 모른다는 것이 중요하다. 즉 우리는 p값도 알 수가 없다

 

그런데 $p$값을 다른 방법으로 추정해볼 수 있는데

 

위와 같은 직사각형 위에서 $N$개의 난수를 뽑는다고 하자.

 

각 난수가 HIT영역에 들어갈 확률은 위에서 구한대로 $p$이고 MISS영역에 들어갈 확률은 $1-p$

 

컴퓨터로 $N$개의 난수를 생성했을 때 HIT영역에 들어간 난수의 수를 (셀 수 있다면) $N _{H}$라고 하자.

 

경험적으로 우리는 난수가 HIT영역에 들어갈 확률이 $\frac{N _{H}}{N}$

 

그렇다면 대수의 법칙(The law of large number)에 의해 경험적확률은 시도 횟수 $N$이 무수히 많아지면 수학적 확률 $p$에 수렴하므로 $$\lim _{N->  \infty } \frac{N _{H}}{N} =p$$

 

이 사실은 $N$을 충분히 크게하면 $p$의 추정값으로 $\frac{N _{H}}{N}$를 쓸 수 있다는 것을 알려준다.

 

그런데 $p= \frac{S}{c(b-a)}$이므로 $${\hat{p}} = \frac{N _{H}}{N}  \approx  \frac{S}{c(b-a)}$$

 

따라서 우리는 빨간색 영역의 넓이 $S$를 다음과 같이 추정할 수 있다.

 

$$S \approx  \frac{N _{H}}{N} c(b-a)= {\hat{S}}$$

 

3. 추정량과 추정치

 

시뮬레이션은 난수를 뽑아서 이루어지기 때문에 사람마다 다른 결과를 얻는다는 점이 중요하다.

 

실험을 한번 수행하면 얻는 $$\frac{N _{H}}{N} c(b-a)= {\hat{S}}$$와 $$\frac{N _{H}}{N} = {\hat{p}}$$은 사람마다 다르다.

 

쉽게 말해 100명의 사람이 위와 같은 실험으로 구해보면

 

$${\hat{S}} _{1} , {\hat{S}} _{2} , {\hat{S}} _{3} , {\hat{S}} _{4} , {\hat{S}} _{5} ,..., {\hat{S}} _{100}$$을 얻는다.

 

근데 서로 다른 이 값 모두 실제 값 $S$에 근사하므로 추정치로 사용할 수 있다.

 

즉 $\frac{N _{H}}{N} c(b-a)= {\hat{S}}$은 하나의 랜덤값이다.

 

이것을 $S$의 추정량(estimator)이라 부르고 하나의 확률분포를 갖는 확률변수가 된다.

 

그리고 어느 한사람이 실험을 통해 계산하여 얻은 하나의 실현값 ${\hat{S}} _{1}$은 추정치(estimate)

 

${\hat{S}} _{1} , {\hat{S}} _{2} , {\hat{S}} _{3} , {\hat{S}} _{4} , {\hat{S}} _{5} ,..., {\hat{S}} _{100}$은 바로 확률변수 $\frac{N _{H}}{N} c(b-a)= {\hat{S}}$의 분포에 따라 랜덤하게 분포한다

 

4. $\hat{S}$의 확률분포

 

$\hat{S}$의 확률분포는 확률변수 $\frac{N_{H}}{N} = \hat{p}$에 의해 결정된다.

 

여기서 N은 일정하므로, $N_{H}$가 확률변수이다.

 

$N_{H}$는 N개의 난수 중 HIT에 들어간 난수의 개수이다.

 

난수 1개를 뽑는 것이 하나의 사건이고, HIT를 하나의 성공, MISS를 하나의 실패라고 한다면, 

 

$N_{H}$는 N번 일어난 사건 중에서 $N_{H}$번 성공이 일어나는 경우로 이항분포(binomial distribution)를 따른다.

 

난수가 HIT에 들어갈 확률, 성공확률은 p이므로 $N_{H} ~ B(N,p)$

 

그러므로 N이 충분히 크면 $\frac{N_{H}}{N}$은 평균이 p이고 분산이 p(1-p)/n인 정규분포를 따른다.

 

 

5. 몬테카를로 시뮬레이션

 

확률변수 $\frac{N _{H}}{N} c(b-a)= {\hat{S}}$의 기댓값은 $$E( {\hat{S}} )=E( {\hat{p}} c(b-a))=pc(b-a)=S$$

 

대수의 법칙(The law of Large number)에 의해 기댓값은 표본평균에 수렴한다.

 

무슨 말이냐면 $\frac{N _{H}}{N} c(b-a)= {\hat{S}}$에서

 

랜덤하게 실현한 추정치 ${\hat{S}} _{1} , {\hat{S}} _{2} , {\hat{S}} _{3} , {\hat{S}} _{4} , {\hat{S}} _{5} ,..., {\hat{S}} _{100}$가 크기가 100인 하나의 표본이고

 

이들의 평균 $\sum _{i=1} ^{100} \frac{\hat{S}_{i}}{100}$는 ${\hat{S}}$의 기댓값 $E( {\hat{S}} )=S$에 수렴한다

 

따라서 참값 $E( {\hat{S}} )=S$의 추정을 위해서는 실험을 1번만 해서는 안된다.

 

실험을 수없이 반복하여 얻은 실현치 ${\hat{S}} _{i}$들의 평균을 구해야한다.

 

추정량을 구하는 시뮬레이션을 할 때 이런 실수를 생각보다 많이하게 된다

 

 

6. 예시 코드

 

 

원주율 $\pi$의 근삿값을 몬테카를로 시뮬레이션으로 구해보자.

 

다음과 같이 제1사분면 위에 한 변의 길이가 1인 정사각형을 올려놓고

 

정사각형안에 반지름의 길이가 1인 원의 1/4을 그려넣는다

 

그림1. 길이가 1인 정사각형과 1/4 원

그러면 $U(0,1)$에서 매 회 2개의 sample을 random하게 추출하는 것이 길이가 1인 정사각형 안에 있는 점 $(x,y)$를 뽑는 과정과 같다.

 

원의 방정식 $x^2+y^2=1$을 이용하여 x좌표를 넣어서 $y \leq \sqrt{1-x^2}$를 만족시키면 $(x,y)$는 원 안에 존재한다.

 

import numpy as np

def random_generator(n):
    
    random_number = []

    for _ in range(n):
        
        random_number.append(np.random.uniform(0,1,2))

    return random_number

def circle_function(x):
    
    return np.sqrt(1-x**2)

 

위와 같이 정사각형 내에 점 $(x,y)$를 뽑는 random_generator와 원의 방정식 $y \leq \sqrt{1-x^2}$을 나타내는 circle_function을 구현

 

in_circle = []

pi_list = []

for _ in range(1000):
    
    random_number = random_generator(100000)

    for x,y in random_number:
        
        if y <= circle_function(x):
            
            in_circle.append((x,y))

    pi = 4*(len(in_circle)/100000)

    pi_list.append(pi)

    in_circle = []

np.mean(pi_list)

3.14163476

 

random_generator로 100000개의 난수를 생성하여 random_number 리스트에 집어넣고

 

원 안에 존재하는 난수만 in_circle에 넣어서 기하학적 확률로 $\pi$를 추정한다

 

이러한 실험을 1000번 반복하여 평균을 낸 값을 구해보면 3.14163476...으로 소수점 셋째자리까지 정확하다.

 

 

내가 실수한 점은 이 실험을 1번만 했다는 것인데 통계적 추정의 원리에 의해 충분히 많이 여러번 실험해서 평균값을 구해야 한다는 점이다.

 

 

 

7. 표준오차(standard error)

 

참값 $S$의 추정량 $\hat{S}$는 $E(\hat{S}) = S$를 만족시킨다. 이를 불편추정량(unbiased estimator)이라고 부른다.

 

참값과 추정값 사이에는 분명한 차이가 있는데, 이를 오차(error)라고 부른다.

 

$$e = \hat{S} - S$$

 

많은 경우에 참값은 알 수 없는 값이니까 e를 구하는 것은 불가능하다.

 

그런데 $\hat{S}$가 확률변수이므로, e도 확률변수이다.

 

따라서, 여러번 해보면 e가 어느정도 나오리라고 기대할 수 있는데 오차 제곱의 기댓값

 

$E((\hat{S} - S)^{2})$을 Mean square for error, MSE라고 부른다.

 

이것을 최소로 하는 추정량 $\hat{S}$를 선택하는 것이 좋은 추정이다.

 

그런데 $\hat{S}$가 불편추정량 $E(\hat{S}) = S$이므로 분산의 정의에 따라 $E((\hat{S} - S)^{2}) = Var(\hat{S})$이다.

 

그러므로 추정량의 오차 $\hat{S}$를 $\hat{S}$의 표준편차로 구할 수 있다.

 

추정량의 표준편차는 표준오차(standard error)라고 부른다.

 

 

8. 구간추정

 

확률변수 $\hat{S}$는 N이 충분히 크면 중심극한정리에 의해 정규분포를 따른다.

 

$\hat{S} = N_{H}(c(b-a)/N)$의 기댓값은 $N_{H}$가 B(N,p)를 따르므로,

 

$E(\hat{S}) = S$이고, $Var(N_{H}/N) = p(1-p)/N$이니까, $Var(\hat{S}) = c^{2}(b-a)^{2}p(1-p)/N$

 

그러므로 N이 충분히 크면 정규분포로 근사시켜서 구간추정을 할 수 있다.

 

예를 들어 95% 신뢰구간은 다음과 같다.

 

 

 

 

9. 원하는 오차 이내를 얻기 위해 시행해야하는 횟수

 

추정량의 오차 $|\hat{S} - S|$가 실험자가 미리 지정한 $\epsilon$이내로 들어올 확률을 $1-\alpha$라고 하자.

 

$$P(|\hat{S} - S| < \epsilon) >= 1-\alpha$$

 

$E(\hat{S}) = S$이므로 표준편차로 나누기만 하면 표준화시킬 수 있다.

 

$$P(\frac{|\hat{S}-S|}{\sqrt{Var(\hat{S})}} < \frac{\epsilon}{\sqrt{Var(\hat{S})}}) >= 1-\alpha$$

 

그런데 $Var(\hat{S}) = c^{2}(b-a)^{2}p(1-p)/N$

 

확률 p는 0이상 1이하이므로, p(1-p)는 1/4이하이다.

 

따라서, $Var(\hat{S}) = c^{2}(b-a)^{2}p(1-p)/N <= c^{2}(b-a)^{2}/4N$

 

그러므로 $k = \frac{\epsilon}{\sqrt{Var(\hat{S})}} >= 4N\epsilon/(c^{2}(b-a)^{2}) = k'$

 

표준정규분포는 0에서 대칭이므로 다음과 같이 변형해보면...

 

$$P(Z < k) = 1-P(Z >= k) >= 1-P(Z >= k') =1-(1-P(Z < k')) = P(Z < k') >= 1-\alpha$$

 

그러므로 다음과 같은 표준정규분포의 누적분포함수 부등식으로 구해진다.

 

$$P(Z < 4N\epsilon/(c^{2}(b-a)^{2})) >= 1-\alpha$$

 

a,b,c, $\epsilon$, $\alpha$은 모두 정해진 상수이고 N만 미지수이다.

 

표준정규분포의 누적확률분포로부터 N을 구할 수 있다.