階層ベイズで作成したパラメータのシミュレーションについて その1

里村先生の『マーケティング・モデル』をベースに、階層ベイズのモデルを作り、シミュレーションを行いました。

マーケティング・モデル 第2版 (Rで学ぶデータサイエンス 13)

• 作者: 里村 卓也
• 出版社/メーカー: 共立出版
• 発売日: 2015/04/09
• メディア: 単行本

コードは、教科書に載っていますし、共立出版のホームページからもダウンロードできます。

以下、里村先生のコードをベースに修正したものを載せていますが、
一部、教科書そのままだと、上手く動かない部分があるので修正しています。
（オリジナルのコードは、上記を参照してください。）

まずは、初期設定から。

## 階層ベイズ2項ロジットモデルの推定
library(bayesm)
library(psych) ## パラメータの散布図で使用
　
set.seed(1234)
　
nvar <- 2　　## ロジットモデルの説明変数の数
hh <- 1000　　## 個人数
nobs <- 10　　## 個人あたりの選択回数
nz <- 2　　　　## 個人属性の説明変数の数
　
Z <- matrix(c(rep(1,hh),runif(hh,min=-1,max=1)),hh,nz)
Delta <- matrix(c(-2,-1,1,2),nz,nvar)
iota <- matrix(1,nvar,1)
Vbeta <- diag(nvar)+.5*iota%*%t(iota)

階層部分のパラメータ（Delta）は、
> Delta
　　[,1] [,2]
[1,]　　-2　　1
[2,]　　-1　　2
で与えています。

続いて、シミュレーションデータの発生とロジットモデルの尤度関数を定義しています。
単純な二項ロジットなので、尤度関数も簡単です。

## シミュレーションデータの発生
hhdata=NULL
　
for (i in 1:hh) {
　beta <- t(Delta)%*%Z[i,]+as.vector(t(chol(Vbeta))%*%rnorm(nvar))
　X <- matrix(runif(nobs*nvar),nobs,nvar)
　prob <- exp(X%*%beta)/(1+exp(X%*%beta))
　unif <- runif(nobs,0,1)
　y <- ifelse(unif < prob,1,0)
　hhdata[[i]] <- list(y=y,X=X,beta=beta)
}
　
## 2項ロジットモデルの対数尤度関数の定義
loglike <- function(y, X, beta) {
　　p1 <- exp(X %*% beta)/(1 + exp(X %*% beta))
　　ll <- y * log(p1) + (1 - y) * log(1 - p1)
　　sum(ll)
}

推定方法ですが、メトロポリス・ヘイスティングで推定しています。
シミュレーション回数は1万回。

## ベイズ推定のための設定
R<- 10000
sbeta<- 1.5
keep<- 1
　
nhh <- length(hhdata)
nz <- ncol(Z)
　
nvar <- ncol(X)
　
## 事前分布のパラメータ
nu<- nvar+3
V<- nu * diag(rep(1,nvar))
ADelta <- 0.01 * diag(nz)
Deltabar <- matrix(rep(0, nz * nvar), nz, nvar)
　
## サンプリング結果の保存スペースの作成
Vbetadraw <- matrix(double(floor(R/keep) * nvar * nvar), floor(R/keep), nvar * nvar)
betadraw <- array(double(floor(R/keep) * nhh * nvar), dim = c(nhh, nvar, floor(R/keep)))
Deltadraw <- matrix(double(floor(R/keep) * nvar * nz), floor(R/keep), nvar * nz)
　
## 棄却率と対数尤度
reject <- array(0, dim = c(R/keep))
llike <- array(0, dim = c(R/keep))
　
## 初期値の設定
oldbetas <- matrix(double(nhh * nvar), nhh, nvar)
oldVbeta <- diag(nvar)
oldVbetai <- diag(nvar)
oldDelta <- matrix(double(nvar * nz), nz, nvar)
betad <- array(0, dim = c(nvar))
betan <- array(0, dim = c(nvar))
　
## 階層ベイズ2項ロジットモデルによる推定
for (iter in 1:R) {
　　rej <- 0
　　logl <- 0
# sV <- sbeta * oldVbeta # 里村先生のオリジナル
　　sV <- sbeta * V
　　root <- t(chol(sV))
　
　　## MH法による個人別betaのサンプリング
　　for (i in 1:nhh) {
　　　　betad <- oldbetas[i, ]
#　　　betan <- betad + root %*% rnorm(nvar) # 里村先生のオリジナル
　　　　betan <- betad + root %*% rnorm(nvar) * 0.1
　　　　lognew <- loglike(hhdata[[i]]$y, hhdata[[i]]$X, betan)
　　　　logold <- loglike(hhdata[[i]]$y, hhdata[[i]]$X, betad)
　　　　logknew <- -0.5 * (t(betan) - Z[i, ] %*% oldDelta) %*% oldVbetai %*% (betan - t(Z[i, ] %*% oldDelta))
　　　　logkold <- -0.5 * (t(betad) - Z[i, ] %*% oldDelta) %*% oldVbetai %*% (betad - t(Z[i, ] %*% oldDelta))
　　　　alpha <- exp(lognew + logknew - logold - logkold)
　　　　if (alpha == "NaN")
　　　　　　alpha = -1
　　　　u <- runif(1)
　　　　if (u < alpha) {
　　　　　　oldbetas[i, ] <- betan
　　　　　　logl <- logl + lognew
　　　　}
　　　　else {
　　　　　　logl <- logl + logold
　　　　　　rej <- rej + 1
　　　　}
　　}
　　
　　## 多変量回帰によるDeltaのギブスサンプリング(bayesmのrmultiregを利用)
　　out <- rmultireg(oldbetas, Z, Deltabar, ADelta, nu, V)
　　oldDelta <- out$B
　　oldVbeta <- out$Sigma
　　oldVbetai <- solve(oldVbeta)
　
　　## 現在までのサンプリング数の表示
　　if ((iter%%100) == 0) {
　　　　cat("繰り返し数", iter, fill = TRUE)
　　}
　　## keep回毎にサンプリング結果を保存
　　mkeep <- iter/keep
　　if (iter%%keep == 0){
　　　　Deltadraw[mkeep, ] <- as.vector(oldDelta)
　　　　Vbetadraw[mkeep, ] <- as.vector(oldVbeta)
　　　　betadraw[, , mkeep] <- oldbetas
　　　　llike[mkeep] <- logl
　　　　reject[mkeep] <- rej/nhh
　　}
}

1万回のシミュレーションなので、数分待てば、計算が終了します。

続けて、出てきた結果を確認。

## log-likelihoodのプロット
plot(llike,type="l",xlab="Iterations/10",ylab="",main="Log Likelihood")
　
## rejectのプロット
plot(reject,type="l",xlab="Iterations/10",ylab="", main="Rejection Rate of MH Algorithm")
　
## 平均パラメータのプロット
Deltas<-as.vector(Delta)
matplot(Deltadraw,type="l",xlab="Iterations/10",ylab="", main="Draw of Delta")
　
## burn-in 期間
R0<-floor(R/2)+1
R1=R
　
## Deltaの統計値
apply(Deltadraw[R0:R1,],2,mean)
rbind(Deltas,apply(Deltadraw[R0:R1,],2,quantile,probs=c(0.025,0.5,0.975)))

ここまでは、オリジナルコードとほとんど同じです。

続けて、階層ベイズで作ったモデルを予測で使う場合、どうなるかです。