몬테 카를로 시뮬레이션(Monte Carlo Simulation) 의 이해 : 원주율값 구하기 (+파이썬 시뮬레이션 코드)

몬테카를로 시뮬레이션(혹은 알고리듬)이란 무엇인지를 정의하기란 매우 어려운데, 그 이유는 이 방법론을 수학과 응용수학 분야에서 매우 광범위하게 사용하기 때문입니다. 

 

(1 )이 시뮬레이션을 이용하여 계산 하려는 대상(값)

(2) 이 시뮬레이션을 이용해야 하는 이유

(2) 이 시뮬레이션의 전체 혹은 일부를 이용하는 방법

 

등에 따라서, 몬테 카를로 시뮬레이션에 대한 정의는 조금씩 다를 텐데요, 그러나 이 방법론을 사용하는 모든 경우를 아우를 수 있는 핵심점인 키워드는 "반복적인 무작위적 샘플링" 입니다. 위키피디아의 몬테 카를로 시뮬레이션을 소개하는 첫 문장은

 

Monte Carlo methods, or Monte Carlo experiments, are a broad class of computational algorithms that rely on repeated random sampling to obtain numerical results. The underlying concept is to use randomness to solve problems that might be deterministic in principle. 

 

입니다. 우리 말로 적당히 번역하면 

 

광범위한 계산 알고리듬의 일종인데, 특히 반복적인 랜덤 샘플링에 의존하여 결과를 얻는 것. 문제를 풀기 위해 (그 문제가 결정론적인 원리라고 하더라도) 무작위성을 이용한다

 

가 중요한 부분이 될 것 입니다. 보통 알고리듬이라고 하면, "주어진 절차에 따라서 기계적인 방식으로 문제를 해결" 하는 것을 떠올리는데 여기에는 무작위성이 파고 들어갈 공간이 없습니다. 그러나 몬테 카를로 시뮬레이션에서는 무작위성이 핵심적인 역할을 하게 됩니다. 

 

이렇게 저렇게 설명을 길게 하는 것 보다는 간단하긴 하지만 핵심적인 시뮬레이션을 하나 해 보면서 몬테 카를로 시뮬레이션이 무엇인지를 "느껴" 보는 것이 더 중요한 것 같습니다. 몬테 카를로 시뮬레이션 설명에서 꼭 나오는 "원주율값 구하기"를 해 보면서 몬테 카를로 시뮬에이션이 어떻게 사용될 수 있고, 기존의 다른 알고리듬 대비 얼마나, 혹은 어떻게 그 효율성을 줄 수 있는지를 알아 보도록 하겠습니다. 

 

원주율값 $\pi$ = 반지름이 1인 원의 넓이

 

라는 것에 착안하여 반지름이 1인 원의 넓이를 시뮬레이션의 방법을 통해 정확히 구하면 원주율 $\pi$를 구할 수 있습니다. 넓이를 구하는 것이니, 적분을 하면 될텐데요. 반지름의 길이가 1인 원의 방정식은 $x^2 + y^2 =1$ 이므로, 원이 넓이는 $x^2 + y^2 <1$ 인 영역의 넓이를 구하는 것과 같습니다. 이 식을 $y$에 대한 식으로 바꾸면 $y = \sqrt{1-x^2}$ 이므로, 원의 넓이 이면서 원주율 값은

$$S = \pi = 4\int_{0}^{1}\sqrt{1-x^2}dx$$

이 됩니다. 1사분면에서 넓이를 계산하고 여기에 곱하기 4를 해 준것 입니다. 이 적분은 역삼각함수를 이용하면 부정적분을 구할 수 있고, 간단하게 값이 $\pi$가 된다는 것을 알 수 있습니다. 

 

따라서, 시뮬레이션적인 방법을 통해서 원의 넓이 (혹은 1사 분면에서의 원의 넓이)를 구하면 원주율 값을 구할 수 있다는 결론을 얻습니다. 즉, 적분을 해야합니다. 

 

무작위로 점을 찍는 방법을 통해 넓이를 구하자

1사분면에 원을 그렸습니다. 원의 내부는 빨간색, 원의 외부는 파란색으로 칠했습니다. 우리의 목표는 빨간색 부분의 넓이를 구하는 것 입니다. 물론 구분구적법을 통해서 구해도 되지만, 이번에는 다른 방법을 통해서 구해보도록 하겠습니다. 

 

가로 세로의 길이가 1인 정사각형안에 한 점을 무작위적으로 고른다면, 그 점이 붉은색이 될 확률과 파란색이 될 확률은 붉은색 영역의 넓이와 푸른색 영역의 넓이의 비율과 같을 것 입니다. 그런데 여기서 전체 넓이는 1이기 때문에, 

 

붉은색 영역이 선정될 확률 = 붉은색 영역의 넓이 = 반지름이 1인 원의 넓이의 $1/4$ = $\pi/4$

 

가 됩니다. 즉, 우리가 무작위적으로 한 점을 선택했을 때, 이 점이 붉은색 영역에 속한 점이 될 확률을 정확히 원주율이 $1/4$이 됩니다. 따라서, 만일 우리가 무작위적으로 점을 선택(찍는) 과정을 반복한다면 확률의 큰수 법칙에 따라서 전체 점의 갯수 중에서 빨간색 영영역의 점의 갯수의 비율은 정확히 $\pi/4$가 될 것 입니다. 이를 수식적으로 정확하게 표현하면, 

 

$$\pi/4 = \frac{\text{선택된점이 빨간색 영역에 속하는 횟수}}{\text{무작위적으로 한 점을 고르는 횟수}}$$

 

가 됩니다. 여기에서의 가정은 무작위적으로 고르는 시행이 가로 세로 길이가 1인 정사각형의 모든 점을 정말로 무작위적이고 고르게 뽑아야 한다는 것 입니다. 특정 영역에 점을 많이 뽑거나 혹은 반대로 특정 영역의 점을 적게 뽑는다면 위 방법을 통해서 계산한 원주율 값은 우리가 알고 있는 원주율값과 많이 다를 것 입니다. 또한, 위 과정을 매우 많이 반복해야만 실제 원주율 값과 비슷한 값을 얻을 수 있을 것 입니다. 

 

프로그래밍을 통해서 무작위로 점을 찍어 보자

파이썬 프로그래밍을 통해서 위 영역에서 100개의 점을 무작위로 선택했고, 이 점이 원이 내부에 있는경우에는 빨간색으로, 원의 외부에 있는 경우에는 파란색으로 표시하였습니다. 완전히 고르게 퍼져 있는 것은 아니지만, 점들이 나름대로는 넓게 잘 펼쳐져 있는 것 같습니다. 위의 구체적인 예시에서는

 

빨간점 : 81개

파란점 : 19개

 

를 얻었고, 따라서 위 식에 대입하면 $\pi = 4 \frac{81}{100} = 3.24$ 값을 얻게 됩니다. 실제 $\pi$값과는 좀 다르지만 그런대로 괜찮은 값을 얻은 것 같습니다. 앞서 언급한대로 시행의 횟수를 키워야 좀 더 정확한 넓이의 비율을 구할 수 있기 때문에, 시행의 횟수N 를 1000, 10000, 100000으로 늘려 보겠습니다. 

영상으로 보면 더 이해하기 쉬운데요,

https://youtu.be/SV41AeVk_Oc

결과는 위와 같습니다. 각 경우 구해진 $\pi$값은

 

N = 1000 : 3.12

N =10000 : 3.1308

N =100000 : 3.137

 

로 N이 증가하면 증가할 수록 원주율값과 차이가 줄어듦을 알 수 있습니다. 

 

위 시뮬레이션에서는 각 시행 횟수 마다 단 한 번의 실험을 한 것 인데요, 이걸 여러번 반복해 보면 좀 더 정확한 결과를 얻을 수 있지 않을까요?

위 그래프는 각 시행횟수 N 마다 100번의 반복된 실험을 한 결과를 보여주고 있습니다. X축은 1, 2, 3, 4, 5 값을 갖는데 각각 시행 횟수 N = 10, 100, 1000, 10000, 10000 를 의미합니다. 10의 거듭 제곱의 값 입니다. Y 축은 각 시행에서 얻어진 원주율 값 입니다. 

 

이 결과를 보면 "N이 커질 수록 각 시행의 값이 실제 원주율값에 가깝에 얻어지고, 각 시행의 편차가 작아진다" 를 확인할 수 있습니다. 이는 확률과 통계에 대한 우리의 기본적인 직관과 일치합니다. 더 많은 시행을 할 수록 주어진 확률에 수렴하는 것 입니다. 

 

실제로 이 시행은 이항 분포를 통해서 설명할 수 있습니다. 한 점을 무작위로 선택했을 때, 빨간점 혹은 파란점이 선택되면서 (우리가 관심이 있는) 빨간색이 될 확률을 $\pi/4$가 됩니다. 그래서 시행횟수를 N이라 했을 때, 빨간점이 나오는 확률은 이항 분포 $B(N, \pi/4)$가 됩니다. 따라서 빨간점의 기대값은 $N \frac{\pi}{4}$ 이고, 표준편차는 $\sqrt{N \frac{\pi}{4} (1-\frac{\pi}{4} )}$가 됩니다. 최종적으로 구하고 싶은 원주율값은 빨간점의 기대값에 4를 곱하고 N으로 나누는 것이니 시뮬레이션을 통해서 구해지는 원주율값의 표준편차는 $\frac{1}{\sqrt{N}}$에 비례하게 됩니다. 각 시행횟수 N 마다 100번씩 반복 실험한 위 결과가 이항 분포의 성질에 부합됨을 알 수 있습니다.

 

요약하면, 위와 같이 시행횟수 N을 매우 크게 하면, 원주율값과 매우 유사한 결과를 얻을 수 있고, 이때의 에러(실제 원주율 값과 시뮬레이션을 통해 얻은 원주율 값이 차이)는 $\frac{1}{\sqrt{N}}$에 비례하게 됩니다. 

 

넓이는 그냥 구분 구적법을 통해서 구하면 되는게 아닌가?

 

위 방법을 통해 원주율 값을 계산한 과정을 다시 한 번 살펴보면, 구분구적법을 통해서 넓이를 구하지 않고 무작위적인 점찍기 방법을 통해서 붉은색 영역의 넓이를 구했다는 것을 발견할 수 있습니다. 앞서 언급했지만, 붉은 영역의 넓이를 정확하게 구하기 위해서는 구분구적법을 하면 됩니다. 구분구적법은 매우 간단한 알고리듬이라서 설명은 생략하고 바로 결과만을 제시하면, 

위 그래프와 같습니다. 구분구적법에서는 [0, 1] 사이를 많은 조각으로 나누는데, 위 그래프의 x축은 몇 조각으로 나누었는지에 대한 값입니다. 앞의 그래프에서와 같이 위 값은 10의 거듭제곱값으로 1, 2, 3, 4, 5는 각각 10, 100, 1000, 10000, 100000 조각으로 나누어서 구분구적법을 계산한 것 입니다. 

 

100000 조각으로 나누어서 계산한 원주율 값은 3.141612로 우리가 알고 있는 원주율값 3.141592와 매우 유사합니다.

 

이처럼 구분구적법을 통해서 계산을 하면, 몬테 카를로 시뮬레이션을 통해 계산을 하는 것에 비해서

 

(1) 프로그래밍 코드를 만들기 쉬움

(2) 같은 시행 횟수를 하더라도 몬테 카를로 시뮬레이션 대비 더 정확한 계산을 할 수 있음

    몬테 카를로 계산으로 100000번 무작위 선정을 한 경우 원주율 값 : 3.137

    구분구적법으로 100000번 나누어 계산한 원주율 값 : 3.142

   

와 같은 장점이 있습니다. (2)를 조금 더 상세하게 설명하면, 구간 [0, 1]를 N 조각 하여 계산하는 경우, 구분구적법과 몬테 카를로 시뮬레이션 모두 N에 비례하는 계산량이 필요합니다. 이때 에러는 구분구적법의 경우 $\frac{1}{N}$에 비례하게 되고, 몬테 카를로 시뮬레의 경우 $\frac{1}{\sqrt{N}}$에 비례 하게 됩니다. N이 커지면 구분구적법의 에러가 훨씬 더 작아지는 것 입니다. 그렇기 때문에 위와 같은 결과가 얻어진 것 입니다. 

 

위 (1), (2)를 보면 구분구적법 대비 몬테 카를로 시뮬레이션의 장점이 없는 것 같은데, 굳이 왜 몬테 카를로 시뮬레이션을 하는 것 일까요?

 

차원이 높아진다면 어떨까?

 

우리가 다루는 원주율값을 구하는 구체적인 문제에서는 적분을 해야하는 차원은 1차원이었습니다. 즉 $y = f(x)$에서 $x \in R^1$, $y \in R^1$ 이었습니다. 만일 적분을 해야하는 차원이 1차원이 아닌 일반적인 $D$ 차원이라면 (즉, $x$가 실수가 아닌 일반적인 $D$ 차원의 벡터) 어떻게 될까요? 동일한 정도의 에러를 갖는 답을 얻기 위해서, 구분구적법과 몬테 카를로 시뮬레이션에서 몇 번의 계산을 해야하는지를 알아 보겠습니다. 

 

구체적인 예를 들어 $1/1000 = 10^{-3}$ 만큼의 에러를 갖는 값을 계산해야 한다고 생각하겠습니다. 몬테 카를로 계산의 경우 N번 무작위 선택 계산을 할 때, 오차가 $\frac{1}{\sqrt{N}}$에 비례 하기 때문에 무작위 선택 계산을 $10^6$ 번 해야합니다. 구분구적법의 경우 각 차원의 영역을 N 등분 할 때, 오차가 $\frac{1}{N}$이 되기 때문에 영역 등분을 $10^3$ 만큼 해야합니다. 따라서 차원이 $D$라면, 전체 공간에서 그리드의 갯수는 $10^{3D}$ 가 되고, 계산량은 이에 비례해 증가 합니다. 정리하면

 

몬테 카를로 : $10^{6}$, 즉 변화가 없음

구분구적법 : $10^{3D}$, 차원 $D$에 지수적으로 증가

 

가 됩니다. 1차원에서는 구분구적법이 유리하고(계산량이 더 적고), 2차원에서는 동일하며, 3차원 이상에서는 몬테 카를로 시뮬레이션이 더 유리합니다. 이 점이 매우 중요합니다! 몬테 카를로 시뮬레이션의 이 같은 장점은 $x$의 차원이 커지면 커질수록 극명해 지는데요, 만일 $D=10$ 이라고 한다면 몬테 카를로 시뮬레이션의 계산량은 $10^{6}$ 로 변함이 없지만, 구분구적법의 경우 $10^{30}$이나 됩니다. 현존하는 슈퍼 컴퓨터를 동원하더라도 계산이 불가능할 정도로 계산량이 늘어 납니다. 물론 좀 더 효율적인 구분구적법 알고리듬을 사용한다면 계산량을 줄일 수 있습니다. 하지만, 차원이 높아짐에 따라서 계산량이 지수적으로 증가하는 현상은 피할 수 없습니다.

 

정리

 

위 구체적인 예시, 그리고 구분구적법과의 비교를 통해서 몬테 카를로 시뮬레이션이 어떻게 수행 되는지, 그리고 다른 시뮬레이션(알고리듬)과의 차이가 무엇인지를 알아보았습니다. 요약하면 몬테 카를로 시뮬레이션은 

 

(1) 무작위 샘플링(위 예시에서는 무작위로 점을 찍는 것)의 반복을 통해서 원하는 값을 얻는다

(2) 구분구적법과 같이 모든 경우를 전부 계산 하는 방법과 비교하면, 공간 차원이 커질 수록 몬테 카를로 계산이 우월해 진다 (동일한 오차의 값을 얻기 위한 계산량이 적어 진다)

 

의 성질을 갖습니다. 

 

추가적인 포스팅을 통해서, 구체적인 활용에서 위 두 성질이 어떻게 이용되는지를 알아 보도록 하겠습니다. 

 

파이썬 시뮬레이션 코드

 

아래는 위 계산에서 사용된 파이썬 시뮬레이션 코드 입니다. 파이썬을 활용할 줄 안다면, 아래 코드를 복-붙 하여 사용할 수 있습니다. 

# 필요한 라이브러리 불러오기 
import numpy as np
from multiprocessing import Pool
from matplotlib import pyplot as plt

if __name__ == '__main__':
    N = 10000 ### 무작위시행 횟수 정의
    x = np.random.random([N, 2])
    distance = np.sum(x ** 2.0, axis=1)
    in_out = distance <= 1.0
    pi = np.sum(in_out)*4/N ### Pi 값은 천제 시행에서 원 안에 있는 점의 갯수로 정해짐
    color = list(map(lambda x: 'red' if x else 'blue', in_out)) ### 원의 안, 밖에 따른 색상 설정

    plt.figure(figsize=(5, 5)) ### 그림 사이즈
    plt.scatter(x[:,0], x[:,1], color = color, s=5, label ='Result : {}'.format(np.round(pi, 4)))
    
    cx = np.cos(np.linspace(0, np.pi/2, 1000))
    cy = np.sin(np.linspace(0, np.pi/2, 1000))
    plt.plot(cx, cy, color = 'black', lw =2) ### 원의 경계를 그려주는 부분
    plt.legend(loc = 'lower right')

    plt.xlim(0, 1)
    plt.ylim(0, 1)
    plt.show()

위 코드는 무작위 시행을 N 번하는 코드 입니다. 

 

다음 포스팅에서는 좀 더 일반적인 문제를 파이썬 프로그래밍을 통해서 풀어 봅니다.

 

https://studyingrabbit.tistory.com/34

파이썬을 이용한 몬테 카를로 적분(Monte Carlo Integration) 구현