실습 · 연습문제 2개

HMC는 위치 $\theta$와 운동량 $r$을 Hamiltonian $H(\theta, r) = -\log \tilde p(\theta) + \tfrac{1}{2} r^\top r$을 따라 움직인다. leapfrog 한 걸음(보폭 $\varepsilon$)은

$$r \leftarrow r + \tfrac{\varepsilon}{2}\nabla\log\tilde p(\theta), \qquad \theta \leftarrow \theta + \varepsilon\, r, \qquad r \leftarrow r + \tfrac{\varepsilon}{2}\nabla\log\tilde p(\theta)$$

이다. 두 함수를 완성한다.

1. leapfrog(theta, r, grad_log_p, eps, L): leapfrog를 $L$걸음 진행한 (theta, r)을 반환한다. 연속한 두 걸음 사이의 반 걸음 두 개는 한 걸음으로 합쳐도 된다(결과는 같다). 입력 배열을 제자리에서 바꾸지 않는다.
2. hamiltonian(theta, r, log_p): $H(\theta, r) = -\log\tilde p(\theta) + \tfrac12 r^\top r$

테스트는 조화 진동자 $\log \tilde p(\theta) = -\tfrac12\theta^2$에서 다음을 확인한다.

• 에너지 오차가 작고, 보폭을 절반으로 줄이면 대략 1/4로 준다($O(\varepsilon^2)$).
• 끝점에서 운동량의 부호를 뒤집어 다시 적분하면 출발점으로 돌아온다(가역성). 이 성질 덕분에 HMC의 수용 확률이 $\min(1, e^{-\Delta H})$로 단순해진다.
• 보폭이 안정 한계($\varepsilon\omega < 2$, 여기서 $\omega = 1$)를 넘으면 에너지가 폭주한다.

체인 $M$개를 각각 앞뒤 절반으로 나누면 길이 $n$인 체인 $2M$개가 된다. 나눈 체인 $j$의 평균을 $\bar\theta_j$, 표본분산(분모 $n-1$)을 $s_j^2$라 하면

$$W = \frac{1}{2M}\sum_j s_j^2, \qquad B = \frac{n}{2M - 1}\sum_j \big(\bar\theta_j - \bar\theta\big)^2, \qquad \widehat{R} = \sqrt{\frac{\frac{n-1}{n}W + \frac{1}{n}B}{W}}$$

이다. $\bar\theta$는 나눈 체인 평균들의 평균이다.

split_rhat(chains)를 완성한다.

• chains는 shape (M, N)인 배열(또는 중첩 리스트)이다. np.asarray로 바꾼다.
• 2차원이 아니거나 $N < 4$이면 ValueError를 낸다.
• $N$이 홀수이면 마지막 표본을 버리고 반으로 나눈다.

테스트는 세 가지를 확인한다. 같은 분포에서 독립으로 뽑은 체인은 $\widehat{R} \approx 1$이다. 서로 다른 곳에 머문 체인은 $\widehat{R}$가 크다. 체인이 하나뿐이어도 앞뒤 절반이 다르면(추세가 있으면) split 덕분에 $\widehat{R}$가 커진다.