Simulations

Motivation

과 같은 적분을 계산하는 것은, $h(x)$가 복잡한 함수이거나 $x$가 다차원 벡터인 경우 쉬운 일이 아니다. 이런 경우, 적분값을 근사적으로 구하기 위해, 여러가지 Simulation 방법론들이 개발되었다.
Simulation은 특히 Bayesian Inference에서 유용하다.

예를 들어, $X^n = (X_1,...,X_n)$이 observed data일 때, posterior density는 다음과 같이 구할 수 있다.

Bayesian inference에서는 posterior의 mean을 이용해 $\theta$를 추정하는 경우가 많은데, 만약 $\theta = (\theta_1, \theta_2....)$가 다차원 벡터라면, $\theta_1$에 대한 추론을 위해서는 $\theta_1$에 대한 주변사전분포, 즉 $\pi(\theta_1|X^n) = \int_{\theta_2}\int_{\theta_3}...\int_{\theta_k} \pi(\theta|X^n)$를 구한 후, 비로소 $\theta_1$에 대한 추론이 가능해 질 것이다. 이런 과정에서 요구되는 다중적분 들은 계산이 힘들거나 연산량이 많은 경우가 빈번하기에, 이러한 경우 우리는 simulation을 선호하게 되는 것이다.

Basic Monte Carlo Integration.

가장 간단한 형태의 simulation은 Monte-carlo 적분이다.

의 적분을 해야 하는 상황을 생각하자.
만약 h(x)가 복잡한 형태라서, I에 대하여 알려진 closed form의 적분형태가 없어서, 직접적으로 적분값을 구할 수 없다. 이러한 상황에서, 우리는 다음의 사실을 발견 할 수 있다.

만약 우리가 f라는 밀도함수에서 observations 들을 generate 할 수 있다면, Strong Law of Large Numbers 에 의해,

임을 알 수 있으므로,
$\frac{1}{n}\sum_i w(x_i)$를 $I$에 대한 Estimate로 활용 할 수 있다.

Example(Bayesian inference for two binomials)

이 때, prior $\pi(p_1, p_2) = \pi(p_1)\pi(p_2) = 1$이라고 설정하자.

Bayesian method를 통해서 $\delta = p_2 - p_1$을 추정해보고자 한다.

수식적으로 풀면 다음과 같다.

$f(p_1, p_2 | X,Y) \propto p_1^X(1-p_1)^{n-X}p_2^Y(1-p_2)^{m-Y}$
$\delta$에 대한 posterior를 구하기 위해선,
$P(\delta \le c|X,Y) = \int_Af(p_1, p_2 | X,Y)$ where $A =\{p_1, p_2| p_2 - p_1 \le c\}$ 를 구하고 이를 $\delta$에 대하여 미분하는 방법을 생각할 수 있으나, 이는 수치적으로 복잡한 형태를 띌 수 있다.

이러한 번거로움을 피하기 위해, 우리는 Monte-Carlo Simulation을 생각 해볼 수 있다.

$\delta$의 posterior mean은 다음과 같이 표현된다.

이 때, $p_1|X \sim Beta(X+1, n-X+1)$, $p_2|Y \sim Beta(Y+1, m-Y+1)$ 이다.

이제

로 샘플링을 하고, $\delta_i = P_2^{(i)} - P_2^{(i)}$로 두면

Monte Carlo 방법에 따라서
Bayes method에 따른 $\delta$의 estimator은

로 구할 수 있다.

Importance Sampling

몬테칼로 적분에서, 우리는 $h(x) = w(x)f(x)$로 표현되는 pdf $f$를 잡아서, $f$를 따르는 분포에서 표본을 샘플링함으로써 구하고자 하는 적분을 근사 할 수 있었다.

만약, 적당한 f를 찾기 힘들다면, 우리는 다음과 같이 생각할 수도 있다. 사실, 몬테칼로 방법과 다른것이 없다.

Sampling이 가능한 pdf $g(x)$를 잡아서,

where $Y = \frac{h(x)f(x)}{g(x)}$.
그렇다면, 몬테칼로와 같은 방법으로,

와 같이 적분의 근사값을 구할 수도 있다.

여기서 하나의 궁금증이 생길 수 있다.

$g(x)$는 무엇을 고르든 상관이 없는가?

$g(x)$를 무엇을 고르든 구해진 추정값이 참값으로의 일치성을 갖는 것은 자명하다.

다만, $g(x)$의 선택에 따라, 추정값의 표준오차가 달라지게 된다. 다음을 보자.

만약, $g(x)$가 너무 얇은 꼬리를 가진다면, w의 2차적률은 $\infty$에 가까워 질 수 있다. 즉, 우리는 꼬리가 f보다 두꺼운 g를 고려할 필요가 있다.
또한, f가 큰 지역에서 작은 g값을 가지는 g를 선택해도, 분산이 커질 우려가 존재한다.
결론적으로 말하면, 우리는 f와 비슷한 개형을 가지는 g를 선택하는 것이 유리하다.

The Metropolis-Hastings Algorithm.

이제 이 장의 Main Theme인 MCMC(Markov Chain Monte Carlo)의 한 종류인 Metropolis-Hastings algorithm을 소개한다.

MCMC의 IDEA

Stationary distribution이 $f$인 Markov Chain $X_1, X_2,...$을 고려하자.

다음이 성립한다.

특정한 조건 하에서

이는 마코프 체인의 대수의 법칙에 의해서 성립하게 된다.

이러한 아이디어를 활용하기 위하여, Staionary distribution이 f인 markov chain을 생성하는 Metropolis-Hastings Algorithm은 다음과 같다.

Metropolis-Hastings Algorithm
0. proposal distribution $q(y|x)$ 를 하나 정하자. q는 pdf를 의미한다.
1. $X_0$를 임의로 고른다
2. $X_0, X_1, ..., X_i$가 주어졌을 때, $X_{i+1}$을 다음과 같이 생성한다.

• proposal $Y \sim q(y|X_i)$을 뽑는다.
• $r = r(X_i, Y) = \min\{\frac{f(y)}{f(x)}\frac{q(x|y)}{q(y|x)}, 1\}$로 두고,
• $X_{i+1} = \begin{cases} Y & \text{if } with \; probability \;r, \\ X_i & \text{if } with\; probability\; 1-r. \end{cases}$

이에 따르면, 우선 $X_n$은 Markov Chain임을 확인 할 수 있다.

Why M-H algorithm works?

그렇다면, 왜 M-H 알고리즘이 stationary distribution f를 가지는 markov chain을 generate 시킬까?

다음의 증명을 보자

Theorem: M-H algorithm generates Markov chain with stationary distribution f

먼저 M-H algorithm으로 생성되는 Markov chain이 stationary distribution을 가지는지 부터 확인해야 하는데, 이는 M-H algorithm이 ergodic 함을 보이면 된다. (적절하게 proposal distribution q를 설정하면 이는 만족된다.)
또한, Markov Chain 이론에 따르면, Detailed Balance Condition, 즉
$p_{ij}\pi_i = p_{ji}\pi_j$
가 만족되면 markov chain은 stationary distribution $\pi$를 가진다고 했다.

WTS:M-H에 의해 생성되는 markov chain은 detailed balance를 만족한다
즉 이 경우, $f(x)p(x,y) = f(y)p(y,x)$임을 보이면 된다.이 때 x=y일때는 이 조건이 성립함이 명백하므로, $x \neq y$인 상황만 고려하자.

$p(x,y) = rq(y|x) = \min\{\frac{f(y)}{f(x)}\frac{q(x|y)}{q(y|x)}, 1\}q(y|x)$ if $x \neq y$

$\therefore f(x)p(x,y) = \min\{f(y)q(x|y), f(x)q(y|x)\} = f(y)p(y,x)$

따라서 M-H 알고리즘에 의해 생성되는 마코프 체인은 정상분포 f를 가진다.

앞의 monte carlo 방법들에 비해서 MCMC가 가지는 이점은, MCMC는 목표하는 분포 f에 가깝게 데이터를 추출할 수 있는 방법을 제시한다는 것이다.

Gibbs Sampling

깁스추출법은 M-H 알고리즘을 다차원으로 확장시킨 것을 볼 수 있다. 즉, 복잡한 다차원 분포에서 표본을 잘 추출할 수 있도록 하는 알고리즘으로, MCMC의 일종이다.

Gibbs Sampling : 2차원 변수의 경우
Iterate until convergence:

• $X_{n+1} \sim f_{X|Y}(x|Y_n)$
• $Y_{n+1} \sim f_{Y|X}(y|X_{n+1})$
  

Gibbs Sampling : 다차원 변수로의 일반화
Iterate until convergence:

• $X_{n+1} \sim f_{X|Others}(x|Others)$
• $Y_{n+1} \sim f_{Y|Others}(y|Others)$
• $Z_{n+1} \sim f_{Z|Others}(z|Others)$
  ....

NOTE: Gibbs Sampling은 M-H 알고리즘의 한 종류이다.

Suppose : proposal = $(X_{n+1},Y_{n+1}, Z, W_n, T_n...)$
current state = $(X_{n+1},Y_{n+1}, Z_n, W_n, T_n...)$
이 때 z의 분포를 $q(z) = f(z|X_{n+1}, Y_{n+1}, W_n, T_n...)$으로 두면, ($=f_Z(Z|others$))
$r((X_{n+1},Y_{n+1}, Z, W_n, T_n...), (X_{n+1},Y_{n+1}, Z_n, W_n, T_n...))$
$= \min\{\frac{f(X_{n+1},Y_{n+1}, Z, W_n, T_n...)}{f(X_{n+1},Y_{n+1}, Z_n, W_n, T_n...)} \frac{q(Z_n|Z)}{q(Z|Z_n)}, 1\}$
$= \min \{\frac{f(X_{n+1},Y_{n+1}, Z, W_n, T_n...)}{f(X_{n+1},Y_{n+1}, Z_n, W_n, T_n...)}\frac{f(Z_n|X_{n+1}, Y_{n+1}, W_n, T_n...)}{f(Z|X_{n+1}, Y_{n+1}, W_n, T_n...)}, 1\} = 1$

Gibbs Sampling의 Stationary Distribution은 f(X,Y,Z..) 인가?

근데 문제가 하나 있다. proposal을 추출하는 분포 $q(x) = f(X|others)$에서 데이터를 추출하기 어려운 상황이라면 어떻게 해야할까?

Metropolis with Gibbs 추출법은 MH algorithm과 Gibbs Sampling을 결합하여 각 변수별로 제안된 샘플을 조건부 분포에 따라 선택하는 Gibbs 샘플링을 수행하면서, Metropolis-Hastings 알고리즘을 통해 각 샘플의 수용 여부를 결정한다.

Metropolis within Gibbs : 2차원 변수

(1a) $Z \sim q(z|X_n)$으로 proposal을 추출한다.
(1b) $r = \min \{\frac{f(Z, Y_n)}{f(X_n, Y_n)}\frac{q(X_n|Z)}{q(Z|X_n)}, 1\}$
(1c) $X_{n+1} = \begin{cases} Z & \text{if } with \; probability \;r, \\ X_n & \text{if } with\; probability\; 1-r. \end{cases}$

(2a) $Z \sim \tilde{q}(z|Y_n)$으로 proposal을 추출한다.
(2b) $r = \min \{\frac{f(X_{n+1}, Z)}{f(X_{n+1}, Y_n)}\frac{\tilde{q}(Y_n|Z)}{\tilde{q}(Z|Y_n)}, 1\}$
(2c) $Y_{n+1} = \begin{cases} Z & \text{if } with \; probability \;r, \\ Y_n & \text{if } with\; probability\; 1-r. \end{cases}$

q는 sampling 가능한 분폳함수로 정하면 된다.

*즉, Gibbs Sampling은 M-H 알고리즘을 다차원 분포에대하여 적용하기 위해, 다차원 분포 샘플링 문제를 여러개의 1차원 분포 샘플링의 문제로 변환시킨 것이다.

Question

• M-H 알고리즘에서 q는 어떻게 선택해야 하는가?