関数解読bayesprobit
Bayesian Computation with R: Second Edition (Use R!)
2値反応回帰モデルの一種であるプロビットモデルの回帰係数の同時事後分布からのサンプリング 回帰係数の事前分布は精度行列がゼロ行列の場合とそうでない場合の両方に対応する。
後者の場合は対数周辺尤度も計算する。
使用法
bayes.probit(y,X,m,prior=list(beta=0,P=0))
引数
y:二値反応（被説明変数）ベクトル
X:説明変数行列
m:シミュレーションの回数
prior:回帰係数に関する事前分布のパラメータ（リスト形式で保存）
beta:平均
P:精度行列（分散共分散行列の逆行列）
Pがゼロ行列のときは事前分布は無情報事前分布
それ以外の時はPは正値定符号行列でなければならない
デフォールトは無情報事前分布:list(beta = 0, P = 0)
値
beta:m行p列（pは説明変数の数）の行列
各行が回帰係数の周辺事後分布からサンプリングした値に対応。
log.marg:シュミレーションによる対数周辺尤度の推定値。（P=0でない場合に出力。）
参考
BCwithR p240-247
Albert and Chib (1993)JASA
http://scholar.google.co.jp/scholar?cluster=9134134856363849481&hl=en&as_sdt=2000
bayes.probit = function (y, X, m, prior = list(beta = 0, P = 0)) #関数bayes.probitの中身 { #切断された分布からのサンプリングのための関数 # LearnBayesに同一内容の関数が含まれるが， # bayes.probitの中でも定義。 rtruncated = function(n, lo, hi, pf, qf, ...){ qf(pf(lo, ...) + runif(n) * (pf(hi, ...) - pf(lo, ...)), ...)} if (sum(prior$P) == 0) log.marg = NULL # y*=Xbeta+e における係数betaの事前分布が無情報分布の場合 # 引数の精度行列Pはゼロ行列。（分散∞） # この場合（無情報事前分布の場合），周辺尤度は計算できないので， # log.margは作成しない。 beta0 = prior$beta #betaの事前分布の期待値をbeta0に代入。 BI = prior$P # 精度行列PをBIに代入。 N = length(y) # データの数をNとする。 fit = glm(y ~ X - 1, family = binomial(link = probit)) # 一般化線形モデルを推定する関数glm（statパッケージ）を利用 # probitモデルを以下の手法（頻度主義）で推定。 # Xに定数項を含むので説明変数行列を"X-1"とする。 # helpによると推定方法は # iteratively reweighted least squares (IWLS) # BCwithR p242にはMLと書いているが。 beta.s = fit$coef # 推定した係数をbeta.sに保存。 p = length(beta.s) # beta.sの長さ（係数の数）をpとする。 beta = array(beta.s, c(p, 1)) #beta.sをp行1列のbetaに保存。 #以下で実行するMCMCの初期値になる beta0 = array(beta0, c(p, 1)) #beta0をp行1列のbeta0に保存。 BI = array(BI, c(p, p)) #BIをp行p列のBIに保存。 Mb = array(0, dim = c(m, p)) # MCMCの結果を保存する入れ物Mb（m行p列）を作成。 # array：2次元なら行列。3次元以上の設定も可能。 lo = c(-Inf, 0) # -∞と0をベクトルloに保存。 hi = c(0, Inf) # 0と∞をベクトルhiに保存。 LO = lo[y + 1] HI = hi[y + 1] #y:被説明変数となる0 or 1の離散変数ベクトル #y+1:yの各要素に1を加えたベクトル。0→1，1→2 #lo[y+1]:ベクトル"y+1"において # 要素が1の場合，その要素を-Infに置き換え， # 要素が2の場合，その要素を0に置き換えたベクトル。 #hi[y+1]:ベクトル"y+1"において # 要素が1の場合，その要素を0に置き換え， # 要素が2の場合，その要素をInfに置き換えたベクトル。 #LOとHIはrtruncatedの引数になる。 # y=0の場合，z<0だから(-Inf,0)の切断正規分布からサンプリング。 # y=1の場合，z<0だから(0,Inf)の切断正規分布からサンプリング。 # #便利そうだが，外形（lo[y+1],hi[y+1]）からは想像も付かないなぁ。 postvar = solve(BI + t(X) %*% X) #行列"BI + t(X) %*% X"の逆行列を計算。 # BI=0（non informative prior）の場合 # (X'X)^(-1) BCwithR p241一番上 # その他（proper prior）の場合 # BCwithR p244 真ん中当たりのV1 aa = chol(postvar) # postvarをコレスキー分解して得た上三角行列 # t(aa) %*% aa = postvar BIbeta0 = BI %*% beta0 # =0 （BI=0のとき） # BCwithR p244 真ん中当たりのbeta1のパーツ post.ord = 0 # BCwithR p244 真ん中当たりのbeta1のパーツ # MCMCの実行。 for (i in 1:m) { z = rtruncated(N, LO, HI, pnorm, qnorm, X %*% beta, 1) # 正規分布 N(mean=X %*% beta,sd=1)を切断した分布からサンプリング。 # yの各要素について，切断の異なる下限と上限LO,HIが設定されている。 mn = solve(BI + t(X) %*% X, BIbeta0 + t(X) %*% z) # solve(A,b) は A %*% x = bの解：x=A^(-1)b # この場合，betaの平均beta1を求めている。 # proper priorの場合，BCwithR p244のbeta1 beta = t(aa) %*% array(rnorm(p), c(p, 1)) + mn # zとdataを所与としたbetaの条件付分布からのサンプリング。 # noninformative BCwithR p241一番上 # proper BCwithR p244真ん中当たり # Y~N(0,1)のときX=U'Y+muはN(mu,U'U)に従うことを利用。 post.ord = post.ord + dmnorm(beta.s, mn, postvar) # 周辺尤度を計算するためのパーツ。 # dmnorm()は多変量正規分布の密度関数の値を計算する関数。 # 平均beta1，分散共分散行列postvarの場合 # beta=beta.sで密度関数を評価 # beta1はzに依存して毎回異なる # 最終的にpost.ordでm×ΣΦ(beta.s;beta1,V1) # が計算できる。BCwithR p246 一番上の式にmを掛けたもの。 Mb[i, ] = beta #p次元ベクトルbetaをMbのi行目に保存。 } #対数周辺尤度の計算。 if (sum(BI) > 0) { # 精度行列がゼロ行列（分散∞）でないproper priorの場合に実行。 log.f = sum(y * log(fit$fitted) + (1 - y) * log(1 - fit$fitted)) # beta=beta.sにおける尤度の値（yの各要素に対応するベクトル） # fit$fittedを上手く活用（推定がMLならOK。他の場合は？） log.g = dmnorm(beta.s, beta0, solve(BI), log = TRUE) # betaの事前分布の対数密度関数をbeta=beta.sで評価した値。 log.marg = log.f + log.g - log(post.ord/m) # 対数周辺尤度の計算。BCwithR p246 のlog[m(y)] } #MCMCの結果Mbと対数周辺尤度の値log.margをリストにして出力。 return(list(beta = Mb, log.marg = log.marg)) }
例
# 被説明変数 response=c(0,1,0,0,0,1,1,1,1,1) # 説明変数 covariate=c(1,2,3,4,5,6,7,8,9,10) # 説明変数行列（切片項を含む） X=cbind(1,covariate) # 事前分布 その1 prior1=list(beta=c(0,0),P=diag(c(.5,10))) # 事前分布 その２ prior2=list(beta=c(0,0),P=diag(c(0,0))) # MCMCの実行 その1 m=10000 s1=bayes.probit(response,X,m,prior1) par(mfrow=c(2,1)) plot(s1$beta[,1],type="l",col="red") plot(s1$beta[,2],type="l",col="red") par(mfrow=c(2,1)) hist(s1$beta[2001:10000,1],freq=F,breaks="Scott",col="grey") abline(v=mean(s1$beta[2001:10000,1]),col="red",lwd=2) abline(v=median(s1$beta[2001:10000,1]),col="green",lwd=2) legend("topright", c("mean", "median"), lty = c(1,1), lwd = c(2,2), col = c("red", "green")) hist(s1$beta[2001:10000,2],freq=F,breaks="Scott",col="grey") abline(v=mean(s1$beta[2001:10000,2]),col="red",lwd=2) abline(v=median(s1$beta[2001:10000,2]),col="green",lwd=2) legend("topright", c("mean", "median"), lty = c(1,1), lwd = c(2,2), col = c("red", "green")) #対数周辺尤度 s1$log.marg # [1] -7.09734
# MCMCの実行 その2 m=10000 s2=bayes.probit(response,X,m,prior2) par(mfrow=c(2,1)) plot(s2$beta[,1],type="l",col="red") plot(s2$beta[,2],type="l",col="red") par(mfrow=c(2,1)) hist(s2$beta[2001:10000,1],freq=F,breaks="Scott",col="grey") abline(v=mean(s2$beta[2001:10000,1]),col="red",lwd=2) abline(v=median(s2$beta[2001:10000,1]),col="green",lwd=2) legend("topright", c("mean", "median"), lty = c(1,1), lwd = c(2,2), col = c("red", "green")) hist(s2$beta[2001:10000,2],freq=F,breaks="Scott",col="grey") abline(v=mean(s2$beta[2001:10000,2]),col="red",lwd=2) abline(v=median(s2$beta[2001:10000,2]),col="green",lwd=2) legend("topright", c("mean", "median"), lty = c(1,1), lwd = c(2,2), col = c("red", "green"))
