몬테카를로 시뮬레이션은 무작위의 난수를 추출해서 함수의 값을 구하는 방법입니다. 대표적으로, 특정 분포의 근사 혹은 수치적분을 할 수 있는데요. 지난 포스팅에서는 이산확률분포를 근사하는 예제를 보여드렸고, 이번 시간에는 \(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\)의 값을 갖고, 이외의 구간에서는 \(0\) 의 값을 가집니다. 그러므로 위의 식을 아래처럼 쓸 수 있습니다.

$$ \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} {\sum_{i=1}^{n} {g(u_{i})}\over {n}} = E(g(u)) = \theta $$

핵심

대수의 법칙에 의해서 표본의 크기가 커질수록 표본평균은 모평균에 근접하기 때문에, 표본을 많이 생성하면 수치적으로 구한 표본평균은 우리가 원하는 \(\theta = E(g(u))\)에 근접하게 됩니다. 즉, 확률변수 \({U}\)를 \({n}\)개 생성한 다음에 이들을 함수 \(g\)에 입력하고 그 결과를 평균하면 됩니다. 이 때, 표준균일분포를 따르는 난수 \(U\)들의 배열을 생성하면 됩니다.

예제

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

# 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]\)로 변환하면서 그에 맞춰 함수를 적절히 변환하고, 위의 테크닉을 사용하면 됩니다. 그 방법은 다음 포스팅에서 작성하도록 하겠습니다. 감사합니다 :)