jags_test -- Sage

# Rのユーティリティ関数を読み込む attach(DATA + "RUtil.py")

# 引用 http://www.singularpoint.org/blog/r/mcmc-jags-1/ r('library(rjags)') r('library(coda)')

 [1] "rjags"     "coda"      "lattice"   "stats"     "graphics" 
"grDevices" "utils"     "datasets" 
 [9] "methods"   "base"

 [1] "rjags"     "coda"      "lattice"   "stats"     "graphics"  "grDevices" "utils"     "datasets" 
 [9] "methods"   "base"

# jagsのモデル print open(DATA+"lm.jags").read()

model {
  for(i in 1:N){
    mu[i] <- alpha + beta*x[i]
    y[i] ~ dnorm(mu[i],1/(tau^2))
  }
  alpha ~ dnorm(0,0.0001)
  beta ~ dnorm(0,0.0001)
  tau ~ dgamma(0.1,0.1)
}

model {
  for(i in 1:N){
    mu[i] <- alpha + beta*x[i]
    y[i] ~ dnorm(mu[i],1/(tau^2))
  }
  alpha ~ dnorm(0,0.0001)
  beta ~ dnorm(0,0.0001)
  tau ~ dgamma(0.1,0.1)
}

# determine the number of samples r('N <- 50') # data generation r('x <- rnorm(N) ') junk = r('y <- 2 + 3*x + rnorm(N,0.1)')

# compiling and initializing # jags.model 関数の第一引数では先程のモデルファイルを指定する # 第二引数 data では JAGS に渡したいRオブジェクトをリスト形式で与える # n.chains で並列パスの数を指定 jags_file = DATA + 'lm.jags' junk = r('lm.mcmc <- jags.model("%s", data=list("x"=x, "y"=y, "N"=N), n.chains=4, n.adapt=100)'%jags_file)

# mcmc sampling # update で burn-in している r('update(lm.mcmc,1000) # burn-in') # coda.samples で指定した数のサンプルを coda 形式で取得する junk = r("posterior <- coda.samples(lm.mcmc, c('alpha','beta'), 2000)")

# 画像の変換に時間がかかる graph = preGraph("fig1.pdf") # plotting the result r('plot(posterior)') postGraph(graph)

jags_test

4735 days ago by takepwave