4.4 Auxiliary Variable MCMC

실전에서 마주치는 대부분의 상황에서 ABC나 HMC 문제를 제외하고는 대부분의 경우 MCMC 문제를 완벽하게 풀어내는 건 불가능. 이때 주어진 variable 말고 보조변수 (Auxiliary Variable)을 추가함으로써 시뮬레이션 품질을 좀 더 높일 수 있지 않을까 하는 것이 논하고자 하는 바.

4.4.1 Introduction

• Difficulties with MH Algorithm. 일반적인 MH 알고리즘으로 풀어낼 수 없는 2가지 상황이 존재:
  1. Local-trap problem: 에너지 계가 울퉁불퉁한 complex system에서 시뮬레이션을 진행했을 때 끝없이 로컬 최적값에서 빠져나오지 못함. 시뮬레이션을 비효율적으로 만듬.
     • density가 높다는 것은 해당 파트의 에너지가 낮다는 것이며, density가 낮은 에너지가 많은 파트에서 high density로 가는 것은 쉽고 자주 일어나도 역은 드뭄. 조밀하면 움직일 여력이 없으니까. 이것이 local trap의 원인
     • 에너지는 이하로 표시 가능: energy function \(= -log \pi(\theta \vert x)\), 즉 negative log posterior, 혹은 negative log density.
  2. Doubly-intractable normalizing constants problem:
     • Inability to sample from distributions with intractable integrals
       • 보통이라면, \(pi(\theta \vert x) \propto \kappa(x) f(x\vert\theta)\pi(\theta)\). \(r= \dfrac{pi(\theta ' \vert x)}{pi(\theta^{(t)} \vert x)} = \dfrac{\kappa(x) f(x\vert\theta ' )\pi(\theta ' )}{\kappa(x) f(x\vert\theta^{(t)})\pi(\theta^{(t)})}\) 과정에서 normarlizing constant \(\kappa\)가 알아서 캔슬되어 MH 돌리는데 문제가 없음.
     • let \(f(x) \propto \kappa(x;\theta) \psi(x)\) 는 알고자 하는 분포. 여기서 \(\kappa(x)\)는 unnormalized density의 함수. 이때 \(\kappa(x)\)는 패러미터의 함수이며 각 이터레이션의 다른 패러미터 추정값마다 변화해버려서 캔슬되지 않음. 그러면 계산하면 되는거 아님? 계산 불가능한 상황 존재 - nearly infinite summation or integration 포함하는 경우. (ex:) 이는 곧 intractable integral. acceptance \(Pr\)이 알 수 없는 비 \(\frac{\kappa(x')}{\kappa(x)}\)를 포함하므로 MH 알고리즘은 사용불가.
       이러한 문제는 bayesian 추론에서 spatial statistical models, random effects models, 그리고 exponential random graph models 등 다양한 통계적 모형에서 부딪히게 된다.
       • ex: Lattice system of areal model (Lattice의 승만큼 연산 필요)
       • e.g., Random Effect Model. 이때는 각 individual별로 Random Effect를 integration 해줘야 하므로 문제터짐
       • ex: Exponential Random Graph model: 네트워크에 사용되는 모델. 얘도 power임.
       • 이러한 상황에서는 대부분의 optimization 알고리즘도 다 먹통됨

• 이러한 2개의 문제점을 극복하기 위해 다양한 진보된 MCMC 방법론이 제시되었음.
  1. Auxiliary variable-based methods
  2. Population-based methods
  3. Importance weight-based methods
  4. Stochastic approximation-based methods

4.4.1.1 Auxiliary Variable MCMC Methods

\(f(x)\)를 가지는 mv 분포에서의 샘플링을 생각해보자. Rao-Blackwellization(#)이 MC 시뮬레이션에의 최우선원칙임은 알려져 있다. 시뮬레이션의 수렴을 좀 더 강화하기 위해 우리는 가능한한 많은 \(x\)의 구성물을 integrate하는 것을 시도해보아야 한다. 하지만 이하의 두가지 경우(이외에도 존재)에 시뮬레이션을 양질로 만들기 위해 우리는 1개 이상의 변수를 추가하는 상황을 고려할 수 있다.

1. 타겟분포 \(f(x)\)가 multimodal. 온도 혹은 아직 관측되지 않은 측정값과 같은 auxiliary variable이 계가 로컬 트랩에서 빠져나올 수 있도록 도움을 줌. multimodal 상황.
2. 타겟분포 \(f(x)\)가 intractable normalizing constant 포함. \(X\)의 auxiliary 실현값이 시뮬레이션에 포함됨으로써 시뮬레이션에서 normalizing constant 를 무력화시킴.

MH 알고리즘 \(\dfrac{ f(\theta ' \vert x )}{f(\theta^{(t)} \vert x )} \dfrac{ g(\theta ' \rightarrow \theta^{(t)} )}{ g(\theta^{(t)} \rightarrow \theta ')}\)은 이하의 2가지 기본적인 부품을 가지고 있다.

1. 타겟분포 (左)
2. proposal 분포 (右)

이에 더해서 auxiliary variable 방법론은 이하의 2가지 방법으로 행해질 수 있다. 타겟과 제안 어느쪽에 변수를 추가하는지에 대한 이야기이다.

1. 타겟분포 augmentation 방법론: Augmenting auxiliary variables to the target distribution
   • auxiliary variable \(u\)와 조건부 분포 \(f(u \rvert x )\)를 정의한다. joint 분포 \(f(x,u) = f(u \rvert x) f(x)\)를 만들기 위해. 이후 MH 알고리즘이나 GS를 사용해 \((x,u)\)를 업데이트. \(f(x)\)의 샘플은 \((X, U)\)의 실현값 \((x_1, u_1), \cdots, (x_N, u_N)\)를 이용해 marginalization이나 프로젝션 등을 이용해 획득될 수 있다.
2. Method of Proposal Distribution Augmentation: Augmenting auxiliary variables to the proposal distribution.
   • proposal 분포 \(T(x', u \rvert x)\)를 특정하고, 이의 reversible version \(T(x, u \rvert x')\)도 특정한다. 즉슨 \(\int T(x', u \vert x)du = T(x' \vert x)\), \(\int T(x, u \vert x')du = T(x \vert x')\)의 관계가 성립한다.
     이제 proposal \(T(x', u \vert x)\) 로부터 후보 (candidate) 샘플 \(x'\)를 생산하고, 이를 with probability \(\min \left\{ 1, r(x, x', u) \right \}\). 이때 \(r(x, x', u) = \dfrac {f(x')} {f(x)} \dfrac {T(x,u \vert x')} {T(x',u \vert x)}\).

실현값 (realizations) \(x_1 , \cdots, x_N\)을 생산할 때까지 이를 반복한다. 이제 \(N\)이 충분히 크다면, 이 실현값들은 근사적으로 \(f(x)\)에 의해 분포되어 있다.

이러한 방법론의 타당성은 이하를 통해 보일 수 있다.

\[ K(x' \vert x) = \int_{\mathcal{u}} s(x, x', u) du + \mathbf{1}(x=x') \left[ 1-\int_{\mathcal{X}}\int_{\mathcal{u}} s(x, x', u) du dx' \right] \]

이는 \(x\)로부터 \(x'\)로의 integrated transition kernel 을 의미하며, 이때 \(s(x, x', u) = T(x', u \rvert x) \ast r(x, x', u)\). Then,

\[ f(x) \int_{\mathcal{u}} s(x, x', u) du = \int_{\mathcal{u}} \min \left[ f(x') T(x, u \vert x'), \; \; f(x)T(x', u \vert x) \right] du \]

이는 \(x\)와 \(x'\) 에 대해 symmetric. 이는 곧 \(f(x)K(x' \vert x) = f(x')K(x \vert x')\) 임을 의미한다.

original density?

4.4.2 Multimodal Target Distribution

4.4.2.1 Simulated Tempering

분포 \(f(x) \propto \exp \left(-H(x) \right), x \in X\) 에서 샘플링하는데에 관심이 있다고 하자. simulated annealing에서 그러했던 것처럼, simulated tempering \(f(x, T) \propto \exp \left( -\dfrac {H(x)} {T} \right)\)로 타겟 분포를 확장시켰다. 이는 auxiliary variable인 temperture \(T\)를 포함함으로서 이루어진다. \(T\)는 사용자가 미리 지정한 값들의 finite set이 된다. \(H(x)\)는 사실상 energy function.

Temperature Transition Matrix \(T = \begin{bmatrix} q_{11} & q_{12} & \cdots & q_{1n} \\ q_{21} & \ddots & & \\ \vdots & & \ddots & \\ q_{n1} & \cdots & & q_{nn} \end{bmatrix}\). 이때 row는 current 온도 \(T_1, \cdots, T_n\), column은 행선지 온도.

---

Parallel Tempering은 인접한 온도로만 이동 가능 (가장 높은 온도에서 가장 낮은 온도로 한단계 한단계씩). 온도 자체를 시뮬레이션한게 아니라 온도의 chain이 주어져 있어 각 온도 간의 움직임을 만드는 것에 그친다. 따라서 이는 multiple chain을 이용하는 population MC 방법론 쪽에 소속됨.

Simulated Tempering과는 이 점에서 차이를 보임. 후자는 어느 온도로든 다 이동. 온도 매트릭스 만들어놓고, \(U(0,1)\) 분포에서 온도 하나 생산하고 이 온도로 이동할 것인지의 여부를 MH 알고리즘으로 결정.

---

\(U(0,1)\)에서 랜덤하게 숫자를 뽑고, \(j\)의 값을 proposal transition matrix \((q_{ij})\)에 따라서 정한다. \(u<q_{11}\)이면 \(T_1 \rightarrow T_1\), \(q_{11}<u<q_{11} + q_{12}\)이면 \(T_1 \rightarrow T_2\), …. - if \(j=i_t\), let \(i_{t+1}=i_t\), and let \(x_{t+1}\)을 MH kernal \(K_{i_t}(x,y)\)에서 뽑는다. 이때 \(K_{i_t}(x,y)\)는 \(f(x, T_{i_t})\)을 invariant distribution로 허용하는 아이이다. 즉 새로운 \(x\)를 생산하면 된다. - if \(j \not= i_t\), let \(x_{t+1}=x_t\)하고 proposal을 이하의 \(Pr\)에 따라 채택한다. 이때 \(Z\)는 \(Z_i\)의 측정값이다. 채택된다면 \(i_{t+1} = j\)이고, 그외의 경우에는 \(i_{t+1} = i_t\)로 한다. 새로운 \(x\)를 생산하는 것이 아니라 들고 있던 \(x\)를 쓰되, 이걸 accept 할건지 안할건지를 체크한다.

\[ \min \left[ 1, \; \; \dfrac {\hat Z_j} {\hat Z_{i_t}} \ast \exp \right \{ -H(x) \left( \dfrac{1}{T_j}-\dfrac{1}{T_{i_t}} \right) \left \} \cdot \dfrac {q_{j,i_t}} {q_{i_t , j}} \right] \]

이때 \(\dfrac {q_{j,i_t}} {q_{i_t , j}}\) 는 proposal distribution이라고 생각할 수 있다. 나머지는 Likelihood part이며, 이때 \(\dfrac {\hat Z_j} {\hat Z_{i_t}}\) 가 normalizing constant의 ratio이다. 온도가 변화하였으므로 두 식의 normalizing constant가 같지 않기 때문이다.

---

Issues on Simulated Tempering:

1. Temperature Ladder를 어떻게 고를 것인가. → 각 chain별로 이동이 원활하게 잡는 것이 핵심.
   • 가장 높은 온도 \(T_1\)은 대부분의 uphill move가 해당 레벨에서 accept 될 수 있도록 설정되어야 한다.
   • 사이의(intermediate) 온도들은 sequential manner로 설정될 수 있다. \(T_1\)에서 시작해서, 점차적으로 다음으로 낮은 온도를 \(Var_i \left\{ H(x) \right\} \ast \delta^2 = O(1)\)을 만족하도록 설정하는 것이다. 이때 \(\delta = \dfrac {1}{T_{i+1}} - \dfrac{1}{T_i}\)이며, \(Var_i(\cdot)\)은 \(H(x)\) (taken with respect to \(f(x, T_i)\)) 의 분산을 의미한다.
     • 이러한 조건들은 \(f(x,T_i), f(x,T_{i+1})\) 사이에 상당히 겹치는 점이 많아야 한다는 것을 의미하기도 한다. 실전에선 \(Var_i \left( H(x) \right)\)는 샘플러를 레벨 \(T_i\)에서 예비적으로(preliminary) 돌려보았던 결과에서 러프하게나마 예측될 수 있다.
       
2. \(Z_i\)를 어떻게 estimate 할 것인가. → accept 여부가 normalizing constant에도 의존해서 이거 이상하게 고르면 효율 떨어짐. 엄청난 단점이라서 요즘은 이 알고리즘 자체를 잘 안씀
   • 이는 simulated tempering의 효율에 직결되는 부분이다. \(Z_i\)들이 잘 estimate 되었다면, simulated tempering은 temperature ladder을 따라 symmetric RW처럼 동작한다. (\(x\)-updating step을 제하고 볼 경우) 그렇지 않다면 이는 특정 temperature 레벨에서 멈춰버린다. 시뮬레이션이 실패함은 물론이다(rendering).
     실전에서 \(Z_i\)들은 stochastic approximation MC 방법론을 사용해서 estimate 가능하다. 혹은 reverse logistic regression 방법론을 사용해서도 \(Z_i\)를 estimate 할 수 있다.

4.4.2.2 Slice Sampler

density \(f(x), \; \; \; x \in \mathcal X\)에서 샘플링하고자 한다. \(x \sim f(x)\)에서 샘플링하는 것은, \(f(x)\) 그래프 이하의 영역에서 uniform하게 샘플링하는 것과 동등하다. 해당 영역은 \(A = \{ (x,u): 0 \le u \le f(x) \}\)이며, 이것이 acceptance-rejection 알고리즘의 기초(basis)였다. 이 목적을 달성하기 위해 우리는 타겟분포 \(f\)를 auxiliary variable \(U\)를 사용하여 확장해볼 수 있다. 이 \(U\)는, \(x\)에 대해서 조건부일 때, 구간 \([0, f(x)]\)에서 uniform하게 분포되어 있다.

따라서, \((X, U)\)의 joint density function은 \(f(x,u)=f(x)f(u \rvert x) \propto I_{(x,u)\in \textit A}\). 후자의 인디케이터는 언급되었던 영역 안에 속한다는 의미. 이는 GS에 의해 이하와 같이 샘플링 가능하다.

1. draw \(u_{t+1} \sim U[0, f(x_t)]\).
2. draw \(x_{t+1}\) uniformly from the region \(\{ x: f(x) \ge u_{t+1} \}\).

위의 샘플러는 slice sampler라고 불림. 이는 multimodal 분포들에 대해 단순 MH 알고리즘보다 더 나을 가능성이 있음. slice 때문에 b/w-mode-transition에 자유롭기 때문. 현재도 핫한 샘플러중 하나임. horseshoe prior 에서의 패러미터 estimate에 대표적으로 이녀석이 쓰인다.

4.4.3 Doubly-intractable Normalizing Constants

Spatial models, e.g., the autologistic model, the Potts model, and the autonormal model (Besag, 1974)는 많은 과학적 문제들을 위한 모델링에 쓰이고 있음. 이러한 모델들에 해당하는 주요한 문제는 normalizing constant가 doubly-intractable하다는데 있음.

for dataset \(X\), 패러미터 \(\theta\), normalizing constant \(\kappa (\theta)\). 이때 \(\kappa (\theta)\)는 \(\theta\)에 의존하나 closed form으로는 만들 수 없음. 이하는 dataset을 생산한 likelihood function.

\[ \begin{align*} X \sim f(x \vert \theta) = \dfrac{1}{\kappa (\theta)} exp \{ -U(x, \theta) \}, &x \in \mathcal{X}, &\theta \in \Theta \end{align*} \]

\(\pi(\theta)\)는 \(\theta\)의 prior. 이 경우 post는 \(f(\theta \vert x) \propto \dfrac{1}{\kappa (\theta)} exp \{ -U(x, \theta) \} \ast \pi(\theta)\).

4.4.3.1 Boltzmann Density

known as Ising Model, 그리고 ~로 확장될 경우 autologistic model.

Consider a 2-D Ising model with the Boltzmann density

\[ f(\pmb x) \propto \exp \left\{ K \sum_{i\sim j} x_i x_j \right\} \]

• spins \(x_i = \pm 1\) (S극이 -1)
• \(K\)는 inverse temperature (measure for interaction : \(x_i\)가 주변에 있는 값과 얼마나 많은 같은 값을 가지는지, 다른 값을 가지는지에 대해 측정해주는 패러미터) 온도가 낮을수록 interaction가 강해지며, 이에 의해 동일값 확률이 높아짐.
• \(i\sim j\)는 lattice 상의 가장 가까운 neighbors.

온도가 높다면, 이 모델은 GS를 사용해 쉽게 시뮬레이션 가능하다. 조건부 분포에 따라 각 spin의 값을 iteratively 초기화한다. 아래의 식에서 \(n(i)\)는 spin \(i\)의 neighbors의 집합 (set). 이하의 수식은 autologistic 과 그 과정이 유사하다.

$$ \[\begin{align*} P(x_i =1 \vert x_j, \; \; j \in n(i)) &= \dfrac {1}{1+ \exp \left \{ -2K \sum_{j \in n(i)} \right\}} \\ P(x_i =-1 \vert x_j, \; \; j \in n(i)) &= \dfrac {\exp \left \{ -2K \sum_{j \in n(i)} \right\}}{1+ \exp \left \{ -2K \sum_{j \in n(i)} \right\}} &= 1- P(x_i =1 \vert x_j, \; \; j \in n(i)) \end{align*}\] $$

하지만, GS는 temperature가 critical temperature로 근접하거나 이하로 내려갈 경우 GS가 빠르게 느려진다. 온도가 낮으면 interaction이 강해져, 주변값과 비슷한 값을 generate 해야만 하기 때문이다. 이렇게 샘플링이 어려워지는 지점, 온도를 critical point라고 부른다. 이는 대략 \(\theta \approx 0.43\). 이것이 소위 critical slowing down 이라고 불리는 현상이다.

4.4.3.1.1 Perfect Sampler

과거 샘플들의 굉장히 많은 조합을 커플링해서 샘플을 생산. previous realization 전체에 대해 (이는 그 이전의 샘플, 아니면 그 이전의 샘플, 혹은 original 데이터에 대해서조차도) independent한 샘플을 생산해내는 sampler. 즉 그 어떤 것에서도 independent한 sample을 생산해낸다. 문제는 이 샘플러는 \(\theta>0.43\)인 순간 바로 작동을 안함. \(\theta>0.32, 0.35\) 정도로 엔간 크기만 해도 드럽게 느림.

4.4.3.1.2 Swendsen-Wang Algorithm

slice sampling에서 Boltzmann 덴시티는 이하의 형으로 다시 쓰여진다. 이때 \(\beta = 2K\). indicator function으로 변형했을 때 저 둘이 어떻게 equation이 성립하는지 유의.

\[ f(\pmb x) \; \propto \; \prod_{i\sim j} \exp \left\{ K(1+x_i x_j) \right\} \; = \; \prod_{i\sim j} \exp \left\{ \beta \ast \mathbf{1}(x_i = x_j) \right\} \]

이때 우리가 auxiliary variable \(\pmb u = (u_{i \sim j})\), where each component \(u_{i \sim j}\), conditional on \(x_i\) and \(x_j\), is uniformly distributed on \(\left[ 0, \; \exp \{\beta \ast \mathbf{1}(x_i = x_j)\} \right]\), then

\[ f(\pmb x, \pmb u) \; \propto \; \prod_{i \sim j} \mathbf{1} \left( 0 \le u_{i \sim j } \le \exp\left\{ \beta \ast \mathbf{1} (x_i = x_j) \right\} \right) \]

이때 \(u_{i \sim j}\) 자체는 bond variable이라고 명명된다. 이는 spin \(i\)와 spin \(j\) 사이의 가장자리에 물리적으로 앉아 있는 변수로서 생각될 수 있다. (i와 j가 묶여져 있는지, 같은 group 안에 존재하는 것인지 아닌지에 대한 indicator가 되는 variable)

• if \(u_{i \sim j}>1\), then \(\exp \left\{ \beta \ast \mathbf{1}(x_i = x_j) \right \}>1\), 따라서 반드시 \(x_i = x_j\).
• if \(u_{i \sim j}<1\), 이 경우 \(x_i, x_j\)에 제약 (constraint) 이 없다.

\(b_{i \sim j}\)가 제약에 대한 indicator variable이라고 정의하자. 즉, \(x_i, x_j\)가 같도록 제약되었다면, \(b_{i \sim j}=1\)이며 이외엔 0이다. for any 2개의 “like-spin” (i.e. 2개의 spin이 같은 값을 가진다) neighbors에 대해서, 이 둘은 with probability \(1-\exp (-\beta)\)를 따라 bonded 될 수 있는 가능성이 있음을 기억하라. \(\pmb u\)의 설정 (configuration)에 따라, “mutual bond” (i.e., \(b_{i \sim j}=1\)) 을 통하여 연결될 수 있는지 없는지 여부에 따라 spin들을 군집 (cluster) 할 수 있다. (위에서 i와 j가 같다고 indicator가 판별했을 경우에만 이런 cluster 로 묶는 것이 가능하다) Then 동일 클러스터 내의 모든 spin은 같은 값을 가질 것이다. 또한 군집 내부의 모든 spin을 동시에 뒤집는 (flip) 것은 \(f(\pmb x , \pmb u)\)의 평형 (equilibrium)을 해치지 않을 것이다.

Proceeds:

1. Update the bond values: check all “like-spin” neighbors, and set \(b_{i \sim j}=1\) with probability \(1-\exp (-\beta)\).
2. Update the spin values: Cluster spins by connecting neighboring sites with a mutual bond, and then flip each cluster with probability \(0.5\).

For the Ising model, the introduction of the auxiliary variable \(\pmb u\) has the dependence between neighboring spins partially decoupled, and the resulting sampler can thus converge substantially faster than the single site updating algorithm. As demonstrated by Swendsen and Wang (1987), this algorithm can eliminate much of the critical slowing down.

같은 값들이 모여있는 cluster를 판별하여 각각을 grouping. grouping을 랜덤으로 하므로 인접해 있는 동일값임에도 그룹에 포함되지 못하는 경우가 존재함. Swendsen-Wang에서는 이렇게 그룹을 만든 후, 해당 그룹을 통채로 toggling. group을 통채로 토글링하기 때문에 dependency가 있는 것들이 통채로 toggling되어서 dependency가 있는 것들은 나머지 것들과 인제 이렇게 independent한 것도 있지만 dependent한 것을 통채로 묶어서 하는 것이므로 좀더 한꺼번에 뒤집으니까 실제로 우리가 업데이트하는 것은 group 내부 말고 group 외부 간들에는 independent하다고 가정될 수 있는 몇몇개의 group들만이 남음. 이 덩어리들을 한꺼번에 업데이트하므로 따라서 샘플러 generate가 상대적으로 쉬움. 하지만 이 만든 덩어리는 매 이터레이션마다 덩어리를 새로 만들어야 함. 매 이터레이션마다 클러스터를 새로 만들고 flip하여 이를 accept할지 말지를 결정하는 이런 형태의 구조를 가짐.

4.4.3.2 Møoller’s Algorithm

auxiliary variable \(y\), 이는 \(x\)와 같은 state space를 공유한다고 정의. 그 경우 이하의 joint pdf \(f\) 를 생각해볼 수 있다. \(f(y \vert \theta , x)\)는 \(y\)의 분포.

\[ f(\theta, y \vert x) = f(x \vert \theta) \ast f(\theta) \ast f(y \vert \theta , x) \]

\(f(\theta, y \vert x)\) 에서 MH 알고리즘을 통해 시뮬레이트하기 위해서는 이하와 같은 제안분포 \(q\) 를 사용해볼 수 있다. 이는 패러미터 벡터 \(\theta \rightarrow \theta'\)의 usual change에 상응하며, 이 후에는 \(q(\cdot \vert \theta ' )\)에서 \(y'\)를 추출하는 exact sampling step이 따른다.

$$ q(’ , y’ , y) = q(‘, y) q(y’ ’) $ $

\(q(y' \vert \theta ' )\)가 \(f(y' \vert \theta)\)로 설정되었다면, MH ratio \(r\)은 이하와 같이 쓰일 수 있다. 이때 unknown normalizing constant \(\kappa(\theta)\)가 상쇄 (cancel) 되었음에 주목하라.

$$ \[\begin{align*} r(\theta, y, \theta', y' \vert x) &= \dfrac {f(x \vert \theta') f(\theta') f(y' \vert \theta' , x) \ast q(\theta\vert \theta' , y') q(y \vert \theta)} {f(x \vert \theta) f(\theta) f(y \vert \theta , x) \ast q(\theta'\vert \theta , y) q(y' \vert \theta')} \\ &= \dfrac {f(\theta', y' \vert x)}{f(\theta, y \vert x)} \ast \dfrac {q(\theta , y \vert \theta' , y'))}{q(\theta' , y' \vert \theta , y)} \\ &= \dfrac { \dfrac {f(\theta', y' \vert x)}{q(\theta' , y' \vert \theta , y)} } { \dfrac{f(\theta, y \vert x) }{q(\theta , y \vert \theta' , y'))} } \end{align*}\] $$

여기서 계산을 간단하게 하기 위해 제안분포 \(q\)와 auxiliary distribution을 이하와 같이 정리하는 것을 생각해볼 수 있다. 이때 \(\hat \theta\)는 \(\theta\)의 estimate로써, 예를 들어 pseudo-likelihood function을 극대화하는 것으로 얻어진 값이다.

$$ \[\begin{align*} q(\theta' \vert \theta , y) &= q(\theta' \vert \theta ) , q(\theta \vert \theta ', y') &= q(\theta \vert \theta ') \\ f(y \vert \theta , x) &= f(y \vert \hat \theta ), f(y' \vert \theta' , x) &= f(y' \vert \hat \theta ) \end{align*}\] $$

분포 \(f(x \vert \theta)\)를 auxiliary variable, 가령 normalizing constant ratio \(\dfrac {\kappa(\theta)} {\kappa(\theta')}\) 등으로 살찌워놓은 것은, 시뮬레이션 진행 과정에서 상쇄시키는 것이 가능하다.

1. generate \(\theta \sim q(\theta' \vert \theta_t)\)
2. generate exact sample \(y' \sim f(y \vert \theta')\)
3. accept \((\theta', y')\) with probability \(\min (1, r)\), \(r=\dfrac {f(x \vert \theta') f(\theta') f(y' \vert \hat \theta') \ast q(\theta_t \vert \theta') q(y \vert \theta_t)} {f(x \vert \theta_t) f(\theta_t) f(y \vert \hat \theta) \ast q(\theta'\vert \theta_t) q(y' \vert \theta')}\).
   • 채택된다면, set \((\theta_{t+1}, y_{t+1}) = (\theta', y' )\).
   • o.w., \((\theta_{t+1}, y_{t+1}) = (\theta_t, y'_t)\).

4.4.3.3 Exchange Algorithm

Møller’s 알고리즘을 parallel tempering 개념을 도입하여 개선.

1. propose candidate point $’ $ proposal distribution \(q(\theta' \vert \theta, x)\).
2. propose auxiliary variable \(y \sim\) perfect sampler \(f(y \vert \theta ' )\).
3. accept $’ $ with probability \(\min \{ 1, r(\theta, \theta' \vert x) \}\).
   • \(r(\theta, \theta' \vert x) = \dfrac{\pi(\theta') }{\pi(\theta) } \ast \dfrac{f(x \vert \theta ') }{f(x \vert \theta)} \ast \dfrac{f(y \vert \theta) }{f(y \vert \theta')} \ast \dfrac{f(\theta \; \vert \theta', x) }{f(\theta' \vert \theta, x)}\)

The exchange algorithm can be viewed as an auxiliary variable MCMC algorithm with the proposal distribution being augmented, for which the proposal distribution can be written as

$

\[\begin{align} T \left( \theta \rightarrow (\theta' , y) \right) &= q(\theta ' \vert \theta) f(y \vert \theta ') \\ T \left( \theta' \rightarrow (\theta , y) \right) &= q(\theta \vert \theta ' ) f(y \vert \theta) \end{align}\] $

This simply validates the algorithm, following the arguments for auxiliary variable Markov chains.

The exchange algorithm generally improves the performance of the Møller algorithm, as it avoids an initial estimation step (for \(\theta\)) that required by the Møller.

Although the Møller’s and exchange algorithms work well for some discrete models, such as the Ising and autologistic models, they cannot be applied to many other models for which perfect sampling is not available.

Even for the Ising and autologistic models, perfect sampling may be very expensive when the temperature is near or below the critical point.

---

4.4.3.4 Adaptive Exchange Algorithm

Object: An adaptive exchange algorithm (AEX) is an adaptive Monte Carlo version of the exchange algorithm, where the auxiliary variables are generated via an importance sampling procedure from a Markov chain running in parallel.

• Advantage
  • Removes the requirement of perfect sampling
  • Overcomes its theoretical difficulty caused by inconvergence of finite MCMC runs

AEX consists of two chains running in parallel.

The first chain is auxiliary, which is run in the space \({\mathcal{x}}\) with an aim to draw samples from a family of distributions \(f(X \vert \theta^{(1)}), \; \; \cdots, \; \; f(X \vert \theta^{(m)})\) for a set of pre-specified parameter values \(\theta^{(1)}, \; \cdots, \; \theta^{(m)}\).

The second chain is the target chain, which is run in the space \(\theta\) with an aim to draw samples from the target posterior \(\pi(\theta \vert y)\). For a candidate point \(\theta'\), the auxiliary variable \(x\) is resampled from the past samples of the auxiliary chain via an importance sampling procedure.

Assume that the neighboring distributions \(f(X \vert \theta^{(i)})\)’s have a reasonable overlap and the set \(\left \{ \theta^{(1)}, \; \cdots, \; \theta^{(m)} \right \}\) has covered the major part of the support of \(\pi (\theta \vert y)\).

• ALGORITHM: PART 1 - (Auxiliary Chain) Auxiliary Sample Collection via SAMC

1. (Sampling) Choose to update \(\vartheta\) or \(\pmb z_t \vert \vartheta\) with pre-specified probabilities, e.g., \(0.75\) for updating \(\vartheta\) and \(0.25\) for updating \(z_t\).

   1.-a Update \(\vartheta_{t}\) : Draw \(\vartheta '\) from the set \(\left \{ \theta^{(1)}, \; \cdots, \; \theta^{(m)} \right \}\) according to a proposal distribution \(T_1 ( \; \cdot \; \vert \vartheta_{t})\), set \((\vartheta_{t+1}, \pmb z_{t+1}) = (\vartheta ' , \pmb z_{t+1} )\) with probability \(\min \left\{ 1, \; \; \dfrac{\omega_t^{J(\vartheta_t)}}{\omega_t^{J(\vartheta ')}} \ast \dfrac {\varphi (\pmb z_{t} \vert \vartheta ')} {\varphi (\pmb z_{t} \vert \vartheta_{t})} \ast \dfrac{T_1 (\vartheta_{t} \vert \vartheta ' )}{T_1 (\vartheta ' \vert \vartheta_{t} )} \right\}\), and set \((\vartheta_{t+1}, \pmb z_{t+1}) = (\vartheta_{t}, \pmb z_t)\) with remaining probability, where \(J(\vartheta_t)\) denotes the index of \(\vartheta_t\), i.e., \(J(\vartheta_t) = j\) if \(\vartheta_t = \theta_i^{(k)}\) and \(\varphi(\pmb z \vert \vartheta)\) is an unnormalized density of \(f(\pmb z \vert \vartheta)\).

   1.-b. Update \(\pmb z_t\) : Draw \(\pmb z '\) according to a proposal distribution \(T_2 ( \; \cdot \; \vert \pmb z_t)\), set \((\pmb z_{t+1} , \vartheta_{t+1}) = (\pmb z ' , \vartheta_{t})\) with probability \(\min \left\{ 1, \; \; \dfrac {\varphi (\pmb z ' \vert \vartheta_{t})} {\varphi (\pmb z_{t} \vert \vartheta_{t})} \ast \dfrac{T_2 (\pmb z_{t} \vert \pmb z ' )}{T_2 (\pmb z ' \vert \pmb z_{t} )} \right\}\), and set \((\pmb z_{t+1} , \vartheta_{t+1}) = (\pmb z_t , \vartheta_{t})\)

2. (Abundance Factor Updating) Set

\[ \log (\omega_{t+0.5}^{(j)}) =\log (\omega_{t}^{(j)}) + a_{t+1} (e_{t+1, \; j} - p_j), \; \; \; \; \; j=1, \cdots, m \]

where \(e_{t+1, \; j} = \begin{cases} 1 & && \vartheta^{t+1} = \theta^{(j)} \\ 0 & && o.w. \end{cases}\).

If \(\omega_{t+0.5}^{(j)} \in \mathcal{K}_{\varsigma_t}\), set \((\omega_{t+1}, \pmb z_{t+1}) = (\omega_{t+0.5}, \pmb z_{t+1})\) and \(\varsigma_{t+1} = \varsigma_t\).

o.w., set \((\omega_{t+1}, \pmb z_{t+1}) = \mathbb{T}(\omega_{t}, \pmb z_{t})\) and \(\varsigma_{t+1} = \varsigma_t + 1\).

3. (Auxiliary Sample Collection) Append the sample \(\left(\pmb z_{t+1} , \vartheta_{t+1}, \omega_{t+1}^{J(\vartheta_{t+1}} \right)\) to the collection \(S_t\). Denote the new collection by \(S_{t+1}\) which is set by \(S_{t+1} = S_t \cup \left\{ \left(\pmb z_{t+1} , \vartheta_{t+1}, \omega_{t+1}^{J(\vartheta_{t+1}} \right) \right\}\).

• ALGORITHM: PART 2 - (Target Chain) Adaptive Exchange Sampler

1. (Proposal) Propose a candidate point \(\theta '\) from a proposal distribution \(q(\theta ' \vert \theta)\)

2. (Resampling for Auxiliary Variables) Resample an auxiliary variable \(\pmb x\) from the collection \(S_{t+1}\) via a dynamic importance sampling procedure;

that is, setting \(\pmb x = \pmb z_k\) with probability

\[ P(\pmb x = \pmb z_k) \dfrac {\sum_{j=1}^{\vert S_{t+1} \vert} \omega_t^{\left( J(\vartheta_j) \right)} \tfrac{\varphi(\pmb z_k \vert \theta ' )} {\varphi(\pmb z_k \vert \vartheta_j ' )} \ast I(\pmb z_j = \pmb z_k )} {\sum_{j=1}^{\vert S_{t+1} \vert} \omega_t^{\left( J(\vartheta_j) \right)} \tfrac{\varphi(\pmb z_k \vert \theta ' )} {\varphi(\pmb z_k \vert \vartheta_j ' )}} \]

• \(\left(\pmb z_{j} , \vartheta_{j}, \omega_{t}^{J(\vartheta_{j}} \right)\) denotes the \(j\)-th element of the set \(S_{t+1}\).
• \(\vert S_{t+1} \vert\)는 \(S_{t+1}\)의 size.

3. (Exchange Algorithm) Set \(\theta_{t+1} = \theta '\) with the probability \(\alpha(\theta_t , \pmb x, \theta' )\), and \(\theta_{t+1} = \theta_{t}\) with probability \(1-\alpha(\theta_t , \pmb x, \theta' )\).

\[ \alpha(\theta_t , \pmb x, \theta' ) = \dfrac {\pi(\theta ' )} {\pi(\theta_t )} \dfrac {\varphi(\pmb y \vert \theta ' )} {\varphi(\pmb y \vert \theta_t )} \dfrac {q(\theta_t \vert \theta ' )} {q(\theta' \vert \theta_t )} \dfrac {\varphi(\pmb x \vert \theta_t )} {\varphi(\pmb x \vert \theta ' )} \]

• Why this algorithm is adaptive?

Since the underlying true proposal distribution for generating auxiliary variables in part II is changing from iteration to iteration, the new algorithm falls into the class of adaptive MCMC algorithms (for which the proposal distribution is changing from iteration to iteration).