몬테카를로 시뮬레이션 : 함수의 적분값 구하기 (1)

지난 포스팅에서 몬테카를로 시뮬레이션을 이용하면 손으로 계산하기 어려운 작업을 쉽게 할 수 있다고 말씀드렸습니다. 대표적으로 어떤 함수를 수치적으로 적분할 때 응용할 수 있기 때문에, 베이지안 통계학에서 많이 사용합니다. 이번 시간에는 $exp(-x^2)$ 함수를 수치적으로 적분하면서 그 원리를 설명하도록 하겠습니다. $exp(-x^2)$함수는 부정적분을 원시함수로 표현할 수 없는 함수 중 하나입니다.

수학적 근거

먼저 어떤 함수 $g(u)$를 구간 $[0, 1]$ 사이에서 적분해보겠습니다. 이 결과를 $\theta$라고 하겠습니다.

$\theta = \int_{0}^{1}{g(u)du}$

이 식을 수치적으로 적분하려면, 기대값과 강대수의 법칙(Strong Law of Large Number)을 사용합니다. 즉, 다음처럼 생각할 수 있습니다.

확률변수 $U$가 표준균일분포를 따른다고 하면 밀도함수 $f_{U}{(u)}$는 $[0, 1]$ 사이에서 $1$입니다. 그러므로 위의 식을 아래처럼 쓸 수 있습니다.

$\theta = \int_{-\infty}^{\infty}{g(u)f_{U}{(u)}du}=\int_{0}^{1}{g(u)f_{U}{(u)}du} = E(g(u))$

그리고 $U$가 확률변수이므로 확률변수들의 함수인 $g(u)$도 확률변수입니다. 따라서 강대수의 법칙에 의해 다음이 성립합니다.

$\lim_{n \to \infty} {{g(u_{1}) + g(u_{2}) + ... + g(u_{n})}\over {n}} = E(g(u))$

핵심

대수의 법칙에 의해서 표본의 크기가 커질수록 표본평균은 모평균에 근접하기 때문에, 표본을 많이 생성할수록 수치적으로 구한 표본평균은 우리가 원하는 $\theta = E(g(u))$에 근접하게 됩니다. 즉, 확률변수 ${U}$를 ${n}$개 생성한 다음에 이들을 함수 $g$에 입력하고 그 결과를 평균하면 됩니다. 그럼 어떤 분포를 따르는 확률변수를 생성해야 할까요? 통상적으로 관심 있는 확률변수를 ${X}$로 표현하는데 계속 ${U}$라고 표현한 것을 찾으신 분들은 표준균일분포를 따르는 난수 $U$들의 배열을 생성하면 된다는 것을 파악하셨을지도 모르겠습니다.

예제

$exp(-x^2)$을 구간 $[0, 1]$에서 직접 구한 결과와 파이썬 커뮤니티에서 사용하는 sympy 패키지에서 동일한 연산을 수행한 결과를 비교하도록 하겠습니다.

# integral of exp(-x^2) from 0 to 1
random_number = np.random.uniform(low = 0, high = 1, size = 100000)
exp_x_square = np.exp(-np.square(random_number))
manual_solution = np.mean(exp_x_square)

scipy_solution = integrate.quad(lambda x : math.exp(-x**2), 0, 1)
print("manual : {}, scipy : {}".format(manual_solution, scipy_solution))
# manual : 0.746285643796329, scipy : (0.7468241328124271, 8.291413475940725e-15)

오차가 0.0006정도 나기는 하지만, 컴퓨터에서는 이런 식으로 적분하는 것을 알 수 있습니다.

Next

구간 $[0, 1]$ 사이에서 임의의 함수를 적분하는 방법을 알아봤습니다. 그렇다면 구간이 $[a, b]$로 나오는 함수들은 어떻게 적분해야 할까요? 하나의 방법은 이미 배운 $[0, 1]$ 사이에서 함수 적분하는 방법을 사용하기 위해 $[a, b]$를 $[0, 1]$로 변환하면서 함수도 적절한 변환을 해 주고 위의 테크닉을 사용하면 됩니다. 그 방법은 다음 포스팅에서 작성하도록 하겠습니다. 감사합니다 :)