몬테카를로 시뮬레이션 (5) - 적분 계산하기

몬테카를로 시뮬레이션 (1) - 파이 계산하기

몬테카를로 시뮬레이션 (2) - 온라인 판매자 예제

몬테카를로 시뮬레이션 (3) - 부트스트랩 알고리즘 (bootstrap)

몬테카를로 시뮬레이션 (4) - 범용 몬테카를로 엔진

---

1차원 적분

몬테카를로 시뮬레이션을 이용해서 충분히 매끈한 곡선을 가지는 피적분 함수에 대해 적분을 계산할 수 있습니다. (충분히 매끈하다는 것은 수학적으로 무한차원의 미분이 가능하다는 것이지만 여기서는 1차원 미분이 가능하다고 가정하겠습니다.) 먼저 $[a,b]$ 사이에서 정의된 $f(x)$ 적분은 다음과 같이 정의됩니다.

$I = \int_a^b f(x) dx$

몬테카를로 적분의 시작은 다음 식들을 만족하는 두 개의 함수 $g(x), p(x)$를 정의하는 것입니다.

$p(x) = 0, x\in [-\infty, a] \cup [b, \infty]$

$\int_{-\infty}^{\infty} p(x) dx = 1$

$g(x) = f(x)/p(x)$

따라서 $p(x)$를 확률질량함수로 해석해서 $g(X)$에 대한 기대값으로 $f(x)$의 적분을 유도할 수 있는데요, 여기서 $X$는 분포가 $p(x)$인 확률 변수이며 $p(x)$는 $[a,b]$ 영역에서만 0이 아닙니다.

$E[g(X)]=\int_{-\infty}^{\infty} g(x)p(x) dx=\int_{a}^{b} f(x) dx = I$

위 수식에서 가장 간단한 방법은 $X$가 $[a, b]$에서 정의된 균등 확률 변수가 되게 하는 방법으로 정의하는 것으로, $[a, b]$에서 $p(x)=1/(b-a)$, $g(x) = (b-a)f(x)$와 같이 됩니다.

이제 $[a,b]$ 영역 안에서 균등 분포인 $N$개의 무작위 점 $x_i$를 만들어내고 피적분 함수 $f$ 값을 각 점에서 계산한 이후에 평균을 낸 다음, 영역의 크기 $(b-a)$를 곱해서 적분을 계산할 수 있습니다.

$I=E[g(X)]=\frac{1}{N}\Sigma_{i} g(x_i)$

이를 범용 몬테카를로 엔진 클래스를 이용해 구현하면 다음과 같습니다.

class MCIntegrator(MCEngine):
    def __init__(self, f, a, b):
        self.f = f
        self.a = a
        self.b = b

    def simulate_once(self):
        a, b, f = self.a, self.b, self.f
        x = a + (b - a) * random.random()
        g = (b - a) * f(x)
        return g


def main():
    s = MCIntegrator(f=lambda x: math.sin(x), a=0, b=1)
    print(s.simulate_many())

main()

 

N 차원 적분

이를 $N$ 차원으로 쉽게 확장할 수 있습니다. $\Omega$는 적분이 정의된 $N$ 차원 영역이며 1차원 적분 떄와 마찬가지로 $p(x), g(x)$를 정의합니다.

$I=\int_{\Omega} f(x_1, ..., x_n) dx_1...dx_n$

$\int p(x_1, ..., x_n)dx_1...dx_n=1$, $p(x_1, ..., x_n)=0 for \notin \Omega$

다음 수식에 대해 몬테카를로 적분을 구할 수 있습니다. (1차원 적분에서 $g(x)$를 구할 때 영역의 넓이 $b-a$를 곱해주었듯이 여기서는 $[0,1]\times[0,1]\times[0,1]\times[0,1]$의 면적인 1을 곱해줍니다.

$I=\int_0^1dx_1\int_0^1dx_2\int_0^1dx_3\int_0^1dx_4 sin(x_1+x_2+x_3+x_4)$

class MCIntegrator(MCEngine):
    def simulate_once(self):
        volume = 1.
        while True:
            x = [random.random() for d in range(4)]
            if sum(xi**2 for xi in x) < 1:  # p(x)의 적분값이 1
                break
        return volume * self.f(x)

s = MCIntegrator()
s.f = lambda x: math.sin(x[0] + x[1] + x[2] + x[3])
print(s.simulate_many())