뷔퐁의 바늘에서 맨해튼 프로젝트까지

뷔퐁의 바늘

몬테카를로 방법의 씨앗은 18세기로 거슬러 올라간다. 1777년, 프랑스의 수학자 조르주 루이 르클레르 드 뷔퐁은 평행선이 일정 간격으로 그어진 바닥에 바늘을 무작위로 던졌을 때 바늘이 선에 걸칠 확률을 계산하는 문제를 제안했다. 이른바 뷔퐁의 바늘Buffon's Needle 실험은 기하학적 확률을 이용해 원주율 $\pi$를 추정할 수 있음을 보인 최초의 사례로, 무작위 시행을 통해 결정론적 상수를 도출할 수 있다는 몬테카를로 방법의 원형이다.

바늘의 모양과 위치는 바늘의 길이 $l$, 그리고 바늘과 평행선이 이루는 각도 $\theta (0 \le \theta \le \pi)$에 의해 결정되는데, 이처럼 사건이 연속적으로 하면 그 개수를 셀 수 없기 때문에 길이나 영역의 넓이를 확률의 계산에 활용한다. 즉, 어떤 특정한 사건이 일어나는 영역의 크기를 사건이 발생할 수 있는 모든 영역의 크기로 나누어 기하학적 확률을 구한다.

맨해튼 프로젝트

체계적인 몬테카를로 방법이 탄생한 곳은 제2차 세계대전 중 미국 로스앨러모스 국립연구소에서 진행된 맨해튼 프로젝트Manhattan Project였다. 과학자들은 핵분열 물질 내에서 중성자가 어떻게 이동하고, 충돌하며, 증식하는지를 예측해야 했는데, 이는 다차원 적분 방정식으로 표현되는 문제로서 기존의 해석적 방법으로는 풀이가 불가능했다.

이때 폴란드 출신의 수학자 스타니스와프 울람Stanislaw Ulam은 질병 회복기 동안 카드 놀이를 하다가 핵심적인 직관을 얻었다. 카드를 성공적으로 배열할 확률을 조합론적으로 계산하려다 실패하자, 카드를 수천 번 무작위로 섞어 게임을 해본 뒤 성공 횟수를 세는 것이 훨씬 효율적이라는 사실을 깨달은 것이다. 울람은 이 아이디어를 존 폰 노이만John von Neumann에게 전했고, 폰 노이만은 이를 당시 개발 중이던 세계 최초의 범용 전자식 컴퓨터 ENIAC에 구현하여 중성자 수송 문제를 시뮬레이션했다.

이 방법론은 극비 군사 프로젝트의 일부였기에 코드명이 필요했다. 울람의 동료였던 니콜라스 메트로폴리스Nicholas Metropolis가 모나코의 유명한 카지노 도시 '몬테카를로'의 이름을 따서 명명했다. 룰렛 휠의 회전처럼 예측 불가능한 무작위성이 알고리즘의 핵심 동력이라는 점에서 적절한 은유였다. 1949년, 울람과 메트로폴리스가 "The Monte Carlo Method"라는 논문을 발표하면서 이 이름은 학계에 공식적으로 정착되었다.

수학적 토대: 무작위성 속의 질서

몬테카를로 방법이 단순한 무작위 추측이 아닌 과학적 정밀성을 갖는 이유는 확률론의 두 가지 핵심 정리, 대수의 법칙과 중심극한정리에 의해 뒷받침되기 때문이다.

대수의 법칙

대수의 법칙(Law of Large Numbers, LLM)은 몬테카를로 시뮬레이션의 수렴성convergence을 보장하는 이론적 근거다. 독립적이고 동일한 확률분포(독립항등분포, i.i.d)를 따르는 표본의 수가 무한히 커질수록, 표본의 평균은 모집단의 실제 기댓값으로 수렴한다는 정리이다.

어떤 다차원 공간 $D$에서 정의된 함수 $f(x)$의 적분값 $I$를 구한다고 하자.

$I = \int_{D} f(x) \, dx$

영역 $D$의 부피를 $V$라 하고, $D$에서 균등 분포를 따르는 확률 변수 $X$가 있다고 할 때 적분값 $I$는 기댓값 $E[f(X)] \times V$로 표현된다. 대수의 법칙에 의해 $N$개의 무작위 표본 $X_1, X_2, \dots, X_N$을 추출하여 계산한 추정값 $I_N$은 $N \to \infty$일 때 $I$로 수렴한다.

$\lim_{N \to \infty} I_N = \lim_{N \to \infty} \frac{V}{N} \sum_{i=1}^{N} f(X_i) = I \quad (\text{거의 확실하게})$

이 정리는 충분히 많은 시뮬레이션을 수행한다면, 반드시 참값에 근접할 수 있다고 보장한다.

중심극한정리와 오차 추정

하지만 수렴한다는 사실만으로는 충분하지 않다. "얼마나 많은 샘플이 필요한가?"와 "현재의 추정값은 얼마나 정확한가?"에 대한 답이 필요하다. 중심극한정리(Central Limit Theorem, CLT)가 바로 이 답을 제공한다.

CLT에 따르면, 표본의 크기 $N$이 충분히 클 때, 추정 오차 $(I_N - I)$의 분포는 정규분포 $\mathcal{N}(0, \sigma^2 / N)$에 근사한다. 여기서 표준 오차Standard Error는 $\sigma / \sqrt{N}$에 비례하므로, 몬테카를로 방법의 수렴 속도는 $O(N^{-1/2})$이다. 이는 표본 수를 100배 늘려야 오차가 10분의 1로 줄어든다는 뜻이므로 언뜻 보기에는 느려 보인다. 그러나 이 수렴 속도가 강력한 이유는 문제의 차원 $d$와 무관하다는 데 있다. 사다리꼴 공식¹이나 심슨 공식² 같은 격자 기반 수치 적분은 차원이 $d$일 때 수렴 속도가 $O(N^{-1/d})$로 떨어진다. 차원이 높아질수록 필요한 계산량이 지수적으로 폭증하는 '차원의 저주'에 걸리는 것이다. 반면 몬테카를로 방법은 100차원이든 1,000차원이든 수렴 속도가 여전히 $O(N^{-1/2})$를 유지한다. 이것이 몬테카를로 방법이 수십 개의 기초자산을 다루는 금융 파생상품이나 수많은 입자의 자유도를 다루는 물리학 문제에서 사실상 유일한 해법으로 통하는 이유다.

핵심 알고리즘

몬테카를로 시뮬레이션은 문제의 성격에 따라 다양한 형태로 구현되지만, 일반적으로 모델 정의 → 반복적인 샘플링 → 결과 집계의 3단계 프로세스를 따른다.

1. 모델 정의: 시뮬레이션을 위해서는 실제 세계를 본뜬 모델이 필요하다. 어떤 시스템을 본뜨기 위해서는 확률 밀도 함수가 필요하다. 어떤 결과물을 원하는지, 그 결과물이 어떻게 쓰일지, 결과물이 얼마나 정확해야 하는지, 입력은 어떻게 정의ㅚ는지 등 확륢 밀도 함수를 결정하기 위한 여러가지 질문이 선행되어야 한다.
2. 반복적인 샘플링: 확률 밀도 함수는 시뮬레이션의 특정 단계에서 가능한 사건의 범위와 확률을 정의하는 함수다. 음수가 아닌 실수 값만을 가지며 가능한 사건의 범위를 넘어서는 구간에 대한 적분은 무조건 1이어야 한다. 샘플링을 위해서는 누적 분포 함수, 즉 확률 분포 함수의 적분을 사용한다.
3. 결과 집계: 몬테카를로 시뮬레이션은 일반적으로 똑같은 과정을 반복하며 일련의 사건들을 만들어낸다. 각 사건은 그들의 속성으로 기록된다.

기본 몬테카를로 적분

가장 기초적인 형태는 난수를 이용해 함수의 면적이나 상수를 추정하는 것이다. 원주율 $\pi$를 추정해보자. 한 변의 길이가 2인 정사각형($[-1, 1] \times [-1, 1]$) 안에 반지름 1인 원을 내접시키고, 정사각형 내부에 점 $(x, y)$를 균등 분포로 무작위 생성한 뒤 $x^2 + y^2 \le 1$을 만족하는 점의 수를 센다. 그러면 전체 점의 수 $N$과 원 내부 점의 수 $N_{\text{in}}$의 비율로 $\pi$를 추정할 수 있다.

$\pi \approx 4 \times \frac{N_{\text{in}}}{N}$

정사각형의 넓이가 4이고 원의 넓이가 $\pi$이므로, 점이 원 안에 떨어질 확률은 $\pi / 4$가 된다. 양변에 4를 곱하면 위와 같은 식이 나온다. 개념적으로 단순하지만, 실제로는 유효숫자 두세 자리를 얻는 데도 수백만 번의 시행이 필요하다. 이 느린 수렴을 개선하기 위한 분산 감소 기법은 아래에서 다룬다.

마르코프 체인 몬테카를로(Markov Chain Monte Carlo, MCMC)

고차원의 복잡한 확률 분포에서 표본을 추출해야 할 때, 기각 샘플링Rejection Sampling 같은 단순한 몬테카를로 방법은 효율이 급격히 떨어진다. 마르코프 체인 몬테카를로(MCMC)는 이 문제를 해결하기 위해 등장한 기법이다. MCMC의 핵심 아이디어는 표본들이 서로 독립적이지 않고, 이전 상태에 의존하는 마르코프 체인Markov Chain을 형성하도록 하는 것이다. 시간이 지남에 따라 이 체인의 분포는 목표로 하는 확률 분포Target Distribution에 수렴하게 된다.

메트로폴리스-헤이스팅스 알고리즘

1953년 메트로폴리스가 제안하고 1970년 헤이스팅스가 일반화한 이 알고리즘은 MCMC의 가장 범용적인 형태다. 절차는 다음과 같다.

1. 현재 상태 $x_t$에서 제안 분포Proposal Distribution $q(x' \mid x_t)$를 통해 후보 상태 $x'$를 생성한다.

2. 수락 확률Acceptance Probability을 계산한다. $\alpha = \min\!\left(1, \; \frac{P(x') \, q(x_t \mid x')}{P(x_t) \, q(x' \mid x_t)}\right)$

3. 균등 분포에서 난수 $u \in [0, 1]$을 뽑아 $u \le \alpha$이면 이동을 수락($x_{t+1} = x'$)하고, 그렇지 않으면 현재 상태에 머문다($x_{t+1} = x_t$).

이 알고리즘의 강력한 점은 정규화 상수를 몰라도 분포의 비정규화된 형태만 알면 샘플링이 가능하다는 것이다. 수락 확률을 계산할 때 분자와 분모에 동일한 정규화 상수가 곱해지므로 자연스럽게 약분되기 때문이다. 이 성질 덕분에 메트로폴리스-헤이스팅스 알고리즘은 사후 분포의 정규화 상수를 구하기 어려운 베이지안 통계학의 부흥을 이끌었다.

깁스 샘플링

깁스 샘플링Gibbs Sampling은 메트로폴리스-헤이스팅스 알고리즘의 특수한 형태로, 다변수 분포에서 변수들을 한 번에 하나씩 조건부 확률 분포Conditional Distribution에 따라 업데이트한다. 제안 분포를 조건부 분포 자체로 사용하므로 수락 확률이 항상 1이 되어 효율적이다. 이미지 처리나 은닉 마르코프 모델(HMM) 학습 등에 널리 쓰인다.

분산 감소 기법

기본 몬테카를로 방법의 느린 수렴 속도 $O(N^{-1/2})$를 개선하기 위해, 추정량의 편향은 건드리지 않으면서 분산 $\sigma^2$만을 줄이는 기법들이 개발되었다. 이를 통해 적은 샘플 수로도 높은 정밀도를 얻을 수 있다.

대조 변수법

대조 변수법Antithetic Variates은 입력 변수 간의 음의 상관관계를 인위적으로 유도하여 오차를 상쇄시키는 기법이다. 난수 $U$를 생성하여 $f(U)$를 계산할 때, 그와 대칭되는 $1 - U$를 이용하여 $f(1 - U)$도 함께 계산한다.

$\hat{I} = \frac{f(U) + f(1 - U)}{2}$

이 추정량의 분산은 다음과 같다.

$\text{Var}(\hat{I}) = \frac{1}{4}\bigl(\text{Var}(f(U)) + \text{Var}(f(1 - U)) + 2\,\text{Cov}(f(U), f(1 - U))\bigr)$

$f$가 단조 함수라면 공분산 $\text{Cov}$가 음수가 되므로 전체 분산이 크게 감소한다.

제어 변수법

제어 변수법Control Variates은 적분하려는 복잡한 함수 $f(x)$와 형태가 유사하면서, 기댓값 $E[g(x)]$를 이미 알고 있는 보조 함수 $g(x)$를 활용하는 기법이다. 새로운 추정량을 다음과 같이 구성한다.

$f^{*}(x) = f(x) - c\bigl(g(x) - E[g(x)]\bigr)$

적절한 상수 $c$를 선택하면 $f^{*}$의 기댓값은 원래의 기댓값과 같으면서 분산은 훨씬 작아진다. $f$와 $g$의 상관계수가 높을수록 분산 감소 효과가 크다.

중요 샘플링

중요 샘플링Importance Sampling은 결과값에 큰 영향을 미치는 '중요한' 영역을 더 자주 샘플링하고, 덜 중요한 영역은 드물게 샘플링하여 효율을 높이는 방법이다. 이 기법은 특히 희귀 사건 시뮬레이션에서 필수적이다. 원자로 사고나 금융 시장 붕괴처럼 발생 확률이 극히 낮은($10^{-6}$ 이하) 사건을 시뮬레이션할 때, 일반적인 샘플링으로는 수억 번을 돌려도 해당 사건이 관측되지 않을 수 있기 때문이다.

원래의 확률 밀도 함수 $P(x)$ 대신, 희귀 사건이 잘 발생하는 새로운 분포 $Q(x)$에서 샘플을 추출한 뒤, 결과에 가중치 $P(x) / Q(x)$(우도비Likelihood Ratio)를 곱하여 편향을 보정한다. 이 기법은 금융 리스크 관리(VaR 계산)나 입자 물리학의 충돌 시뮬레이션에서 표준적으로 사용된다.

---

출처

• Metropolis, N. and Ulam, S. (1949) 'The Monte Carlo method', Journal of the American Statistical Association, 44(247), pp. 335–341.
• Metropolis, N. et al. (1953) 'Equation of state calculations by fast computing machines', The Journal of Chemical Physics, 21(6), pp. 1087–1092.
• Hastings, W.K. (1970) 'Monte Carlo sampling methods using Markov chains and their applications', Biometrika, 57(1), pp. 97–109.
• Robert, C.P. and Casella, G. (2004) Monte Carlo statistical methods. 2nd edn. New York: Springer.
• Kroese, D.P., Taimre, T. and Botev, Z.I. (2011) Handbook of Monte Carlo methods. Hoboken, NJ: Wiley.
• Glasserman, P. (2003) Monte Carlo methods in financial engineering. New York: Springer.
• 중심 극한 정리(CLT, Central Limit Theorem)
• Introduction To Monte Carlo Simulation

\ F: \text{파생상품의 가격}, S: \text{기초자산의 가격}, r: \text{무위험 이자율}, t: \text{시간}, σ: \text{변동성} \end{array}$$

Footnotes

1. 닫힌구간 $\displaystyle [t_{0},t_{N}]$위의 적분 가능 함수 $f\colon [t_{0},t_{N}]\to \mathbb {R}$ 및 수열 $t_{0}\leq t_{1}\leq \dotsc \leq t_{N-1}\leq t_{N}$이 주어졌을 때, 이에 대한 적분

   $$F=\int _{t_{0}}^{t_{N}}f(x)\,\mathrm {d} x$$

   의 사다리꼴 공식 근사는 다음과 같다.

   $$\tilde {F} = \sum _{i=0}^{N-1}{\frac {(t_{i+1}-t_{i})(f(t_{i+1})+f(t_{i}))}{2}}$$ ↩

2. 심슨 공식은 $P(x)$라는 이차방정식을 이용해 $f(x)$의 근사값을 구한다. 이때 $P(x)$는 $a$, $b$, 그리고 둘의 중간값 $m=\frac {a+b}{2}$에서 $f(x)$와 같은 값을 갖는 근사식이다. 라그랑주의 다항식 보강법을 사용해서 $P(x)$를 구하면 다음을 얻는다.

   $P(x)=f(a){\frac {(x-m)(x-b)}{(a-m)(a-b)}}+f(m){\frac {(x-a)(x-b)}{(m-a)(m-b)}}+f(b){\frac {(x-a)(x-m)}{(b-a)(b-m)}}$

   이 식을 전개하면 심프슨 공식으로 알려진 다음 공식을 구할 수 있다.

   $\int _{a}^{b}f(x)\,dx\approx \int _{a}^{b}P(x)\,dx={\frac {b-a}{6}}\left[f(a)+4f\left({\frac {a+b}{2}}\right)+f(b)\right]$ ↩