ポアソン回帰モデルをMCMCに書き直してみる

MCMCでポアソン回帰モデル

一般化線形モデルの１つであるポアソン回帰モデルをMCMCで書いてみました。
ただし、大事な収束判定のところとか、その辺りについては全然入っていません。
簡単なmcmcのアルゴリズムを自分で書いてみたいという人の参考になればと。

使っているのは、MHアルゴリズム（メトロポリス・ヘイスティングス）です。中での計算は、簡略化するために一度（正規化していない）事後分布を紙で計算して、必要なところだけに絞っています。

あと、回帰係数の提案分布は正規分布なので、結局MHの採択・非採択のところの計算で見事に消えてくれるので、気にしないでも良いことになっています。なので、MHの特殊な場合である「ランダムウォークメトロポリス」となっています。

僕は、昔授業でMHつくってという課題で、ネット上でRのプログラムを探して見つけられずに困ったので、そんな人のためも想って。

注）終わってから気づいたのですが、、、ベータの事前分布がimproperになっているので、必要な人はそこのところだけ、分散が非常に大きい正規分布等で書き直してください。

#poisson regression
set.seed(20150505)

#simulation
x = data.frame(x1 = rnorm(100,0,1), x2 = rnorm(100,0,1))
x = as.matrix(x)
beta = c(0.8,0.6)
lambda = exp(x%*%beta)
y = rpois(100,lambda)

par(mfcol=c(1,2))
plot(x[,1],y,xlab="",ylab="")
plot(x[,2],y,xlab="",ylab="")

prop = function(beta_before,sigma){
	n = length(beta_before)
	mu = beta_before
	Cov = matrix(c(sigma,0,0,sigma),2,2)
	sample = mvrnorm(1, mu =mu, Sigma = Cov)
	return(sample)
}


log_post = function(y,x,beta){
	log_lambda = x%*%beta
	result = sum(log_lambda*y) - sum(exp(log_lambda))
	return(result)
}

mcmc_pois = function(y,x,intercept = F,iter=10000,burn_in = 3000,step_width=0.01){
	initial = rnorm(dim(x)[2])
	beta.r = matrix(0,nrow=dim(x)[2],ncol=iter)
	log_likeli = rep(0,iter)
	accept = 0
	beta.r[,1] = initial
	log_likeli[1] = log_post(y,x,initial)
	alpha.r = rep(0,iter)
	for(i in 2:iter){
		#betaの候補を生成
		candidate = prop(beta.r[,i-1],step_width)
		#この値に基づくlog-postを計算する
		log_likeli_candidate = log_post(y,x,candidate)
		#棄却率alphaの計算
		alpha = exp(log_likeli_candidate - log_likeli[i-1])
		#alphaに基づく条件分岐
		if(runif(1) < alpha){
			beta.r[,i] = candidate
			log_likeli[i] = log_likeli_candidate
			accept = accept + 1
		}
		else{
			beta.r[,i] = beta.r[,i-1]
			log_likeli[i] = log_likeli[i-1]
		}
	}
	print(paste("accept rate is ",accept,sep=""))
	return(list(sample = beta.r,burn_in_sample = beta.r[,-c(1:3000)]))
}

result = mcmc_pois(y,x,iter=10000,burn_in=3000)

par(mfcol=c(2,2))
plot(result$burn_in_sample[1,])
hist(result$burn_in_sample[1,])
plot(result$burn_in_sample[2,])
hist(result$burn_in_sample[2,])

#分布の平均値/標準誤差(monte carlo approximation)
mean(result$burn_in_sample[1,]) #発生元は0.8
#0.7905004
mean(result$burn_in_sample[2,]) #発生元は0.6
#0.6605761
sd(result$burn_in_sample[1,])
#0.05299724
sd(result$burn_in_sample[2,])
#0.06061163


#glmの結果と比較
glm(y~.,data=data.frame(x),family="poisson")

#Coefficients:
#            Estimate Std. Error z value Pr(>|z|)    
#(Intercept) -0.03028    0.11725  -0.258    0.796    
#x1           0.80767    0.07370  10.959   <2e-16 ***
#x2           0.67309    0.07200   9.349   <2e-16 ***
#---

補足

ちなみに、こんなことしなくても、ポアソン回帰モデルの頻度論的な結果を得たい場合には「glm」コマンドで一発です。そのほか、glmmMLのパッケージを使うこともできますし、BUGS, Stanなどでももちろん実行可能です。mcmc系のアルゴリズムを実行する環境は非常に充実しているんですが、やっぱり僕は、自分で書いていろいろ試したいんですよね。最近は、adaptive mcmcをつくろうとしています（笑）

上のコードをGLMM（一般化線形混合モデル）に拡張したコードはこちらです。<a href="http://tomoshige-n.hatenablog.com/entry/2015/05/06/041539">一般化線形混合モデル(GLMM) -Poisson/NormalをMCMCで書き直す - Data Science by R and Python</a>tomoshige-n.hatenablog.com