Monte Carlo Methods몬테카를로 방법

기본 아이디어:
1. 무작위 샘플 N개 추출
2. 함수값 계산
3. 평균으로 기댓값 추정

오차: O(1/√N)  ← 차원에 독립적 (차원의 저주 극복)

import numpy as np

def estimate_pi(n_samples: int) -> float:
    """몬테카를로로 π 추정"""
    # 단위 정사각형 [-1,1]x[-1,1]에 무작위 점 생성
    x = np.random.uniform(-1, 1, n_samples)
    y = np.random.uniform(-1, 1, n_samples)

    # 단위 원 내부 점 개수
    inside = np.sum(x**2 + y**2 <= 1)

    # π ≈ 4 × (원 내부 점수 / 전체 점수)
    return 4 * inside / n_samples

for n in [1_000, 10_000, 100_000, 1_000_000]:
    pi_est = estimate_pi(n)
    print(f"N={n:>10,}: π ≈ {pi_est:.6f} (오차: {abs(pi_est - np.pi):.6f})")

def monte_carlo_integrate(f, a, b, n=100_000):
    """1차원 몬테카를로 적분"""
    x = np.random.uniform(a, b, n)
    return (b - a) * np.mean(f(x))

# sin(x) 적분 [0, π] = 2
result = monte_carlo_integrate(np.sin, 0, np.pi)
print(f"∫sin(x)dx from 0 to π ≈ {result:.4f} (정답: 2.0)")

# 다차원 적분 (5차원 구의 부피)
def is_in_sphere(x):
    return np.sum(x**2, axis=1) <= 1

def hypersphere_volume(d, n=1_000_000):
    points = np.random.uniform(-1, 1, (n, d))
    fraction = np.mean(is_in_sphere(points))
    return (2**d) * fraction  # 초정육면체 부피 × 비율

for d in range(2, 7):
    vol = hypersphere_volume(d)
    print(f"{d}차원 구의 부피 ≈ {vol:.4f}")

import numpy as np

def metropolis_hastings(target_pdf, n_samples, x_init=0.0, step_size=0.5):
    """Metropolis-Hastings 알고리즘"""
    samples = [x_init]
    x_current = x_init

    for _ in range(n_samples):
        # 제안 분포 (가우시안)
        x_proposed = x_current + np.random.normal(0, step_size)

        # 수용 확률
        accept_ratio = target_pdf(x_proposed) / target_pdf(x_current)
        if np.random.random() < min(1, accept_ratio):
            x_current = x_proposed
        samples.append(x_current)

    return np.array(samples)

# 정규분포 샘플링 (검증)
target = lambda x: np.exp(-x**2 / 2)
samples = metropolis_hastings(target, n_samples=10_000)
print(f"MCMC 평균: {samples.mean():.3f}, 표준편차: {samples.std():.3f}")