■
前回の続き.
そもそも分布を導出する理論モデルのパラメータを階層ベイズMCMCで推定できるのか,という確認をやってみる.
まずデータとして,理論モデルからのサンプルをつくる.
今回はJAGS
をR
から呼び出すので, コードは全部R
で済むように書き換える.
data生成用コードはざっと以下のような感じ
# Generate a person's income under the given condition: set.seed(47406) gmi <- function(y1, b0, p0, n0){ y = y1 for(i in 1:n0){ if(runif(1) < p0){ y <- y + (y*b0)} else{y <- y - (y*b0)} }#for end y } #generate income data for m people m=2000; y1=10; b0=0.2; p0=0.5; n0=15 dat0 <- c() for(i in 1:m){ dat0[i] <- gmi(y1, b0, p0, n0)}
ここで定義しているパラメータが真の値に相当する
次にJAGSに渡すモデルの記述
modelstring0 = " model{ for( i in 1 : N ) { y[i] ~ dlnorm(m, 1/(s^2)) } m <- log(y0)+n*log(1-b)+log((1+b)/(1-b))*n*p s <- sqrt(n*p*(1-p)*(log((1+b)/(1-b)))^2) n ~ dpois(a)#nをポアソン分布で近似 a ~ dunif(1,20) y0 ~ dpois(c)# c ~ dunif(1,20) p ~ dbeta(1,1)#p,bの事前分布にβ分布 b ~ dbeta(1,1)}"
実行には例によってJohn K. Kruschke.さんのラッパー関数DBDA2E-utilities.Rをお借りする. (MCMCの設定はデフォルト)
さて,実行してみると,chainは収束してるみたいだけど,
なんかうまく推定できていない気がする.
よくよく考えたら,n
だのy_0
だの理論的に定数であるべきパラメータを(非負だからといって)
ポアソン分布に従うとか雑に仮定したのがよくないのでは?
というわけで修正した統計モデルがこちら
modelstring1 = " model{ for( i in 1 : N ) { y[i] ~ dlnorm(m, 1/(s^2)) } m <- log(y1)+n0*log(1-b)+log((1+b)/(1-b))*n0*p s <- sqrt(n0*p*(1-p)*(log((1+b)/(1-b)))^2) #prior for hyper parameters p ~ dbeta(1,1) b ~ dbeta(1,1)}"
p, b
以外は定数扱いとして,事前に定義しておく.推定しているのは実質的には
p, b
のみである.具体的にはp
が投資成功確率でb
が投資する資本割合である.
で,結果は
こんな感じでpの事後分布のモードは真値に近いけど,b
の方はあまりよくない.
中心極限定理の観点から考えるとp
の分散は0でなくてもいいけど(リアプノフの条件から直感的にそう思う),b
のほうが理論モデルにあっていないのかもしれない.つまり,そもそもb
が確率分布していると,対数正規分布をGMIから導出できないのかもしれない.これでは元も子もない.
ここまでの疑問を整理すると
- ハイパーパラメータはそもそも何個指定できるのか?
- 実際に所得データからパラメータを推定する際に「足りない情報」をどの変数で補うのか?
- 分布Aのパラメータが階層ベイズの仮定として(Aと共役でない)確率分布Bに従う場合,データがAに従うという仮定は妥当か?
といったところ.
続く