読者です 読者をやめる 読者になる 読者になる

チラシの裏の落書き日記

統計とか,研究とか,日常についての備忘録的なもの。

jagsUIなるパッケージ

R

MCMCサンプラーのうちの1つであるJAGSをRから利用する新しいパッケージができたようだ。基本はrjagsのラッパーのようだが,並列計算の指定などがやりやすそうだ。今後,これを中心に使っていくかもしれないのでヘルプにあった例をコピペしておく。初期値の指定が関数になっているが,いいのだろうか。詳細はのちほど検証のこと。

マニュアル
https://cran.r-project.org/web/packages/jagsUI/jagsUI.pdf
作者のgithub
GitHub - kenkellner/jagsUI: Run JAGS (Just Another Gibbs Sampler) analyses from within R.


#install.packages("jagsUI")
library("jagsUI")

#Analyze Longley economic data in JAGS
#Number employed as a function of GNP
######################################
##   1. Collect and Package Data    ##
######################################
#Load data (built into R)
data(longley)
head(longley)
#Separate data objects
gnp <- longley$GNP
employed <- longley$Employed
n <- length(employed)
#Input data objects must be numeric, and must be
#scalars, vectors, matrices, or arrays.
#Package together: several possible ways

#1. A named list of the objects
data <- list(gnp=gnp,employed=employed,n=n)
#2. A character vector of the names of the objects
data <- c( gnp , employed , n )
#3. A list of names of the objects
data <- list( gnp , employed , n )

######################################
##      2. Write model file         ##
######################################
#Write a model in the BUGS language
#Generate model file directly in R
#(could also read in existing model file)
writeLines("
           model{
           #Likelihood
           for (i in 1:n){
           employed[i] ~ dnorm(mu[i], tau)
           mu[i] <- alpha + beta*gnp[i]
           }
           #Priors
           alpha ~ dnorm(0, 0.00001)
           beta ~ dnorm(0, 0.00001)
           sigma ~ dunif(0,1000)
           tau <- pow(sigma,-2)
           }
           ", con="model.txt")
#Identify filepath of model file;
#in this case in the working directory
modfile <-  "model.txt"

######################################
##    3. Initialize Parameters      ##
######################################
#Best to generate initial values using function
inits <- function(){
  list(alpha=rnorm(1,0,1),beta=rnorm(1,0,1),sigma=runif(1,0,3))
}
inits
#In many cases, JAGS can pick initial values automatically;
#you can leave argument inits=NULL to allow this.

######################################
##  4. Set parameters to monitor    ##
######################################
#Choose parameters you want to save output for
#Only parameters in this list will appear in output object
#(deviance is added automatically if DIC=TRUE)
#List must be specified as a character vector
params <- c('alpha' , 'beta' , 'sigma' )

######################################
##        5. Run Analysis           ##
######################################
#Call jags function; specify number of chains, number of adaptive iterations,
#the length of the burn-in period, total iterations, and the thin rate.
out <- jags(data = data,
            inits = inits,
            parameters.to.save = params,
            model.file = modfile,
            n.chains = 3,
            n.adapt = 100,
            n.iter = 1000,
            n.burnin = 500,
            n.thin = 2,
            parallel = T,
            n.cores=3)
#Arguments will be passed to JAGS; you will see progress bars
#and other information
#Examine output summary
out
#Look at output object elements
names(out)
#Plot traces and posterior densities
plot(out)
#Plot traces
traceplot(out)
#Update model another 1000 iterations
out <- update(out,n.iter = 1000)

追記
JAGSの並列計算を楽にしてくれる,runjagsパッケージもあるようだ。
https://cran.r-project.org/web/packages/runjags/runjags.pdf
このへんの使い方も理解しておきたい。
MCMCをするときの計算時間が尋常じゃないので,こういうものを使って少しでも楽になればいい。並列計算のやりかたを調べたり何だりしているときの虚しさたるや。。。それはもう悲惨なので。。。