1 개요
지금까지 확률 표본의 이론적 성질(표본분포, 수렴 등)을 다루었다. 이 포스트에서는 방향을 뒤집어, 주어진 분포에서 확률 표본을 컴퓨터로 생성하는 방법을 다룬다.
난수 생성은 순수한 이론 문제가 아니다. Monte Carlo 시뮬레이션, 부트스트랩, MCMC, 베이지안 계산 — 현대 통계학의 핵심 도구들이 모두 “특정 분포에서 표본을 뽑는 능력”에 의존한다 (Casella & Berger, 2002, Ch.5).
출발점은 단순하다: 균일분포 \(U \sim \text{Uniform}(0, 1)\) 난수를 생성할 수 있다고 가정한다. 문제는 이 균일 난수를 원하는 분포의 난수로 변환하는 것이다.
2 방법론의 전체 구조
균일 난수 U ~ Uniform(0,1)
↓
┌─── 직접 방법 (Direct Methods) ───────────────────────┐
│ │
│ 역변환법: Y = F^{-1}(U) │
│ → CDF 역함수가 존재하면 정확 │
│ │
│ Box-Muller: 균일 2개 → 정규 2개 │
│ → 정규분포 전용, 매우 효율적 │
│ │
│ 분포 간 관계 이용: 지수 → 감마, 카이제곱, 베타 │
│ → 지수분포를 기반으로 다양한 분포 생성 │
│ │
├─── 간접 방법 (Indirect Methods) ─────────────────────┤
│ │
│ 수락-기각 알고리즘: 봉투 분포에서 뽑고 조건부 수락 │
│ → 범용적, M < ∞ 필요 │
│ │
│ Metropolis 알고리즘: 마르코프 체인으로 수렴 │
│ → M < ∞ 불필요, MCMC의 기초 │
│ │
└───────────────────────────────────────────────────────┘3 직접 방법 1: 역변환법 (Inverse Transform Method)
3.1 원리
확률적분변환(Probability Integral Transform, Theorem 2.1.10)의 역이다:
\(Y\) 가 CDF \(F_Y\) 를 가진 연속 확률변수이면:
\[ Y = F_Y^{-1}(U), \quad U \sim \text{Uniform}(0, 1) \]
는 분포 \(F_Y\) 를 따른다.
증명: \(P(Y \leq y) = P(F_Y^{-1}(U) \leq y) = P(U \leq F_Y(y)) = F_Y(y)\) . \(\square\)
직관: CDF \(F_Y\) 는 \([0, 1]\) 로의 단조증가 함수이다. 균일 난수 \(U\) 를 \(y\) -축에 놓고 수평선을 그어 CDF와 만나는 점을 찾으면, 그 \(y\) 좌표가 원하는 분포의 난수이다.
3.2 예시: 지수분포
\(Y \sim \text{Exponential}(\lambda)\) 이면 \(F_Y(y) = 1 - e^{-y/\lambda}\) 이므로:
\[ F_Y^{-1}(u) = -\lambda \log(1 - u) \]
\(U\) 와 \(1 - U\) 가 같은 분포이므로 \(Y = -\lambda \log(U)\) 로 단순화할 수 있다.
3.3 지수분포에서 파생되는 분포
지수분포를 기반으로 다양한 분포를 생성할 수 있다. \(U_j\) 가 iid \(\text{Uniform}(0, 1)\) 이면 \(Y_j = -\lambda\log(U_j)\) 는 iid \(\text{Exponential}(\lambda)\) 이고:
| 변환 | 결과 분포 |
|---|---|
| \(Y = -2\sum_{j=1}^{\nu} \log(U_j)\) | \(\chi^2_{2\nu}\) |
| \(Y = -\beta\sum_{j=1}^{a} \log(U_j)\) | \(\text{Gamma}(a, \beta)\) (정수 \(a\) ) |
| \(Y = \frac{\sum_{j=1}^{a} \log(U_j)}{\sum_{j=1}^{a+b} \log(U_j)}\) | \(\text{Beta}(a, b)\) (정수 \(a, b\) ) |
이 변환들은 모두 지수분포-균일분포 변환에서 출발하며, 감마분포가 지수분포의 합이라는 사실과 베타분포와 감마분포의 관계에 기반한다.
3.4 예시: 와이블 분포
\(Y \sim \text{Weibull}(k, \lambda)\) 이면 \(F_Y(y) = 1 - e^{-(y/\lambda)^k}\) 이므로:
\[ F_Y^{-1}(u) = \lambda[-\log(1 - u)]^{1/k} \]
와이블 분포는 생존분석과 신뢰성 공학에서 핵심적으로 사용되며, 역변환법이 깔끔하게 적용되는 좋은 예이다. \(k = 1\) 이면 지수분포로 환원된다.
3.5 예시: 로그정규 분포
\(Y \sim \text{LogNormal}(\mu, \sigma^2)\) 는 \(\log Y \sim N(\mu, \sigma^2)\) 로 정의된다. 따라서:
- Box-Muller로 \(Z \sim N(0, 1)\) 을 생성한다
- \(Y = e^{\mu + \sigma Z}\) 로 변환한다
이것은 역변환법과 Box-Muller의 조합이다. 정규분포의 역함수가 닫힌 형태가 없으므로 역변환법 단독으로는 불가능하지만, Box-Muller를 경유하면 간단하다.
3.6 직접 방법의 적용 가능성 판단
| CDF 역함수 | 적용 가능한 방법 | 예시 |
|---|---|---|
| 닫힌 형태 존재 | 역변환법 직접 적용 | 지수, 와이블, 파레토, 코시 |
| 닫힌 형태 없음 + 특수 구조 | 분포 간 관계 이용 | \(\chi^2_{2\nu}\) (지수합), 로그정규 (Box-Muller + 변환) |
| 닫힌 형태 없음 + 구조 없음 | 간접 방법 필요 | \(\chi^2_1\) , 일반 베타, 절단 정규 |
3.7 역변환법의 한계
CDF의 역함수가 닫힌 형태(closed form)로 존재하지 않으면 각 난수 생성마다 적분 방정식을 풀어야 한다. 예를 들어 \(\chi^2_1\) (자유도 1인 카이제곱)의 CDF 역함수는 닫힌 형태가 없으므로 역변환법이 비실용적이다. 이런 경우 다른 방법이 필요하다.
3.8 이산분포에 대한 역변환법
이산 확률변수 \(Y\) 가 값 \(y_1 < y_2 < \cdots\) 를 가질 때:
- \(U \sim \text{Uniform}(0, 1)\) 을 생성한다
- \(F_Y(y_i) < U \leq F_Y(y_{i+1})\) 이면 \(Y = y_{i+1}\) 로 설정한다
예시: \(Y \sim \text{Binomial}(4, 5/8)\) 을 생성하려면:
| \(y\) | \(P(Y = y)\) | \(F_Y(y)\) | \(U\) 의 범위 |
|---|---|---|---|
| 0 | 0.020 | 0.020 | \((0, 0.020]\) |
| 1 | 0.132 | 0.152 | \((0.020, 0.152]\) |
| 2 | 0.329 | 0.481 | \((0.152, 0.481]\) |
| 3 | 0.366 | 0.847 | \((0.481, 0.847]\) |
| 4 | 0.153 | 1.000 | \((0.847, 1.000]\) |
이 방법은 범위가 무한인 분포(포아송, 음이항)에도 적용 가능하다. 실무적으로는 평균 근처부터 탐색을 시작하면 효율이 높아진다.
4 직접 방법 2: Box-Muller 알고리즘
역변환법으로는 표준정규분포 \(N(0, 1)\) 을 직접 생성하기 어렵다 (CDF 역함수가 닫힌 형태 없음). Box-Muller 알고리즘은 이 문제를 우아하게 해결한다.
\(U_1, U_2 \overset{iid}{\sim} \text{Uniform}(0, 1)\) 에 대해:
\[ R = \sqrt{-2\log U_2}, \quad \Theta = 2\pi U_1 \]
로 놓으면:
\[ X = R\cos\Theta, \quad Y = R\sin\Theta \]
은 독립인 \(N(0, 1)\) 확률변수이다 (Casella & Berger, 2002, Ch.5).
직관: 2차원 표준정규분포 \((X, Y)\) 를 극좌표 \((R, \Theta)\) 로 변환하면, \(R^2 \sim \text{Exponential}(1/2)\) 이고 \(\Theta \sim \text{Uniform}(0, 2\pi)\) 이며 둘은 독립이다. Box-Muller는 이 변환의 역이다.
단일 정규 난수를 위한 직접 변환은 없지만, 2개를 동시에 생성하는 변환은 존재한다. 이것이 Box-Muller의 핵심 아이디어이다.
5 간접 방법 1: 수락-기각 알고리즘 (Accept/Reject Algorithm)
직접 변환이 없을 때 사용하는 범용적 방법이다.
5.1 기본 아이디어
목표 분포 \(f_Y\) 에서 직접 표본을 뽑을 수 없지만, 후보 분포 \(f_V\) 에서는 뽑을 수 있다고 하자. \(f_V\) 에서 뽑은 후보 \(V\) 가 \(f_Y\) 와 “닮았는지” 검정하여 수락하거나 기각한다.
목표 밀도 \(f_Y\) 와 후보 밀도 \(f_V\) 가 같은 지지(support)를 가지고:
\[ M = \sup_y \frac{f_Y(y)}{f_V(y)} < \infty \]
이면, 다음 알고리즘은 \(Y \sim f_Y\) 를 생성한다:
- \(U \sim \text{Uniform}(0, 1)\) , \(V \sim f_V\) 를 독립으로 생성한다
- \(U < \frac{1}{M}\frac{f_Y(V)}{f_V(V)}\) 이면 \(Y = V\) 로 설정 (수락); 아니면 1로 돌아감 (기각)
(Casella & Berger, 2002, Ch.5)
5.2 증명
생성된 \(Y\) 의 CDF는:
\[ \begin{aligned} P(Y \leq y) &= P(V \leq y \mid \text{수락}) \\ &= \frac{P(V \leq y, U < \frac{1}{M}f_Y(V)/f_V(V))}{P(U < \frac{1}{M}f_Y(V)/f_V(V))} \\ &= \frac{\int_{-\infty}^{y} \frac{1}{M}\frac{f_Y(v)}{f_V(v)} f_V(v) \, dv}{\frac{1}{M}} \\ &= \int_{-\infty}^{y} f_Y(v) \, dv = F_Y(y) \quad \square \end{aligned} \]
5.3 효율성
수락까지 필요한 시행 횟수 \(N \sim \text{Geometric}(1/M)\) 이므로 \(EN = M\) 이다. \(M\) 이 작을수록(후보 분포가 목표 분포에 가까울수록) 효율적이다.
| 후보 분포 선택 | \(M\) | 효율 |
|---|---|---|
| \(f_V = \text{Uniform}\) , \(f_Y = \text{Beta}(2.7, 6.3)\) | 2.67 | 37.5% 수락률 |
| \(f_V = \text{Beta}(2, 6)\) , \(f_Y = \text{Beta}(2.7, 6.3)\) | 1.67 | 59.9% 수락률 |
후보 분포를 목표 분포에 가깝게 선택하면 \(M\) 이 줄어든다. 그러나 후보 분포에서 표본을 생성하는 비용도 고려해야 한다.
\(M = \sup f_Y/f_V < \infty\) 는 후보 분포의 꼬리가 목표 분포의 꼬리보다 두꺼워야 함을 의미한다. 코시 분포(두꺼운 꼬리)에서 정규분포(목표)를 생성할 때는 \(M < \infty\) 이지만, 정규분포에서 코시분포를 생성하려면 \(M = \infty\) 여서 수락-기각이 불가능하다. 이런 경우 Metropolis 알고리즘이 필요하다.
6 간접 방법 2: Metropolis 알고리즘
수락-기각 알고리즘의 \(M < \infty\) 조건이 만족되지 않을 때 사용한다. MCMC(Markov Chain Monte Carlo)의 기초이다.
목표 밀도 \(f_Y\) , 후보 밀도 \(f_V\) (같은 지지):
- \(V \sim f_V\) 를 생성하여 \(Z_0 = V\) 로 설정
\(i = 1, 2, \ldots\) 에 대해:
- \(U_i \sim \text{Uniform}(0, 1)\) , \(V_i \sim f_V\) 를 생성하고:
\[ \rho_i = \min\left\{\frac{f_Y(V_i)}{f_V(V_i)} \cdot \frac{f_V(Z_{i-1})}{f_Y(Z_{i-1})}, \; 1\right\} \]
\(\rho_i\) 는 두 비율의 곱이다. 각각의 역할을 따로 보면:
- \(\dfrac{f_Y(V_i)}{f_V(V_i)}\): 후보 \(V_i\) 가 “목표 분포 대비 후보 분포에서 얼마나 과다 추출되었는가” — 목표 분포에서 밀도가 낮고 후보 분포에서 밀도가 높은 곳(과다 추출된 곳)은 이 비율이 작아 수락이 어렵다.
- \(\dfrac{f_V(Z_{i-1})}{f_Y(Z_{i-1})}\): 현재 위치 \(Z_{i-1}\) 이 “후보 분포 대비 목표 분포에서 얼마나 과다 추출되었는가” — 현재 위치가 목표 분포에서 과다 추출된 곳(목표 밀도가 높은 곳)이라면 이동을 더 쉽게 허용한다.
두 비율의 곱, 즉 \(\frac{f_Y(V_i)/f_Y(Z_{i-1})}{f_V(V_i)/f_V(Z_{i-1})}\) 은 “후보 쪽에서 목표 분포로의 이동이 얼마나 적절한가”를 측정한다. 이 비율이 1보다 크면 항상 이동, 작으면 그 비율 확률로 이동하는 것이 마르코프 체인을 \(f_Y\) 의 정상분포(stationary distribution)로 수렴시킨다.
- \(U_i \leq \rho_i\) 이면 \(Z_i = V_i\) (이동), 아니면 \(Z_i = Z_{i-1}\) (유지)
\(i \to \infty\) 이면 \(Z_i \xrightarrow{d} Y\) (Casella & Berger, 2002, Ch.5).
6.1 수락-기각과의 차이
| 수락-기각 | Metropolis | |
|---|---|---|
| \(M < \infty\) 필요 | 필요 | 불필요 |
| 출력 | 독립 표본 | 상관된 마르코프 체인 |
| 수렴 | 즉시 (정확한 분포) | 점근적 (분포수렴) |
| 초기값 의존 | 없음 | 있음 (burn-in 필요) |
Metropolis의 핵심 아이디어: 현재 위치 \(Z_{i-1}\) 에서 후보 \(V_i\) 로 “이동할지”를 \(\rho_i\) 로 결정한다. \(f_Y(V_i)/f_V(V_i)\) 가 크면(목표 분포에서 밀도가 높은 곳) 이동 확률이 높고, 작으면 현재 위치에 머문다. 이 과정이 반복되면 체인은 \(f_Y\) 에 수렴한다.
6.2 Metropolis-Hastings 확장
Metropolis 알고리즘은 대칭 후보 분포 \(f_V(v|z) = f_V(z|v)\) 를 가정한다. Hastings (1970)는 비대칭 후보 분포로 일반화했다:
\[ \rho_i = \min\left\{\frac{f_Y(V_i)}{f_Y(Z_{i-1})} \cdot \frac{f_V(Z_{i-1} | V_i)}{f_V(V_i | Z_{i-1})}, \; 1\right\} \]
대칭 후보( \(f_V(v|z) = f_V(z|v)\) )이면 두 번째 비율이 1이 되어 원래 Metropolis로 환원된다. 실무에서 가장 흔한 선택은 Random Walk Metropolis: \(V_i = Z_{i-1} + \epsilon\) , \(\epsilon \sim N(0, \sigma^2_{\text{step}})\) 이다. 이 경우 \(\sigma_{\text{step}}\) (걸음 크기)의 선택이 수렴 속도를 결정한다:
| \(\sigma_{\text{step}}\) | 행동 | 문제 |
|---|---|---|
| 너무 작음 | 거의 항상 수락, 느리게 이동 | 공간 탐색 느림 |
| 너무 큼 | 자주 기각, 제자리 머무름 | 수락률 낮음 |
| 적절함 | 수락률 \(\approx\) 23-44% | 최적 탐색 |
Roberts et al. (1997)은 목표 분포가 다변량 정규에 가까울 때 최적 수락률이 약 23.4%임을 보였다.
7 간접 방법 3: Gibbs Sampler
다변량 분포에서 결합분포 \(f(x_1, \ldots, x_k)\) 에서 직접 표본을 뽑기 어렵지만, 조건부 분포 \(f(x_i | x_{-i})\) 에서는 뽑을 수 있는 경우에 사용한다.
결합분포 \(f(x_1, \ldots, x_k)\) 에서 표본을 생성하려면:
- 초기값 \((x_1^{(0)}, \ldots, x_k^{(0)})\) 를 설정한다
각 반복 \(t = 1, 2, \ldots\) 에서:
- \(x_1^{(t)} \sim f(x_1 | x_2^{(t-1)}, x_3^{(t-1)}, \ldots, x_k^{(t-1)})\)
- \(x_2^{(t)} \sim f(x_2 | x_1^{(t)}, x_3^{(t-1)}, \ldots, x_k^{(t-1)})\)
- \(\vdots\)
- \(x_k^{(t)} \sim f(x_k | x_1^{(t)}, x_2^{(t)}, \ldots, x_{k-1}^{(t)})\)
\(t \to \infty\) 이면 \((x_1^{(t)}, \ldots, x_k^{(t)}) \xrightarrow{d} (X_1, \ldots, X_k) \sim f\) .
Gibbs Sampler는 Metropolis-Hastings의 특수한 경우로, 수락 확률이 항상 1이다 (조건부 분포에서 직접 추출하므로 기각이 없다). 베이지안 추론에서 사후분포의 표본 추출에 핵심적으로 사용된다.
7.1 예시: 이변량 정규분포
\((X, Y)\) 가 상관계수 \(\rho\) 인 이변량 정규분포를 따를 때, 조건부 분포는:
\[ X | Y = y \sim N(\rho y, 1 - \rho^2), \quad Y | X = x \sim N(\rho x, 1 - \rho^2) \]
Gibbs Sampler는 이 두 조건부 정규분포를 번갈아 추출한다. 물론 이변량 정규는 직접 생성할 수 있으므로 Gibbs가 불필요하지만, 원리를 보여주는 교과서적 예시이다.
8 간접 방법 4: Importance Sampling
Importance Sampling은 난수 생성이라기보다 기댓값 추정의 방법이지만, 간접 표본 추출과 밀접하게 관련된다.
\(Eh(Y)\) 를 계산하고 싶지만 \(f_Y\) 에서 직접 추출이 어려울 때, 추출 가능한 \(f_V\) 를 사용한다:
\[ Eh(Y) = \int h(y) f_Y(y) \, dy = \int h(y) \frac{f_Y(y)}{f_V(y)} f_V(y) \, dy = E_V\left[h(V) \frac{f_Y(V)}{f_V(V)}\right] \]
따라서 \(V_1, \ldots, V_n \overset{iid}{\sim} f_V\) 를 추출하고:
\[ \widehat{Eh(Y)} = \frac{1}{n}\sum_{i=1}^{n} h(V_i) \frac{f_Y(V_i)}{f_V(V_i)} \]
비율 \(w_i = f_Y(V_i)/f_V(V_i)\) 를 중요도 가중치(importance weight)라 한다.
| 수락-기각 | Importance Sampling | |
|---|---|---|
| 출력 | 목표 분포의 표본 | 가중 평균 (기댓값 추정) |
| 표본 낭비 | 있음 (기각된 표본) | 없음 (모든 표본 사용) |
| \(M < \infty\) 필요 | 필요 | 불필요 (단, 분산이 커질 수 있음) |
| 용도 | 분포에서 표본 필요 | 기댓값/확률 추정만 필요 |
9 후보 분포 선택 전략
간접 방법의 성능은 후보 분포 \(f_V\) 의 선택에 크게 의존한다. 실무적 지침을 정리한다.
| 원칙 | 이유 | 예시 |
|---|---|---|
| 꼬리가 목표보다 두꺼워야 | \(M < \infty\) 보장 | 정규 목표 → \(t\) 또는 코시 후보 |
| 모양이 비슷해야 | \(M\) 이 작아 효율적 | Beta(2.7, 6.3) 목표 → Beta(2, 6) 후보 |
| 생성이 쉬워야 | 전체 계산 비용 최소화 | 가능하면 정규, 감마, 균일 |
| 밀도 비율 계산이 쉬워야 | \(f_Y/f_V\) 평가 비용 | 정규화 상수를 알 필요 없는 경우 유리 |
베이지안 추론에서 사후분포 \(f_Y(y) \propto \mathcal{L}(\theta|x)\pi(\theta)\) 는 정규화 상수를 모르는 경우가 대부분이다. Metropolis-Hastings와 Gibbs Sampler는 밀도의 비율만 사용하므로 정규화 상수가 소거된다. 이것이 MCMC가 베이지안 추론에서 핵심적인 이유이다.
10 시뮬레이션 검증과 대수의 법칙
난수 생성의 정당성은 대수의 법칙에 의해 보장된다. \(Y_1, \ldots, Y_n\) 이 iid이면:
\[ \frac{1}{n}\sum_{i=1}^{n} h(Y_i) \xrightarrow{P} Eh(Y) \quad (n \to \infty) \]
따라서 기댓값, 분산, 확률 등 \(Eh(Y)\) 형태로 표현 가능한 모든 양을 시뮬레이션으로 근사할 수 있다. 이것이 Monte Carlo Simulation의 수학적 기반이다.
11 코드 예시
11.1 Step 1: 순수 Python 구현 (역변환법 + Box-Muller)
import math
import random
random.seed(42)
n = 10000
# --- 1. 역변환법: 지수분포 ---
lam = 2.0
exp_samples = [-lam * math.log(random.random()) for _ in range(n)]
mean_exp = sum(exp_samples) / n
var_exp = sum((x - mean_exp)**2 for x in exp_samples) / (n - 1)
print(f"=== 역변환법: 지수분포(lambda={lam}) ===")
print(f" E[Y] = {mean_exp:.4f} (이론: {lam})")
print(f" Var[Y] = {var_exp:.4f} (이론: {lam**2})")
# --- 2. Box-Muller: 정규분포 ---
normal_samples = []
for _ in range(n // 2):
u1 = random.random()
u2 = random.random()
r = math.sqrt(-2 * math.log(u2))
theta = 2 * math.pi * u1
normal_samples.append(r * math.cos(theta))
normal_samples.append(r * math.sin(theta))
mean_norm = sum(normal_samples) / len(normal_samples)
var_norm = sum((x - mean_norm)**2 for x in normal_samples) / (len(normal_samples) - 1)
print(f"\n=== Box-Muller: 정규분포(0,1) ===")
print(f" E[X] = {mean_norm:.4f} (이론: 0)")
print(f" Var[X] = {var_norm:.4f} (이론: 1)")
# --- 3. 지수 → 감마 변환 ---
alpha = 5 # 정수
beta = 2.0
gamma_samples = []
for _ in range(n):
y = -beta * sum(math.log(random.random()) for _ in range(alpha))
gamma_samples.append(y)
mean_gam = sum(gamma_samples) / n
var_gam = sum((x - mean_gam)**2 for x in gamma_samples) / (n - 1)
print(f"\n=== 지수합 → 감마({alpha},{beta}) ===")
print(f" E[Y] = {mean_gam:.4f} (이론: {alpha * beta})")
print(f" Var[Y] = {var_gam:.4f} (이론: {alpha * beta**2})")
# --- 4. 이산분포 역변환: 포아송 ---
lam_pois = 7.0
pois_samples = []
for _ in range(n):
u = random.random()
y = 0
cumprob = math.exp(-lam_pois)
while u > cumprob:
y += 1
cumprob += math.exp(-lam_pois) * lam_pois**y / math.factorial(y)
pois_samples.append(y)
mean_pois = sum(pois_samples) / n
var_pois = sum((x - mean_pois)**2 for x in pois_samples) / (n - 1)
print(f"\n=== 이산 역변환: 포아송({lam_pois}) ===")
print(f" E[Y] = {mean_pois:.4f} (이론: {lam_pois})")
print(f" Var[Y] = {var_pois:.4f} (이론: {lam_pois})")11.2 Step 2: scipy 구현 (수락-기각 + Metropolis + 검증)
import numpy as np
from scipy import stats
from scipy.special import gamma as gamma_func
np.random.seed(42)
n_sim = 50000
# --- 1. 수락-기각: Beta(2.7, 6.3) 생성 ---
print("=== 수락-기각: Beta(2.7, 6.3) ===")
a, b = 2.7, 6.3
target = stats.beta(a, b)
# 후보: Uniform(0,1) → M = max f_Y(y)
M_uniform = target.pdf(target.ppf(0.5)) # 근사
# 정확한 M = max of beta pdf
from scipy.optimize import minimize_scalar
res = minimize_scalar(lambda x: -target.pdf(x), bounds=(0, 1), method='bounded')
M_uniform = -res.fun
accepted = []
total_trials = 0
while len(accepted) < n_sim:
v = np.random.uniform(0, 1)
u = np.random.uniform(0, 1)
total_trials += 1
if u < target.pdf(v) / M_uniform:
accepted.append(v)
accepted = np.array(accepted)
accept_rate = n_sim / total_trials
print(f" M = {M_uniform:.3f}, 수락률 = {accept_rate:.3f} (이론: {1/M_uniform:.3f})")
ks_stat, p_val = stats.kstest(accepted, 'beta', args=(a, b))
print(f" KS test: stat={ks_stat:.4f}, p={p_val:.4f}")
print(f" E[Y] = {accepted.mean():.4f} (이론: {a/(a+b):.4f})")
# --- 2. Metropolis: N(0,1) from Cauchy 후보 ---
print(f"\n=== Metropolis: N(0,1), 후보=Cauchy ===")
n_chain = 100000
burn_in = 5000
chain = np.zeros(n_chain)
chain[0] = 0.0 # 초기값
for i in range(1, n_chain):
# 후보: 현재 위치 + Cauchy 증분
v = np.random.standard_cauchy()
rho = min(
stats.norm.pdf(v) / stats.cauchy.pdf(v) *
stats.cauchy.pdf(chain[i-1]) / stats.norm.pdf(chain[i-1]),
1.0
)
if np.random.uniform() <= rho:
chain[i] = v
else:
chain[i] = chain[i-1]
# burn-in 제거
samples = chain[burn_in:]
print(f" E[Z] = {samples.mean():.4f} (이론: 0)")
print(f" Var[Z] = {samples.var():.4f} (이론: 1)")
ks_stat, p_val = stats.kstest(samples, 'norm')
print(f" KS test: stat={ks_stat:.4f}, p={p_val:.4f}")
# --- 3. Box-Muller 정규성 검증 ---
print(f"\n=== Box-Muller 정규성 검증 ===")
u1 = np.random.uniform(0, 1, n_sim)
u2 = np.random.uniform(0, 1, n_sim)
r = np.sqrt(-2 * np.log(u2))
theta = 2 * np.pi * u1
x_bm = r * np.cos(theta)
y_bm = r * np.sin(theta)
# X와 Y의 독립성 확인
corr = np.corrcoef(x_bm, y_bm)[0, 1]
ks_x, p_x = stats.kstest(x_bm, 'norm')
ks_y, p_y = stats.kstest(y_bm, 'norm')
print(f" X: KS={ks_x:.4f}, p={p_x:.4f}")
print(f" Y: KS={ks_y:.4f}, p={p_y:.4f}")
print(f" Corr(X,Y) = {corr:.4f} (이론: 0, 독립)")
# --- 4. 역변환법 vs scipy 비교: 지수분포 ---
print(f"\n=== 역변환법 vs scipy: 지수(3) ===")
lam = 3.0
inv_samples = -lam * np.log(np.random.uniform(0, 1, n_sim))
scipy_samples = np.random.exponential(lam, n_sim)
print(f" 역변환: E={inv_samples.mean():.4f}, Var={inv_samples.var():.4f}")
print(f" scipy: E={scipy_samples.mean():.4f}, Var={scipy_samples.var():.4f}")
print(f" 이론: E={lam:.4f}, Var={lam**2:.4f}")이 코드는 네 가지를 검증한다: (1) 수락-기각 알고리즘으로 Beta(2.7, 6.3) 생성 + KS 검정, (2) Metropolis 알고리즘으로 정규분포 생성 (Cauchy 후보) + burn-in 처리, (3) Box-Muller로 생성한 두 정규 변수의 정규성과 독립성, (4) 역변환법과 scipy의 지수분포 생성 결과 비교.
12 응용 분야
| 분야 | 난수 생성의 역할 | 구체적 예시 |
|---|---|---|
| Monte Carlo 적분 | 고차원 적분의 시뮬레이션 근사 | 옵션 가격 산정 (블랙-숄즈 시뮬레이션) |
| 부트스트랩 | 경험적 분포에서 재표집 | 추정량의 표준오차와 신뢰구간 |
| 베이지안 추론 | 사후분포에서 표본 추출 (MCMC) | Gibbs Sampler, Hamiltonian MC |
| 가설 검정 | 순열 검정, 시뮬레이션 기반 p-값 | 포아송 분산 검정 (포스트 내 예시) |
| 최적화 | 확률적 탐색 (Simulated Annealing) | 조합 최적화 문제 |
| 기계학습 | SGD, 드롭아웃, 데이터 증강 | 미니배치 표본 추출 |
| 임상시험 | 시뮬레이션 기반 표본 크기 결정 | 복잡한 설계의 검정력 시뮬레이션 |
13 관련 주제
선행 지식
- 확률 표본의 기본 개념 (Basic Concepts of Random Samples)
- 확률 표본의 성질 개요 (Properties of a Random Sample: Overview)
- 수렴 개념 (Convergence Concepts) — WLLN/SLLN이 시뮬레이션을 정당화
후속 주제
관련 개념
- 순서통계량 (Order Statistics) — 균일분포 순서통계량과 베타분포
- 정규 모집단에서의 표본분포 (Sampling from the Normal Distribution) — 카이제곱, t, F 분포 생성