The One with ...

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

JAGSによる対数正規分布のパラメータ推定

Kruschke氏のblogでは,対数正規分布のデータを生成するのに,一度正規乱数作ってから基準化して指数化するというやや煩雑な手順を用いている.

doingbayesiandataanalysis.blogspot.jp

最初にコードを見たとき,BUGSのクセとあいまって,何をやっているのか分からなかった. 単純化のため,対数正規分布を使って,直接データを生成してみる.

set.seed(47405)
y = rlnorm( n=500, 5.0, 0.5)  

# Package the data for shipping to JAGS:
dataList = list(
  y = y ,
  N = length(y) ,
  meanOfX = mean(log(y)) ,#事前分布のパラメータに使う
  sdOfX = sd(log(y)))#事前分布のパラメータに使う

JAGSに渡すコードは以下の通り

modelstring = "
  model {
    for( i in 1 : N ) {
      y[i] ~ dlnorm( mu, 1/sigma^2)#JAGS codeの定義より
    }
  sigma ~ dunif( 0.001*sdOfX , 1000*sdOfX )
  mu ~ dunif( 0.001*meanOfX , 1000*meanOfX )
  meanY <- exp(mu+(sigma^2)/2)
  modeOfY <- exp(mu-sigma^2)
  sigmaOfY <- sqrt(exp(2*mu+sigma^2)*(exp(sigma^2)-1))
}
" # close quote for modelstring
writeLines(modelstring,con="model.txt")

一応これでもMCMC推定はうまくいく.

(今思いついたが,事前分布のパラメータの範囲を指定するのに正規分布を1回経由したほうが便利だから,ああいうコードになっていたのかもしれん)

理論モデルから生成したデータに基づく推定

次の手順として,対数正規分布のパラメータを y_0, b, p, nの関数とみなして,ダイレクトに対数正規分布確率密度関数を使って推定してみる.

要するに理論モデル(GMI)に基づいてデータを人工的に発生させて,理論上のパラメータ(潜在変数)をMCMC推定できるのか? という実験である.

まずデータは次のように生成する

#直接対数正規分布のパラメータを指定する
set.seed(47406)

y0=20; b=0.2; p=0.65; n=10;
mu=log(y0)+n*log(1-b)+log((1+b)/(1-b))*n*p
sigma= sqrt(n*p*(1-p)*(log((1+b)/(1-b)))^2) 

y = rlnorm( n=100, mu, sigma)  

JAGS用モデルは次の通り(これをモデルと呼ぶことに依然として抵抗を感じるが今は無視する).

#modelstring0はパラメータをn,y0,p,b,を全て推定
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(13,18)
y0 ~ dpois(c)#
c ~ dunif(10,18) 
p ~ dbeta(1,1)#p,bの事前分布にβ分布
b ~ dbeta(1,1)}"
# close quote for modelstring

この場合真値をHDIに含む推定をしてくれる.ただ関数gmiは人工データとしては,あまりうまく機能しないようだ.

(自由度に問題がないのか,若干気になる.直感的に言って自由度が高すぎてMCMCが収束しない場合はp以外のパラメータを決定論的に 指定する必要がある.)