The code BART 1.m or BART 1.R applies the model to data from a single subject, known as George, provided in the file GeorgeSober.txt. (우리 실습에서는 Jags 코드를 씀)
In [1]:
# 이 부분은 파이썬 커널로 변경해서 실행할 것
!head GeorgeSober.txt
In [15]:
# BART Model of Risky Decision-Making
BART_1_s = "
model{
# Optimal Number of Pumps
omega <- -gplus/log(1-p)
# Choice Data
for (j in 1:ntrials){
for (k in 1:options[j]){
theta[j,k] <- 1-(1/(1+max(-15,min(15,exp(beta*(k-omega))))))
d[j,k] ~ dbern(theta[j,k])
}
}
# Priors
gplus ~ dunif(0,10)
beta ~ dunif(0,10)
}
"
In [2]:
library(R2jags)
In [4]:
p <- .15 # (Belief of) bursting probability
ntrials <- 90 # Number of trials for the BART
In [5]:
Data <- matrix (data = as.numeric (as.matrix (read.table ("GeorgeSober.txt"))[-1,]), ntrials, 8)
head(Data)
Out[5]:
In [8]:
d <- matrix (, ntrials, 30) # Data in binary format
head(d)
Out[8]:
In [9]:
cash <-(Data[,7]!=0)*1 # Cash or burst?
npumps <- Data[,6] # Nr. of pumps
options <- cash + npumps # Nr. of decision possibilities
In [10]:
for (j in 1:ntrials)
{
if (npumps[j]>0) {d[j, 1:npumps[j]] <- rep (0, npumps[j])}
if (cash[j]==1) {d[j, (npumps[j]+1)] <- 1}
}
In [11]:
data <- list("ntrials", "p", "options", "d") # to be passed on to JAGS
In [12]:
myinits <- list(
list(gplus = 1.2, beta = 0.5))
In [13]:
# parameters to be monitored:
parameters <- c("gplus", "beta")
In [16]:
# The following command calls JAGS with specific options.
# For a detailed description see the R2jags documentation.
samples <- jags(data, inits=myinits, parameters,
textConnection(BART_1_s), #model.file = "BART_1.txt",
n.chains=1, n.iter=5000, n.burnin=2000, n.thin=1)
In [17]:
# Now the values for the monitored parameters are in the "samples" object,
# ready for inspection.
gplus <- samples$BUGSoutput$sims.list$gplus
beta <- samples$BUGSoutput$sims.list$beta
In [18]:
#################### PLOT RESULTS
par (cex.main = 2.5, cex.lab = 2, cex.axis = 1.5, mar = c(5, 5, 4, 0), las = 1)
layout (matrix (1:3, 1, 3))
hist(npumps, main = " ", xlab="Number of Pumps", ylab="Frequency", breaks=c(0:7),
col="lightblue", axes=F)
axis(2)
axis(1, at=seq(.5,6.5,by=1), labels = c("1", "2", "3", "4", "5", "6", "7"))
plot(density(gplus), xlab = expression (gamma^'+'),
main = " ", bty = 'n', lwd=2, ylab="Posterior Density")
plot(density(beta), xlab = expression (beta),
main = " ", bty = 'n', lwd=2, lab = c(5,3,5), ylab="Posterior Density")
In [19]:
library(rstan)
model <- "
// BART Model of Risky Decision-Making
data {
int<lower=1> ntrials;
real<lower=0,upper=1> p;
int<lower=1> options[ntrials];
int d[ntrials,30];
}
parameters {
// Priors
real<lower=0,upper=10> gplus;
real<lower=0,upper=10> beta;
}
transformed parameters {
real<lower=0> omega;
// Optimal Number of Pumps
omega <- -gplus / log1m(p); // log1m(p) equals log(1 - p), but faster
}
model {
// Choice Data
for (j in 1:ntrials) {
for (k in 1:options[j]) {
real theta;
theta <- 1 - inv_logit(-beta * (k - omega));
d[j,k] ~ bernoulli(theta);
}
}
}"
p <- .15 # (Belief of) bursting probability
ntrials <- 90 # Number of trials for the BART
Data <- matrix(data=as.numeric(as.matrix(read.table("GeorgeSober.txt"))[-1, ]),
ntrials, 8)
d <- matrix(-99, ntrials, 30) # Data in binary format
cash <-(Data[, 7] != 0) * 1 # Cash or burst?
npumps <- Data[, 6] # Nr. of pumps
options <- cash + npumps # Nr. of decision possibilities
for (j in 1:ntrials) {
if (npumps[j] > 0) {d[j, 1:npumps[j]] <- rep(0, npumps[j])}
if (cash[j] == 1) {d[j, (npumps[j] + 1)] <- 1}
}
# To be passed on to Stan:
data <- list(ntrials=ntrials, p=p, options=options, d=d)
myinits <- list(
list(gplus=1.2, beta=.5))
parameters <- c("gplus", "beta") # Parameters to be monitored
# The following command calls Stan with specific options.
# For a detailed description type "?rstan".
samples <- stan(model_code=model,
data=data,
init=myinits, # If not specified, gives random inits
pars=parameters,
iter=5000,
chains=1,
thin=1,
warmup=2000, # Stands for burn-in; Default = iter/2
# seed=123 # Setting seed; Default is random seed
)
# Now the values for the monitored parameters are in the "samples" object,
# ready for inspection.
gplus <- extract(samples)$gplus
beta <- extract(samples)$beta
#################### PLOT RESULTS
par (cex.main = 2.5, cex.lab = 2, cex.axis = 1.5, mar = c(5, 5, 4, 0), las = 1)
layout (matrix (1:3, 1, 3))
hist(npumps, main = " ", xlab="Number of Pumps", ylab="Frequency", breaks=c(0:7),
col="lightblue", axes=F)
axis(2)
axis(1, at=seq(.5,6.5,by=1), labels = c("1", "2", "3", "4", "5", "6", "7"))
plot(density(gplus), xlab = expression (gamma^'+'),
main = " ", bty = 'n', lwd=2, ylab="Posterior Density")
plot(density(beta), xlab = expression (beta),
main = " ", bty = 'n', lwd=2, lab = c(5,3,5), ylab="Posterior Density")
In [1]:
#### 입력 파일을 살펴보자. 이 부분은 파이썬 커널로 바꿔서 실행할 것.
!head GeorgeSober.txt GeorgeTipsy.txt GeorgeDrunk.txt
In [20]:
# Hierarchical BART Model of Risky Decision-Making
BART_2_s = "
model{
# Choice Data
for (i in 1:nconds){
gplus[i] ~ dnorm(mug,lambdag)T(0,)
beta[i] ~ dnorm(mub,lambdab)T(0,)
omega[i] <- -gplus[i]/log(1-p)
for (j in 1:ntrials){
for (k in 1:options[i,j]){
theta[i,j,k] <- 1-(1/(1+max(-15,min(15,exp(beta[i]*(k-omega[i]))))))
d[i,j,k] ~ dbern(theta[i,j,k])
}
}
}
# Priors
mug ~ dunif(0,10)
sigmag ~ dunif(0,10)
mub ~ dunif(0,10)
sigmab ~ dunif(0,10)
lambdag <- 1/pow(sigmag,2)
lambdab <- 1/pow(sigmab,2)
}
"
In [21]:
library(R2jags)
p <- .15 # (Belief of) bursting probability
ntrials <- 90 # Number of trials for the BART
In [22]:
################### READ IN THE DATA
Data <- list(
matrix (data = as.numeric (as.matrix (read.table ("GeorgeSober.txt"))[-1,]), ntrials, 8),
matrix (data = as.numeric (as.matrix (read.table ("GeorgeTipsy.txt"))[-1,]), ntrials, 8),
matrix (data = as.numeric (as.matrix (read.table ("GeorgeDrunk.txt"))[-1,]), ntrials, 8)
)
In [24]:
str(Data)
In [23]:
head(Data)
Out[23]:
In [26]:
nconds <- length(Data)
cash <- npumps <- matrix (, nconds, ntrials) # Cashes and nr. of pumps
d <- array (, c(nconds, ntrials, 30)) # Data in binary format
In [27]:
for (i in 1:nconds)
{
cash[i,] <- (Data[[i]][,7]!=0)*1 # Cash or burst?
npumps[i,] <- Data[[i]][,6] # Nr. of pumps
for (j in 1:ntrials)
{
if (npumps[i,j]>0) {d[i, j, 1:npumps[i,j]] <- rep(0, npumps[i,j])}
if (cash[i,j]==1) {d[i, j, (npumps[i,j]+1)] <- 1}
}
}
In [28]:
options <- cash + npumps # Nr. of decision possibilities
In [29]:
data <- list("nconds", "ntrials", "p", "options", "d") # to be passed on to JAGS
In [31]:
myinits <- list(
list(mug = 1.2, sigmag = 0.1, mub = 0.8, sigmab = 0.8),
list(mug = 1.5, sigmag = 0.2, mub = 1.0, sigmab = 1.2))
In [32]:
# parameters to be monitored:
parameters <- c("gplus", "beta", "mug", "sigmag", "mub", "sigmab")
In [33]:
# The following command calls JAGS with specific options.
# For a detailed description see the R2jags documentation.
samples <- jags(data, inits=myinits, parameters,
textConnection(BART_2_s), #model.file = "BART_2_jags.txt",
n.chains=2, n.iter=5000, n.burnin=2000, n.thin=1)
In [34]:
# Now the values for the monitored parameters are in the "samples" object,
# ready for inspection.
gplus.sober <- samples$BUGSoutput$sims.list$gplus[,1]
gplus.tipsy <- samples$BUGSoutput$sims.list$gplus[,2]
gplus.drunk <- samples$BUGSoutput$sims.list$gplus[,3]
beta.sober <- samples$BUGSoutput$sims.list$beta[,1]
beta.tipsy <- samples$BUGSoutput$sims.list$beta[,2]
beta.drunk <- samples$BUGSoutput$sims.list$beta[,3]
In [35]:
#################### PLOT SOME RESULTS
par (cex.main = 2.5, cex.lab = 2, cex.axis = 1.5, mar = c(5, 5, 4, 0), las = 1)
layout (matrix (1:9, 3, 3, byrow = T))
hist(npumps[1,], xlab=" ", main = "Sober: # pumps", breaks=c(0:max(npumps[1,])), xlim = range(c(0:9)), col="lightblue", axes=F)
axis(2)
axis(1, at=seq(.5,8.5,by=1), labels = c("1", "2", "3", "4", "5", "6", "7", "8", "9"))
plot(density(gplus.sober), xlab = expression (gamma^'+'),
main = expression (paste ("Sober: Posterior ", gamma^'+')), xlim = c(0.6,1.8), bty = 'n')
plot (density(beta.sober), xlab = expression (beta),
main = expression (paste ("Sober: Posterior ", beta)), xlim = c(0.2,1.4), bty = 'n')
hist(npumps[2,], xlab=" ", main = "Tipsy: # pumps", breaks=c(0:max(npumps[2,])), xlim = range(c(0:9)), col="lightblue", axes=F)
axis(2)
axis(1, at=seq(.5,8.5,by=1), labels = c("1", "2", "3", "4", "5", "6", "7", "8", "9"))
plot(density (gplus.tipsy), xlab = expression (gamma^'+'),
main = expression (paste ("Tipsy: Posterior ", gamma^'+')), xlim = c(0.6,1.8), bty = 'n')
plot(density (beta.tipsy), xlab = expression (beta),
main = expression (paste ("Tipsy: Posterior ", beta)), xlim = c(0.2,1.4), bty = 'n')
hist(npumps[3,], xlab="Number of Pumps", main = "Drunk: # pumps", breaks=c(0:max(npumps[3,])), xlim = range(c(0:9)), col="lightblue", axes=F)
axis(2)
axis(1, at=seq(.5,8.5,by=1), labels = c("1", "2", "3", "4", "5", "6", "7", "8", "9"))
plot(density(gplus.drunk), xlab = expression (gamma^'+'),
main = expression(paste ("Drunk: Posterior ", gamma^'+')), xlim = c(0.6,1.8), bty = 'n')
plot(density(beta.drunk), xlab = expression (beta),
main = expression(paste ("Drunk: Posterior ", beta)), xlim = c(0.2,1.4), bty = 'n')
In [36]:
library(rstan)
model <- "
// BART Model of Risky Decision-Making
data {
int<lower=1> nconds;
int<lower=1> ntrials;
real<lower=0,upper=1> p;
int<lower=1> options[nconds,ntrials];
int d[nconds,ntrials,30];
}
parameters {
vector<lower=0>[nconds] gplus;
vector<lower=0>[nconds] beta;
// Priors
real<lower=0,upper=10> mug;
real<lower=0,upper=10> sigmag;
real<lower=0,upper=10> mub;
real<lower=0,upper=10> sigmab;
}
transformed parameters {
vector<lower=0>[nconds] omega;
// Optimal Number of Pumps
omega <- -gplus / log1m(p);
}
model {
for (i in 1:nconds) {
gplus[i] ~ normal(mug, sigmag)T[0,];
beta[i] ~ normal(mub, sigmab)T[0,];
}
// Choice Data
for (i in 1:nconds) {
for (j in 1:ntrials) {
for (k in 1:options[i,j]) {
real theta;
theta <- 1 - inv_logit(-beta[i] * (k - omega[i]));
d[i,j,k] ~ bernoulli(theta);
}
}
}
}"
p <- .15 # (Belief of) bursting probability
ntrials <- 90 # Number of trials for the BART
################### READ IN THE DATA
Data <- list(
matrix(data=as.numeric(as.matrix(read.table("GeorgeSober.txt"))[-1, ]),
ntrials, 8),
matrix(data=as.numeric(as.matrix(read.table("GeorgeTipsy.txt"))[-1, ]),
ntrials, 8),
matrix(data=as.numeric(as.matrix(read.table("GeorgeDrunk.txt"))[-1, ]),
ntrials, 8)
)
nconds <- length(Data)
cash <- npumps <- matrix (, nconds, ntrials) # Cashes and nr. of pumps
d <- array (-99, c(nconds, ntrials, 30)) # Data in binary format
for (i in 1:nconds) {
cash[i, ] <- (Data[[i]][, 7] != 0) * 1 # Cash or burst?
npumps[i, ] <- Data[[i]][, 6] # Nr. of pumps
for (j in 1:ntrials) {
if (npumps[i, j] > 0) {d[i, j, 1:npumps[i, j]] <- rep(0, npumps[i, j])}
if (cash[i, j] == 1) {d[i, j, (npumps[i, j] + 1)] <- 1}
}
}
options <- cash + npumps # Nr. of decision possibilities
# To be passed on to Stan:
data <- list(nconds=nconds, ntrials=ntrials, p=p, options=options, d=d)
myinits <- list(
list(gplus=rep(1.2, 3), beta=rep(.5, 3), mug=1.2, sigmag=.1, mub=.8, sigmab=.8),
list(gplus=rep(1.2, 3), beta=rep(.5, 3), mug=1.5, sigmag=.2, mub=1, sigmab=1.2))
# Parameters to be monitored:
parameters <- c("beta", "gplus", "mub", "mug", "sigmab", "sigmag")
# The following command calls Stan with specific options.
# For a detailed description type "?rstan".
samples <- stan(model_code=model,
data=data,
init=myinits, # If not specified, gives random inits
pars=parameters,
iter=1500,
chains=2,
thin=1,
warmup=500, # Stands for burn-in; Default = iter/2
seed=1234 # Setting seed; Default is random seed
)
# Now the values for the monitored parameters are in the "samples" object,
# ready for inspection.
gplus.sober <- extract(samples)$gplus[,1]
gplus.tipsy <- extract(samples)$gplus[,2]
gplus.drunk <- extract(samples)$gplus[,3]
beta.sober <- extract(samples)$beta[,1]
beta.tipsy <- extract(samples)$beta[,2]
beta.drunk <- extract(samples)$beta[,3]
#################### PLOT SOME RESULTS
windows()
par (cex.main = 2.5, cex.lab = 2, cex.axis = 1.5, mar = c(5, 5, 4, 0), las = 1)
layout (matrix (1:9, 3, 3, byrow = T))
hist(npumps[1,], xlab=" ", main = "Sober: # pumps", breaks=c(0:max(npumps[1,])),
xlim = range(c(0:9)), col="lightblue", axes=F)
axis(2)
axis(1, at=seq(.5,8.5,by=1),
labels = c("1", "2", "3", "4", "5", "6", "7", "8", "9"))
plot(density(gplus.sober), xlab = expression (gamma^'+'), xlim = c(0.6,1.8),
main = expression (paste ("Sober: Posterior ", gamma^'+')), bty = 'n')
plot (density(beta.sober), xlab = expression (beta), bty = 'n',
main = expression (paste ("Sober: Posterior ", beta)), xlim = c(0.2,1.4))
hist(npumps[2,], xlab=" ", main = "Tipsy: # pumps", breaks=c(0:max(npumps[2,])),
xlim = range(c(0:9)), col="lightblue", axes=F)
axis(2)
axis(1, at=seq(.5,8.5,by=1),
labels = c("1", "2", "3", "4", "5", "6", "7", "8", "9"))
plot(density (gplus.tipsy), xlab = expression (gamma^'+'), xlim = c(0.6,1.8),
main = expression (paste ("Tipsy: Posterior ", gamma^'+')), bty = 'n')
plot(density (beta.tipsy), xlab = expression (beta), xlim = c(0.2,1.4),
main = expression (paste ("Tipsy: Posterior ", beta)), bty = 'n')
hist(npumps[3,], xlab="Number of Pumps", main = "Drunk: # pumps",
breaks=c(0:max(npumps[3,])), xlim = range(c(0:9)), col="lightblue", axes=F)
axis(2)
axis(1, at=seq(.5,8.5,by=1),
labels = c("1", "2", "3", "4", "5", "6", "7", "8", "9"))
plot(density(gplus.drunk), xlab = expression (gamma^'+'), xlim = c(0.6,1.8),
main = expression(paste ("Drunk: Posterior ", gamma^'+')), bty = 'n')
plot(density(beta.drunk), xlab = expression (beta), xlim = c(0.2,1.4),
main = expression(paste ("Drunk: Posterior ", beta)), bty = 'n')