마르코프 연쇄 몬테카를로

Markov chain Monte Carlo; MCMC; 마르코프 연쇄 몬테카를로

원하는 확률 분포를 정상 분포로 갖는 마르코프 연쇄를 만들어 오래 돌리고, 그 경로를 표본으로 삼아 기댓값과 분포를 근사하는 표본 추출 방법의 총칭

마르코프 연쇄 몬테카를로(Markov chain Monte Carlo, MCMC)는 직접 표본을 뽑기 어려운 확률 분포 $\pi$에서 표본을 얻는 방법이다. $\pi$를 정상 분포로 갖는 마르코프 연쇄를 설계해 충분히 오래 돌리면 연쇄가 머무는 위치의 분포가 $\pi$에 가까워지므로, 그 경로를 (서로 상관된) 표본으로 쓴다. $\pi$의 정규화 상수를 몰라도 된다는 점 때문에 베이즈 통계의 사후 분포 계산, 통계 물리, 확률적 프로그래밍의 중심 도구가 되었다.

마르코프 연쇄의 정상 분포·세부 균형·에르고딕 정리와 메트로폴리스-헤이스팅스의 수락 확률 식, 1953년과 1970년의 출발점은 마르코프 연쇄 문서의 해당 절에 정리되어 있다. 이 문서는 그 위에서 알고리즘별 작동 방식, 튜닝, 수렴 진단, 구현을 다룬다.

왜 필요한가

베이즈 추론에서 매개변수 $\theta$의 사후 분포는

$p(\theta \mid D)=\frac{p(D\mid \theta )p(\theta )}{\int p(D\mid {\theta }^{\prime })p({\theta }^{\prime })d{\theta }^{\prime }}$

인데, 분자는 쉽게 계산되지만 분모의 적분(증거)은 차원이 조금만 높아도 계산할 수 없다. 우리가 원하는 것은 대개 사후 평균, 신용 구간, 예측 분포 같은 기댓값 ${𝔼}_{\pi }[f(\theta )]$이다. 독립 표본 ${\theta }^{(1)},\dots ,{\theta }^{(n)}\sim \pi$가 있으면 평균 $\frac{1}{n}\sum _{t}f({\theta }^{(t)})$로 근사하면 되지만(몬테카를로 방법), 그 독립 표본을 뽑는 것부터 어렵다. 기각 샘플링이나 중요도 샘플링은 제안 분포가 목표와 조금만 달라도 차원이 커질수록 효율이 기하급수로 떨어진다.

MCMC는 독립 표본을 포기하고, 현재 위치 근처에서 조금씩 움직이며 $\pi$가 큰 영역을 오래 돌아다니는 연쇄를 만든다. 비율 $\pi ({\theta }^{\prime })/\pi (\theta )$만 쓰므로 분모의 적분이 약분되어 사라진다.

원리

• 불변성: 전이 커널 $T(x^{\prime }\mid x)$가 $\int \pi (x)T(x^{\prime }\mid x)dx=\pi (x^{\prime })$를 만족해야 한다. 세부 균형 $\pi (x)T(x^{\prime }\mid x)=\pi (x^{\prime })T(x\mid x^{\prime })$는 이를 보장하는 충분조건이고, 대부분의 알고리즘이 이 조건으로 설계된다.
• 에르고딕성: 연쇄가 기약이고 비주기적이면 초기값과 무관하게 분포가 $\pi$로 수렴하고, 시간 평균이 기댓값으로 수렴한다(에르고딕 정리).
• 몬테카를로 오차: 표본이 자기상관되어 있으므로 표본 $n$개의 정보량은 독립 표본보다 적다. 지연 $k$의 자기상관을 ${\rho }_k$라 하면 유효 표본 크기(effective sample size, ESS)는

$\mathrm{E}\mathrm{S}\mathrm{S}=\frac{n}{1+2\sum _{k=1}^{\mathrm{\infty }}{\rho }_k}$

이고, 추정값의 몬테카를로 표준 오차(MCSE)는 대략 ${\mathrm{s}\mathrm{d}}_{\pi }(f)/\sqrt{\mathrm{E}\mathrm{S}\mathrm{S}}$다. 알고리즘의 효율은 결국 "계산 시간당 ESS"로 비교한다.

주요 알고리즘

메트로폴리스-헤이스팅스

현재 상태 $x$에서 제안 분포 $q(x^{\prime }\mid x)$로 후보를 뽑고 확률 $\alpha =\min (1,\frac{\pi (x^{\prime })q(x\mid x^{\prime })}{\pi (x)q(x^{\prime }\mid x)})$로 받아들인다. 거절하면 현재 상태를 한 번 더 표본으로 기록한다(버리면 안 된다).

• 랜덤 워크 메트로폴리스: $x^{\prime }=x+ϵ$, $ϵ\sim 𝒩(0,s^{2}I)$. 대칭이라 $q$ 비율이 약분된다. 보폭 $s$가 너무 작으면 거의 다 수락되지만 조금씩만 움직이고, 너무 크면 대부분 거절되어 제자리에 머문다. 목표가 독립 성분의 곱이고 차원이 클 때 최적 수락률은 약 0.234라는 결과가 있어 튜닝 기준으로 쓴다. 이 경우 $d$차원에서 독립 표본 하나를 얻는 데 대략 $O(d)$번의 단계가 든다.
• 독립 제안: $q(x^{\prime }\mid x)=q(x^{\prime })$. 제안 분포가 목표와 비슷하면 효율적이지만, 꼬리가 목표보다 얇으면 매우 나빠진다.
• 실무에서는 예열(warm-up) 구간에서 수락률을 보며 보폭과 제안 공분산을 조정하고, 그 구간 표본은 버린다.

깁스 샘플링

변수를 $x=(x_1,\dots ,x_d)$로 나눠, 하나씩 나머지가 주어졌을 때의 완전 조건부 분포 $\pi (x_j\mid x_{-j})$에서 다시 뽑는다. 이는 그 조건부 분포를 제안 분포로 쓰는 메트로폴리스-헤이스팅스이며 수락 확률이 항상 1이다. 1984년 제먼 형제(Stuart Geman, Donald Geman)가 영상 복원에서 깁스 분포의 표본을 뽑는 방법으로 이름 붙였고, 1990년 겔펀드(Alan Gelfand)와 스미스(Adrian Smith)가 베이즈 주변 분포 계산에 쓸 수 있음을 보이면서 통계학 전반에 퍼졌다.

• 켤레 사전 분포를 쓰는 모형에서는 완전 조건부 분포가 알려진 분포족이 되어 구현이 쉽다. 베이지안 네트워크에서는 마르코프 블랭킷만 보고 뽑으면 된다. BUGS, JAGS 같은 도구가 이 방식을 자동화했다.
• 블록 깁스: 강하게 얽힌 변수들을 묶어 한꺼번에 뽑는다.
• 붕괴(collapsed) 깁스: 일부 변수를 해석적으로 적분해 없애고 남은 변수만 뽑는다. 잠재 디리클레 할당의 표준 추론법이다.
• 약점: 변수 사이 상관이 강하면 좌표축 방향으로만 조금씩 움직여 느려진다. 상관계수 $\rho$인 2차원 정규 분포에서 한 좌표의 연쇄는 지연 1 자기상관이 ${\rho }^2$인 자기회귀 과정이 되어, $\rho =0.9$면 ESS가 표본 수의 약 $(1-{\rho }^2)/(1+{\rho }^2)\approx 10\mathrm{\%}$로 떨어진다.

해밀토니언 몬테카를로

해밀토니언 몬테카를로(Hamiltonian Monte Carlo, HMC)는 목표 분포의 그래디언트를 써서 멀리 떨어진 곳까지 한 번에 움직이면서도 수락률을 높게 유지한다. 1987년 격자 양자 색역학 계산을 위해 하이브리드 몬테카를로라는 이름으로 나왔고, 래드퍼드 닐(Radford Neal)이 통계와 신경망 학습에 소개했다.

위치 $q$(매개변수)에 운동량 $p$를 덧붙이고, 위치 에너지 $U(q)=-\log \pi (q)$, 운동 에너지 $K(p)=\frac{1}{2}{p}^{\mathrm{\top }}{M}^{-1}p$로 해밀토니언 $H(q,p)=U(q)+K(p)$를 정의한다. 결합 분포 $\propto e^{-H}$의 $q$ 주변 분포가 $\pi$다. 한 번의 이동은 다음과 같다.

• 운동량을 새로 뽑는다: $p\sim 𝒩(0,M)$
• 해밀턴 방정식을 립프로그(leapfrog) 적분으로 보폭 $\epsilon$, $L$단계 푼다: $p\leftarrow p-\frac{\epsilon }{2}\mathrm{\nabla }U(q)$, $q\leftarrow q+\epsilon M^{-1}p$, $p\leftarrow p-\frac{\epsilon }{2}\mathrm{\nabla }U(q)$를 반복
• 수치 적분 오차 때문에 $H$가 조금 변하므로 확률 $\min (1,e^{H(q,p)-H(q^{\prime },p^{\prime })})$로 수락한다. 립프로그는 부피를 보존하고 시간 가역적이라 이 보정만으로 정확한 MCMC가 된다.

정확한 해밀턴 역학은 $H$를 보존하므로 오차가 작으면 수락률이 1에 가깝고, 등고선을 따라 멀리 이동하므로 자기상관이 작다. 차원이 커질 때의 비용 증가도 랜덤 워크보다 완만하다. 대신 $\mathrm{\nabla }\log \pi$가 필요하므로 연속 매개변수에만 쓰고(이산 변수는 적분하거나 다른 방법과 섞는다), 성능이 $\epsilon$과 $L$에 매우 민감하다. $\epsilon$이 크면 에너지 오차가 커져 거절이 늘고, 곡률이 급한 곳에서 궤적이 폭주하는 발산(divergence)이 생긴다. $L$이 작으면 랜덤 워크처럼 되고, 크면 궤적이 되돌아와 계산을 낭비한다.

NUTS

NUTS(No-U-Turn Sampler)는 HMC의 $L$을 자동으로 정한다. 궤적을 앞뒤로 두 배씩 늘려 가다가 궤적의 양 끝이 서로 가까워지기 시작하면(U턴) 멈추고, 만들어진 점들 가운데 세부 균형을 지키는 방식으로 다음 상태를 고른다. 보폭 $\epsilon$은 예열 구간에서 목표 평균 수락률에 맞추도록 이중 평균(dual averaging)으로 조정한다. 호프만(Matthew Hoffman)과 겔먼(Andrew Gelman)은 NUTS가 잘 튜닝한 HMC와 같거나 더 효율적이면서 사람이 손으로 튜닝할 필요가 없다고 보고했다. Stan과 PyMC의 기본 샘플러다. PyMC 5의 NUTS는 목표 수락률 target_accept=0.8, 최대 트리 깊이 max_treedepth=10(한 번에 최대 $2^{10}$ 립프로그 단계)이 기본값이다. 발산이 보고되면 target_accept를 0.9~0.99로 올려 보폭을 줄이거나 모형을 다시 매개변수화(중심화 제거 등)한다.

그 밖의 방법

방법	핵심	쓰임
슬라이스 샘플링	밀도 그래프 아래 영역에서 균등하게 뽑는 보조 변수 방법. 보폭 튜닝이 덜 민감하다	1차원 조건부 분포, 깁스 안의 한 단계
랑주뱅 방법(MALA)	그래디언트 방향으로 한 걸음 + 잡음, 메트로폴리스 보정	HMC의 $L=1$ 특수형
확률적 그래디언트 랑주뱅 동역학(SGLD)	미니배치 그래디언트와 잡음으로 갱신, 보폭을 줄여 가며 메트로폴리스 보정 생략	대규모 데이터의 베이지안 심층 신경망
병렬 템퍼링	온도가 다른 여러 연쇄를 돌리며 상태를 교환	봉우리가 여럿인 분포
가역 점프 MCMC	차원이 다른 모형 사이를 오가는 제안	모형 선택, 성분 수가 정해지지 않은 혼합 모형

알고리즘 비교

알고리즘	필요한 것	장점	약점	언제 쓰나
랜덤 워크 MH	$\log \pi$ 값	구현이 가장 쉽다, 그래디언트 불필요	고차원·강한 상관에서 매우 느림, 보폭 튜닝 필요	저차원, 그래디언트를 못 구하는 모형
깁스	완전 조건부 분포에서 뽑는 방법	튜닝 불필요, 수락률 1, 이산 변수 가능	상관 강하면 느림, 조건부 분포 유도 필요	켤레 모형, 이산 잠재 변수(LDA 등)
HMC	$\mathrm{\nabla }\log \pi$	고차원 연속 공간에서 효율적	$\epsilon ,L$ 튜닝, 이산 변수 불가	매끄러운 연속 사후 분포
NUTS	$\mathrm{\nabla }\log \pi$	HMC 효율 + 자동 튜닝, 발산으로 문제 진단	단계당 계산이 많고, 곡률이 급변하는 분포(깔때기 모양)에서 고전	확률적 프로그래밍의 기본값

수렴 진단

MCMC 결과가 맞는지 확실히 증명할 방법은 없고, 문제가 있다는 신호를 찾을 수 있을 뿐이다.

• 예열과 번인: 초기값의 영향이 남은 앞부분을 버린다. 적응형 샘플러는 이 구간에 보폭 등을 튜닝한다(PyMC의 tune, 기본 1000).
• 여러 연쇄: 서로 흩어진 초기값에서 연쇄를 4개 이상 돌린다. 한 연쇄만 보면 한 봉우리에 갇힌 것을 알아채기 어렵다.
• R-hat: 연쇄 사이 분산과 연쇄 안 분산을 비교한다. $W$를 연쇄 안 분산의 평균, $B/n$을 연쇄 평균들의 분산이라 하면 $\hat{R}=\sqrt{\frac{\frac{n-1}{n}W+\frac{B}{n}}{W}}$이고, 모든 연쇄가 같은 분포를 훑었다면 1에 가깝다. 겔먼과 루빈(1992)이 제안했고, 베흐타리(Aki Vehtari) 등은 연쇄를 반으로 나누고 순위로 정규화한 개선판을 내면서 $\hat{R}<1.01$과 순위 정규화 ESS 400 이상을 기준으로 권했다. 겔먼과 루빈이 처음 권한 것보다 훨씬 엄격한 기준이다.
• ESS: 분포 중심부 추정용(bulk)과 꼬리 분위수용(tail)을 따로 본다.
• 추적 그림(trace plot): 연쇄들이 같은 영역에서 고르게 섞이는 "털뭉치" 모양이어야 한다. 한쪽으로 흐르거나 연쇄마다 다른 높이에 머물면 문제다.
• HMC 전용 진단: 발산 전이 수(0이어야 한다), 최대 트리 깊이 도달 비율, 에너지 분포(E-BFMI).
• 솎아내기(thinning): 몇 개마다 하나만 남기는 것은 저장 공간을 줄일 뿐 추정 정확도는 오히려 떨어진다. 메모리가 문제일 때만 쓴다.

구현

NumPy로 세 알고리즘 비교

상관계수 0.9인 2차원 정규 분포를 목표로 랜덤 워크 MH, 깁스, HMC를 각각 네 연쇄씩 돌려 수락률, R-hat, ESS를 비교한다(NumPy 2.1에서 실행).

import numpy as np

rho = 0.9
S_inv = np.linalg.inv(np.array(((1.0, rho), (rho, 1.0))))
log_p = lambda x: -0.5 * x @ S_inv @ x          # 정규화 상수 없는 로그 밀도
grad_log_p = lambda x: -S_inv @ x

def metropolis(n, step, rng, x0):
    x, out, acc = x0.copy(), np.empty((n, 2)), 0
    for t in range(n):
        prop = x + step * rng.normal(size=2)        # 대칭 제안: q 비율이 약분된다
        if np.log(rng.uniform()) < log_p(prop) - log_p(x):
            x, acc = prop, acc + 1
        out[t] = x                                  # 거절돼도 현재 상태를 기록
    return out, acc / n

def gibbs(n, rng, x0):
    x, out = x0.copy(), np.empty((n, 2))
    sd = np.sqrt(1 - rho**2)
    for t in range(n):
        x[0] = rng.normal(rho * x[1], sd)           # p(x1 | x2)
        x[1] = rng.normal(rho * x[0], sd)           # p(x2 | x1)
        out[t] = x
    return out

def hmc(n, eps, L, rng, x0):
    x, out, acc = x0.copy(), np.empty((n, 2)), 0
    for t in range(n):
        p0 = rng.normal(size=2)                     # 운동량 새로 뽑기
        q, p = x.copy(), p0 + 0.5 * eps * grad_log_p(x)
        for i in range(L):                          # 립프로그 적분
            q = q + eps * p
            if i < L - 1:
                p = p + eps * grad_log_p(q)
        p = p + 0.5 * eps * grad_log_p(q)
        H0 = -log_p(x) + 0.5 * p0 @ p0
        H1 = -log_p(q) + 0.5 * p @ p
        if np.log(rng.uniform()) < H0 - H1:         # 이산화 오차를 MH로 보정
            x, acc = q, acc + 1
        out[t] = x
    return out, acc / n

def rhat(chains):
    """chains: (m, n). 겔만-루빈 R-hat(분할 없는 기본형)."""
    m, n = chains.shape
    W = chains.var(axis=1, ddof=1).mean()
    B = n * chains.mean(axis=1).var(ddof=1)
    return np.sqrt(((n - 1) / n * W + B / n) / W)

def ess(x, max_lag=1000):
    """자기상관이 처음 음수가 되기 전까지 더하는 단순 추정."""
    x = x - x.mean()
    acf = np.correlate(x, x, "full")[len(x) - 1:] / (x @ x)
    s = 0.0
    for k in range(1, max_lag):
        if acf[k] < 0:
            break
        s += acf[k]
    return len(x) / (1 + 2 * s)

rng = np.random.default_rng(0)
starts = [np.array(s) for s in ([-4.0, 4.0], [4.0, -4.0], [4.0, 4.0], [-4.0, -4.0])]
burn, n = 1000, 5000
for name, run in [
    ("랜덤 워크 MH", lambda x0: metropolis(burn + n, 0.5, rng, x0)),
    ("깁스", lambda x0: (gibbs(burn + n, rng, x0), 1.0)),
    ("HMC", lambda x0: hmc(burn + n, 0.25, 8, rng, x0)),
]:
    res = [run(x0) for x0 in starts]
    ch = np.array([r[0][burn:, 0] for r in res])   # 첫 좌표만 진단
    print(f"{name:10s} 수락률={np.mean([r[1] for r in res]):.2f} "
          f"R-hat={rhat(ch):.3f} ESS(연쇄 1개)={ess(ch[0]):.0f}/{n} "
          f"상관 추정={np.corrcoef(np.concatenate([r[0][burn:] for r in res]).T)[0, 1]:.3f}")

랜덤 워크 MH   수락률=0.54 R-hat=1.010 ESS(연쇄 1개)=114/5000 상관 추정=0.897
깁스         수락률=1.00 R-hat=1.001 ESS(연쇄 1개)=574/5000 상관 추정=0.900
HMC        수락률=0.99 R-hat=1.000 ESS(연쇄 1개)=2302/5000 상관 추정=0.900

세 방법 모두 상관 0.9를 복원하지만 같은 표본 수에서 ESS는 20배 가까이 차이 난다. 깁스의 ESS(574)는 위에서 말한 이론값(약 525)과 비슷하다. 다만 HMC는 그래디언트를 단계당 8번 계산했으므로 계산량까지 따지면 차이가 줄어든다. 같은 코드에서 HMC를 $\epsilon =0.2,L=10$으로 바꾸면 ESS가 876으로 떨어지는데, 이런 민감성이 NUTS가 나온 이유다. 실무에서는 이 ESS 함수 대신 ArviZ의 az.ess, az.rhat(순위 정규화 분할 R-hat)을 쓴다.

PyMC

확률적 프로그래밍 라이브러리에서는 모형만 선언하면 그래디언트 계산, NUTS, 보폭 적응, 진단을 자동으로 해 준다. 아래는 베이즈 선형 회귀다(PyMC 5.25, ArviZ 0.23에서 실행).

import numpy as np
import pymc as pm
import arviz as az

def main():
    rng = np.random.default_rng(42)
    x = rng.normal(size=100)
    y = 1.0 + 2.5 * x + rng.normal(scale=0.8, size=100)   # 참값: a=1, b=2.5, sigma=0.8

    with pm.Model() as model:
        a = pm.Normal("a", mu=0, sigma=10)
        b = pm.Normal("b", mu=0, sigma=10)
        sigma = pm.HalfNormal("sigma", sigma=5)
        pm.Normal("y", mu=a + b * x, sigma=sigma, observed=y)
        idata = pm.sample(draws=1000, tune=1000, chains=4, cores=1,
                          target_accept=0.9, random_seed=0, progressbar=False)

    cols = ["mean", "sd", "hdi_3%", "hdi_97%", "ess_bulk", "r_hat"]
    summary = az.summary(idata, var_names=["a", "b", "sigma"], kind="all", round_to=2)
    print(summary[cols])
    print("발산 전이 수:", int(idata.sample_stats["diverging"].sum()))

if __name__ == "__main__":
    main()

       mean    sd  hdi_3%  hdi_97%  ess_bulk  r_hat
a      1.00  0.08    0.84     1.14   3763.00    1.0
b      2.60  0.11    2.39     2.79   4799.19    1.0
sigma  0.79  0.06    0.69     0.90   3450.48    1.0
발산 전이 수: 0

R-hat이 1.0이고 ESS가 4000개 표본 가운데 3400 이상이며 발산이 없으므로 진단상 문제가 없다. hdi_3%~hdi_97%는 94% 최고 밀도 구간이다. chains를 비워 두면 cores와 2 가운데 큰 값이 된다.

인공지능에서의 쓰임

• 확률적 프로그래밍과 베이즈 모형: Stan, PyMC, NumPyro 같은 도구로 계층 모형, 시계열, A/B 시험 분석의 사후 분포를 구한다.
• 토픽 모델: 잠재 디리클레 할당의 붕괴 깁스 샘플링.
• 베이즈 신경망: 닐은 HMC로 신경망 가중치의 사후 분포에서 표본을 뽑는 연구를 했다. 규모가 큰 모형에는 SGLD 같은 확률적 그래디언트 MCMC나 변분 추론을 쓴다(베이지안 심층 신경망).
• 에너지 기반 모델과 볼츠만 머신: 제한 볼츠만 머신의 대조 발산(contrastive divergence) 학습은 깁스 샘플링을 몇 단계만 돌려 그래디언트를 근사한다. 에너지 기반 모델도 랑주뱅 MCMC로 표본을 만든다.
• 점수 기반 생성 모델: 송양(Yang Song)과 에르몬(Stefano Ermon)은 데이터 분포의 로그 밀도 그래디언트(점수)를 신경망으로 학습하고 잡음 수준을 줄여 가는 랑주뱅 동역학으로 이미지를 생성했다. 이는 확산 모델과 같은 계열로 이어졌다.

변분 추론과 비교

구분	MCMC	변분 추론
방식	표본 추출	다루기 쉬운 분포족 안에서 최적화
정확성	오래 돌리면 점근적으로 정확	분포족의 한계만큼 편향이 남음(분산 과소 추정 경향)
속도·규모	느림, 대규모 데이터에 불리	빠름, 미니배치·GPU와 잘 맞음
진단	R-hat, ESS 등 성숙한 진단	근사 품질 판단이 어려움
쓰임	중소 규모 정밀 추론, 기준 답 계산	대규모 모형, 변분 오토인코더, 온라인 학습

블레이 등의 리뷰는 정밀함이 중요하고 계산 여유가 있으면 MCMC, 큰 데이터나 빠른 모형 탐색이 필요하면 변분 추론이 맞는다고 정리했다.

한계

• 수렴했다는 보장을 얻을 수 없다. 진단은 실패를 드러낼 뿐이다.
• 봉우리가 여럿이고 서로 멀리 떨어진 분포에서는 연쇄가 한 봉우리에 갇혀, 모든 진단이 정상이어도 답이 틀릴 수 있다.
• 매 단계 전체 데이터의 우도를 계산하므로 데이터가 크면 느리다.
• 깔때기 모양(계층 모형의 분산 매개변수 근처)처럼 곡률이 크게 변하는 분포는 HMC·NUTS도 발산을 일으킨다. 비중심화 매개변수화 같은 모형 재작성이 필요하다.

같이 보기

• 마르코프 연쇄
• 몬테카를로 시뮬레이션
• 몬테 카를로 알고리즘
• 베이즈 정리
• 베이지안 네트워크
• 깁스 분포
• 잠재 디리클레 할당
• EM 알고리즘
• 변분 추론
• 베이지안 심층 신경망
• 에너지 기반 모델
• 확산 모델

참고 문헌

• Steve Brooks, Andrew Gelman, Galin L. Jones, Xiao-Li Meng 엮음, 《Handbook of Markov Chain Monte Carlo》, Chapman & Hall/CRC, 2011.
• Andrew Gelman 외, 《Bayesian Data Analysis》 3판, CRC Press, 2013. 저자 공개본
• Michael Betancourt, 〈A Conceptual Introduction to Hamiltonian Monte Carlo〉, 2017. arXiv:1701.02434