KDE, Mixture 밀도 추정 - 기니피그의 치아길이
Kong, Seokkyu
2022-09-06
Summary
본 내용은 코세라 BAYES3 과정의 Lesson 6.2 의 코드를 이용해서, 기니피그의 치아 길이에 대한 밀도 추정을 진행합니다. 밀도 추정은 KDE 와 Mixture 모델로 진행하며, Mixture 모델은 3개의 location normal distribution 으로 구성합니다.
Mixture 모델의 파라미터 추정은 EM, MCMC 를 이용합니다.
목표: 치아 길이의 분포가 multi modal 인지 확인하고, 얼마나 많은 mode 가 존재하는지 이해하고, 밀도 추정량을 서로 비교합니다.
데이터셋 설명
?ToothGrowth
반응은 60마리의 기니피그에서 상아모세포(치아 성장을 담당하는 세포)의 길이입니다. 각 동물은 오렌지 주스 또는 아스코르브산(비타민 C의 한 형태이며 VC로 코딩됨)의 두 가지 전달 방법 중 하나로 비타민 C의 세 가지 용량 수준(0.5, 1, 2mg/day) 중 하나를 받았습니다.
분석하고자 하는 대상은 len 반응 변수이며, 3가지 용량(dose) 에 따라 Mixture 모델은 3개의 성분으로 구성합니다.
데이터 셋 확인
dose 에 따라 len 의 분포를 boxplot 으로 확인합니다.
head(ToothGrowth)## len supp dose
## 1 4.2 VC 0.5
## 2 11.5 VC 0.5
## 3 7.3 VC 0.5
## 4 5.8 VC 0.5
## 5 6.4 VC 0.5
## 6 10.0 VC 0.5str(ToothGrowth)## 'data.frame': 60 obs. of 3 variables:
## $ len : num 4.2 11.5 7.3 5.8 6.4 10 11.2 11.2 5.2 7 ...
## $ supp: Factor w/ 2 levels "OJ","VC": 2 2 2 2 2 2 2 2 2 2 ...
## $ dose: num 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...boxplot(len ~ dose, data = ToothGrowth)
밀도 추정
KDE 의 density 구하는 함수를 기본으로만 사용하면, mode 를 파악하기 어렵다. bandwidth 를 조절해서 mode 가 캡쳐되도록 조정했다.
# rm(list=ls())
### Loading data and setting up global variables
# 데이터 로딩 부분
library(MASS)
library(MCMCpack)
# 혼합모델 내의 성분 개수
KK = 3 # 문서 설명 및 데이터의 사전 탐색을 통해서 3개의 mode 가 있을 것으로 추정함
x = ToothGrowth$len
n = length(x)
set.seed(2022)
# ==============================================================================
### First, compute the "Maximum Likelihood" density estimate associated
# with a location mixture of 6 Gaussian distributions using the EM algorithm
# ==============================================================================
## Initialize the parameters
# 가중치, 평균, 표준편차에 대한 초기 추측값 제공
w = rep(1,KK)/KK
mu = rnorm(KK, mean(x), sd(x))
sigma = sd(x)/KK
# 알고리즘 수렴을 제어하는 파라미터 초기화화
epsilon = 0.000001
s = 0
sw = FALSE
KL = -Inf
KL.out = NULL
while(!sw){
## E step
# ------------------
v = array(0, dim=c(n,KK))
for(k in 1:KK){
v[,k] = log(w[k]) + dnorm(x, mu[k], sigma,log=TRUE)
}
for(i in 1:n){
v[i,] = exp(v[i,] - max(v[i,]))/sum(exp(v[i,] - max(v[i,])))
}
## M step
# ------------------
# Weights
w = apply(v,2,mean)
mu = rep(0, KK)
for(k in 1:KK){
for(i in 1:n){
mu[k] = mu[k] + v[i,k]*x[i]
}
mu[k] = mu[k]/sum(v[,k])
}
# Standard deviations
sigma = 0
for(i in 1:n){
for(k in 1:KK){
sigma = sigma + v[i,k]*(x[i] - mu[k])^2
}
}
sigma = sqrt(sigma/sum(v))
##Check convergence
KLn = 0
for(i in 1:n){
for(k in 1:KK){
KLn = KLn + v[i,k]*(log(w[k]) + dnorm(x[i], mu[k], sigma, log=TRUE))
}
}
if(abs(KLn-KL)/abs(KLn)<epsilon){
sw=TRUE
}
KL = KLn
KL.out = c(KL.out, KL)
s = s + 1
#print(paste(s, KLn))
}
## 서로 다른 location 에서, 밀도 추정량 density estimator 를 계산한다.
## XXX: xx 범위 계산이 중요하다.
xx = seq(1, 40,length=40)
nxx = length(xx)
density.EM = rep(0, nxx)
for(s in 1:nxx){
for(k in 1:KK){
density.EM[s] = density.EM[s] + w[k]*dnorm(xx[s], mu[k], sigma)
}
}
# ==============================================================================
# 두 번째 방식으로 베이지안 커널 사용
# ==============================================================================
## Priors set up using an "empirical Bayes" approach
aa = rep(1,KK) # 가중치에 대한 사전분포로 uniform 을 사용한다.
eta = mean(x) # 평균들에 대한 평균으로 데이터의 평균을 사용한다. (경험적 기반하에서)
# 평균들에 대한 분산으로, 데이터의 분산을 사용한다.
tau = sqrt(var(x))
# 커널의 분산에 대한 사전분포 제공. 역감마 분포를 사전분포로 제공한다.
# shape 파라미터 2 를 사용한다는 것은, 많은 집중도를 갖지 않는다.
# 거의 무한한 분산을 가지는거지.
dd = 2
qq = var(x)/KK # 중심화 작업
## Initialize the parameters
# 사후분포 샘플들을 저장할 변수들 지정
w = rep(1,KK)/KK
mu = rnorm(KK, mean(x), sd(x))
sigma = sd(x)/KK
cc = sample(1:KK, n, replace=T, prob=w)
## Number of iterations of the sampler
rrr = 12000 # 전체 샘플 개수
burn = 2000 # 번인으로 버리는 샘플
## Storing the samples
cc.out = array(0, dim=c(rrr, n))
w.out = array(0, dim=c(rrr, KK))
mu.out = array(0, dim=c(rrr, KK))
sigma.out = rep(0, rrr)
logpost = rep(0, rrr)
for(s in 1:rrr){
# Sample the indicators
# 지시자 샘플링
for(i in 1:n){
v = rep(0,KK)
for(k in 1:KK){
v[k] = log(w[k]) + dnorm(x[i], mu[k], sigma, log=TRUE) #Compute the log of the weights
}
v = exp(v - max(v))/sum(exp(v - max(v)))
cc[i] = sample(1:KK, 1, replace=TRUE, prob=v)
}
# Sample the weights
# 가중치 새플링
# tabulate: 0 보다 큰 정수에 대해서 각 개수가 몇개가 나타나는지 계산한다.
# tabulate(c(2,3,5))
# a <- rdirichlet(1, c(1,2,3,2,1)); sum(a)
w = as.vector(rdirichlet(1, aa + tabulate(cc, nbins=KK)))
# Sample the means
# 평균 샘플링
for(k in 1:KK){
nk = sum(cc==k)
xsumk = sum(x[cc==k])
tau2.hat = 1/(nk/sigma^2 + 1/tau^2)
mu.hat = tau2.hat*(xsumk/sigma^2 + eta/tau^2)
mu[k] = rnorm(1, mu.hat, sqrt(tau2.hat))
}
# Sample the variances
# 분산 샘플링
dd.star = dd + n/2
qq.star = qq + sum((x - mu[cc])^2)/2
sigma = sqrt(1/rgamma(1, dd.star, qq.star))
# Store samples
# 로그 사후분포 계산
cc.out[s,] = cc
w.out[s,] = w
mu.out[s,] = mu
sigma.out[s] = sigma
for(i in 1:n){
logpost[s] = logpost[s] + log(w[cc[i]]) + dnorm(x[i], mu[cc[i]], sigma, log=TRUE)
}
logpost[s] = logpost[s] + log(ddirichlet(w, aa))
for(k in 1:KK){
logpost[s] = logpost[s] + dnorm(mu[k], eta, tau, log=TRUE)
}
logpost[s] = logpost[s] + dgamma(1/sigma^2, dd, qq, log=TRUE) - 4*log(sigma)
if(s/500==floor(s/500)){
#print(paste("s =",s))
}
}
## Compute the samples of the density over a dense grid
# 샘플링 이후에 MCMC 각 반복에 대해서 밀도 계산
density.mcmc = array(0, dim=c(rrr-burn,length(xx)))
for(s in 1:(rrr-burn)){
for(k in 1:KK){
density.mcmc[s,] = density.mcmc[s,] + w.out[s+burn,k]*dnorm(xx,mu.out[s+burn,k],sigma.out[s+burn])
}
}
# 점 추정으로써, 사후분포 평균을 계산한다.
density.mcmc.m = apply(density.mcmc , 2, mean)
# -----------------------------------------------------
## 커널 밀도 추정과 비교하기 위해 플롯한다.
## Plot Bayesian estimate with pointwise credible bands along with kernel density estimate and frequentist point estimate
colscale = c("black", "blue", "red")
yy = density(x, bw=1.2) # XXX: KDE 를 위한 density 구함
# yy = density(x) # KDE 를 위한 density 구함
# polygon 플롯
density.mcmc.lq = apply(density.mcmc, 2, quantile, 0.025)
density.mcmc.uq = apply(density.mcmc, 2, quantile, 0.975)
plot(xx, density.mcmc.m, type="n",ylim=c(0,max(density.mcmc.uq)),xlab="Dose", ylab="Len")
polygon(c(xx,rev(xx)), c(density.mcmc.lq, rev(density.mcmc.uq)), col="grey", border="grey")
lines(xx, density.mcmc.m, col=colscale[1], lwd=2) # MCMC
lines(xx, density.EM, col=colscale[2], lty=2, lwd=2) # EM
lines(yy, col=colscale[3], lty=3, lwd=2) # KDE
points(x, rep(0,n)) # 데이터 표시
legend('topright', c("KDE","EM","MCMC"), col=colscale[c(3,2,1)], lty=c(3,2,1), lwd=2, bty="n")
참고로, KDE 에서 bandwidth 를 조금씩 변경했을 때의 플롯을 보자.
par(mar = c(2,1,2,1))
par(mfrow = c(5, 2))
for (i in seq(0.1, 3, length.out = 10)) {
plot(density(x, bw=i), main = paste0('bw=', round(i, 2)))
}