JW Blog
Gibbs Sampler

Gibbs Sampler

Gibbs Sampler

๊ธฐ๋ณธ ์›๋ฆฌ

(์ถ”ํ›„ ์—…๋ฐ์ดํŠธ)

์˜ˆ์‹œ

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

## Data load 
data = dget('https://www2.stat.duke.edu/~pdh10/FCBS/Inline/yX.o2uptake')
y = data[,1]
X = data[,-1]
inv = solve

## set prior
g = length(y)
nu0 = 1
s20 = summary(lm(y~-1+X))$sigma^2
n = length(y)
p = ncol(X)

## setup
S = 1000
set.seed(2021)
BETA = matrix(NA, nrow=S, ncol=p)
sigma2 = matrix(NA, nrow=S, ncol=1)
BETA[1,] = inv(t(X) %*% X) %*% t(X) %*% y
sigma2[1,] = s20

## gibbs sampling
nun = nu0 + n
betan = (g/(g+1)) * inv(t(X) %*% X) %*% t(X) %*% y
for(s in 2:S){
  s2n = nu0*s20 + t(y-X%*%BETA[s-1,]) %*% (y-X%*%BETA[s-1,])
  sigma2[s,] = 1/rgamma(1, shape=nun/2, rate=s2n/2)
  
  Sigman = (g/(g+1)) * sigma2[s,] * inv(t(X) %*% X)
  BETA[s,] = MASS::mvrnorm(n=1, betan, Sigman)
}

## graph
colnames(BETA) = colnames(X)
gather(as.data.frame(BETA)) %>% 
  ggplot(aes(y=value, fill=key)) + geom_histogram() +
  coord_flip() + facet_wrap(~key, scales='free_x') +
  ggtitle('Posterior samples of Beta') +
  theme(legend.position = 'None')

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

์ข…๋ฅ˜

1. Systematic Sweep Gibbs Sampler

๋ชจ๋“  ๋ชจ์ˆ˜์— ๋Œ€ํ•ด ์ƒ˜ํ”Œ๋ง์„ ๋ฐ˜๋ณตํ•˜์—ฌ ์—…๋ฐ์ดํŠธํ•˜๋Š” ๋ฐฉ๋ฒ•

2. Random Sweep Gibbs Sampler

๋ฌด์ž‘์œ„๋กœ ๋ชจ์ˆ˜๋ฅผ ๋ฝ‘์•„์„œ ์—…๋ฐ์ดํŠธํ•˜๋Š” ๋ฐฉ๋ฒ•

3. Grouped Gibbs Sampler

์—ฌ๋Ÿฌ ๊ฐœ์˜ ๋ชจ์ˆ˜๋ฅผ ํ•œ๋ฒˆ์— ๋ฝ‘์•„์„œ ์—…๋ฐ์ดํŠธํ•˜๋Š” ๋ฐฉ๋ฒ•
์ด๋•Œ ๋ชจ์ˆ˜ ๊ฐ„์˜ ์—ฐ๊ด€์„ฑ์ด ์ƒ๊ธด๋‹ค.

4. Collapsed Gibbs Sampler

์ ๋ถ„์„ ํ†ตํ•ด์„œ ํŠน์ • ๋ชจ์ˆ˜์™€ ๋…๋ฆฝ์ ์ธ ๋ถ„ํฌ๋ฅผ ๊ณ„์‚ฐํ•œ ํ›„, ์ƒ˜ํ”Œ๋งํ•˜์—ฌ ์—…๋ฐ์ดํŠธํ•˜๋Š” ๋ฐฉ๋ฒ•

์ฐธ๊ณ ์‚ฌ์ดํŠธ

[1] https://niceguy1575.tistory.com/entry/%EB%B2%A0%EC%9D%B4%EC%A7%80%EC%95%88-%ED%86%B5%EA%B3%84%ED%95%99-4-Gibbs-Sampler%EA%B9%81%EC%8A%A4-%EC%83%98%ED%94%8C%EB%9F%AC-%EC%9D%98-%EC%A2%85%EB%A5%98%EC%99%80-%EC%84%B1%EC%A7%88

๋ชฉ์ฐจ