■

前回の続き．

そもそも分布を導出する理論モデルのパラメータを階層ベイズ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に従うという仮定は妥当か？

といったところ．

続く