階層ベイズで作成したパラメータのシミュレーションについて その1 [階層ベイズ]
里村先生の『マーケティング・モデル』をベースに、階層ベイズのモデルを作り、シミュレーションを行いました。
コードは、教科書に載っていますし、共立出版のホームページからもダウンロードできます。
以下、里村先生のコードをベースに修正したものを載せていますが、
一部、教科書そのままだと、上手く動かない部分があるので修正しています。
(オリジナルのコードは、上記を参照してください。)
まずは、初期設定から。
階層部分のパラメータ(Delta)は、
> Delta
[,1] [,2]
[1,] -2 1
[2,] -1 2
で与えています。
続いて、シミュレーションデータの発生とロジットモデルの尤度関数を定義しています。
単純な二項ロジットなので、尤度関数も簡単です。
推定方法ですが、メトロポリス・ヘイスティングで推定しています。
シミュレーション回数は1万回。
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)))
ここまでは、オリジナルコードとほとんど同じです。
続けて、階層ベイズで作ったモデルを予測で使う場合、どうなるかです。