4.2 Markov Chain Monte Carlo

\(f\) 가 측정은 되는데 샘플화가 안되면, MC를 통해 유사한 샘플을 만들어낼 수 있었다. 이를 넘어서 MCMC는 $ f $ 의 모사함수에서 샘플링하는 게 가능하지만, 이 이상으로 이는 임의의 함수 \(p\)에 대해 \(E[p(X)]\)가 신뢰도 높게 측정되는 경우에만 샘플링 가능한 별개의 방법론으로 보는 게 정확하다.

MC MCMC numerical integration approach
iterative nature
can be customized to very diverse &
difficult problem
무관하며 implementation이 complex하지도 않음 문제가 고차원이면 수렴이 느려짐

시퀀스 \(\{\textbf X^{(t)}\}\)는 MC, \(t = 0, 1, 2, ….\). \(\textbf X^{(t)} = (X_1^{t} , \cdots, X_p^{(t)})\) 와 state space 는 양쪽 모두 연속이거나 discrete.

For the types of Markov chains, \(\{ \textbf X^{(t)} \}\)의 분포는 체인의 limiting stationary distribution으로 수렴한다. 체인이 irreducible, aperiodic 할 때.

MCMC의 샘플링 전략은 irreducible, aperiodic MC를 만드는 것. stationary distribution이 목표분포 \(f\) 와 일치하는.

t가 충분히 크다면 이 체인으로부터의 \(\textbf X^{(t)}\)의 실현값은 근사적으로 마지널 분포 \(f\) 를 갖는다.

이런 MCMC의 특성은 베이지안 추론에 크게 도움이 되며 자주 쓰인다.




Markov Chain 자체는 어떤 상태에서 다른 상태로 넘어갈 때, 바로 전 단계의 상태에만 영향을 받는 (Markov Property) 확률 과정을 의미한다.

  • 보통 사람들은 전날 먹은 식사와 유사한 식사를 하지 않으려는 경향이 있다.
  • 가령, 어제 짜장면을 먹었다면 오늘은 면종류의 음식을 먹으려고 하지 않는다.






4.2.1 MH Algorithm

MCMC 중 가장 유명한 적용법은 MH 알고리즘. \(t=0\)에서 시작. 시작 distribution \(g\)에서 추출한, \(f(\textbf x^{(0)} )> 0\) 을 만족하는 \(\textbf x^{(0)}\)를 \(\textbf X^{(0)} = \textbf x^{(0)}\) 로 잡고 개시한다.

이때 제안분포 \(g\) 에서 후보 \(\textbf X^{( \ast )}\) 를 하나 만들고, MH ratio \(R (\textbf {x}^{(t)}, \textbf X^{\ast} )\) 는

\[ R (\textbf {x}^{(t)}, \textbf X^{\ast} ) = \dfrac {f(\textbf X^{( \ast )}) g(\textbf x^{( t )} | \textbf X^{( \ast )})} {f(\textbf x^{( t )}) g(\textbf X^{( \ast )} | \textbf x^{( t )}) } =\dfrac {\dfrac {f(\textbf X^{( \ast )})} {g(\textbf X^{( \ast )} | \textbf x^{( t )})} } {\dfrac {f(\textbf x^{( t )})} {g(\textbf x^{( t )} | \textbf X^{( \ast )})} } =\dfrac {\dfrac {f(\textbf X^{( t+1 )})} {g(\textbf X^{( t+1 )} | \textbf x^{( t )})} } {\dfrac {f(\textbf x^{( t )})} {g(\textbf x^{( t )} | \textbf X^{( t+1 )})} } \]

warning

여기서 단순 Metropolis 알고리즘은 단순히 $ {f(x_0)}$ > $1 $ 이기만 하면 새로운 샘플을 수용한다. 이인즉 \(g\)로 표준화해주는 것의 가장 주요한 요점은 둘의 시작 높이, 즉 쥐고 나온 수저가 다를 수 있으므로 이를 표준화해준다는 것이다. 아웃풋이 높은 \(x_i\)를 선택하는 것은 MLE 관점에 기반한다.

단 언제나 그렇듯 이렇게 샘플을 다쳐내면 오히려 음질의 결과가 나온다. 따라서 샘플의 풀을 넓히기 위해 탈락할 녀석들도 확률적으로 살려서 합류시킨다. 이게 고정된 기준점으로 샘플을 쳐내는 것이 아닌, \(\dfrac {f(x_1)} {f(x_0)}\) > \(u \sim U {(0,1)}\) 을 기준으로 삼아 샘플을 생존시키는 것이다. 이 조건을 실패하면 생산해두었던 샘플 \(\textbf X^{(t)}\) 는 버려지고 새로운 샘플을 \(t+1\)으로 설정해 재진행한다.

이후

\[ \textbf {X}^{(t+1)} = \left\{\begin{array}{@{}lr@{}} \textbf {X}^{\ast}, & \text{with probability } min \left\{ R \left( \textbf {x}^{(t+1)}, \textbf {X}^{\ast} \right), 1 \right\} \\ \textbf {x}^{(t)}, & o.w. \end{array}\right\} \]

이러한 MH 알고리즘에 의해 생성된 MC가 aperiodic & irreducible 이라면, 해당 체인은 정적분포로 수렴.

우리는 이러한 MH 체인에 의해 생성된 정적분포의 실현값들을 평균함으로써 rv의 함수의 기댓값을 구할 수 있다.

\(E \left[ h \left( \textbf {X} \right) \right] \approx \dfrac {1} {n} \sum_{i=1}^n {h \left( \textbf {x}^{(i)} \right)},\)

\(E { \left\{ h \left( \textbf {X} \right) - E \left[ h \left( \textbf {X} \right) \right] \right\} }^2\),

\(E \{ I_{h ( \textbf {X} \le q )} \}\)

시퀀스 \(\{ \textbf x^{(\inf)} \}\) 는 state space의 몇몇 포인트들의 multiple copies를 가질 수 있다는 것을 명심. 이는 \(\textbf {X}^{(t+1)}\)가 제안값 \(\textbf {x}^{(\ast )}\)가 아니라 \(\textbf {x}^{(t)}\)를 채택했을 때 발생.

  • Burn-in Period: 체인의 초기값에 대한 persistent dependence는 이의 성능을 심각하게 낮출 수 있다. 이는 샘플 평균을 계산할 때 체인의 초기 실현값 일부를 제하는 것으로 보정될 수 있다.

consistent 결과들을 관측하기 위해 MCMC를 여러 시작점에서 각각 돌려본다.

잘 골라진 제안분포 $ g $가 생산하는 후보값들은 stationary 분포의 서포트를 합리적인 반복 안에서 전부 커버하고, 내놓는 후보값들이 지나치게 여러번 accepted되거나 rejected 되지 않는다.

  • proposal 분포 $ g $ 가 지나치게 퍼져있으면, 후보값들은 자주 reject되고 체인은 타겟분포 $ f $ 의 space를 탐색하기 위해 많은 반복을 요구하게 된다.
  • proposal 분포 $ g $ 가 지나치게 모여있으면, 체인은 다회의 반복동안 타겟분포 $ f $의 한 작은 구역에 모여있게 된다. 따라서 타겟분포의 다른 영역은 적절하게 탐색되지 못한다.






4.2.1.1 Independent Chains

acceptance 여부 결정시에 ratio 자체는 MH ratio를 사용한다. 이 MC ratio에는 과거 실현값(\(x^{(t)}\))이 들어있다. 따라서 이는 MCMC 방법론에 해당한다. 하지만 새로운 value \(g(x')\)을 생산할 때, 이 자체는 \(g(x'\rvert x^{(t)})=g(x' \rvert \cdot )\)을 따르게, 즉 \(g(x'\rvert x^{(t)})=g(x' )\) 마냥 과거의 실현값에 dependent 하지 않게 새로운 값을 생산해내는 방법론. 즉 새로이 제시되는 candidate value가 과거의 실현값들과 independent 하므로 명칭이 저러한 것이다.

즉 MH ratio 자체는 과거의 샘플을 이용해서 MCMC 범주에 들어가나 샘플 자체는 과거의 샘플과 independent하게 생산.

MH 알고리즘의 제안분포는 고정된 덴시티 \(g\)에 대해서 $ g ( ^{} ^{(t)} )$ 따위로 생성. 이는 independent chain 이라고 불리며, 이에 사용되는 각 후보값들은 과거에서 독립적으로 추출되었다. MH ratio는

\[ R (\textbf {x}^{(t)}, \textbf X^{\ast} ) = \dfrac {f(\textbf X^{( \ast )}) g(\textbf x^{( t )} | \textbf X^{( \ast )})} {f(\textbf x^{( t )}) g(\textbf X^{( \ast )} | \textbf x^{( t )}) } = \dfrac {f(\textbf X^{( \ast )}) g(\textbf x^{( t )})} {f(\textbf x^{( t )}) g(\textbf X^{( \ast )}) } =\dfrac {\dfrac {f(\textbf X^{( \ast )})} {g(\textbf x^{( t )})} } {\dfrac {f(\textbf x^{( t )})} {g(\textbf X^{( \ast )})} } =\dfrac {\dfrac {f(\textbf X^{( \ast )})} {g(\textbf X^{( \ast )})} } {\dfrac {f(\textbf x^{( t )})} {g(\textbf x^{( t )})} } = \dfrac {w^{\ast}} {w^{(t)}} \tag{1} \]

  • importance ratio of \(X'\), importance ratio of \(X^{(t)}\). This can be seen as weight also.

이때 가장 우측의 등식은 importance ratio (\(r\))의 등식으로 재표현된 것이며, 이때 \(f\) 는 타겟분포, \(g\) 는 그것의 envelope로 본 것이다.






4.2.1.1.1 Example: Bayesian Inference, Mixture Distribution

MCMC 알고리즘은 \(p(\theta \rvert y ) = c \cdot p(\theta) L(\theta \rvert y)\) 로 표현될 수 있는 베이지안 추론에서 특히 강력하다. 베이지안 추론에서 \(c\)의 계산이 드럽게 어렵다는 것이 이것 이상의 다른 추론전략을 방해하기 때문이다.

보유하는 stationary 분포가 타겟 post인 MCMC에서 생산된 샘플들은 post 모먼트, tail 확률, 그리고 다른 유용한 quantity 계산에 쓰일 수 있다.

independent 체인에서 prior를 proposal 분포(\(g\))로 쓰자. 즉 \(f\)가 post, \(g\)가 prior다. 이런다면

\[ R ({\theta}^{(t)}, \theta^{\ast} ) =\dfrac {\dfrac {f(\theta^{( \ast )})} {g(\theta^{( \ast )})} } {\dfrac {f(\theta^{( t )})} {g(\theta^{( t )})} } =\dfrac {c \; \cdot \; \pi( \theta^{\ast} ) L( \theta^{\ast} \rvert y )} {c \; \cdot \; \pi( \theta^{(t)} ) L( \theta^{(t)} \rvert y )} \cdot \dfrac { \pi (\theta^{(t)} )}{ \pi (\theta^{\ast})} =\dfrac {\dfrac {\pi (\theta^{\ast} \rvert y)} {\pi (\theta^{\ast})} } {\dfrac {\pi (\theta^(t) \rvert y)} {\pi (\theta^{(t)})} } = \dfrac {L \left( \theta^{\ast} \rvert y \right)} {L \left( \theta^{(t)} \rvert y \right)} \]

  • Mixing Properties:
    • Good Mixing: 첫번째 그림의 MC는 시작점에서 빠르게 멀어지며 \(\delta\)에 대한 post에서의 모든 서포트에 해당하는 패러미터 space의 모든 부분을 훑으면서 샘플을 뽑아내는게 쉬워보인다.
    • Bad Mixing: 두번째 그림은 starting value에서 멀어지는 것도 느리고, posterior support의 영역을 탐색하는 게 시원찮아보인다.
  • Burn-in 이후의 실현값들을 히스토그램으로 만들어서 살펴보면 \(BETA(1,1)\) proposal 덴시티를 쓴 MCMC만이 \(\delta\)의 참값을 잘 모사하는듯.










4.2.2 Random Walk Chains (Most Widely Used)

MH method 에 해당.

let \(\textbf {X}^{\ast} = \textbf {X}^{(t)} + \epsilon , \epsilon \sim h(\epsilon )\). 이때 \(h\)는 임의의 덴시티.

  • \(h\) 로 자주 선택되는건 \(U, \textbf{N}, \text {Student's } t\).

이러면 proposal 덴시티 \(g\)는 어떻게 되는가? \(x^{'} g \sim N( \cdot \rvert x^{(t)}, \sigma^2 )\). 가장 빈번하게 쓰이는게 \(N\)이므로 \(N\)으로 설명. 여기서 \(x^{(t)}\)는 평균으로 사용되었고, \(\sigma^2\)은 Jumping Rule에 해당한다.

  • proposal can be
    • too diffused: Jumping rule is too big
    • too focused: Jumping rule is too small
  • Random Walk
    1. generate \(x^{'} g \sim N( \cdot \rvert x^{(t)}, \sigma^2 )\)
    2. \(u \sim U(0,1)\).
    3. calculate MH ratio: \(r = \dfrac {f(x')}{f(x^{(t)}} \dfrac {g(x' \rightarrow x^{(t)}} {g(x^{(t)} \rightarrow x'} = \dfrac {f(x')}{f(x^{(t)}}\), cause it’s \(N\).
    4. if \(u<r\), \(x^{(t+1)} = x'\). o.w., \(x^{(t+1)} = x^(t)\).






4.2.2.1 Example: Mixture Distribution

let proposal \(\delta^{(t+1)} = \delta^{(t)} + U(-a, a)\)1. 이인즉 몇몇 proposal들은 \([0, 1]\) 이외에서 생산된다.

  • note that \(\forall \delta \notin [0,1]\), post is zero,
  • Reparameterize the problem by letting \(U = \log{\dfrac {\delta}{1-\delta}}\)(logit). 왜? Probabilty Space is bounded: \(0 \le P(\cdot) \le 1\).

Run a random walk chain on \(U\) by adding \(U(-b,b)\).

Two ways of Reparamaterization:

  • Run MCMC in \(\delta\)-space. \(u\) 값을 다시 T해와서 \(\delta\)-space에서 돌림.
  • Run MCMC in \(u\)-space. \(u\)-space에서 돌리고 마지막에 모델 샘플들을 전부 \(\delta\)-space로 환원.


4.2.2.1.1 \(\delta\)-space에서 돌리는 방법

개략적으로는 \(T^{-1}: \delta' \leftarrow u' \sim g( \cdot \rvert u^{(t)})\) 와 같은 형을 띤다. 조건부 proposal \(g\)에서 생산된 u를 하나하나마다 \(\delta\)로 역변환해서 그 하나하나의 역변환 값으로 MH 알고리즘을 돌린다. proposal 덴시티 \(g( \cdot \rvert u^{(t)} )\) 는 \(\delta\)-space에서의 proposal 덴시티로 변환되어야 한다. 이경우 MH ratio는

\[ R (\textbf {x}^{(t)}, \textbf X^{\ast} ) =\dfrac {\dfrac {f(\delta^{\ast})} {g \left( logit(\delta^{\ast}) \rvert logit(\delta^{(t)}) \right) \cdot \left| J(\delta^{\ast})\right|} } {\dfrac {f(\delta^{(t)})} {g \left( logit(\delta^{(t)})|logit(\delta^{\ast}) \right) \cdot \left| J(\delta^{(t)})\right|} } \]

$J(^{(t)}) $ 는 $T: u $에 대한 \(J\)를 \(\delta^{(t)}\)에서 측정한 값. 주의해야 할 것이 해당 방법론에서는 \(T\)한 value를 사용하였으므로 \(g\)에 대한 \(J\)를 구해야 한다.


4.2.2.1.2 \(u\)-space에서 돌리는 방법

이 상황에서 쓰이는 proposal은 \(\dfrac {g(u^{(t)} \rightarrow u')} {g(u' \rightarrow u^{(t)})}\). \(\delta\)에 대한 타겟 덴시티는 \(u\)에 대한 덴시티로 변형되어야 한다. 이때 \(\delta = logit^{-1}(U) = \dfrac{\exp(U)}{1+\exp(U)}\) 였으므로, \(U^{\ast} = u^{\ast}\) 로 두었을 때 생산되는 MH ratio는



$$ R (^{(t)}, ^{} )

= { {f(logit^{-1} {(u^{})})} {g (u{}|u{(t)}) ) * | J(u^{(t)})|} } { {f(logit^{-1} {(u^{(t)})})} {g (u{(t)}|u{()}) ) * | J(u^{})|} } $$



우리가 transform value를 사용한데가 \(f\) 덴시티이므로 \(f\) 덴시티에 대한 야코비안을 붙여줘야 하는데 그 덴시티에 대한 야코비안은 \(u\)에 대한 야코비안이므로 쓰인 야코비안은 위와 같다. 이때 \(\lvert J(u^{\ast})\ \rvert = \dfrac {1} {\lvert J(\delta^{\ast}) \rvert}\) 이므로 위와 아래에서 만들어지는 MH ratio는 같다. 따라서 두 관점은 equivalent한 체인을 생산한다.

  • Sample paths for \(\delta\) from RW chains in Ex. 7.3, run in \(u\)-space iwth b=1 (top) and b=0.01 (bottom).

4.2.2.2 Example: Autocorrelation Plot (ACF)

배우지 않은 MCMC 방법론 중 하나.

reminder. thinning을 하더라도 거의 줄지 않아서 MCMC 샘플로서 거의 가치가 없는 케이스가 존재한다.

no










4.2.3 Basic Gibbs Sampler

  • need to derive the conditional density for all, 모든 joint density에 대해 coefficient를 1개씩 제한 상황의 모든 conditional density를 구한 후 GS 제작이 가능
    • if conditional densities are not available, we can use MH algorithm when updating \(x_i\) -> 이런식으로 접근할 경우 이를 MH-within-Gibbs라고 부른다.
      • 1PL IRT HW

let \(\textbf {X} = (X_1 , \cdots, X_p )^{'}\) , \(\textbf {X}_{-i} = (X_1 , \cdots, X_{i-1}, X_{i+1}, \cdots, X_p )^{'}\).

시작값 $^{(0)} $를 잡고, \(t=0\)으로 설정한다. 이후 각각을 \(t+1\) 단계의 시퀀스의 구성요소 각각을 \(X_i^{(t+1)} \vert \; \; \cdot \sim f \left( x_1 \rvert x_2^{(t)}, \cdots, x_p^{(t)} \right)\) 에 따라서 생산한다.

Gibbs Sampler의 stationary 분포는 \(f\).

\(X_i^{(t)}\)의 limiting 마지널 분포는 \(i\)번째 coordinate에 따른 타겟분포의 단변량 마지널化와 같다.

MH 알고리즘과 마찬가지로, 우리는 \(X\)의 임의의 함수 \(g(X)\)에 대해 \(E \left[ g(X) \right]\) 를 추정하기 위해 체인에서의 실현값을 사용할 수 있다.






4.2.3.1 Example: Fur Seal Pup Capture-Recapture Model

1800년대 후반 (by the late 1800s) 뉴질랜드 물범은 거의 전멸했다가 요즘 들어 폭증함 (abundance). 물범의 고름 숫자를 capture-recapture 사용해서 해보자.

사이즈 불명인 모집단의 크기 파악 위에 반복 연구 실행. 각 연구마다 포획었던 개체는 표식 새기고 풀어줌. 후속 연구에서 또 포획되면 재포획으로 표기. 높은 재포획 비율은 참 모집단 사이즈값이 포획되었던 개체들의 총량을 크게 넘지 않을 것임을 암시.

  • \(N\): 불명인 모집단 사이즈. \(l\)회의 조사 통해 얻어진 각 회의 총 포획 숫자는 각각 \(c=(c_1, \cdots, c_l)\)로 저장. 모집단 사이즈는 샘플링 동안에는 변동 없다(죽음, 출생, 이주 없음 inconsequential)고 가정한다.
  • \(r\): 연구 동안에 포획되었던 이질 동물들의 총 숫자.
  • 각 연구 시도에서 상응하는 구분되고 알려지지 않은 포획 확률은 \(\alpha = (\alpha_1 , \cdots, \alpha_l )\). 이 모델은 모든 동물들이 각 1회의 포획 발생에서 잡힐 가능성 자체는 각각의 동물에 대해서 동일하나, 이 被포획 확률은 시간이 지남에 따라 변할 수 있다는 것을 말함.

이 모델의 likelihood는

\[ L \left( N, \alpha \rvert c, r \right) \propto \dfrac {N!}{(n-r)!} \prod_{=1}^{l} \alpha_i^{c_i} \] Fur Seal Data for Seven Studies in One Season on the Otago Peninsula가 주어졌으며, prior은 \(\pi(N) \propto 1\), \(\pi (\alpha_i ) = BETA(\theta_1 , \theta_2)\) 이다. 계산하라.

해당 모델의 conditional posterior distribution에서 시뮬레이트하는 것으로 Gibbs Sampler를 제작할 수 있다.

\[ N^{(t+1)}-r \rvert \cdot \sim Negative Binomial \left( r+1, 1- \prod_{i=1}^7 \left( 1- \alpha_i^{(t)} \right) \right) \]

\[ \alpha_i^{(t+1)} \rvert \cdot \sim BETA \left( c_i + \dfrac{1}{2}, N^{(t+1)} - c_i + \dfrac{1}{2} \right) \]

for \(i= 1, \cdots, 7\), \(r = \sum_{i=1}^7 {m_i} =84\). 이는 unique fur seals were observed during the sampling period.

  • Split boxplots of \(\bar {\alpha}^{(t)}\) against $N^{(t)} $ for the seal pup example.
  • Estimated marginal posterior probabilities for \(N\) for the seal pup example.






4.2.3.2 MH-within-Gibbs Sampler

실제 implementation에 무지막지 유용하다. 이는 각각의 사이클을 GS의 사이클로 만들어놓고, conditional density는 MH 알고리즘으로 획득하는 것이다. Gibbs Sampler는 MH Sampler의 특별한 경우라고 볼 수 있다. MH 알고리즘의 proposal 분포를 시간에 따라 변화하도록 함으로써, GS와 MH 알고리즘 사이에 연결고리가 생긴다. 각 Gibbs 사이클은 \(p\) 개의 MH 스텝으로 구성되어 있다. 사이클 내에서의 \(i\)번째 Gibbs 스텝은, 체인의 현 상태 \((x_{1}^{(t+1)}, \cdots, x_{i-1}^{(t+1)}, x_{i+1}^{(\underline{t})}, \cdots, x_{p}^{(\underline{t})})\) 가 주어졌을 때, 효과적으로 후보 벡터 \((x_{1}^{(t+1)}, \cdots, x_{i-1}^{(t+1)}, {\underline{X_i^{*}}}, x_{i+1}^{(\underline{t})}, \cdots, x_{p}^{(\underline{t})})\) 를 생산한다. 밑줄 차이점에 주목.

\(i\)번째 단변량 Gibbs 업데이트는 이하와 같이 MH 스텝 drawing으로 볼 수 있다.

\[ {\underline{X_i^{*}}} \rvert (x_{1}^{(t+1)}, \cdots, x_{i-1}^{(t+1)}, x_{i+1}^{(\underline{t})}, \cdots, x_{p}^{(\underline{t})}) \sim g_i \left( \cdot \rvert (x_{1}^{(t+1)}, \cdots, x_{i-1}^{(t+1)}, x_{i+1}^{(\underline{t})}, \cdots, x_{p}^{(\underline{t})}) \right) \]

where

\[ g_i \left( \cdot \rvert (x_{1}^{(t+1)}, \cdots, x_{i-1}^{(t+1)}, x_{i+1}^{(\underline{t})}, \cdots, x_{p}^{(\underline{t})})\right) = \left\{ \begin{array}{@{}lr@{}} f(x_i^{*} \rvert x_i^{(t)} ), & \text{if } \; \; \; X_i^{*} = x_i^{(t)} \\ 0 & o.w. \end{array}\right\} \]

이 경우, MH ratio는 1과 같아진다. 즉슨 모든 후보들은 언제나 accept 된다. 즉슨 GS는 MH 알고리즘에서 acceptance ratio가 항상 1인 경우에 해당한다. 따라서 conditional density을 구할 수만 있으면 GS를 사용하는게 좋다. 샘플을 버릴 필요가 없고, 버려지는 샘플이 없기 때문.






4.2.3.3 Update Ordering

Random Scan Gibbs Sampling: 기본 GS에서 $ $ 에 가해지는 업데이트의 순서는 한 사이클에서 다음 사이클로 넘어갈 때마다 바뀔 수 있다. 패러미터들이 높은 수준에서 상호연관되어있을 경우, 각 사이클을 랜덤하게 순서배치하는 것은 효과적일 수 있다.2 특정 모델에 대한 전문화된 지식이 없다면, 한 이터레이션에서 다음으로 넘어갈 때 패러미터끼리 높이 상호연관되어있다면 deterministic 한 방법과 RSGS 양쪽 모두를 시도해보는 것이 권장된다.

4.2.3.4 Blocking

with \(p=4\), e.g., 각 사이클을 다음 절차를 따르면서 업데이트:

  1. \(X_1^{(t+1)}\rvert \cdot \sim f \left(x_1 \rvert x_2^{(t)}, x_3^{(t)}, x_4^{(t)} \right)\).
  2. \(X_2^{(t+1)}, X_3^{(t+1)}\rvert \cdot \sim f \left(x_2, x_3 \rvert x_1^{(t+1)}, x_4^{(t)} \right)\).
  3. \(X_4^{(t+1)}\rvert \cdot \sim f \left(x_4 \rvert x_1^{(t+1)}, x_2^{(t+1)}, x_3^{(t+1)} \right)\).

Blocing은 \(X\)의 구성요소들이 서로 상관관계가 있을 때 \(X_i\) 내부가 상관이라는 건가, 아니면 \(X_t, X_{t+1}\) 이 상관이라는 건가? 유용함. 해당 알고리즘을 통해 더욱 상관된 구성요소끼리는 한 블럭 안에서 샘플링됨.

4.2.3.5 Hybrid Gibbs Sampling

하나 이상의 \(X\)에 대한 조건부 분포는 대부분 closed form으로 만들 수 없음. 깁스 샘플러의 주어진 스텝에서, 적절한 조건부 분포에서 샘플링하기 위해 MH 알고리즘이 쓰인다면 Hybrid MCMC 알고리즘이 완성됨.

with \(p=5\), e.g., 하이브리드 MCMC 알고리즘은 다음 절차를 따르면서 업데이트:

  1. Update \(X_1^{(t+1)} \rvert \left( x_2^{(t)}, x_3^{(t)}, x_4^{(t)}, x_5^{(t)} \right)\) with 깁스 스텝.
  2. Update ( x_2^{(t+1)}, x_3^{(t+1)} ) $ ( x_1^{(t+1)}, x_4^{(t)}, x_5^{(t)} )$ with MH 스텝.
  3. Update \(X_4^{(t+1)} \rvert \left( x_1^{(t+1)}, x_2^{(t+1)}, x_3^{(t+1)}, x_5^{(t)} \right)\) with 랜덤워크.
  4. Update \(X_5^{(t+1)} \rvert \left( x_1^{(t+1)}, x_2^{(t+1)}, x_3^{(t+1)}, x_4^{(t+1)} \right)\) with 깁스 스텝.
4.2.3.5.1 Example: Fur Seal Pup Capture-Recapture Study





4.2.4 Implementation

MCMC의 목적은 타겟분포 \(f\)의 특징들을 알아내는 것. 모든 MCMC는 정답인 limiting 정적분포를 가지고 있음. 실전에는 체인을 얼마나 충분히 오래 돌릴지를 결정하는 게 중요함. \(X\)의 dimensionality(차원)이 높다면 수렴이 엄청 느려서 엄청 긴 run을 요할 수도 있음. 이하의 요건들을 생각해서 long run을 결정해야 함.

  • Has the chain run long enough?
  • Is the first portion of the chain highly influenced by the starting value?
  • Should the chain be run from several different starting values?
  • Has the chain traversed all portions of the region of support of \(g\)?
  • Are the sampled values approximate draws from \(f\)?
  • How shall the chain output be used to produce estimates and assess their precision?



4.2.4.1 Ensuring Good Mixing and Convergence

MCMC 알고리즘이 대상 문제에 얼마나 쓸만한 정보를 주는지 고민해야 함. 이는 곧

  1. 체인이 얼마나 빠르게 체인의 starting value를 까먹는가
  2. 얼마나 빠르게 체인이 타겟분포 \(f\)의 모든 서포트를 훑는가
  3. 체인이 그것의 정적분포에 근사적으로나마 닿는가를 고민.
  4. There is substantial overlap between the goals of diagnosing convergence to the stationary distribution and investigating the mixing properties of the chain.



4.2.4.2 Simple Graphical Diagnostics

트레이스 플롯 (sample path)은 이터레이션 횟수 \(t\)와 \(X^{(t)}\)의 실현값 간의 플롯이다.

  • 체인의 mixing이 구리면 이는 장기간의 이터레이션동안 동일값 근처에서 머무르게 됨.
  • 체인의 mixing이 좋으면 시작값에서 빠르게 떠나서 \(f\)의 서포트에 해당하는 영역을 열심히 훑음.



autocorrelation 플롯은 \(X^{(t)}\)의 시퀀스에서 다른 이터레이션 래그에서의 상관관계를 서술한다. 래그 \(i\)에서의 autocorrelation은 \(i\) 이터레이션만큼 떨어진 이터레이트 간의 상관관계이다. 구린 mixing 프로퍼티를 가지는 체인은 이터레이션 간의 래그가 증가하더라도 autocorrelation의 부식(decay)이 느림.

MCMC 체인에서의 첫 \(D\) 개의 값은 보통 burn-in period라고 해서 버려짐. 체인의 시작점에 대한 의존이 강하게 남아있을지도 모르기 때문. 어느정도가 적절한 번인 피리어드인가는 Gelman-Rubin diagnostics에 의해 결정된다. 결정된 번인 피리어드가 제대로 된 값을 못내면 \(D\)가 늘어나던가 \(L\)이 늘어나던가 둘다 늘리던가 해야함. - Motivated by an Analysis of Variance - 사슬간 분산 (between-chain variance)이 사슬내 분산 (within-chain variance)보다 유의하게 크면 번인 피리어드나 MCMC 길이가 늘어나야 함 - Difficulties - multimodal \(f\)에서 적절한 스타팅 밸류를 찾는건 어렵고, 체인이 로컬 region이나 mode에서 갇힐수 있음 - 이의 단일차원성 (uni-dimensionality) 때문에 타겟분포 \(f\)가 멀티디멘션이면 타겟분포의 수렴에 대한 정보에 대한 잘못된 직관을 줘버릴 수 있음



proposal \(g\) 선택할 때 고려해야할 요소는? mixing은 proposal 분포 \(g\)의 특질에 큰 영향을 받으며 특히 이의 스프레드(spread)에 큰 영향을 받음. 타겟분포 \(f\)와 \(g\)간의 닮음에 있어서, 프로포절 \(g\)의 tail behavior의 닮음은 high density의 닮음보다 훨씬 중요함. \(f/g\)가 bounded 라면, MC의 정적분포로의 수렴은 대체로(overall) 빠름. 실전에선 정보를 줄 수 있는 이터레이티브 프로세스를 통해 proposal 분포의 분산이 택해질 수 있음. (In practice, the variance of the proposal distribution can be selected through an informal iterative process.) 20%~50% 사이의 acceptance rate가 선호되어야 함.



  • Reparameterization

모델의 reparameterization은 MCMC 알고리즘의 mixing behavior에 상당한 기능향상을 가져올 수 있음. reparameterization은 dependence를 낮추기 위해 가장 우선되어야 할 전략 중 하나임. 서로 다른 모델들은 서로 다른 reparameterization 전략을 적용해야 함. reparameterization 접근법은 보통 특정 모델에 대한 원오프로서 채택되므로 일반화된 조언을 하기에는 어려운 부분이 있음.



  • Comparing Chains

MCMC 실현값이 크게 상관관계있다면, MCMC의 각 이터레이션에서 주어지는 정보는 run length에서 주어지는 정보 대비 보잘것 없다(will be less than suggested by the run length). 여기서 reduced information은 effective sample size라고 칭해지는 더 작은 iid 샘플에 담겨있는 정보와 동등하다. 여기서 샘플의 총 숫자와 effective sample size 사이의 차이는 잃게 된 효율을 의미한다. 관심있는 변량을 확인하기 위해 우리가 MC 체인에서의 correlated 샘플들을 사용했을 때 잃게된 효율.



  • Number of Chains

모델 진단에 있어서 가장 진단을 어렵게 하는 부분은 체인이 타겟분포 \(f\)의 1개 이상의 모드에 걸리냐 안걸리냐 하는 부분. 이 경우 모든 수렴진단은 체인이 수렴하긴 한다는 결론을 내리게 된다. 정작 체인이 타겟분포 \(f\)의 모든 특성을 나타내지 못하는데도. 그래서 여러번 돌려보게 되는 것. 최소한 그 여러번의 run들 중 체인 1개에서라도 타겟분포 \(f\)의 모든 흥미로운 특질들이 드러났으면 하니까. 이러한 특질을을 찾아내는데 실패한 개별 체인들의 경우에는 mixing을 더 좋게 만들기 위해 체인 길이가 늘어나거나 문제가 reparameterize 되어야 함.