극한 정리와 시뮬레이션
웨이퍼 한 장의 문턱 전압 산포를 히스토그램으로 그리면 거의 언제나 종 모양이 나온다. 공정 엔지니어 누구도 "정규 분포가 되게 하라"고 지시한 적이 없는데도 그렇다. 이유는 하나의 정리에 있다. 서로 독립인 작은 요인이 많이 더해지면, 요인 하나하나가 어떤 모양이든 합은 정규 분포에 가까워진다. 이 장은 그 중심 극한 정리와 짝인 큰 수의 법칙, 분포를 몰라도 쓸 수 있는 부등식, 그리고 이 정리들을 계산 도구로 바꾼 몬테카를로 방법을 다룬다. 아래에서 모집단 막대를 끌어 아무리 이상한 모양을 만들어도, 표본 평균은 결국 종 모양으로 쌓인다.
- 마르코프 부등식 \(P(X\ge a)\le E[X]/a\)와 체비쇼프 부등식 \(P(|X-\mu|\ge k\sigma)\le 1/k^2\)을 유도하고, 분포를 모를 때의 꼬리 확률 상한으로 쓰며, 그 상한이 얼마나 느슨한지 판단한다.
- 표본 평균의 분산 \(\sigma^2/n\)에서 약한 큰 수의 법칙을 증명하고, 확률 수렴과 거의 확실한 수렴의 뜻, 평균이 없는 코시 분포에서 법칙이 깨지는 이유를 설명한다.
- 중심 극한 정리 \((\bar X-\mu)/(\sigma/\sqrt n)\to N(0,1)\)을 진술하고, 표본 평균·합의 확률을 근사 계산하며, 모집단의 왜도가 수렴 속도를 정한다는 것(\(\gamma/\sqrt n\))을 안다.
- 갈톤 보드와 공정 산포의 예로 "많은 작은 독립 요인의 합 → 정규"를 이해하고, 한 요인이 지배하면 정규성이 깨지는 경우를 알아본다.
- 이항·포아송 분포를 정규 분포로 근사하고, 연속성 수정 \(\pm0.5\)를 써서 \(P(X\le k)\), \(P(X=k)\) 등을 계산한다.
- 몬테카를로 적분의 추정량과 표준 오차 \(\sigma/\sqrt N\)을 구하고, 원하는 정밀도에 필요한 표본 수를 계산하며, 오차가 차원과 무관하다는 장점을 안다.
- 선형 합동 생성기와 헐–도벨 정리로 의사 난수의 주기와 격자 구조를 이해하고, 시드를 고정해 시뮬레이션을 재현한다.
마르코프·체비쇼프 부등식
공정 데이터에서 평균과 표준 편차는 쉽게 얻지만, 분포 모양까지 확신하기는 어렵다. 그럴 때도 "규격 밖으로 나갈 확률이 많아야 얼마"라는 말은 할 수 있다. 평균과 분산 두 숫자만으로 꼬리 확률의 상한을 주는 두 부등식이 이 장의 모든 극한 정리의 출발점이다.
\(X\ge 0\)인 확률 변수와 \(a>0\)에 대해 \(\;P(X\ge a)\le \dfrac{E[X]}{a}.\)
마르코프 부등식(Markov's inequality)의 증명은 한 줄이다. 사건 \(\{X\ge a\}\)의 지시 함수 \(\mathbf 1\{X\ge a\}\)는 \(X\ge a\)이면 1, 아니면 0이다. \(X\ge0\)이므로 어디서나 \(a\,\mathbf 1\{X\ge a\}\le X\)이고(아래 그림), 양변의 기댓값을 취하면 \(a\,P(X\ge a)\le E[X]\)이다. 평균이 작은 음이 아닌 양은 큰 값을 자주 낼 수 없다는 뜻이다. 평균이 1인 양이 10 이상일 확률은 10%를 넘지 못한다 — 넘으면 그것만으로 평균이 1을 넘기 때문이다.
이제 \(X\) 대신 음이 아닌 양 \((X-\mu)^2\)에 마르코프 부등식을 쓰고 \(a = k^2\sigma^2\)으로 두면, \(P\big((X-\mu)^2\ge k^2\sigma^2\big)\le E[(X-\mu)^2]/k^2\sigma^2 = 1/k^2\)이다. 이것이 체비쇼프 부등식(Chebyshev's inequality)이다.
| k | 체비쇼프 상한 \(1/k^2\) | 정규 \(2(1-\Phi(k))\) | 균일 분포 | 지수 분포 \(e^{-(1+k)}\) |
|---|---|---|---|---|
| 1.5 | 0.444 | 0.134 | 0.134 | 0.0821 |
| 2 | 0.250 | 0.0455 | 0 | 0.0498 |
| 3 | 0.111 | 0.00270 | 0 | 0.0183 |
균일 분포는 \(\mu\pm\sqrt3\,\sigma\) 밖에 값이 없으므로 \(k\ge\sqrt3\)이면 꼬리 확률이 0이다. 지수 분포(\(\mu=\sigma=1/\lambda\))는 왼쪽 꼬리가 없고 오른쪽 꼬리 \(P(X\ge\mu+k\sigma)=e^{-(1+k)}\)만 남는다. 어느 쪽이든 체비쇼프 상한은 실제보다 몇 배에서 수십 배 크다. 아래 시뮬레이터로 분포를 바꿔 가며 확인해 보자.
정규 분포에서 \(\mu\pm2\sigma\) 밖으로 나갈 실제 확률은 체비쇼프 상한 0.25와 비교해 어느 정도일까?
꼬리 확률과 부등식 상한
그렇다면 체비쇼프 상한은 쓸모없이 느슨한가? 그렇지 않다. 평균과 분산만 알 때는 이보다 좋은 상한이 있을 수 없다. 세 점 \(\mu-k\sigma,\ \mu,\ \mu+k\sigma\)에 확률 \(\tfrac{1}{2k^2},\ 1-\tfrac1{k^2},\ \tfrac1{2k^2}\)을 두면 평균이 \(\mu\), 분산이 \(\sigma^2\)이고 꼬리 확률이 정확히 \(1/k^2\)이기 때문이다.
막두께가 평균 100, 표준 편차 2.0인 공정에서 규격을 평균 ± 5.0로 잡으면 \(k = c/\sigma\) = 이다. 분포 모양을 전혀 모르더라도 규격 밖 비율은 많아야 이다. 분포가 정규라면 실제로는 로 .
규격이 한쪽뿐(누설 전류 상한 등)이면 체비쇼프를 그대로 쓰는 것보다 칸텔리 부등식 \(P(X-\mu\ge k\sigma)\le \dfrac{1}{1+k^2}\)이 낫다. \(k=2\)이면 0.2(체비쇼프 양쪽 0.25). 증명은 \(P(X-\mu\ge k\sigma)\le P\big((X-\mu+t)^2\ge(k\sigma+t)^2\big)\)에 마르코프를 쓰고 \(t=\sigma/k\)로 고르면 된다.
(a) 다이 한 개의 결함 수 \(X\)는 평균 2개다. 분포를 모를 때 \(P(X\ge10)\)의 상한은? 실제로 \(X\sim\mathrm{Poisson}(2)\)라면? (b) 막두께가 \(\mu = 100\) nm, \(\sigma = 2\) nm다. 분포를 모를 때 \(|X-100|\ge5\) nm일 확률의 상한은? 정규라면?
풀이 보기
큰 수의 법칙
같은 측정을 여러 번 반복해 평균을 내면 값이 안정된다는 것은 누구나 경험으로 안다. 2장에서는 동전을 많이 던지면 앞면의 상대 도수가 1/2에 다가간다는 직관으로 확률을 소개했다. 이제 그 직관을 정리로 바꾼다. \(X_1, X_2, \dots\)가 평균 \(\mu\), 분산 \(\sigma^2\)인 독립 동일 분포(iid)라 하고 표본 평균을 \(\bar X_n = (X_1+\cdots+X_n)/n\)이라 하자. 6장의 선형 결합 공식에서
여기에 체비쇼프 부등식을 쓰면 곧바로 다음 정리가 나온다.
임의의 \(\varepsilon>0\)에 대해 \(\;P\big(|\bar X_n-\mu|\ge\varepsilon\big)\le\dfrac{\sigma^2}{n\varepsilon^2}\;\longrightarrow\;0\quad(n\to\infty).\)
이처럼 "벗어날 확률이 0으로 가는" 수렴을 확률 수렴(convergence in probability)이라 하고 \(\bar X_n\xrightarrow{P}\mu\)로 쓴다. 이 정리를 큰 수의 법칙(law of large numbers)이라 한다. 아래 그림은 증명의 흐름을 보여 준다. \(n\)이 커질수록 \(\bar X_n\)의 분포가 좁아져 \(\pm\varepsilon\) 밖의 면적이 줄어든다.
- \(n=1\): 측정 하나는 \(\pm\varepsilon\) 밖으로 자주 나간다(빨강 면적 0.62). 체비쇼프 상한 \(\sigma^2/\varepsilon^2 = 4\)는 1보다 커서 아무것도 말하지 않는다.
- \(n=4\): 분산이 \(\sigma^2/4\)로 줄어 분포가 반으로 좁아졌다. 상한 \(\sigma^2/(n\varepsilon^2) = 1\).
- \(n=16\): 상한 0.25. 실제 값(0.046)은 훨씬 작지만, 상한만으로도 0으로 가는 것이 보이기 시작한다.
- \(n=64\): 상한 \(4/n\)은 \(n\)이 커지면 반드시 0으로 간다. 분포의 모양을 몰라도 분산만 유한하면 된다 — 약한 큰 수의 법칙의 증명 끝.
더 강한 진술도 있다. 강한 큰 수의 법칙은 \(P\big(\lim_{n\to\infty}\bar X_n = \mu\big) = 1\), 곧 표본 평균의 경로 하나하나가 (확률 1로) 결국 \(\mu\)에 붙어 버린다는 것이다. 이런 수렴을 거의 확실한 수렴(almost sure convergence)이라 한다. 약한 법칙은 "큰 \(n\) 하나를 고정하면 벗어난 경로가 드물다"는 말이고, 강한 법칙은 "각 경로가 어느 순간부터 영원히 \(\pm\varepsilon\) 안에 머문다"는 말이다. 강한 법칙은 평균 \(E|X|\)만 유한하면 성립한다(분산이 필요 없다). 아래 시뮬레이터의 스파게티 경로가 두 법칙을 함께 보여 준다.
표본 평균의 경로와 깔때기
큰 수의 법칙을 "동전 앞면이 많이 나왔으니 이제 뒷면이 나올 차례"라는 뜻으로 오해하면 안 된다. 동전에는 기억이 없다. 앞면 비율이 1/2로 가는 것은 앞선 치우침이 보상되어서가 아니라, 이후의 수많은 공정한 던지기에 묽어지기 때문이다.
동전을 10번 던졌을 때보다 10 000번 던졌을 때, "앞면 개수 − 던진 횟수의 절반"의 크기는 보통 어떻게 될까?
큰 수의 법칙에는 전제가 있다. 평균이 존재해야 한다. 2장에서 만난 코시 분포 \(f(x) = \frac{1}{\pi(1+x^2)}\)은 꼬리가 너무 두꺼워 \(E|X| = \infty\)이고, 놀랍게도 코시 표본 \(n\)개의 평균은 다시 같은 코시 분포를 따른다. 표본을 아무리 많이 모아도 평균이 한 개의 측정보다 나아지지 않는다.
불량률 \(p\)를 표본 비율 \(\hat p\)로 추정해 오차 ± 0.020 안에 들 확률을 95.0 이상으로 만들고 싶다. \(p(1-p)\le 1/4\)이므로 체비쇼프 부등식으로는 \(n\ge \dfrac{1/4}{(1-\text{신뢰})\,\varepsilon^2}\) = 개면 충분하다. 다음 절의 중심 극한 정리를 쓰면 \(n\approx (z/2\varepsilon)^2\) = 개로 .
새 검사 장비의 과검출률 \(p\)를 추정하려 한다. 추정 오차가 0.02 이내일 확률이 적어도 0.95가 되려면 체비쇼프 부등식으로 몇 개를 검사해야 하는가? 중심 극한 정리(다음 절)를 쓰면?
풀이 보기
중심 극한 정리
큰 수의 법칙은 \(\bar X_n\)이 \(\mu\)로 간다고 말할 뿐, 가는 도중에 \(\mu\) 주위에서 어떤 모양으로 흔들리는지는 말하지 않는다. \(\bar X_n-\mu\)는 \(1/\sqrt n\) 크기로 줄어드니, 그만큼 확대해 \(\sqrt n\) 배를 해 보자. 놀랍게도 확대한 모양은 모집단이 무엇이든 언제나 같은 곡선 — 정규 분포 — 로 간다.
\(X_1,\dots,X_n\)이 평균 \(\mu\), 분산 \(0<\sigma^2<\infty\)인 iid이면, 모든 \(z\)에서 $$\lim_{n\to\infty}P\!\left(\frac{\bar X_n-\mu}{\sigma/\sqrt n}\le z\right) = \lim_{n\to\infty}P\!\left(\frac{S_n-n\mu}{\sigma\sqrt n}\le z\right) = \Phi(z).$$
이것이 중심 극한 정리(central limit theorem, CLT)다. 누적 분포 함수가 수렴하는 이런 수렴을 분포 수렴(convergence in distribution)이라 하고 \(Z_n\xrightarrow{d}N(0,1)\)로 쓴다. 실무에서는 큰 \(n\)에 대해 다음 근사로 쓴다.
오른쪽으로 크게 치우친 지수형 모집단에서 크기 \(n=5\)인 표본의 평균을 수천 번 뽑아 히스토그램을 그리면?
어떤 모집단이든 표본 평균은 정규로
얼마나 큰 \(n\)이면 "충분히 큰가"는 모집단에 달렸다. 모집단이 대칭이면 \(n = 5\)–10만으로도 정규 근사가 훌륭하고, 많이 치우치면 수십, 수백 개가 필요하다. 수렴 속도를 정하는 것은 왜도 \(\gamma = E[(X-\mu)^3]/\sigma^3\)와 첨도다.
표준화한 \(Y_i = (X_i-\mu)/\sigma\)의 적률 생성 함수(4장)를 \(M(t) = 1 + \tfrac{t^2}{2} + \tfrac{\gamma t^3}{6} + \cdots\)라 하자(평균 0, 분산 1). \(Z_n = \sum Y_i/\sqrt n\)이므로 \(M_{Z_n}(t) = M(t/\sqrt n)^n = \big(1 + \tfrac{t^2}{2n} + O(n^{-3/2})\big)^n \to e^{t^2/2}\)이고, 이는 표준 정규의 적률 생성 함수다. 모집단 정보는 \(t^3\) 이상의 항에만 들어 있는데 그 항들이 \(n\)이 커지며 사라진다. 그래서 극한은 모집단과 무관하다. 오차의 크기는 베리–에센 정리가 준다: \(\sup_z|P(Z_n\le z)-\Phi(z)|\le \dfrac{0.4748\,E|X-\mu|^3}{\sigma^3\sqrt n}\).
CLT가 왜 성립하는지는 6장의 합성곱으로도 볼 수 있다. 독립인 두 확률 변수의 합의 밀도는 두 밀도의 합성곱이고, 합성곱은 모양을 "뭉개는" 연산이다. 균일 분포를 거듭 합성곱하면 사각형이 삼각형으로, 삼각형이 둥근 언덕으로 바뀐다.
- \(U(0,1)\) 하나를 표준화한 밀도는 높이 0.289인 사각형이다. 정규 곡선과 전혀 다르다.
- 두 개의 합은 사각형과 사각형의 합성곱인 삼각형이다. 벌써 가운데가 솟았다.
- 네 개의 합은 매끈한 언덕이다. 꼬리 끝(±3.46 밖)만 잘려 있다.
- 열두 개의 합은 정규 곡선과 거의 겹친다. 예전 계산기는 \(U_1+\cdots+U_{12}-6\)을 표준 정규 난수로 썼다.
위 시뮬레이터에서 그린 모집단으로 같은 일을 정확히 계산해 볼 수 있다. 막대 16개를 값 0, 1, …, 15를 갖는 이산 분포로 보면, 합의 확률 질량 함수는 합성곱으로 정확히 구해진다(표본 추출 없이).
어떤 소자의 수명은 평균 1000 h인 지수 분포를 따른다. 무작위로 고른 36개의 평균 수명이 1200 h를 넘을 확률을 CLT로 구하라. 정확한 값(평균이 감마 분포를 따름)과 비교하라.
풀이 보기
갈톤 보드와 수많은 작은 요인
1870년대 프랜시스 갈톤은 못을 삼각형으로 촘촘히 박은 판 위에서 구슬을 떨어뜨려 정규 분포를 "기계로" 만들었다. 이것이 갈톤 보드(Galton board)다. 구슬은 못에 부딪힐 때마다 오른쪽(확률 \(p\)) 또는 왼쪽으로 반 칸 비켜난다. \(R\)단을 지나면 오른쪽으로 간 횟수 \(K\)가 칸 번호가 되므로 \(K\sim\mathrm{Bin}(R,p)\)다. 그리고 \(K\)는 베르누이 \(R\)개의 합이니 CLT에 따라 \(N(Rp,\ Rp(1-p))\)에 가까워진다.
12단 보드(\(p = 1/2\))에 구슬 1000개를 떨어뜨리면 가운데 칸(\(k = 6\))에는 대략 몇 개가 들어갈까?
갈톤 보드
갈톤 보드의 교훈은 반도체 공정으로 그대로 옮겨진다. 트랜지스터 하나의 문턱 전압 \(V_t\)는 수많은 작은 요인의 합으로 흔들린다. 어느 장비·챔버에서 처리됐는지, 웨이퍼의 중심인지 가장자리인지, 채널에 도펀트 원자가 몇 개 들어갔는지, 게이트 가장자리가 얼마나 거친지, 산화막이 몇 Å 두꺼운지, 금속 게이트의 결정립 방향이 어떤지. 각 요인의 분포는 정규가 아니어도(두 챔버면 두 값, 위치면 균일에 가깝다) 독립이고 고만고만한 크기라면 합은 정규 분포에 가깝다.
작은 요인 여덟 개가 만드는 Vt 산포
공정 산포를 예측할 때 자주 쓰는 계산이 공차 누적이다. 층 \(n\)개가 각각 \(\pm t\) 안에서 균일하게 흔들리면, 모든 층이 한쪽 끝에 몰린 최악의 경우 합은 \(\pm nt\)까지 간다. 그러나 독립인 층들이 동시에 같은 쪽 끝에 몰릴 확률은 거의 없다. CLT가 주는 "통계적 공차"는 훨씬 좁다.
두께가 각각 ± 1.0 안에서 균일하게 흔들리는 막 10층을 쌓으면, 최악의 경우 합은 ± 까지 가지만, 합의 표준 편차는 \(t\sqrt{n/3}\) = 이고 CLT로 본 99.73% 구간(±3σ)은 ± 로 최악값의 이다. 합이 ± 3.0 안에 들 확률은 약 이다.
각 층의 두께 오차가 \(-1\) nm에서 \(+1\) nm 사이에 균일하게 분포하는 독립인 막 10층을 쌓는다. 전체 두께 오차가 ±3 nm 안에 들 확률을 CLT로 구하라.
풀이 보기
요인이 더해지지 않고 곱해지면(예: 여러 단계의 수율, 누설 전류의 지수 의존성 \(I\propto e^{-V_t/nV_T}\)) 로그를 취한 값이 합이 되므로 CLT는 \(\ln X\)에 적용된다. 그 결과 \(X\)는 정규가 아니라 로그 정규 분포(5장)를 따른다. 누설 전류·결함 크기·수명 분포가 오른쪽으로 긴 꼬리를 갖는 이유다.
이항·포아송의 정규 근사와 연속성 수정
갈톤 보드가 보여 주었듯 이항 확률 변수 \(X\sim\mathrm{Bin}(n,p)\)는 베르누이 \(n\)개의 합이므로 CLT가 그대로 적용된다. 컴퓨터가 없던 시절에는 \(\binom{400}{25}\) 같은 계산을 피하는 유일한 방법이었고, 지금도 손 계산과 표본 크기 설계, 신뢰 구간(10장)의 바탕이다.
그런데 이산 분포를 연속 분포로 근사하면 한 가지 문제가 생긴다. 이항 분포의 확률은 정수 \(k\)에만 있고, 히스토그램에서 \(P(X = k)\)는 \(k-0.5\)부터 \(k+0.5\)까지의 막대 넓이다. 정규 곡선 아래의 넓이로 그 막대를 덮으려면 경계를 반 칸씩 넓혀야 한다. 이것이 연속성 수정(continuity correction)이다.
- \(\mathrm{Bin}(10, 0.5)\)의 확률을 폭 1인 막대로 그리면 막대 \(k\)는 \(k-0.5\)에서 \(k+0.5\)까지를 차지한다.
- 평균 5, 분산 2.5인 정규 곡선이 막대 꼭대기를 지난다.
- \(P(X\le6)\)을 정규 곡선 아래 \(x\le6\)의 넓이로 계산하면 막대 6의 오른쪽 절반이 빠져 0.736이 된다. 정확한 값은 0.828이다.
- 경계를 6.5로 옮기면 막대 0–6을 모두 덮어 0.829가 된다. 오차가 0.09에서 0.0005로 줄었다.
| 구하는 확률 | 연속성 수정을 한 정규 근사 |
|---|---|
| \(P(X\le k)\) | \(\Phi\!\left(\dfrac{k+0.5-np}{\sqrt{np(1-p)}}\right)\) |
| \(P(X< k) = P(X\le k-1)\) | \(\Phi\!\left(\dfrac{k-0.5-np}{\sqrt{np(1-p)}}\right)\) |
| \(P(X\ge k)\) | \(1-\Phi\!\left(\dfrac{k-0.5-np}{\sqrt{np(1-p)}}\right)\) |
| \(P(a\le X\le b)\) | \(\Phi\!\left(\dfrac{b+0.5-np}{\sigma}\right)-\Phi\!\left(\dfrac{a-0.5-np}{\sigma}\right)\) |
| \(P(X = k)\) | \(\Phi\!\left(\dfrac{k+0.5-np}{\sigma}\right)-\Phi\!\left(\dfrac{k-0.5-np}{\sigma}\right)\) |
요령은 "구하는 정수 막대를 모두 덮도록 반 칸 넓힌다"이다. \(<\)와 \(\le\)를 먼저 정수 경계로 바꾼 뒤 수정하면 헷갈리지 않는다.
\(X\sim\mathrm{Bin}(100, 0.5)\)에서 \(P(X = 50)\)을 연속성 수정 없이 정규 근사로 계산하면?
이항 분포의 정규 근사 · 구간 과녁
\(X\sim\mathrm{Bin}(\)400, 0.05\()\)에서 \(P(X\le\)25\()\)는 정확히 이다. 정규 근사는 연속성 수정을 하면 (z = ), 하지 않으면 (z = )이다.
불량률 5%인 공정에서 칩 400개를 검사한다. 불량이 25개 이하일 확률을 정규 근사(연속성 수정 포함·제외)로 구하고 정확한 이항 확률과 비교하라.
풀이 보기
포아송 분포의 정규 근사
\(\mathrm{Poisson}(\lambda)\)는 독립인 \(\mathrm{Poisson}(1)\) \(\lambda\)개의 합으로 볼 수 있으므로(정수 \(\lambda\)일 때; 일반 \(\lambda\)도 같은 결론) \(\lambda\)가 크면 \(N(\lambda,\lambda)\)에 가까워진다. 왜도는 \(1/\sqrt\lambda\)로 줄어든다. 보통 \(\lambda\ge10\)–20이면 연속성 수정과 함께 쓴다. 상대 산포 \(\sigma/\mu = 1/\sqrt\lambda\)는 광자·전자·입자를 세는 모든 측정의 근본 잡음, 곧 샷 잡음이다.
EUV 노광에서 어떤 컨택 홀 영역에 흡수되는 광자 수 \(N\)이 평균 100인 포아송 분포를 따른다. (a) 광자가 90개 이하(10% 이상 부족)일 확률, (b) 80개 이하일 확률을 연속성 수정한 정규 근사로 구하고 정확한 값과 비교하라.
풀이 보기
몬테카를로 적분과 오차
큰 수의 법칙과 CLT를 거꾸로 쓰면 강력한 계산 도구가 된다. 구하고 싶은 양을 어떤 확률 변수의 기댓값 \(\theta = E[g(X)]\)로 쓸 수 있으면, \(X\)를 컴퓨터로 \(N\)번 뽑아 평균을 내면 된다. 큰 수의 법칙이 수렴을, CLT가 오차의 크기를 보장한다. 이것이 몬테카를로 방법(Monte Carlo method)이다. 1946–47년 로스앨러모스 연구소의 울람과 폰 노이만이 중성자 확산 계산에 쓰기 시작했고, 동료 메트로폴리스가 모나코의 카지노 이름을 붙였다.
가장 유명한 예는 다트로 \(\pi\)를 구하는 것이다. 단위 정사각형에 다트를 균일하게 던지면 사분원 안에 들어갈 확률은 넓이 비 \(\pi/4\)이다. 지시 변수 \(g = \mathbf 1\{U_1^2+U_2^2\le1\}\)의 평균이 \(\pi/4\)이므로 \(\hat\pi = 4\times(\text{안에 든 수})/N\)이고, 표준 오차는 \(4\sqrt{p(1-p)/N} = 1.642/\sqrt N\)이다(\(p=\pi/4\)).
다트 1만 개로 얻은 \(\hat\pi\)의 오차가 대략 0.016이다. 오차를 0.0016(10분의 1)로 줄이려면 다트가 몇 개 필요할까?
다트로 π 구하기 · 소수 둘째 자리 도전
함수의 적분
적분도 기댓값이다. \(U\sim U(a,b)\)이면 \(E[f(U)] = \frac{1}{b-a}\int_a^b f(x)\,dx\)이므로
손으로 그린 함수의 몬테카를로 적분
그렇다면 왜 정확한 수치 적분(사다리꼴·심프슨) 대신 몬테카를로를 쓸까? 1차원에서는 격자법이 압도적으로 낫다. 하지만 \(d\)차원 적분에서 점 \(N\)개를 격자로 깔면 축마다 \(N^{1/d}\)개밖에 못 놓아 중점법의 오차가 \(N^{-2/d}\)로 느려진다. 몬테카를로의 오차 \(N^{-1/2}\)는 차원과 무관하다. 공정 변수 수십 개가 동시에 흔들리는 수율 계산이 몬테카를로의 무대인 이유다. 격자법이 차원에 따라 급격히 나빠지는 현상을 차원의 저주(curse of dimensionality)라 한다.
시뮬레이션 한 번의 출력 표준 편차가 \(\sigma\) = 0.20일 때, 표본 4.0개로 얻은 평균의 표준 오차는 , 95% 오차 한계는 ± 이다. 오차 한계를 ± 0.0010로 만들려면 \(N = (1.96\sigma/e)^2\) = 개가 필요하다.
\(I = \int_0^1 e^{-x^2}\,dx\)를 표본 평균법으로 추정한다. \(f(U) = e^{-U^2}\)의 표준 편차가 \(\sigma_f = 0.2010\)일 때 (a) \(N = 10^4\)이면 표준 오차와 95% 오차 한계는? (b) 95% 오차 한계를 0.001로 만들려면 \(N\)은? (c) 명중–실패법(\([0,1]^2\)에 점을 뿌림)이면 같은 \(N=10^4\)에서 표준 오차는?
풀이 보기
아날로그 회로의 공정 변동을 넣은 몬테카를로 SPICE를 2000번 돌렸더니 37번이 사양을 벗어났다. 수율의 추정값, 표준 오차, 95% 신뢰 구간(정규 근사)을 구하라.
풀이 보기
SE \(=\sigma_g/\sqrt N\)에서 \(N\)을 늘리는 대신 \(\sigma_g\)를 줄이는 방법들이 있다. 층화 추출은 구간을 나눠 각 층에서 정해진 수만큼 뽑는다(라틴 하이퍼큐브 표집이 다차원판이다). 대조 변수는 \(U\)와 \(1-U\)를 짝지어 음의 상관으로 오차를 상쇄한다. 중요도 표집은 \(g\)가 큰 곳(희귀 불량 영역)에서 더 자주 뽑고 가중치 \(f/q\)로 보정한다. 몬테카를로 SPICE의 "고수율 추정" 옵션들이 이 원리다.
의사 난수와 재현성
몬테카를로의 재료는 난수인데, 컴퓨터는 결정론적 기계다. 실제로 쓰는 것은 정해진 규칙으로 계산하지만 통계적으로 무작위처럼 보이는 수열, 곧 의사 난수(pseudo-random number)다. 생성기는 내부 상태 \(s_k\)를 갱신 규칙 \(s_{k+1} = T(s_k)\)로 바꾸고 출력 \(u_k = G(s_k)\in[0,1)\)을 낸다. 처음 상태를 정하는 값이 시드(seed)다. 상태의 개수가 유한하므로 수열은 언젠가 반드시 되풀이되고, 되풀이되기까지의 길이를 주기(period)라 한다.
가장 오래되고 단순한 생성기가 선형 합동 생성기(linear congruential generator, LCG)다.
- \(m = 16\)이면 상태는 0–15의 16개뿐이다. 시드 7(초록)에서 출발한다.
- 규칙을 한 번씩 적용하면 7 → 6 → 1 → 8 → 11로 원 위를 건너뛴다.
- \(a = 5,\ c = 3\)은 헐–도벨 조건을 모두 만족해 16개 상태를 한 번씩 다 지난 뒤 7로 돌아온다(주기 16).
- \(a = 3\)으로 바꾸면 \(a-1 = 2\)가 4로 나뉘지 않아 주기가 8로 줄고, 상태 공간이 두 고리로 쪼개진다. 어느 고리를 도는지는 시드가 정한다.
\(m = 64\), \(a = 13\), \(c = 2\)인 LCG는 시드에 상관없이 주기가 64일까?
선형 합동 생성기의 고리와 격자
격자 구조의 악명 높은 예가 1960–70년대 IBM 컴퓨터에 널리 쓰인 RANDU(\(a = 65539,\ c = 0,\ m = 2^{31}\))다. \(a = 2^{16}+3\)이라서 \(x_{k+2} = 6x_{k+1}-9x_k \pmod{2^{31}}\)이 성립하고, 연속한 세 수 \((u_k, u_{k+1}, u_{k+2})\)는 정육면체 안의 평면 15장 위에만 놓인다. 이 생성기로 계산한 수많은 3차원 시뮬레이션 결과가 나중에 의심받았다.
현대의 범용 생성기는 메르센 트위스터(주기 \(2^{19937}-1\), 파이썬·R·MATLAB의 기본), PCG, xoshiro 계열이다. 이 책의 시뮬레이터는 짧고 빠른 mulberry32를 쓴다(`Stat.rng(seed)`). 32비트 상태 하나를 더하고 비트를 섞는 방식으로 주기는 \(2^{32}\approx4.3\times10^9\)이다. 교과서 실험에는 충분하지만, 10억 개 넘는 난수를 쓰는 대규모 시뮬레이션에는 짧다.
상태가 32비트인 생성기의 주기는 많아야 \(2^{b}\) = 이다. 초당 9개씩 난수를 쓰면 만에 수열이 한 바퀴 돈다. 경험적으로는 주기의 제곱근(개)쯤부터 생일 문제 같은 통계 검정에서 흠이 드러나므로, 이 생성기로 안심하고 쓸 수 있는 난수는 그 정도다.
의사 난수의 가장 큰 장점은 재현성(reproducibility)이다. 같은 시드로 시작하면 같은 수열이 나오므로, 몬테카를로 결과가 이상할 때 그 실행을 그대로 되살려 원인을 찾을 수 있다. 시뮬레이션 보고서에는 생성기 이름과 시드를 기록하고, 비교 실험(설계 A 대 B)에서는 같은 시드를 써서 난수 차이가 아닌 설계 차이만 남게 한다(공통 난수법).
① 표본마다 현재 시각으로 시드를 다시 정하기: 같은 밀리초 안의 시드가 같아져 똑같은 "무작위" 표본이 반복된다. 시드는 실험 시작에 한 번만 정한다. ② 병렬 실행 100개에 시드 1, 2, …, 100을 주기: 단순한 생성기는 이웃한 시드의 수열이 서로 상관될 수 있다. 병렬 흐름을 지원하는 생성기(PCG의 스트림, 파이썬 `SeedSequence.spawn`)를 쓴다. ③ 결과가 마음에 들 때까지 시드를 바꾸기: 그것은 표본을 고르는 것이고, 11장의 다중 비교 문제와 같다.
균일 난수 \(U\)만 있으면 다른 분포는 변환으로 만든다. 5장의 역변환법 \(X = F^{-1}(U)\)이 가장 일반적이고(`Stat.dist.*.sample`이 이 방식), 정규 난수는 두 균일 난수를 쓰는 극 좌표법(마르사글리아, `Stat.rng().normal`)이나 박스–뮐러 변환 \(Z = \sqrt{-2\ln U_1}\cos(2\pi U_2)\)으로 만든다. 갈톤 보드와 표본 평균 실험도 결국 이렇게 만든 균일 난수 위에 서 있다.
핵심 정리
- 마르코프 \(P(X\ge a)\le E[X]/a\)(\(X\ge0\)), 체비쇼프 \(P(|X-\mu|\ge k\sigma)\le1/k^2\). 분포를 몰라도 성립하는 대신 느슨하다(정규 ±2σ: 실제 0.046, 상한 0.25). 세 점 분포가 등호를 만들므로 평균·분산만으로는 더 좋은 상한이 없다.
- \(E[\bar X_n] = \mu\), \(\operatorname{Var}(\bar X_n) = \sigma^2/n\). 약한 큰 수의 법칙 \(P(|\bar X_n-\mu|\ge\varepsilon)\le\sigma^2/(n\varepsilon^2)\to0\)(확률 수렴), 강한 법칙은 경로마다 \(\bar X_n\to\mu\)(거의 확실한 수렴). 평균이 없는 코시 분포에서는 성립하지 않는다. 비율은 수렴해도 개수의 차이 \(|S_n-n\mu|\)는 \(\sqrt n\)으로 자란다.
- 중심 극한 정리: 유한 분산 iid이면 \((\bar X_n-\mu)/(\sigma/\sqrt n)\xrightarrow{d}N(0,1)\). 실무 근사 \(\bar X\;\dot\sim\;N(\mu,\sigma^2/n)\), \(S_n\;\dot\sim\;N(n\mu,n\sigma^2)\). 수렴 속도는 왜도 \(\gamma/\sqrt n\)이 정한다 — "n ≥ 30"은 보장이 아니다.
- 갈톤 보드의 칸 번호는 \(\mathrm{Bin}(R,p)\)이고 정규로 간다. 많은 독립 소요인의 합(\(V_t\), CD, 두께)은 정규에 가깝고 \(\sigma = \sqrt{\sum\sigma_i^2}\)(RSS). 한 요인이 지배하면 정규성이 깨지고, 곱으로 쌓이면 로그 정규가 된다.
- \(\mathrm{Bin}(n,p)\approx N(np,np(1-p))\)(\(np, n(1-p)\ge5\)), \(\mathrm{Poisson}(\lambda)\approx N(\lambda,\lambda)\)(\(\lambda\gtrsim10\)). 연속성 수정: 구하는 정수 막대를 모두 덮도록 경계를 ±0.5 옮긴다. \(P(X = k)\)도 \([k-0.5,k+0.5]\)의 넓이로 구한다. 먼 꼬리에서는 치우침 때문에 근사가 나빠진다.
- 몬테카를로: \(\theta = E[g(X)]\approx\frac1N\sum g(X_i)\), \(\operatorname{SE} = \sigma_g/\sqrt N\), 필요한 표본 수 \(N = (z\sigma_g/e)^2\). 오차를 1/10로 줄이려면 100배. 적분 \(\int_a^bf = (b-a)E[f(U)]\). 오차가 차원과 무관해 고차원·수율 계산에 쓴다. 희귀 사건은 \(N\approx(1-p)/(p\,r^2)\)이 필요하므로 분산 감소 기법을 쓴다.
- 의사 난수는 시드로 정해지는 결정론적 수열이다. LCG \(x_{k+1} = (ax_k+c)\bmod m\)는 헐–도벨 조건에서 주기 \(m\)이지만 연속한 수가 격자 위에 놓인다(RANDU의 평면 15장). 시드를 기록하면 실험을 재현할 수 있고, 비교 실험에는 공통 난수를 쓴다.
연습문제
표준 정규 분포 값은 부록 A의 계산기나 표를 쓴다. 따로 말이 없으면 답이 ±1% 이내면 정답이다(표를 읽는 문제는 ±2%).
세정 후 웨이퍼 한 장의 파티클 수는 평균 4개다. 분포를 모를 때 한 장에 파티클이 20개 이상일 확률의 상한은?
마르코프 부등식: \(P(X\ge20)\le E[X]/20 = 4/20 = 0.2\).
어떤 분포든 평균 ±3σ 안에 들어가는 비율은 적어도 얼마인가?
\(P(|X-\mu|<3\sigma)\ge 1-1/3^2 = 8/9 = 0.889\). 정규 분포라면 0.9973이다.
표준 편차 4 nm인 측정을 100번 반복해 평균을 낸다. 표본 평균이 참값에서 1 nm 이상 벗어날 확률의 체비쇼프 상한과, CLT로 근사한 값을 구하라.
\(\operatorname{Var}(\bar X) = 16/100 = 0.16\). 체비쇼프: \(P(|\bar X-\mu|\ge1)\le0.16/1^2 = 0.16\). CLT: \(z = 1/0.4 = 2.5\), \(2(1-\Phi(2.5)) = 2(0.00621) = 0.0124\).
룰렛에서 빨강이 10번 연속 나왔다. 큰 수의 법칙에 비추어 옳은 말은?
독립 시행에는 기억이 없다. 큰 수의 법칙은 비율이 수렴한다는 것이지 앞선 치우침을 "되갚는다"는 것이 아니다. 실제로 빨강 개수와 기대 개수의 차이는 판이 늘수록 \(\sqrt n\)처럼 커지는 것이 보통이다(도박사의 오류).
평균 50, 표준 편차 10인 모집단에서 크기 25인 표본을 뽑는다. 표본 평균이 53보다 클 확률은?
\(\operatorname{SD}(\bar X) = 10/5 = 2\), \(z = (53-50)/2 = 1.5\), \(P = 1-\Phi(1.5) = 0.0668\).
공정한 주사위 48개를 던진 눈의 합 \(S\)가 180 이상일 확률을 연속성 수정한 정규 근사로 구하라. (주사위 하나의 평균 3.5, 분산 35/12)
\(\sigma_S = \sqrt{140} = 11.83\), \(z = (179.5-168)/11.83 = 0.972\), \(P\approx1-\Phi(0.972) = 0.166\). 합성곱으로 구한 정확한 값은 0.166이다.
CD 측정의 표준 편차가 5 nm다. 표본 평균이 참값에서 1 nm 이내에 들 확률을 0.95 이상으로 하려면 최소 몇 개를 재야 하는가? (\(z_{0.025} = 1.96\))
\(1.96\times5/\sqrt n\le1\)에서 \(n\ge(9.8)^2 = 96.04\). 정수로 올려 \(n = 97\).
\(X\sim\mathrm{Bin}(100, 0.1)\)에서 \(P(X\le15)\)를 연속성 수정한 정규 근사로 구하라.
\(\mu = 10\), \(\sigma = \sqrt{9} = 3\). \(z = (15.5-10)/3 = 1.833\), \(\Phi(1.833) = 0.967\). 정확한 값은 0.960이다(\(np = 10\)이라 치우침이 조금 남는다).
\(X\sim\mathrm{Bin}(100, 0.5)\)에서 \(P(X = 50)\)을 연속성 수정한 정규 근사로 구하라.
\(\mu = 50\), \(\sigma = 5\). \(P(49.5\le Y\le50.5) = \Phi(0.1)-\Phi(-0.1) = 0.5398-0.4602 = 0.0797\). 정확한 값 \(\binom{100}{50}/2^{100} = 0.0796\).
한 픽셀에 흡수되는 광자 수가 평균 64인 포아송 분포를 따른다. 광자가 50개 이하일 확률을 연속성 수정한 정규 근사로 구하라.
\(\sigma = 8\), \(z = (50.5-64)/8 = -1.6875\), \(\Phi(-1.69) = 0.0458\). 정확한 포아송 값은 0.0417로, 왼쪽 꼬리에서 정규 근사가 10% 과대평가한다.
다트 \(N = 10\,000\)개로 \(\pi\)를 추정할 때 \(\hat\pi = 4\times(\text{사분원 안의 비율})\)의 표준 오차는?
\(p = \pi/4 = 0.7854\), \(\operatorname{SE}(\hat\pi) = 4\sqrt{p(1-p)/N} = 4\sqrt{0.7854\times0.2146}/100 = 4(0.4105)/100 = 0.0164\).
회로 수율이 약 90%로 예상된다. 몬테카를로 시뮬레이션으로 수율을 95% 신뢰도에서 ±0.5%p 이내로 추정하려면 실행이 몇 번 필요한가? (\(z = 1.96\))
\(N = z^2p(1-p)/e^2 = 1.96^2(0.9)(0.1)/0.005^2 = 13\,830\).
SRAM 셀의 불량 확률이 \(p = 10^{-4}\)쯤이다. 단순 몬테카를로로 \(\hat p\)의 상대 표준 오차 \(\operatorname{SE}(\hat p)/p\)를 10%로 만들려면 실행이 몇 번 필요한가? 이 결과가 6σ(\(p\approx10^{-9}\)) 셀에 대해 뜻하는 것은?
\(\sqrt{p(1-p)/N} = 0.1p\)에서 \(N = (1-p)/(0.01p) = 0.9999/10^{-6} \approx 1.00\times10^6\). 불량 사건을 약 100개 보아야 한다는 뜻이다(\(Np = 100\)). \(p = 10^{-9}\)이면 \(10^{11}\)번이 필요해 단순 몬테카를로로는 불가능하고, 중요도 표집 같은 분산 감소 기법을 써야 한다.
LCG \(x_{k+1} = (5x_k+3)\bmod16\)을 시드 \(x_0 = 7\)로 시작한다. \(x_3\)과 이 수열의 주기를 구하라. 헐–도벨 조건으로 주기를 확인하라.
\(x_1 = 38\bmod16 = 6\), \(x_2 = 33\bmod16 = 1\), \(x_3 = 8\). 계속하면 11, 10, 5, 12, 15, 14, 9, 0, 3, 2, 13, 4, 7로 16번째에 시드로 돌아온다. 헐–도벨: \(\gcd(3,16) = 1\), \(16\)의 소인수 2가 \(a-1 = 4\)를 나눔, \(4\mid16\)이고 \(4\mid4\) — 세 조건 모두 만족하므로 주기는 \(m = 16\)이다.
10단 갈톤 보드(\(p = 0.5\))에서 구슬이 가운데 칸(\(k = 5\))에 들어갈 확률을 이항 분포로 정확히, 그리고 연속성 수정한 정규 근사로 구하라.
정확: \(\binom{10}{5}/2^{10} = 252/1024 = 0.246\). 근사: \(\mu = 5\), \(\sigma = \sqrt{2.5} = 1.581\), \(\Phi(0.5/1.581)-\Phi(-0.5/1.581) = 2\Phi(0.316)-1 = 0.248\).
표준 코시 분포를 따르는 측정 1000개의 평균에 대해 옳은 것은?
코시 분포는 평균과 분산이 없어 큰 수의 법칙과 CLT의 전제가 깨진다. 독립 코시 \(n\)개의 평균은 다시 표준 코시다(특성 함수 \(e^{-|t|}\)가 \((e^{-|t|/n})^n = e^{-|t|}\)로 보존). 극단값 하나가 평균 전체를 끌고 간다. 이런 자료에는 중앙값처럼 강건한 통계량을 쓴다.
개념 확인 퀴즈
중심 극한 정리에 대한 설명으로 옳은 것은?
\(X\sim\mathrm{Bin}(50, 0.4)\)에서 \(P(X<25)\)를 연속성 수정한 정규 근사로 구할 때 쓰는 경계는?
몬테카를로 추정의 표준 오차를 절반으로 줄이려면?
평균과 분산만 알려진 분포에서 \(P(|X-\mu|\ge2\sigma)\)에 대해 확실히 말할 수 있는 것은?
몬테카를로 시뮬레이션의 결과를 동료가 똑같이 재현하려면 무엇을 알려 주어야 하는가?