読者です 読者をやめる 読者になる 読者になる

Data Science by R and Python

統計学を、広く、深く、わかりやすく。

パッケージを使わないで、一般化線形混合モデルのMCMCアルゴリズムを1から作る.

こんにちは.
金曜日の夜になり、激しめの睡魔に襲われております.

先日のこちらの記事で公開したスライドの後半にあるシミュレーションで、地域差を考慮したPoisson - Normalモデルを構築しているのですが、そのコードを載せておきます。tomoshige-n.hatenablog.com

MCMCを自分で設計

最近だと、MCMCはStanだったり、BUGSだったり、色々と便利なツールもできてきていて、パッケージを用いて推定している人も多いみたいです。が、僕は、なんというか、中で何してるのかが気になって、提案分布の設計とか実際のところ、どうしてるのかとか、ステップ幅をどうやって決めてるのかとか、自動的にやられてしまうのはどうも気持ち悪いので、自分で書いてみましょうということになってしまいがちなんですよね。

ということで、ポアソン分布で、個体差を考慮できるようなMCMCアルゴリズムを書いてみました.ちなみに、大変申し訳ないのですが、こちらの記事のMCMCのコードには、いくつか間違いがありますので(そのまま載せてはおきますが)、、ご注意ください.

tomoshige-n.hatenablog.com

書いたRのコード

まぁ、提案分布のステップ幅などは、きちんと調節しなくては行けないですし、収束の判定をするようには作っていないので、、、申し訳ないですが、その辺が必要な方はコードを書き加えていただければ...(追記)「パッケージ使わないで、、、」とかいいつつ、truncated normalとかはパッケージ使って計算してます.許してください.

#poisson regression

set.seed(20150513) #wed seminar date
library(dummies)
library(truncnorm)
library(MASS)
#simulation
a = 400
b = 400
c = 400
x = data.frame(
	x1 = c(runif(a,-2,2), runif(b,-1.8,2.2),runif(c,-3,2.4)),
	#x2 = c(rnorm(a,-0.5,1), rnorm(b,0,1),rnorm(c,0,1))
	x2 = c(rep(1,a), rep(0,b), rep(0,c))
)
x = as.matrix(x)
z_fac = c(rep("A",a),rep("B",b),rep("C",c))
z_fac1 = dummy(as.factor(z_fac))
z = z_fac1[,-1]
beta = c(0.8,-0.5)
u = c(rnorm(a,0,0.2),rnorm(b,0,0.3),rnorm(c,0,0.3))
lambda = exp(x%*%beta + u) #z%*%u は正規分布からの乱数と考えればよい。
y = rpois(a+b+c,lambda)

par(mfcol=c(1,2))
plot(x[,1],y,xlab="",ylab="",col=as.factor(z_fac))
plot(x[,2],y,xlab="",ylab="",col=as.factor(z_fac))

#個体差なし
y_norandom = rpois(a+b+c,exp(x%*%beta))
par(mfcol=c(1,2))
plot(x[,1],y_norandom,xlab="",ylab="",col=as.factor(z_fac))
plot(x[,2],y_norandom,xlab="",ylab="",col=as.factor(z_fac))


############################
#mcmc calculation from here#
############################

#1.proposal distribution

## regression coefficient parameter
betaprop = function(beta_before,step_width,step_width2,sigma,p,q){
	n = length(beta_before)
	#mu = c(beta_before[seq(1,p,1)],rep(0,q))
	mu = beta_before
	Cov = diag(c(rep(step_width,p),rep(step_width2,q)))
	sample = mvrnorm(1, mu = mu, Sigma = Cov)
	return(sample)
}

## area dispersion parameter
sigmaprop = function(sigma,sd_trunc,q){
	mu = sigma
	sample = rep(0,q)
	for(i in 1:q){
		sample[i] = rtruncnorm(1, a=0.05, b=Inf, mean = mu[i], sd = sd_trunc)
	}
	return(sample)
}

## to calculate proposal distribution probability (q(x_{n} | x_{n-1}))
prop_lograte = function(beta_after,beta_before,sigma_after,sigma_before,p,q,sd_trunc,step_width2){
	#beta_log = sum( - (step_width2/2)*((beta_after[seq(p+1,p+q,1)])**2))
	sigma_log = sum(
		-(1/2)*((sigma_after-sigma_before)/sd_trunc)**2 - log(1 - pnorm((0.05-sigma_before)/sd_trunc))
	)
	return(sigma_log)
}

#2. posterior

prior_logsigma = function(sigma,theta,k){
	result = sum((k-1)*log(sigma) - (1/theta)*(sigma))
	return(result)
}

#mixed effectの所のlogは計算する必要がある。
prior_logmixbeta = function(beta,sigma){
	#betaとsigmaは次元分あるので、最後にsumを取る
	result = sum(-(1/2)*((beta**2)*sigma) + (1/2)*log(sigma))
	return(result)
}

#improperな事前分布は収束しないこともあるので、betaのpriorを分散1000の正規分布にしてみた.
prior_beta = function(beta){
	#betaとsigmaは次元分あるので、最後にsumを取る
	result = sum(-(1/(2*1000))*((beta**2)) - (1/2)*log(10))
	return(result)
}

#正規化定数なしの事後分布の確率のlog
log_post = function(y,x,z,beta,sigma,theta,k){
	p = dim(x)[2]
	q = dim(z)[2]
	beta1 = beta[seq(1,p,1)]
	beta2 = beta[seq(p+1,p+q,1)]
	log_lambda = x%*%beta1 + z%*%beta2
	likeli = sum(log_lambda*y) - sum(exp(log_lambda))
	#prior calculation
	#beta given sigma
	beta_given_sigma_prior = prior_logmixbeta(beta2,sigma)
	#sigma
	prior_logsigma = prior_logsigma(sigma,theta,k)
	#beta
	#improperな事前分布を避ける
	beta_prior = prior_beta(beta1)
	result = likeli + beta_given_sigma_prior + prior_logsigma + beta_prior
	return(result)
}

#アルゴリズム
mcmc_pois = function(y,x,z,intercept = F,iter=1000,burn_in = 500,step_width=0.01,step_width2=0.01,theta=0.01,sd_trunc=0.01,k=0.01){
	p = dim(x)[2]
	q = dim(z)[2]
	initial_beta = c(rnorm(p),rep(0,q))
	initial_sigma = runif(q,0.5,1)
	beta.r = matrix(0,nrow=p+q,ncol=iter)
	sigma.r = matrix(0,nrow=q,ncol=iter) 
	log_likeli = rep(0,iter)
	accept = 0
	#初期値の設定
	beta.r[,1] = initial_beta
	sigma.r[,1] = initial_sigma
	log_likeli[1] = log_post(y,x,z,initial_beta,initial_sigma,theta,k)
	accept.r = rep(0,iter)
	alpha.r = rep(0,iter)
	for(i in 2:iter){
		#sigmaの候補を生成
		candidate_sigma = sigmaprop(sigma.r[,i-1],sd_trunc,q)
		#betaの候補を生成
		candidate_beta = betaprop(beta.r[,i-1],step_width,step_width2,candidate_sigma,p,q)
		#この値に基づくlog-postを計算する
		log_likeli_candidate = log_post(y,x,z,candidate_beta,candidate_sigma,theta,k)
		#棄却率alphaの計算
		alpha = exp(
			log_likeli_candidate - 
			log_likeli[i-1] - 
			prop_lograte(candidate_beta,beta.r[,i-1],candidate_sigma,sigma.r[,i-1],p,q,sd_trunc,step_width2) +
			prop_lograte(beta.r[,i-1],candidate_beta,sigma.r[,i-1],candidate_sigma,p,q,sd_trunc,step_width2)
		)
		#alphaに基づく条件分岐
		alpha.r[i] = alpha
		if(runif(1) < alpha){
			beta.r[,i] = candidate_beta
			sigma.r[,i] = candidate_sigma
			log_likeli[i] = log_likeli_candidate
			accept.r[i] = 1
		}
		else{
			beta.r[,i] = beta.r[,i-1]
			sigma.r[,i] = sigma.r[,i-1]
			log_likeli[i] = log_likeli[i-1]
			accept.r[i] = 0
		}
	}
	print_obj = sum(accept.r)*100/iter 
	print(paste("accept rate is ",print_obj,"%",sep=""))
	return(list(
		sample_beta = beta.r,
		burn_in_beta = beta.r[,-c(1:burn_in)],
		sample_sigma = sigma.r,
		burn_in_sigma = sigma.r[,-c(1:burn_in)],
		accept=accept.r,
		alpha = alpha.r,
		ll = log_likeli))
}

result = mcmc_pois(y,x,z,intercept = F,iter=1000000,burn_in = 200000,step_width=0.001,step_width2=0.001,k=0.01,theta=(1/0.01),sd_trunc=0.7)

あとは、自己相関とか、収束先の分布などをプロットするために、以下のコードを使えばいいと思います。

#サンプルの自己相関の確認
par(mfcol=c(2,1))
acf(result$burn_in_beta[1,], type = "correlation", plot = TRUE, lag.max=400)
acf(result$burn_in_beta[2,], type = "correlation", plot = TRUE, lag.max=400)
par(mfcol=c(2,1))
acf(result$burn_in_beta[3,], type = "correlation", plot = TRUE, lag.max=400)
acf(result$burn_in_beta[4,], type = "correlation", plot = TRUE, lag.max=400)
par(mfcol=c(2,1))
acf(1/result$burn_in_sigma[1,], type = "correlation", plot = TRUE, lag.max=800)
acf(1/result$burn_in_sigma[2,], type = "correlation", plot = TRUE, lag.max=800)

#自己相関が小さくなる適切な間隔毎にデータを引っこ抜く
sample_thin = function(x,lag){
	n = dim(x)[2]
	num = seq(1,n,lag)
	return(x[,num])
}

#上記の操作で得られるMCMCのサンプル
mcmc_sample_beta = sample_thin(result$burn_in_beta,400)
mcmc_sample_sigma = sample_thin(result$burn_in_sigma,400)

#事後分布
par(mfrow=c(2,2))
hist(mcmc_sample_beta[1,],xlab="",ylab="",main="beta1",breaks=15) #0.8
hist(mcmc_sample_beta[2,],xlab="",ylab="",main="beta2",breaks=15) #-0.5
hist(sqrt(1/mcmc_sample_sigma[1,]),breaks=50,main="sigmaB") #発生元は0.3
hist(sqrt(1/mcmc_sample_sigma[2,]),breaks=50,main="sigmaC") #発生元は0.3

あと、やはり予測値というか、値の信頼区間内にデータがどの程度入るのかが気になるので、予測値の90%信頼区間書いてみます。ここは、もっと効率よくできるはずですが、、、やる気スイッチが切れたのでこれで勘弁してください。

#上のMCMCサンプルデータを一つのマトリックスにまとめる
prdt_para = data.frame(
	beta1 = mcmc_sample_beta[1,],
	beta2 = mcmc_sample_beta[2,],
	sigma1 = 1/mcmc_sample_sigma[1,],
	sigma2 = 1/mcmc_sample_sigma[2,]
)

#推定値を保存するところ
prdt_para = as.matrix(prdt_para)
p = dim(x)[2]
q = dim(z)[2]
n = length(y)
r = dim(prdt_para)[1]
result2 = matrix(0,n,r)
result3 = matrix(0,n,r)

#各点の結果を予測する
for(i in 1:n){
	result_ind = rep(0,r) #各地点の結果を保存する
	for(j in 1:r){
		fix_eff = prdt_para[j,seq(1,p,1)]
		#rand_eff = mvrnorm(1,rep(0,q),diag(prdt_para[j,seq(p+1,p+q,1)],q))
		rand_eff = mvrnorm(1,rep(0,q),diag(c(0.26,0.28),q))
		mu = exp(t(x[i,]) %*% fix_eff + t(z[i,]) %*% rand_eff)
		result2[i,j] = mu
		result3[i,j] = rpois(1,mu)
	}	
}

#上側97.5%点
up95 = function(x){
	quantile(x,0.975)
}

#下側97.5%点
bo05 = function(x){
	quantile(x,0.025)
}

#ここから、図を描く作業
upper_line = apply(result3,1,up95)
bottom_line = apply(result3,1,bo05)
predicted_mu = apply(result2,1,mean)
mu_ord = order(predicted_mu)

plot(predicted_mu,y,col = as.factor(z_fac),pch=16,ylim=c(0,25))
abline(a=0,b=1)
points(predicted_mu[mu_ord],bottom_line[mu_ord],type="l",lwd=0.7)
points(predicted_mu[mu_ord],upper_line[mu_ord],type="l",lwd=0.7)
#upper line
loess_up = loess(upper_line[mu_ord] ~ predicted_mu[mu_ord])
xl = seq(0,20,0.01)
lines(xl, predict(loess_up,xl), col='red', lwd=2)
#bottom line
loess_bo = loess(bottom_line[mu_ord] ~ predicted_mu[mu_ord])
xl = seq(0,20,0.01)
lines(xl, predict(loess_bo,xl), col='red', lwd=2)

ちなみに、信頼区間的なものを書いておりますが、ポアソン分布は「左右非対称」なので、Highest posterior region(もっとも事後分布の値が大きい95%の区間)を求めるのは、難しいので、上記のような処理に落ち着けました...

結果でこんな図が書けます.

事後分布

f:id:tomoshige_n:20150515194554p:plain

信頼区間とか

f:id:tomoshige_n:20150515195036p:plain

まとめ

意外に自分でMCMCを書くのは楽しかったです.特に、計算を速くさせるための工夫とか、やっぱり大事だなと。特に、積はlogとれば和の計算というのは、数値計算では有効だなーと思いました。特にガンマ関数、階乗の計算はさせるべきじゃありません...MCMCのRのコードを書きながら、昔、MCMCについて素人の頃に、全然MCMCのコードを書いてる記事が見つからなくて困ったのを想い出しました。そんな人の少しでも助けになればと思いながら、まとめとさせてもらいます。

おすすめ文献

参考にした文献が今回はないので、読んでいる本を載せておきます.僕の研究分野の本です.

Hierarchical Modeling and Analysis for Spatial Data, Second Edition (Chapman & Hall/CRC Monographs on Statistics & Applied Probability)

Hierarchical Modeling and Analysis for Spatial Data, Second Edition (Chapman & Hall/CRC Monographs on Statistics & Applied Probability)