The One with ...

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

研究会用メモ

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

というわけで

  1. gmcから直接対数正規分布のパラメータを指定して,乱数を発生
  2. そのデータからパラメータを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

modelstring1n,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)

モデル1の推定結果

f:id:hamada7418:20161027223220p:plain

やはり理論的に固定すべき変数を真値で与えたデータの場合,未知のp,bの推定が 正確である.

問題はnt,y0tが現実のデータから推測できないことである.さてどうしよう...