The One with ...

研究用のメモなどを公開しています

前回の続き.

そもそも分布を導出する理論モデルのパラメータを階層ベイズMCMCで推定できるのか,という確認をやってみる.

まずデータとして,理論モデルからのサンプルをつくる.

今回はJAGSRから呼び出すので, コードは全部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が投資する資本割合である.

で,結果は

f:id:hamada7418:20161002163636p:plain:w600

こんな感じでpの事後分布のモードは真値に近いけど,bの方はあまりよくない.

中心極限定理の観点から考えるとpの分散は0でなくてもいいけど(リアプノフの条件から直感的にそう思う),bのほうが理論モデルにあっていないのかもしれない.つまり,そもそもbが確率分布していると,対数正規分布をGMIから導出できないのかもしれない.これでは元も子もない.

ここまでの疑問を整理すると

  • ハイパーパラメータはそもそも何個指定できるのか?
  • 実際に所得データからパラメータを推定する際に「足りない情報」をどの変数で補うのか?
  • 分布Aのパラメータが階層ベイズの仮定として(Aと共役でない)確率分布Bに従う場合,データがAに従うという仮定は妥当か?

といったところ.

続く