パッケージを使わないで、一般化線形混合モデルのMCMCアルゴリズムを1から作る.
こんにちは.
金曜日の夜になり、激しめの睡魔に襲われております.
先日のこちらの記事で公開したスライドの後半にあるシミュレーションで、地域差を考慮したPoisson - Normalモデルを構築しているのですが、そのコードを載せておきます。tomoshige-n.hatenablog.com
MCMCを自分で設計
最近だと、MCMCはStanだったり、BUGSだったり、色々と便利なツールもできてきていて、パッケージを用いて推定している人も多いみたいです。が、僕は、なんというか、中で何してるのかが気になって、提案分布の設計とか実際のところ、どうしてるのかとか、ステップ幅をどうやって決めてるのかとか、自動的にやられてしまうのはどうも気持ち悪いので、自分で書いてみましょうということになってしまいがちなんですよね。
ということで、ポアソン分布で、個体差を考慮できるようなMCMCのアルゴリズムを書いてみました.ちなみに、大変申し訳ないのですが、こちらの記事のMCMCのコードには、いくつか間違いがありますので(そのまま載せてはおきますが)、、ご注意ください.
書いた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%の区間)を求めるのは、難しいので、上記のような処理に落ち着けました...
結果でこんな図が書けます.
事後分布
信頼区間とか
まとめ
意外に自分で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)](http://ecx.images-amazon.com/images/I/51BBabyaIEL._SL160_.jpg)
- 作者: Sudipto Banerjee,Bradley P. Carlin,Alan E. Gelfand
- 出版社/メーカー: Chapman and Hall/CRC
- 発売日: 2014/07/28
- メディア: ハードカバー
- この商品を含むブログを見る