研究会用メモ - The One with ...

案の定，2週間ほど前にやった計算をすっかり忘れて，同じコードを書くところであった．ブログにメモを残しておいたおかげで，思い出した．

というわけで

gmcから直接対数正規分布のパラメータを指定して，乱数を発生
そのデータからパラメータをMCMC推定する

という実験の続きを再開する（10/11のエントリの続きです）．参照したコードは，例によって Kruschke氏のコードである．

準備はこれだけ

source("DBDA2E-utilities.R")
fileNameRoot="Jags-lognormal-direct-sim"
graphFileType="png"

まずR内で人工データを発生させる

# "t" means "true value"
y0t=10; bt=0.1; pt=0.3; nt=25;
mu=log(y0t)+nt*log(1-bt)+log((1+bt)/(1-bt))*nt*pt
sigma= sqrt(nt*pt*(1-pt)*(log((1+bt)/(1-bt)))^2) 
gmc = rlnorm( n=500, mu, sigma)

pt, bt, nt が推定すべき真のパラメータである． mu, sigmaがこのような決定論的関数で表される理由はgmcで証明が与えられる（ここ重要）

次にJAGSにわたすためのデータを定義する

dataList = list(
  y = gmc ,
  N = length(gmc),
  nt=nt,
  y0t=y0t)#JAGSに渡すパラメータはlist内で再定義する
# nt, y0tはR内で定義されているがJAGSにはdataを介して
# 渡している．

nt,y0tはスカラである．dataListはRオブジェクトだが，modelstringで参照するために，このようなリストを作っておく

JAGSにわたすモデルを複数定義する．今回は

modelstring0

modelstring1

modelstring2

の三つ．

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(mean.n)#nをポアソン分布で近似
mean.n ~ dunif(1,2*nt)
y0 ~ dpois(mean.y0)#
mean.y0 ~ dunif(1,2*y0t) # y0t=10; bt=0.2; pt=0.65; nt=25;
p ~ dbeta(1,1)#p,bの事前分布にβ分布
b ~ dbeta(1,1)}"
# close quote for modelstring

modelstring0は推定するパラメタが一番多い．

つぎは，より理論モデルに忠実な設定．

#modelstring1は理論的に固定すべきパラメータn,y0を定数にする

modelstring1 = "
model{
for( i in 1 : N ) {
y[i] ~ dlnorm(m, 1/(s^2))
}
y0 <- y0t
n <- nt
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)
#prior for hyper parameters
p ~ dbeta(1,1)
b ~ dbeta(1,1)
}"
# close quote for modelstring

modelstring1はn,y0を定数として固定する．

ここで，先ほど定義したリスト内のスカラを割り当てる（最初はRのグローバル変数を参照していたのだが,modelstringの中からは参照できないことが分かった．まあ，よく考えればこっちはJAGSが読むコードなので，当然である）

最後のモデル．

# dilectly assign mean of poisson dist. of n,y0
modelstring2 = "
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(nt)#nをポアソン分布で近似
y0 ~ dpois(y0t)#
p ~ dbeta(1,1)#p,bの事前分布にβ分布
b ~ dbeta(1,1)}"

最後に，ハイパーパラメータだけ，真の値を代入してその他を推定するモデル

モデルの選択はこんな風にswitchで定義する．（文字型しかダメなのだろうか？）

a<-"1"

modelstring<-switch(a,
"0"=modelstring0,
"1"=modelstring1,
"2"=modelstring2)

モデルによって推定するパラメータは直に指定する

parameters0 = c("p","b","n","mean.n","mean.y0","y0","m","s")
parameters1 = c("p","b","m","s")
parameters2 = c("p","b","n","y0","m","s")

parameters<-switch(a,
  "0"=parameters0,
  "1"=parameters1,
  "2"=parameters2)

あとは細かい話だが，plotPostに渡すパラメタがモデルによって異なるので， if文でモデル毎に指定しておく．

分析の度に，モデル選択変数aを指定することに注意する．

    #全model共通で推定するパラメタ
    postInfo = plotPost( mcmcChain[,"p"] , xlab="success p" )
    postInfo = plotPost( mcmcChain[,"b"] , xlab="interest b" )

  #以下モデルに合わせてplotするパラメタを選択  
    if(a=="0") postInfo = plotPost( mcmcChain[,"n"] , xlab="time n" )
    if(a=="0") postInfo = plotPost( mcmcChain[,"mean.n"] , xlab="mean of n" )
    if(a=="0") postInfo = plotPost( mcmcChain[,"mean.y0"] , xlab="mean of y0" )  
    
    if(a=="1") postInfo = plotPost( mcmcChain[,"m"] , xlab="mu of y")
    if(a=="1") postInfo = plotPost( mcmcChain[,"s"] , xlab="sigma of y")
    
    if(a=="2") postInfo = plotPost( mcmcChain[,"n"] , xlab="time n" )
    if(a=="2") postInfo = plotPost( mcmcChain[,"y0"] , xlab="y0" )
    if(a=="2") postInfo = plotPost( mcmcChain[,"m"] , xlab="mu of y")
    if(a=="2") postInfo = plotPost( mcmcChain[,"s"] , xlab="sigma of y")

    saveGraph(file=fileNameRoot,type=graphFileType)