> [!abstract] Introduction > 해석적으로 풀 수 없는 문제 앞에서, 과학자들은 역설적이게도 '무작위성'에서 답을 찾았다. 몬테카를로 방법*Monte Carlo Method*은 난수를 이용한 반복 시뮬레이션으로 함수의 값을 근사하거나 복잡한 시스템의 거동을 모사하는 계산 기법을 통칭한다. 핵무기 개발의 부산물로 탄생한 이 방법론은 오늘날 금융 파생상품 가격 결정, 방사선 치료 계획, 신약 개발, 자율주행, 그리고 [[Monte Carlo Tree Search(MCTS)|몬테카를로 트리 탐색(MCTS)]]을 통한 AI 추론에 이르기까지, 불확실성을 다루는 모든 분야의 표준 도구로 자리 잡았다. 이 글에서는 몬테카를로 방법의 역사적 기원, 수학적 토대, 핵심 알고리즘, 그리고 다양한 응용 사례를 살펴본다. # 뷔퐁의 바늘에서 맨해튼 프로젝트까지 ## 뷔퐁의 바늘 몬테카를로 방법의 씨앗은 18세기로 거슬러 올라간다. 1777년, 프랑스의 수학자 조르주 루이 르클레르 드 뷔퐁은 평행선이 일정 간격으로 그어진 바닥에 바늘을 무작위로 던졌을 때 바늘이 선에 걸칠 확률을 계산하는 문제를 제안했다. 이른바 뷔퐁의 바늘*Buffon's Needle* 실험은 기하학적 확률을 이용해 원주율 $\pi$를 추정할 수 있음을 보인 최초의 사례로, 무작위 시행을 통해 결정론적 상수를 도출할 수 있다는 몬테카를로 방법의 원형이다. ![[BuffonsNeedle.webp]] 바늘의 모양과 위치는 바늘의 길이 $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)](https://gguguk.github.io/posts/CLT/) - [Introduction To Monte Carlo Simulation](https://pmc.ncbi.nlm.nih.gov/articles/PMC2924739/) [^사다리꼴공식]: 닫힌구간 $\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}} $ [^심슨공식]: 심슨 공식은 $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]$ [^블랙-슐즈공식]: $\begin{array}{l} {\partial F \over \partial t}+{1 \over 2}\sigma ^{2}S^{2}{\partial ^{2}F \over \partial S^{2}}+rS{\partial F \over \partial S}=rF \\ F: \text{파생상품의 가격}, S: \text{기초자산의 가격}, r: \text{무위험 이자율}, t: \text{시간}, σ: \text{변동성} \end{array}$