Commit 1ec5ef9e authored by linushof's avatar linushof
Browse files

Loop for MCMC samples

parent c8e55515
model
{
# Define priors
# define priors
alpha.phi ~ dnorm(0,1)
gamma.phi ~ dnorm(0,1)
......@@ -15,27 +15,36 @@ model
# Trial loop
for (i in min_i:max_i)
for (i in start:stop)
{
v.x.a[i] <- pow(x_A[i],alpha)
v.y.a[i] <- pow(y_A[i],alpha)
v.x.b[i] <- pow(x_B[i],alpha)
v.y.b[i] <- pow(y_B[i],alpha)
w.x.a[i] <- (delta * (pow(px_A[i],gamma))) / (delta * (pow(px_A[i], gamma)) + pow(py_A[i],gamma))
w.y.a[i] <- 1-w.x.a[i]
# value function
w.x.b[i] <- (delta * (pow(px_B[i],gamma))) / (delta * (pow(px_B[i], gamma)) + pow(py_B[i],gamma))
w.y.b[i] <- 1-w.x.b[i]
v.a.o1[i] <- pow(a_o1[i], alpha)
v.a.o2[i] <- pow(a_o2[i], alpha)
v.b.o1[i] <- pow(b_o1[i], alpha)
v.b.o2[i] <- pow(b_o2[i], alpha)
Vf.a[i] <- w.x.a[i] * v.x.a[i] + w.y.a[i] * v.y.a[i]
Vf.b[i] <- w.x.b[i] * v.x.b[i] + w.y.b[i] * v.y.b[i]
# weighting function
w.a.p2[i] <- (delta * (pow(a_p2[i], gamma))) / (delta * (pow(a_p2[i], gamma)) + pow(a_p1[i], gamma))
w.a.p1[i] <- 1-w.a.p2[i]
w.b.p1[i] <- (delta * (pow(b_p1[i], gamma))) / (delta * (pow(b_p1[i], gamma)) + pow(b_p2[i], gamma))
w.b.p2[i] <- 1-w.b.p1[i]
# valuation
Vf.a[i] <- w.a.p1[i] * v.a.o1[i] + w.a.p2[i] * v.a.o2[i]
Vf.b[i] <- w.b.p1[i] * v.b.o1[i] + w.b.p2[i] * v.b.o2[i]
## retransform
# retransform
Vf.a.re[i] = pow(Vf.a[i], (1/alpha))
Vf.b.re[i] = pow(Vf.b[i], (1/alpha))
# choice rule
binval[i] <- (1)/(1+exp((-1*rho)*(Vf.a.re[i]-Vf.b.re[i])))
resp[i] ~ dbern(binval[i])
choice[i] ~ dbern(binval[i])
}
}
pacman::p_load(tidyverse,
R2jags)
# read choice data
cols <- list(.default = col_double(),
strategy = col_factor(),
boundary = col_factor(),
gamble = col_factor(),
rare = col_factor(),
agent = col_factor(),
choice = col_factor())
choices <- read_csv("data/choices/choices.csv", col_types = cols)
# prepare data for JAGS
choices_cpt <- choices %>%
filter(!(is.na(a_ev_exp) | is.na(b_ev_exp))) %>% # remove choices where prospects were not attended
mutate(choice_A = if_else(choice == "A", 1, 0), # for logit choice rule
i = row_number()) # assign trial numbers
## prepare computation of CPT estimates for each distinct strategy-parameter combination
params_sim <- choices_cpt %>% distinct(strategy, s, boundary, a) # group trials of distinct strategy-parameter combinations
choices_grouped <- vector("list", nrow(params_sim))
for(set in seq_len(nrow(params_sim))){
choices_grouped[[set]] <- choices_cpt %>%
filter(strategy == theta[[set, "strategy"]] & s == theta[[set, "s"]] & boundary == theta[[set, "boundary"]] & a == theta[[set, "a"]])
}
estimates_cpt <- vector("list", nrow(params_sim)) # allocate space for JAGS output
# compute estimates of CPT parameters
params_cpt <- c("alpha", # outcome sensitivity
"gamma", # curvature of w(p)
"delta", # elevation of w(p)
"rho") # choice sensitivity
n_chains <- 2 # number of markov chains
for(set in seq_len(nrow(params_sim))){
# get trial data of current strategy parameter combination
trials <- list(choice = choices_grouped[[set]]$choice_A,
a_o1 = choices_grouped[[set]]$a_o1,
a_o2 = choices_grouped[[set]]$a_o2, # higher risky outcome (x)
b_o1 = choices_grouped[[set]]$b_o1, # safe outcome (x)
b_o2 = choices_grouped[[set]]$b_o2,
a_p1 = choices_grouped[[set]]$a_p1,
a_p2 = choices_grouped[[set]]$a_p2, # probability of higher risky outcome (x)
b_p1 = choices_grouped[[set]]$b_p1,
b_p2 = choices_grouped[[set]]$b_p2, # probability of safe outcome
start = min(choices_grouped[[set]]$i),
stop = max(choices_grouped[[set]]$i))
# MCMC samples
MCMC_samples <- jags.parallel(data = trials,
inits = NULL,
parameters.to.save = params_cpt,
model.file = "JAGS/cpt_trial_level.txt",
n.chains = n_chains,
n.iter = 10000,
n.burnin = 5000,
n.thin = 1,
n.cluster = n_chains,
DIC = FALSE,
jags.seed = 8362)
# to get posterior estimates, credibility intervals, and MCMC diagnostics
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment