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

이 포스팅에서는 지난 시간에 설명한 몬테 카를로 시뮬레이션 방법을 프로그램 언어 파이썬으로 구현해 보도록 하겠습니다. 적분은 몬테 카를로 시뮬레이션을 이용해서 할 수 있는 것의 가장 기본적인 어플리케이션입니다. 사실, 몬테 카를로 시뮬레이션은 "어려운 적분을 하기 위한 계산 방법론 이다" 라고 말해도 과언이 아닐 정도로 몬테 카를로 시뮬레이션의 적분에 대한 응용은 기본적이면서도 활용 범위가 넓습니다. 

 

몬테 카를로 적분을 통한 계산 값이 정확하기 위해서는 몬테 카를로 시뮬레이션에서 하는 랜덤 샘플링의 횟수가 매우 커야 하는데요, 따라서 컴퓨터로 계산할 수 밖에 없습니다. 컴퓨터 프로그램 언어를 이용해서 원하는 알고리듬을 코딩 언어로 구현해야하는데, 프로그램 언어 중에서 배우기 쉽고, 활용 범위가 넓은 파이썬을 통해서 몬테 카를로 적분을 구현해 보도록 하겠습니다. 

 

몬테 카를로 적분이란?

 

지난 포스팅에서 몬테 카를로 "시뮬레이션"의 일반에 대해서 설명하긴 했지만, 사실 다양한 몬테 카를로 시뮬레이션 중에서 가장 간단하면서도 핵심적인 몬테 카를로 "적분"을 설명하였습니다. 원주율값 $\pi$를 계산하기 위해 반지름이 1인 원의 넓이를 계산하였습니다. 넓이를 계산한 것이니 당연히 적분입니다. 지난 포스팅을 요약하면 몬테 카를로 적분의 과정은

 

(P-1) 원하는 적분 영역을 포함하는 사각형 $R$을 정의하고, 원하는 영역의 안 $S$과 밖 $S_{out}$을 정의한다.

(P-2) 이 사각형에 균일하게 랜덤 샘플링을 반복한다.

(P-3) 원하는 영역의 넓이 $A$는 $A = R \frac{N_{in}}{N}$ 으로 주어진다. 여기서 $R$은 사각형의 넓이, $N$은 랜덤 샘플링의 갯수, $N_{in}$은 랜덤 샘플링 중, 원하는 영역안에 속하는 샘플링의 갯수

 

입니다. 샘플링이 삼각형 영역 안에서 균일하고, 랜덤 샘플링 횟수 $N$이 커감에 따라서 몬테 카를로 적분을 통해 얻은 값과 실제 적분 값의 차이는 $\frac{1}{\sqrt{N}}$에 비례해 줄어들게 됩니다. 

 

구체적인 예를 들어서 설명해 보자

우리가 계산하고자 하는 영역은 위 그래프에서 파란색으로 색칠한 영역입니다. 

 

(1) $y = f(x) = x^2 + \sin(\pi x)$

(2) $y = 0$

(3) $x = 1$

 

로 둘러쌓인 영역입니다. 간단히 적분식으로 표현하면 $A = \int_{0}^{1} x^2 + \sin(\pi x) dx$ 가 됩니다. 기본적인 다항함수와 삼각함수의 합이기 때문에 해석적으로 적분할 수 있고, 그 답은 $A = \frac{1}{3} + \frac{2}{\pi} = 0.96995$ 가 됩니다. 

 

이 문제에 대해서 위에서 몬테 카를로 적분의 과정 P-1, P-2, P-3 을 차례로 해 보도록 하겠습니다. 

 

P-1 : 원하는 적분 영역을 포함하는 사각형을 정의하고, 원하는 영역의 안과 밖을 정의한다.

 

위 영역을 포함하는 직사각형을 잡으면 됩니다. 원리적으로는 아무런 모양의 직사각형을 잡아도 되지만, 계산의 효율성을 위해서는 해당 영역을 포함하는 최소 넓이의 직사각형을 선택하는것이 좋습니다. 직사각형의 가로의 길이는 자명하게 1이 되고, 세로의 길이는 $f(x)$의 최대값으로 하면 됩니다. $f(x)$의 최대값을 정확하게 구하기는 어려우니 (사실은 어렵다기 보다는 귀찮으니), 최대값 보다 크면서 적당히 작은 1.4를 세로의 길이로 하도록 하겠습니다. 가로의 길이 $W= 1$, 세로의 길이 $H = 1.4$, 직사각형의 넓이 $R = W H = 1.4$가 됩니다. 

그래프로 표현하면 이와 같습니다. 

이때, 적분을 하려는 영역 $S$(하늘색)과 적분을 하려는 영역의 외부(빨간색) $S_{out}$으로 나눌 수 있습니다.  

 

P-2, P-3 : 이 사각형에 균일하게 랜덤 샘플링을 반복한다. $A = R \frac{N_{in}}{N}$

 

과정 P-2와 과정 P-3을 두 과정으로 분리하긴 좀 애매해서 한 번에 하도록 하겠습니다. 

 

랜덤 샘플링을 한 점의 좌표를 $(x, y)$라고 하겠습니다. 이 점이 $S$에 속하는지 $_{out}$에 속하는지는 $y$값과 $f(x)$값의 비교를 통해 알 수 있습니다. $S$는 $y=f(x)$의 아래 부분이니,

 

(1) 만일 $y<f(x)$ 이면, $(x, y) \in S$, 

(2) 그렇지 않으면, $(x, y) \in S_{out}$

 

가 됩니다. 

 

위에서 설명한 것과 같이 적분하려는 영역을 포함하는 넓은 영역에서 균일하게 샘플링을 하고, 선택된 점 중에서 원하는 영역에 속해 있는 점을 고르는 과정을 기각 샘플링(rejection sample)이라고 합니다. 기각 샘플링을 이용하여 보다 복잡한 분포를 갖는 확률 샘플링을 할 수 있는데, 이 포스팅에서 사용한 기각 샘플링은 그 중에서 가장 간단한 형태의 기각 샘플링 입니다. 

 

이번 포스팅의 목표는 파이썬 프로그래밍을 통해서 실제로 위 값을 계산하는 것이니, 이 부분에 해당하는 파이썬 코딩을 해 보도록 하겠습니다. 다양한 방법으로 구현할 수 있지만, 나름 간단하게 구현하는 방법은 아래와 같습니다. 

import numpy as np
from  matplotlib import pyplot as plt

def function(x): ### 함수 정의
    return x**2 + np.sin(np.pi*x)

if __name__ == '__main__':
    N = int(1e+4)  ### 랜덤 샘플링 시행 횟수 정의
    W, H = 1, 1.4  ### 랜덤 샘플링을 할 사각형 R의 가로 세로 정의

    X = np.random.random(N)   ### 각 점의 x 좌표 랜덤 샘플링
    Y = H*np.random.random(N) ### 각 점의 y 좌표 랜덤 샘플링

    F = function(X) ### 각 랜덤 샘플링에 대한 f(x)를 계산
    in_out = Y < F  ### y와 f(x) 값 비교
    
    A = H * W * np.sum(in_out) / N ### 영역 S의 넓이
    print("A = ", A)

매우 짧은 코딩으로 몬테 카를로 적분을 구현하였습니다. 조금 더 효율적으로 짠다면 더 짧게 짤 수도 있지만 설명의 편의를 위해서 여러줄에 걸쳐 표현하였습니다. 

 

파이썬의 numpy 라이브러리를 이용하여 구현하였습니다. np.random.random()을 활용하여 N 개의 랜덤 샘플 포인트 집합 $\{(x_i, y_i) | i = 1, 2, 3, ..., N\}$ 을 생성하였습니다. F = function(X)를 통해서 $y_i$ 값을 계산해 주고, $y_i, f(x_i)$의대소 관계를 얻습니다. 우리가 원하는 영역은 $y < f(x)$ 이니, 이 경우에 in_out값은 True가 됩니다. 

 

마지막 줄은 영역 $S$의 넓이 $A$를 구하는 부분입니다. 사각형의 넓이와 전체 샘플링 포인트 중에서 원하는 영역에 속하는 샘플링 포인트의 비율의 곱이 됩니다. 간단한 파이썬 코드이니 실제로 구현을 해 보고 싶으신 분은 위 코드를 복-붙 하신 다음에 파이썬을 실행하면 됩니다. 

 

여기서 계산을 끝내도 되지만, 시각화를 할 수 있는 이미지를 출력하면 더 좋을 것 같습니다. 파이썬에서 자주 활용하는 그래프 그리는 라이브러리인 matplotlib을 통해서 그래프를 그려 보도록 하겠습니다. 

    color = list(map(lambda x: 'dodgerblue' if x else 'red', in_out))
               ### 샘플링 포인트의 색깔 설정
               
    plt.scatter(X, Y, color=color, s=5, label='A = {}'.format(np.round(S, 4)))
               ### 샘플링 포인트 그리기 
               
    plt.plot([0, W], [0, 0], color='black')  ### 사각형 영역 그리기
    plt.plot([W, W], [0, H], color='black')
    plt.plot([0, W], [H, H], color='black')
    plt.plot([0, 0], [0, H], color='black')

    px = np.linspace(0, W, 1000)
    py = function(px)
    plt.plot(px, py, color = 'black') ### 함수 y = f(x) 그리기

    plt.legend(loc = 'lower right') ### 그래프 그리는 옵션 추가
    plt.xlabel('x')
    plt.ylabel('y')
    plt.grid()
    plt.show()

in_out 값이 True이면 하늘색으로, False이면 빨간색으로 점을 색칠할 것 입니다. 그 이후 부분은 matplotlib의 일반적인 사용법이라 사실 몬테 카를로 적분의 구현과는 큰 상관이 없습니다. 

 

위 전체 프로그램을 실행하면 $N= 4000$번의 샘플링을 통해서 영역 $S$의 넓이 $A$를 구하게 됩니다. 

결과는 위와 같습니다. 1,000번 정도 샘플링을 하면 선택된 점은 전체 사각형 영역에 고르게 분포하게 됩니다. 이 때 적분값은 0.9796으로 실제 정확한 값 0.96700 대비 오차는 0.0126, 비율로 하면 1.3% 가 됩니다. 

 

샘플링 횟수 $N$의 변화에 따른 적분값의 변화를 구해보도록 하겠습니다. 파이썬의 문법 For문을 사용하여 구현할 수 있는데요, 

if __name__ == '__main__':
    N = int(1e+4)
    W, H = 1, 1.4

    As = []
    for d in [1, 2, 3, 4, 5, 6]: ### 10^{1, 2, 3, 4, 5, 6}에 대한 계산 반복
        N = 10**d
        X = np.random.random(N)
        Y = H*np.random.random(N)

        F = function(X)
        in_out = Y < F
        A = H * np.sum(in_out) / N ### 넓이를 계산
        As.append(A)  ### 계산된 넓이를 As에 저장

    plt.scatter([1, 2, 3, 4, 5, 6], As) ### d에 따른 A값을 그리기
    plt.xlabel('Exponent')
    plt.ylabel('A')
    plt.axhline(0.967, color='gray', ls='--') ###참값 표시
    plt.grid()
    plt.show()

와 같이 간단히 구현할 수 있습니다. 앞에서 설명한 부분과 파이썬의 기본 문법인 For와 list 변수의 append를 사용했습니다. 

결과는 위와 같습니다. 랜덤 샘플링의 횟수가 증가할 수록 몬테 카를로 적분값은 참값 (가로 점선)에 가까워 집니다. 무작위성이 들어가는 계산이기 때문에, 계산을 수행할 때 마다 얻어지는 값은 다릅니다. X축은 선형이 아니라 로그 스케일입니다. 

 

모든 것은 일반화 되어야 한다. 

 

몬테 카를로 적분 방법을 이용하여 1차원 함수 $y = f(x)$를 적분하는 계산을 해 보았습니다. 당연히 좀 더 복잡한 형태의 적분도 가능한데요, 예를들어서 $z = f(x, y)$가 주어졌을 때, $A = \int_{0}^{1} \int_{0}^{1} f(x, y) dxdy$ 와 같이 2차원 적분도 가능합니다. 이 경우에도 동일한 방법으로 부피 적분을 할 수 있는데요, 위에서 설명한 방식에서 "사각형 $\rightarrow$ 직육면체"로만 바꾸어 생각하면 됩니다. 

 

다음 포스팅에서는 효율적인 샘플링을 할 수 있는 방법 중에 하나인 메트로폴리스-헤이스팅스 알고리듬에 대해서 알아 봅니다. 

 

https://studyingrabbit.tistory.com/35

메트로폴리스-헤이스팅스 알고리듬을 활용한 분포 생성 및 이를 이용한 몬테 카를로 적분