R Help -
I am trying to use a grid search for a 2 free parameter reinforcement
learning model and the grid search is incredibly slow. I've used optimx but
can't seem to get reasonable answers. Is there a way to speed up this grid
search dramatically?
dat <- structure(list(choice = c(0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1,
1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1,
0, 1, 0, 1, 0, 1, 0,
0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1,
0, 0, 1, 0, 0, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0,
0, 1, 0, 0, 0, 0, 1,
1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1,
1, 0, 0, 0, 0, 0, 0,
1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1,
1, 0, 0, 0, 0, 0, 1,
1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1,
1, 0, 0, 1, 1, 0, 0,
0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1,
0, 0, 1, 0, 0, 0, 0,
1, 0, 1, 1, 1, 0), reward = c(0L, 0L, 0L,
0L, 1L, 1L, 0L, 0L,
1L, 0L, 0L,
0L, 0L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 1L,
1L, 0L, 1L,
0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 1L, 0L, 1L,
1L, 0L, 1L,
0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 1L,
0L, 0L, 1L,
1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 1L,
1L, 1L, 0L,
0L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L,
0L, 0L, 0L,
0L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 0L,
1L, 0L, 0L,
1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L,
0L, 1L, 0L,
0L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L,
0L, 1L, 0L,
1L, 1L, 0L, 1L, 1L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L,
0L, 0L, 1L,
0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 1L, 1L), RepNum = c(1L,
1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L)), .Names
= c("choice", "reward", "RepNum"), row.names = c(NA,
165L), class =
"data.frame")
CNTRACSID <- 0; subjectFit <- 0;
pLlist <- 0; pRlist <- 0; logLikelihood <- 0; trialProb <- 0;
hmmFunc <- function(delta, temperature){
pLlist = 1
pRlist = 1
block = 0
for (i in 1:length(dat$choice))
{
if (dat$RepNum[i] != block)
{
pL = 0.5
pR = 0.5
block = dat$RepNum[i]
}
# Markov Transitions
pL <- pL*(1-delta) + pR*delta
pR <- 1-pL
# Apply feedback
#denom <- p(F|L,C) * p(L) + p(F|R,C) * p(R)
pflc <- ifelse(dat$choice[i] == dat$reward[i], .8, .2)
pfrc <- 1 - pflc
denom <- pflc * pL + pfrc * pR
# What's the new belief given observation
posteriorL <- pflc * pL/denom
posteriorR <- 1-posteriorL
pL <- posteriorL
pR <- posteriorR
pL <- (1/(1 + exp(-temperature * (pL-.5))))
pR <- (1/(1 + exp(-temperature * (pR-.5))))
pLlist[i] = pL
pRlist[i] = pR
if(i > 1){
if(dat$choice[i] == 1){
trialProb[i] <- pLlist[i-1]
} else
{
trialProb[i] <- 1-pLlist[i-1]
}
}
else {
trialProb[1] <- .5
}
}
trialProb2 <- sum(log(trialProb))
subFit <- exp(trialProb2/length(dat$choice))
hmmOutput <- list("logLikelihood" = trialProb2, "subjectFit" = subFit,
"probabilities" = pLlist)
# print(hmmOutput$logLikelihood)
return(hmmOutput)
}
subjectFits <- 0; subLogLike <- 0; bestTemp <- 0; bestDelta= 0;
min = 0.001; max = .5; inc = 0.001;
deltaList = seq(min, max, inc)
mina = 0; maxa = 5; inca = .01
amList = seq(mina, maxa, inca)
maxLogValue <- -1000
for(delta in deltaList){
for(temp in amList){
probabilities <- hmmFunc(delta, temp)
if(probabilities$logLikelihood > maxLogValue){
pList <- probabilities$probabilities
maxLogValue <- probabilities$logLikelihood
subLogLike <- probabilities$logLikelihood
subjectFits <- probabilities$subjectFit
bestTemp <- temp
bestDelta <- delta
}
}
}
--
Edward H Patzelt | Clinical Science PhD Student
Psychology | Harvard University
Systems Neuroscience of Psychopathology Laboratory
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
