研究会用メモ
案の定,2週間ほど前にやった計算をすっかり忘れて,同じコードを書くところであった. ブログにメモを残しておいたおかげで,思い出した.
というわけで
という実験の続きを再開する(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)
モデル1の推定結果
やはり理論的に固定すべき変数を真値で与えたデータの場合,未知のp,b
の推定が
正確である.
問題はnt,y0t
が現実のデータから推測できないことである.さてどうしよう...