본문으로 바로가기



Simulation with R - 2. 깁스 표집(Gibbs sampling)을 사용한 다변량 정규분포 난수 생성

category Statistics/Simulation 2017.07.25 10:53
Simulation with R - 2. 깁스 표집(Gibbs sampling)을 사용한 다변량 정규분포 난수 생성

이전 포스팅에서는 메트로폴리스-해스팅스(Metropolis-Hastings, 이하 MH) 알고리즘에 대하여 알아보았다. 오늘은 MH 알고리즘의 특별 케이스인 깁스 표집기(Gibbs sampler)을 이용하여 다변량 정규분포에서 표본을 추출해 보도록하자. 이렇게 깁스 표집을 이용하여 생성한 난수들은 타겟분포인 다변량 정규분포를 정상분포(stationary distribution)로 갖는 마르코프 체인을 형성하며, 이러한 마르코프 체인은 각각이 독립인 표본이 아닌 서로 상관성이 있는 표본이라는 점을 잊지말자.

Gibbs sampler 알고리즘

Josiah Willard Gibbs (사진출처: 위키피디아)

인터넷을 찾아보니 깁스 표집 방법은 Josiah Willard Gibbs가 생전에 발표한 것이아닌, 그의 사후 1984년에 Stuart와 Donald Geman 형제가 Gibbs distribution을 사용하여 아이디어를 처음 제안했다고 한다. 두 형제가 발표한 이 논문은 학계에 큰 영향을 미쳤는데, 이 영향력은 현재(2017년 7월)까지의 citation 수 20,485번에서도 느낄수 있다. 깁스 표집 알고리즘은 보통 표본을 추출한 목표 분포가 다변량 분포일때 주로 사용된다. 이유는 깁스 표집을 사용하면 다변량 함수에서의 난수를 각 변수에 대한 조건부 분포을 이용하여 뽑을수 있게 되는데, 이러한 각 변수에 대한 조건부 분포는 단변량(univariate) 조건부 분포이기 때문에 계산이 훨씬 쉬워진다.

만약 우리가 난수를 추출할 분포가 d개의 확률변수를 가진 분포라고 하자.

\[ \mathbf{X}=\left(X_{1},...,X_{d}\right)^{T}\in\mathbb{R}^{d} \]

그리고 이 확률벡터의 결합확률밀도 함수는 \(f\left(\mathbf{x}\right)=f\left(x_{1},...,x_{d}\right)\)이라고 하자. 또한 앞에서 정의한 확률벡터 중 i번째 변수만을 제외한 나머지 확률 변수들을 벡터로 만든 확률벡터를 다음과 같이 정의도록 하자.

\[ \mathbf{X}_{\left(-i\right)}=\left(X_{1},...,X_{i-1},X_{i+1},...,X_{d}\right)^{T}\in\mathbb{R}^{d-1} \]

따라서 i번째 확률변수의 나머지 확률벡터, \(\mathbf{X}_{\left(-i\right)}\),값이 모두 주어진 조건부 확률밀도함수는 1변량 함수가되고, \(f\left(x_{i}|\mathbf{x}_{\left(-i\right)}\right)\)로 나타낼 수 있게 된다. 이 조건부 확률밀도함수를 이용하여 깁스 표본은 다음의 과정을 거쳐 마르코프 체인을 생성한다.

  1. 초기값 벡터, \(\mathbf{X}(0)\),를 설정한다. 이 값은 임으로 정해도 된다. 왜냐하면 초기 체인의 값들의 경우 수렴하기까지 시간이 걸리므로 보통 나중에 사용할때 제외시켜(burn in이라는 용어를 사용한다.) 버리기 때문이다. 하지만 수렴분포와 멀리 떨어진 초기값의 경우 수렴까지 걸리는 시간이 길어지므로 너무 동떨어진 값을 설정하는 것은 피해야한다.

  2. 현재의 체인의 인덱스 값을 i라고 할 경우, \(\mathbf{X}\left(i\right)=\mathbf{x}\left(i\right)\), i+1의 체인을 발생시키기 위해서 다음의 과정을 반복한다.

  • 각 변수 \(j=1,...,d\) 마다 조건부 분포를 이용하여 난수를 발생시킨후, 발생된 값을 새로운 값으로 업데이트 시킨다. 따라서 업데이트된 값은 그 다음 변수를 생성할때 영향을 미치게된다. 이것을 좀 더 시각적으로 나타내어보면 다음과 같다.

\[ \begin{align*} X_{1}^{\left(1\right)} & \sim f\left(x_{1}|\mathbf{X}_{\left(-1\right)}^{\left(0\right)}=\left(x_{2}^{\left(0\right)},x_{3}^{\left(0\right)},...,x_{d}^{\left(0\right)}\right)\right)\\ X_{2}^{\left(1\right)} & \sim f\left(x_{1}|\mathbf{X}_{\left(-2\right)}=\left(x_{1}^{\left(1\right)},x_{3}^{\left(0\right)},...,x_{d}^{\left(0\right)}\right)\right)\\ X_{3}^{\left(1\right)} & \sim f\left(x_{1}|\mathbf{X}_{\left(-3\right)}=\left(x_{1}^{\left(1\right)},x_{2}^{\left(1\right)},X_{4}^{\left(0\right)},...,x_{d}^{\left(0\right)}\right)\right)\\ & \vdots\\ X_{d}^{\left(1\right)} & \sim f\left(x_{1}|\mathbf{X}_{\left(-3\right)}=\left(x_{1}^{\left(1\right)},x_{2}^{\left(1\right)},...,x_{d-1}^{\left(1\right)}\right)\right) \end{align*} \]

위와 같이 \(X_1\)에서 \(X_d\)까지를 차례로 발생시킨다. 조건부 확률밀도함수가 업데이트 되는것을 잘 보도록 하자.

  1. 2 단계를 모두 마치면 i+1 인덱스의 마르코프 체인을 발생시키는 것을 완료하게 된다. 체인이 수렴할 때까지 단계 2를 계속해서 반복한다.

다변량 정규분포 난수발생에 적용하기

위에서 설명한 깁스 표집 알고리즘을 사용하여 다변량 정규분포를 정상분포로 갖는 마르코프 체인을 발생시켜 보도록 하자. 깁스 표집기를 사용하기 위해서는 먼저 각 변수별 조건부 분포를 유도해야 할 수 있어야 하므로, 먼저 다변량 정규분포의 조건부 분포에 대하여 간단히 짚고 넘어가도록 하자. (다변량 정규분포는 통계에서 가장 중요한 분포라고 해도 무방할 정도이므로, 추후 이 내용에 관한 포스팅을 해보도록 하겠다.)

다변량 정규분포 조건부 확률분포에 대하여

\(\mathbf{X}\in \mathbb{R}^n\)이 다음과 같이 \(n\)변량 정규분포를 따른다고 하자.

\[ \mathbf{X}=\begin{bmatrix}\mathbf{X}_{1}\\ \mathbf{X}_{2} \end{bmatrix}\text{ with sizes }\begin{bmatrix}q\times1\\ (n-q)\times1 \end{bmatrix}\sim\mathcal{N}\left(\boldsymbol{\mu},\boldsymbol{\Sigma}\right) \]

위의 정규분포의 모수는 다음과 같다.

\[ \boldsymbol{\mu}=\begin{bmatrix}\boldsymbol{\mu}_{1}\\ \boldsymbol{\mu}_{2} \end{bmatrix}\text{ with sizes }\begin{bmatrix}q\times1\\ (n-q)\times1 \end{bmatrix} \]

\[ \boldsymbol{\Sigma}=\begin{bmatrix}\boldsymbol{\Sigma}_{11} & \boldsymbol{\Sigma}_{12}\\ \boldsymbol{\Sigma}_{21} & \boldsymbol{\Sigma}_{22} \end{bmatrix}\text{ with sizes }\begin{bmatrix}q\times q & q\times(n-q)\\ (n-q)\times q & (n-q)\times(n-q) \end{bmatrix} \]

즉 n의 확률벡터로 이루어진 n변량 정규분포를 2부분으로 쪼개서 나타낸 것이다. 이유는 한부분이 주어졌을때 다른 한 부분의 분포를 알아보기 위해서인데, 모든 정규분포의 조건부 분포는 다시 정규분포를 이룬다. 이런 면에서 정규분포는 이론적을 참으로 아름다운 분포이지 않나 싶다. 각설하고, \(\mathbf{X}_{2}\)가 주어졌을때, \(\mathbf{X}_{1}\)의 분포는 다음과 같다.

\[ \mathbf{X}_{1}|\mathbf{X}_{2}\sim\mathcal{n}\left(\boldsymbol{\mu^{\star}},\boldsymbol{\Sigma}^{\star}\right) \]

위의 식에서 \(\mu^\star\)\(\Sigma^\star\)는 다음과 같다.

\[ \begin{align*} \boldsymbol{\mu^{\star}} & =\boldsymbol{\mu}_{1}+\boldsymbol{\Sigma}_{12}\boldsymbol{\Sigma}_{22}^{-1}\left(\mathbf{\mathbf{x}_{2}}-\boldsymbol{\mu}_{2}\right)\\ \boldsymbol{\Sigma}^{\star} & =\boldsymbol{\Sigma_{11}}-\boldsymbol{\Sigma}_{12}\boldsymbol{\Sigma}_{22}^{-1}\boldsymbol{\Sigma}_{21} \end{align*} \]

위의 식은 2변량의 경우 다음과 같아진다. 아래의 경우 \(\rho\)는 상관계수를 의미한다.

\[ X_1\mid X_2=x_2 \ \sim\ \mathcal{N}\left(\mu_1+\frac{\sigma_1}{\sigma_2}\rho( x_2 - \mu_2),\, (1-\rho^2)\sigma_1^2\right) \]

3변량 정규분포에 적용하기

위의 식을 3변량 정규분포에 적용시켜보도록 하겠다. 보통의 수리통계 책에서는 바로 위에서 언급한 변수가 2개인 경우를 가지고 설명을 하고 있지만, 앞선 행렬식을 사용한 공식을 완전히 이해하기 위해서는 3개의 변수를 가지고 놀아봐야만 제대로 이해할 수 있기 때문이다.

우리가 3변량 정규분포에서 깁스 표집을 통하여 난수를 뽑고 싶다고 가정하자.

\[ \left(\begin{array}{c} X_{1}\\ X_{2}\\ X_{3} \end{array}\right)\sim\mathcal{N}\left(\left(\begin{array}{c} 1\\ 6\\ 7 \end{array}\right),\left(\begin{array}{ccc} 4 & 2 & 3\\ 2 & 5 & 0\\ 3 & 0 & 25 \end{array}\right)\right) \]

이 경우 \((X_2,X_3)=(x_1,x_2)\)으로 주어졌을 때의 \(X_1|(X_2,X_3)\)의 조건부 분포는 다음과 같은 평균과 분산을 갖는 정규분포를 따른다.

  • 평균

\[ \begin{align*} \boldsymbol{\mu^{\star}} & =1-\left(\begin{array}{cc} 2 & 3\end{array}\right)\left(\begin{array}{cc} 5 & 0\\ 0 & 25 \end{array}\right)^{-1}\left(\left(\begin{array}{c} x_{2}\\ x_{3} \end{array}\right)-\left(\begin{array}{c} 6\\ 7 \end{array}\right)\right)\\ & =1-\left(\begin{array}{cc} 0.4 & 0.12\end{array}\right)\left(\begin{array}{c} x_{2}-6\\ x_{3}-7 \end{array}\right)\\ & =1-0.4\left(x_{2}-6\right)-0.12\left(x_{3}-7\right) \end{align*} \]

  • 분산

\[ \begin{align*} \boldsymbol{\Sigma}^{\star} & =4-\left(\begin{array}{cc} 2 & 3\end{array}\right)\left(\begin{array}{cc} 5 & 0\\ 0 & 25 \end{array}\right)^{-1}\left(\begin{array}{c} 2\\ 3 \end{array}\right)\\ & =2.72 \end{align*} \]

자, 그렇다면 나머지 두 가지 경우, 즉, \(X_2\)\(X_3\)의 조건부 분포는 어떻게 구할까? 우리가 배운 공식은 확률벡터의 아래쪽 부분의 값이 주어졌을 경우, 위쪽 부분의 조건부 분포였는데 말이다. 이 부분은 정규분포를 바꿔서 다시 써보면 쉽다. 만약 우리가 \(X_3|(X_1, X_2)\)의 조건부 분포를 구하고 싶다면 우리가 배운 공식을 적용할 수 있도록 다음과 같이 바꿔써준다.

\[ \left(\begin{array}{c} X_{3}\\ X_{1}\\ X_{2} \end{array}\right)\sim\mathcal{N}\left(\left(\begin{array}{c} 7\\ 1\\ 6 \end{array}\right),\left(\begin{array}{ccc} 25 & 3 & 0\\ 3 & 4 & 2\\ 0 & 2 & 5 \end{array}\right)\right) \]

같은 방식으로 \(X_2|(X_1, X_3)\)의 조건부 분포를 구하기 위해서는 다음과 같은 정규분포를 이용하면 된다.

\[ \left(\begin{array}{c} X_{2}\\ X_{1}\\ X_{3} \end{array}\right)\sim\mathcal{N}\left(\left(\begin{array}{c} 6\\ 1\\ 7 \end{array}\right),\left(\begin{array}{ccc} 5 & 2 & 0\\ 2 & 4 & 3\\ 0 & 3 & 25 \end{array}\right)\right) \]

깁스 표집 적용하기

위의 방법을 사용하면 깁스 표집 방법의 준비물인 세 개의 조건부 확률분포를 구할 수 있다.

\[ \begin{align*} X_{1}|\left(X_{2},X_{3}\right) & =\left(x_{2},x_{3}\right)\\ \sim & \mathcal{N}\left(1-0.4\left(x_{2}-6\right)-0.12\left(x_{3}-7\right),2.72\right)\\ X_{2}|\left(X_{1},X_{3}\right) & =\left(x_{1},x_{3}\right)\\ \sim & \mathcal{N}\left(6+0.5495\left(x_{1}-1\right)-0.06593\left(x_{3}-7\right),3.901\right)\\ X_{3}|\left(X_{1},X_{2}\right) & =\left(x_{1},x_{2}\right)\\ \sim & \mathcal{N}\left(7+0.9375\left(x_{1}-1\right)-0.375\left(x_{2}-6\right),22.1875\right) \end{align*} \]

R을 사용한 깁스 Sampler 구현하기

먼저 초기값을 설정해준다.

set.seed(2017)

# length of chain
N <- 10000

# chain will be N by 3 matrix
mchain <- matrix(0, nrow = N, ncol = 3)

# coefficients
mu_vec <- c(1, 6, 7)
cond_var <- c(2.72, 3.901, 22.1875)
coef_mat <- matrix(c(1, -0.4, -0.12,
                     6, 0.5459, -0.06593,
                     7, 0.9375, -0.375),
                   ncol = 3, byrow = T)

# init chain
mchain[1,] <- mu_vec

초기값을 설정하였으므로, 깁스 표집의 가장 핵심 부분을 수행하도록 하자. i는 체인에 대한 인덱스이며, j는 각 변수에 대한 인덱스로 사용하였다.

# Generate chain with Gibbs sampler
for (i in 2:N){
  given_value <- mchain[i-1,]
  for (j in 1:3){
    given_value[j] <-
      rnorm(1, mean = t(coef_mat[j,]) %*%
                 c(1, given_value[-j] - mu_vec[-j]),
               sd = sqrt(cond_var[j]))
  }
  mchain[i,] <- given_value
}

# Burn the beginning period
X <- mchain[2000:N, ]

발생시킨 체인이 수렴하는 것을 고려하여 처음 1000개를 버렸다. 이제 우리가 뽑은 표본들이 정말 우리가 목표로 설정한 정규분포를 따르는지 각 변수별 평균과 공분산 행렬을 구해보도록 하자.

colMeans(X)
## [1] 1.011891 5.995341 7.054795
cov(X)
##          [,1]      [,2]       [,3]
## [1,] 3.970512 2.4325182  2.7190872
## [2,] 2.432518 5.4524394  0.2173198
## [3,] 2.719087 0.2173198 24.5006443

Visualization

아래는 우리가 구한 체인의 마지막 3000개를 3차원 상에 표시한 것이다. \(X_2\)\(X_3\)의 경우 서로의 독립이라서 넓게 퍼져있는 반면 다른 변수들간에는 상관성이 존재하여 세밀하게 표시되는 것을 확인할 수 있다.(컴퓨터의 경우 마우스와 휠을 사용해 돌려볼 수 있습니다.)

마치며

이번 포스팅에 사용된 다변량 정규분포의 난수 생성 코드는 3변수에서, 또는 필자가 사용한 정규분포에만 적용되는 코드임을 알아야한다. 이것은 어디까지나 깁스 표집을 좀 더 잘 이해하고, 다변량 정규분포의 조건부 분포를 이해하기 위한 연습이며, 3변수를 넘어서 일반화된 평균과 분산이 주어진 다변량 정규분포의 난수 생성은 좀 더 일반화된 접근이 필요할 것이다. 아울러 좀 더 나은 접근방법을 아시는 독자분들께서는 댓글을 통하여 필자와 소통하였으면 좋겠다. 항상 새로운 더 나은 내용을 위한 제안과 코멘트는 언제나 환영이다.

Reference

[1] Rizzo, Maria L. Statistical computing with R. CRC Press, 2007.


SHARE TO



티스토리 툴바