関数解読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"))
■
長らくアップロードしてなかったけど，一段落ついたので今後はぼちぼちやります。
関数解読discint
Bayesian Computation with R: Second Edition (Use R!)
離散確率分布に対して最高確率区間を計算する。
最高確率区間：こちらから与える確率を含む最短区間
確率分布が単峰でない場合は，2つ以上の区間に分かれる。
確率分布が左右対称でない場合は，区間の位置が分布中心から左右のどちらかにずれる。
使用法
discint(dist, prob)
引数
dist:1列目に確率変数の値，2列目に対応する確率を含む行列
prob:最高確率区間を求めるためにこちらから与える確率の値
値
prob:下記のsetに含まれる確率変数のどれかが実現する正確な確率（>=prob）
set:最高確率区間に含まれる確率変数の値の集合
discint = function (dist, prob) #関数discintの中身 { x = dist[, 1] p = dist[, 2] # distの中身をx,pに分けて保存。 n = length(x) # xの長さをnに保存。 sp = sort(p, index.return = TRUE) # pを昇順で並べ替えて，対応する要素の位置も出力し， # spにリストとして保存。 ps = sp$x # 昇順に並べ替えたpをpsに保存。後で使わず。 i = sp$ix[seq(n, 1, -1)] # seq(n, 1, -1) # nから1まで1ずつ減少する数列 # sp$ix[seq(n, 1, -1)] # pの昇順に対応する添字を降順に変換しiに保存。 # 結局，iにはpの要素を降順で並べ替えた場合の # pの要素の位置が保存されている。 ps = p[i] # pの要素を降順で並べ替えてpsに保存。 # ps = sp$x[seq(n,1,-1)]でも良さそう。 xs = x[i] # xをpの要素を降順に対応するように並べ替えてxsに保存。 cp = cumsum(ps) # psの累積和を計算。 # cp[i]=sum(ps[1:i]) ii = 1:n # 別の添字1〜nの作成。 j = ii[cp >= prob] # cp>=probを満たすcpの要素について，iiの位置でTRUEと表す。 j = j[1] # cp>=probを最初に満たしたcpの要素の位置をjに保存。 # cp[j]>=prob かつ cp[j-1]<prob が成立。 eprob = cp[j] # cp[j]をeprobに保存。 set = sort(xs[1:j]) # xs[1:j]を昇順に並べ替え，setに保存。 v = list(prob = eprob, set = set) # eprobとsetをリスト形式にしてvに保存。 return(v) # vを出力。 }
例
x=0:10 probs=dbinom(x,size=10,prob=.3) # 図示 plot(x,probs,type="h") dist=cbind(x,probs) pcontent=.8 (hpi = discint(dist,pcontent)) $prob [1] 0.8214841 $set [1] 1 2 3 4
関数解読histprior
Bayesian Computation with R: Second Edition (Use R!)
（ヒストグラムのように）均等な区間ごとに確率が定義された確率分布の密度を計算する
使用法
histprior(p,midpts,prob)
引数
p:連続型確率変数において，対応する密度を計算する確率変数の値をまとめたベクトル
以下の例のように，できるだけ細かく設定する。
midpts:各区間の中点をまとめたベクトル
prob:各区間の確率をまとめたベクトル
値
各pに対応する確率密度をまとめたベクトル
histprior = function (p, midpts, prob) #関数histpriorの中身 { binwidth = midpts[2] - midpts[1] #ベクトルmidptsの第2要素と第1要素の差から， #1つの区間の長さを求めて，binwidthに代入。 browser() lo = round(10000 * (midpts - binwidth/2))/10000 # 10000 * (midpts - binwidth/2)の小数点以下を四捨五入 # その後，10000で割り，loに代入。 # # 結局，各区間の下限を求めている。 val = 0 * p # ベクトルpと同じ長さの0ベクトルを作成（空箱作り）。 for (i in 1:length(p)) { val[i] = prob[sum(p[i] >= lo)] # p[i] >= lo # p[i]が各区間の下限以上かどうかを判定 # p[i]が第1区間に含まれるなら # TRUE FALSE FALSE.....となる。 # p[i]が第2区間に含まれるなら # TRUE TRUE FALSE.....となる。 # sum(p[i] >= lo) # p[i] >= loが成立する区間の数を集計 # p[i]が第1区間に含まれるなら，TRUEは1個だから1 # p[i]が第2区間に含まれるなら，TRUEは2個だから2 # prob[sum(p[i] >= lo)] # 引数のprobの第"sum(p[i] >= lo)"要素を出力。 # p[i]が第1区間に含まれるなら，prob[1] # p[i]が第2区間に含まれるなら，prob[2] # # 結局，pの各要素（密度を計算する各点）について， # その点を含む区間の確率を密度と定めるということ。 } return(val)#valを出力する。 }
例
#各区間の中点を設定。 midpts=c(.1,.3,.5,.7,.9) #各区間の確率を設定。 prob=c(.2,.2,.4,.1,.1) #密度を計算する各点を設定。 p=seq(.01,.99,by=.0001) #密度を保存。 histden = histprior(p,midpts,prob) #(p,histden)を図示。 plot(p,histden,type="l")
関数解読beta.select
Bayesian Computation with R: Second Edition (Use R!)
百分位点と対応する事前累積確率に関する2組の情報に基づき，それにマッチするベータ分布のパラメータを求める。
使用法
beta.select(quantile1, quantile2)
引数
quantile1
ある100分位点x1（例 50%点）と対応する事前累積確率p1のリスト
quantile2
ある100分位点x2（例 90%点）と対応する事前累積確率p2のリスト
x1 > x2となるようにquantile1,quantile2を設定しても良い。
値
引数の組に対応するベータ分布のパラメータ(a,b)
#関数beta.selectの中身 beta.select = function (quantile1, quantile2) { betaprior1 = function(K, x, p) { #beta.select内で使う関数betaprior1の定義 # p = pbeta(x, K * m0, K * (1 - m0)) # を満たすm0を求めるための関数 m.lo = 0 m.hi = 1 flag = 0 while (flag == 0) { m0 = (m.lo + m.hi)/2 p0 = pbeta(x, K * m0, K * (1 - m0)) #ベータ分布のパラメータを # (a,b)=(K*m0, K*(1-m0)) と置き， # K>0と0<m0<1を満たすパラメータと対応させる。 if (p0 < p) m.hi = m0 else m.lo = m0 if (abs(p0 - p) < 1e-04) # 許容する誤差 1e-04 = 0.0001 flag = 1 } return(m0) }#ここまで関数betaprior1 # リストquantile1，quantile2から各要素を取り出す。 p1 = quantile1$p x1 = quantile1$x p2 = quantile2$p x2 = quantile2$x logK = seq(-3, 8, length = 100) K = exp(logK) #パラメータKの候補を100個作る。 # なぜlogKを-3から8の範囲の限定したのだろうか？ m = sapply(K, betaprior1, x1, p1) # パラメータ(K[i],x1,p1) i=1,2,,..,100 # に対して関数betaprior1を適用。 # m[i] i=1,2,,..,100を求める。 prob2 = pbeta(x2, K * m, K * (1 - m)) # パラメータ(K[i]*m[i],K[i]*(1-m[i]) )のベータ分布 # のx=x2における累積確率分布の値（ベクトル） # を求めてprob2に代入。 ind = ((prob2 > 0) & (prob2 < 1)) # 0 < prob2 <1 を満たすprob2の要素の位置をTRUEとする。 app = approx(prob2[ind], logK[ind], p2) # 関数approxで線形補間 # 点(x,y)=(prob2[ind],logK[ind])を # つなぐ線上でx=p2に対応するy=logKをapprox$yとして出力。 K0 = exp(app$y) m0 = betaprior1(K0, x1, p1) # 上で求めたK0について， # 点x1においてペアになるパラメータm0を求める。 return(round(K0 * c(m0, (1 - m0)), 2)) # (K0*m0,K0*(1-mo)を小数点第3位を四捨五入して， # 小数点第2位まで出力。 # 部分部分は一応理解できるが， # 全体像にまだスッキリ分かった感がない。 # パラメータが2つあるから， # 累積分布関数の異なる2点の情報があれば， # そのパラメータを求めることができるということか。 }
例
q1=list(p=.5,x=0.25) q2=list(p=.9,x=0.45) # 中央値の事前累積確率は0.25 # 90%分位点の事前累積確率は0.45 beta.select(q1,q2) #関数の作り方から考えて，引数を入れ替えてもOK。 beta.select(q2,q1) #結果は一致。 beta_para = beta.select(q1,q2) [1] 2.67 7.37 pbeta(c(0.25,0.45),beta_para[1],beta_para[2]) # a=beta_para[1]=2.67 # b=beta_para[2]=7.37 を使ってβ分布の累積確率を計算。 # q1,q2に（ほぼ）対応する累積確率を得た。 [1] 0.5001631 0.9001404
関数解読pdiscp
Bayesian Computation with R: Second Edition (Use R!)
比率に関する離散分布を用いた場合の将来の二項分布による実験の成功数に関する予測分布の計算
使用法
pdiscp(p, probs, n, s)
引数
p:比率の候補となる値で構成されるベクトル
probs:各比率の確率で構成されるベクトル
n:将来の二項分布による実験のサンプル数
s:将来の二項分布による実験における成功数
seq(0,n)のうち，関心を持つ要素のみに限定しても良い。
値
成功数の予測確率（ベクトル）
確率probsが事前確率なら，事前予測確率
確率probsが事後確率なら，事後予測確率（BCwithR p29）
#関数pdiscpの中身 pdiscp = function (p, probs, n, s) { pred = 0 * s # ベクトルsと同じ長さのベクトルpredを作成（空箱作り）。 for (i in 1:length(p)) { pred = pred + probs[i] * dbinom(s, n, p[i]) # sの各要素に対して予測確率を計算。 # BCwithR p29一番下の数式 } return(pred) #predを出力。 }
例
#比率と対応する確率（その1） p1=c(.1,.2,.3,.4,.5,.6,.7,.8,.9) prob1=c(0.05,0.10,0.10,0.15,0.20,0.15,0.10,0.10,0.05) #比率と対応する確率（その2） p2=c(.1,.2,.3,.4,.5,.6,.7,.8,.9) prob2=c(0.2,0.2,0.2,0.2,0.2,0,0,0,0)
par(mfrow=(c(2,1))) plot(p1,prob1,type="h") plot(p2,prob2,type="h") #将来の試行回数 n=10 #予測確率を計算する成功数 s=0:10 #予測確率の計算 yf1 = pdiscp(p1,prob1,n,s) yf2 = pdiscp(p2,prob2,n,s) par(mfrow=(c(2,1))) plot(s,yf1,type="h") plot(s,yf2,type="h")
関数解読pdisc
Bayesian Computation with R: Second Edition (Use R!)
離散事前分布を用いた場合の比率に関する事後分布の計算
使用法
pdisc(p, prior, data)
引数
p:（推測する）比率の候補となる値で構成されるベクトル
prior:候補となる各比率の事前確率で構成されるベクトル
data:成功と失敗の数で構成されるベクトル
確率pで成功が起こるベルヌーイ分布でモデル化
値
候補となる各比率の事前確率で構成されるベクトル
#関数pdiscの中身 pdisc = function (p, prior, data) { s = data[1] f = data[2] p1 = p + 0.5 * (p == 0) - 0.5 * (p == 1) #ベクトルpの要素に0があれば，その要素に0.5を加える。 #ベクトルpの要素に1があれば，その要素から0.5を引く。 # 上の操作は何を意味する？ #log(0)は計算できないからか。 #ベクトルpの要素に0も1もなければ，p1=pでこれが通常。 like = s * log(p1) + f * log(1 - p1) #対数尤度(A) （定数項を除く） #確率p1で成功がs回発生 #確率1-p1で失敗がf回発生 like = like * (p > 0) * (p < 1) - 999 * ((p == 0) * (s > 0) + (p == 1) * (f > 0)) #通常は第1項目のみでOK。 # p>0かつp<1の場合は対数尤度(A)をそのまま代入。 # p==0かつs>0 （成功の確率は0なのに成功した結果あり） # p==1かつf>0 （失敗の確率は0なのに失敗した結果あり） # この場合，log(0)=-Infの代わりに-999を代入。 like = exp(like - max(like)) #対数尤度を尤度に戻す。 # max(like)を引く意味は？ # 後でsum(post)=1となるよう基準化するので， # -max(like)はなくても同じ。 # exp(like - max(like))=exp(like)*exp(-max(like)) # exp()の値を発散させないよう念のため付けてる？ product = like * prior #事後分布のカーネル post = product/sum(product) #事後確率 #sum(post)=1がとなるようにproductを基準化 return(post)#事後確率postを出力 }
例
#比率（成功確率）と対応する事前確率 p=c(.0,.25,.3,1) prior=c(.25,.25,.25,.25) #成功と失敗の数 data=c(5,10) #比率（成功確率）の事後確率 pdisc(p,prior,data)
