int<lower=0> N // N은 lower bound가 0인, integer 형 변수라는 의미.
modelString = "
data {
int<lower=0> N ;
int y[N] ; // y is a length-N vector of itegers
}
parameters {
real<lower=0, uppper=1> theata ;
}
model {
theta ~ beta(1, 1) ;
y ~ bernoulli(theta) ;
}
" # close quote for modelString
y ~ bernoulli(theta) ;
for ( i in 1:N ) {
y[i] ~ dbern(theta)
}
library(rstan)
# translate the model into C++ code
# and compile the C++ code into an executable dynamic shared object(DSO)
stanDso = stan_model(model_code=modelSTring)
# Create some fictitious data :
N = 50; z = 10; y = c(rep(1,z), rep(0,N-z))
dataList = list(y=y, N=N)
stanFit = sampling(object=stanDso, data=dataList,
chains=3, iter=1000, warmup=200, thin=1)
# Load rjags, coda, and DBDA2E functions :
source("DBDAE-utlities.R")
# Convert stan format to coda format :
mcmcCoda = mcmc.list(lapply(1:ncol(stanFit),
function(x) { mcmc(as.array(stanFit)[,x,]) } ) )
# Graph chin diagnostics using DBDA2E function :
diagMCMC( mcmcCoda, parName=c("theta"))
data {
... declarations ...
}
transforemd data {
... declarations ... statements ...
}
parameters {
... declarations ...
}
transforemd parameters {
... declarations ... statements ...
}
model {
... declarations ... statements ...
}
generated quantities {
... declarations ... statements ...
}
y ~ normal(mu, sigma)
increment_log_prob(normal_log(y, mu, sigma))
modelString = "
data {
int<lower=0> N ;
int y[N] ; // y is a length-N vector of itegers
}
parameters {
real<lower=0, uppper=1> theata ;
}
model {
theta ~ beta(1, 1) ;
// y ~ bernoulli(theta) ;
}
" # close quote for modelString
# 책에는 없는 코드. 실습을 위해 data 폴더에 넣은 코드를 실행시키기 위한 작업
cur_dir = getwd()
setwd(sprintf("%s/%s", cur_dir, 'data'))
# Example for Stan-Ydich-Xnom1subj-MbernBeta.R
#-------------------------------------------------------------------------------
# Optional generic preliminaries:
#graphics.off() # This closes all of R's graphics windows.
#rm(list=ls()) # Careful! This clears all of R's memory!
#-------------------------------------------------------------------------------
# Load The data
myData = read.csv("z15N50.csv")
myData
y | |
---|---|
1 | 0 |
2 | 1 |
3 | 0 |
4 | 0 |
5 | 0 |
6 | 0 |
7 | 0 |
8 | 0 |
9 | 0 |
10 | 0 |
11 | 1 |
12 | 0 |
13 | 0 |
14 | 0 |
15 | 1 |
16 | 1 |
17 | 1 |
18 | 0 |
19 | 0 |
20 | 1 |
21 | 0 |
22 | 0 |
23 | 0 |
24 | 0 |
25 | 1 |
26 | 1 |
27 | 0 |
28 | 0 |
29 | 0 |
30 | 0 |
31 | 0 |
32 | 1 |
33 | 0 |
34 | 1 |
35 | 0 |
36 | 1 |
37 | 0 |
38 | 0 |
39 | 0 |
40 | 1 |
41 | 0 |
42 | 0 |
43 | 1 |
44 | 1 |
45 | 0 |
46 | 1 |
47 | 0 |
48 | 0 |
49 | 0 |
50 | 0 |
#-------------------------------------------------------------------------------
# Load the functions genMCMC, smryMCMC, and plotMCMC:
source("Stan-Ydich-Xnom1subj-MbernBeta.R")
********************************************************************* Kruschke, J. K. (2015). Doing Bayesian Data Analysis, Second Edition: A Tutorial with R, JAGS, and Stan. Academic Press / Elsevier. *********************************************************************
#-------------------------------------------------------------------------------
# Optional: Specify filename root and graphical format for saving output.
# Otherwise specify as NULL or leave saveName and saveType arguments
# out of function calls.
fileNameRoot = "Stan-Ydich-Xnom1subj-MbernBeta-"
graphFileType = "png" #"eps"
#-------------------------------------------------------------------------------
# Generate the MCMC chain:
mcmcCoda = genMCMC( data=myData , numSavedSteps=10000 , saveName=fileNameRoot )
TRANSLATING MODEL 'modelString' FROM Stan CODE TO C++ CODE NOW. COMPILING THE C++ CODE FOR MODEL 'modelString' NOW. In file included from file33140a60ec6.cpp:8: In file included from /Library/Frameworks/R.framework/Versions/3.2/Resources/library/rstan/include//stansrc/stan/model/model_header.hpp:17: In file included from /Library/Frameworks/R.framework/Versions/3.2/Resources/library/rstan/include//stansrc/stan/agrad/rev.hpp:5: /Library/Frameworks/R.framework/Versions/3.2/Resources/library/rstan/include//stansrc/stan/agrad/rev/chainable.hpp:87:17: warning: 'static' function 'set_zero_all_adjoints' declared in header file should be declared 'static inline' [-Wunneeded-internal-declaration] static void set_zero_all_adjoints() { ^ In file included from file33140a60ec6.cpp:8: In file included from /Library/Frameworks/R.framework/Versions/3.2/Resources/library/rstan/include//stansrc/stan/model/model_header.hpp:21: /Library/Frameworks/R.framework/Versions/3.2/Resources/library/rstan/include//stansrc/stan/io/dump.hpp:26:14: warning: function 'product' is not needed and will not be emitted [-Wunneeded-internal-declaration] size_t product(std::vector<size_t> dims) { ^ 2 warnings generated. SAMPLING FOR MODEL 'modelString' NOW (CHAIN 1). Iteration: 1 / 10500 [ 0%] (Warmup) Iteration: 501 / 10500 [ 4%] (Sampling) Iteration: 1550 / 10500 [ 14%] (Sampling) Iteration: 2600 / 10500 [ 24%] (Sampling) Iteration: 3650 / 10500 [ 34%] (Sampling) Iteration: 4700 / 10500 [ 44%] (Sampling) Iteration: 5750 / 10500 [ 54%] (Sampling) Iteration: 6800 / 10500 [ 64%] (Sampling) Iteration: 7850 / 10500 [ 74%] (Sampling) Iteration: 8900 / 10500 [ 84%] (Sampling) Iteration: 9950 / 10500 [ 94%] (Sampling) Iteration: 10500 / 10500 [100%] (Sampling) # Elapsed Time: 0.006246 seconds (Warm-up) # 0.128609 seconds (Sampling) # 0.134855 seconds (Total) SAMPLING FOR MODEL 'modelString' NOW (CHAIN 2). Iteration: 1 / 10500 [ 0%] (Warmup) Iteration: 501 / 10500 [ 4%] (Sampling) Iteration: 1550 / 10500 [ 14%] (Sampling) Iteration: 2600 / 10500 [ 24%] (Sampling) Iteration: 3650 / 10500 [ 34%] (Sampling) Iteration: 4700 / 10500 [ 44%] (Sampling) Iteration: 5750 / 10500 [ 54%] (Sampling) Iteration: 6800 / 10500 [ 64%] (Sampling) Iteration: 7850 / 10500 [ 74%] (Sampling) Iteration: 8900 / 10500 [ 84%] (Sampling) Iteration: 9950 / 10500 [ 94%] (Sampling) Iteration: 10500 / 10500 [100%] (Sampling) # Elapsed Time: 0.009537 seconds (Warm-up) # 0.125926 seconds (Sampling) # 0.135463 seconds (Total) SAMPLING FOR MODEL 'modelString' NOW (CHAIN 3). Iteration: 1 / 10500 [ 0%] (Warmup) Iteration: 501 / 10500 [ 4%] (Sampling) Iteration: 1550 / 10500 [ 14%] (Sampling) Iteration: 2600 / 10500 [ 24%] (Sampling) Iteration: 3650 / 10500 [ 34%] (Sampling) Iteration: 4700 / 10500 [ 44%] (Sampling) Iteration: 5750 / 10500 [ 54%] (Sampling) Iteration: 6800 / 10500 [ 64%] (Sampling) Iteration: 7850 / 10500 [ 74%] (Sampling) Iteration: 8900 / 10500 [ 84%] (Sampling) Iteration: 9950 / 10500 [ 94%] (Sampling) Iteration: 10500 / 10500 [100%] (Sampling) # Elapsed Time: 0.006508 seconds (Warm-up) # 0.141121 seconds (Sampling) # 0.147629 seconds (Total) SAMPLING FOR MODEL 'modelString' NOW (CHAIN 4). Iteration: 1 / 10500 [ 0%] (Warmup) Iteration: 501 / 10500 [ 4%] (Sampling) Iteration: 1550 / 10500 [ 14%] (Sampling) Iteration: 2600 / 10500 [ 24%] (Sampling) Iteration: 3650 / 10500 [ 34%] (Sampling) Iteration: 4700 / 10500 [ 44%] (Sampling) Iteration: 5750 / 10500 [ 54%] (Sampling) Iteration: 6800 / 10500 [ 64%] (Sampling) Iteration: 7850 / 10500 [ 74%] (Sampling) Iteration: 8900 / 10500 [ 84%] (Sampling) Iteration: 9950 / 10500 [ 94%] (Sampling) Iteration: 10500 / 10500 [100%] (Sampling) # Elapsed Time: 0.010422 seconds (Warm-up) # 0.142794 seconds (Sampling) # 0.153216 seconds (Total)
#-------------------------------------------------------------------------------
# Display diagnostics of chain, for specified parameters:
parameterNames = varnames(mcmcCoda) # get all parameter names
for ( parName in parameterNames ) {
diagMCMC( mcmcCoda , parName=parName ,
saveName=fileNameRoot , saveType=graphFileType )
}
# 창으로 뜨는 것들을 없앤다.
graphics.off()
#-------------------------------------------------------------------------------
# Get summary statistics of chain:
summaryInfo = smryMCMC( mcmcCoda , compVal=0.5 , rope=c(0.45,0.55) ,
saveName=fileNameRoot )
# Display posterior information:
plotMCMC( mcmcCoda , data=myData , # compVal=0.5 , rope=c(0.45,0.55) ,
saveName=fileNameRoot , saveType=graphFileType )
Mean Median Mode ESS HDImass HDIlow HDIhigh CompVal theta 0.3078927 0.3053939 0.3109765 8995.5 0.95 0.1830674 0.4277429 0.5 PcntGtCompVal ROPElow ROPEhigh PcntLtROPE PcntInROPE PcntGtROPE theta 0.19 0.45 0.55 98.35 1.64 0.01
graphics.off()
#-------------------------------------------------------------------------------
# Use Stan display functions instead of DBDA2E functions:
# Load the stanFit object that was saved by genMCMC:
load("Stan-Ydich-Xnom1subj-MbernBeta-StanFit.Rdata")
# Display information:
show(stanFit)
openGraph()
traceplot(stanFit,pars=c("theta"))
openGraph()
plot(stanFit,pars=c("theta"))
Inference for Stan model: modelString. 4 chains, each with iter=10500; warmup=500; thin=4; post-warmup draws per chain=2500, total post-warmup draws=10000. mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat theta 0.31 0.00 0.06 0.19 0.27 0.31 0.35 0.44 9014 1 lp__ -32.59 0.01 0.71 -34.64 -32.73 -32.31 -32.14 -32.10 8875 1 Samples were drawn using NUTS(diag_e) at Wed Jul 8 18:40:05 2015. For each parameter, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat=1).
graphics.off()