diff --git a/app.R b/app.R
new file mode 100644
index 0000000..2ceb602
--- /dev/null
+++ b/app.R
@@ -0,0 +1,199 @@
+#
+# This is a Shiny web application. You can run the application by clicking
+# the 'Run App' button above.
+#
+# Find out more about building applications with Shiny here:
+#
+# http://shiny.rstudio.com/
+#
+
+library(shiny)
+library(tidyverse)
+library(wesanderson)
+
+# Define UI for application that draws a histogram
+ui <- fluidPage(
+
+ # Application title
+ titlePanel("Associative Learning Simulation"),
+
+ withMathJax(),
+
+ # Sidebar with a slider input for number of bins
+ sidebarLayout(
+ sidebarPanel(
+ h3('Set parameter levels'),
+ sliderInput("xi",
+ HTML("\\(\\xi\\) (noise)"),
+ min = -3,
+ max = 3,
+ value = 0,
+ step=.05),
+ sliderInput("ep",
+ HTML("\\(\\epsilon\\) (learning rate)"),
+ min = -3,
+ max = 3,
+ value = 0,
+ step=.05),
+ # sliderInput("b",
+ # "Bias",
+ # min = -2,
+ # max = 2,
+ # value = 0,
+ # step=.05),
+ sliderInput("rho",
+ HTML("\\(\\rho\\) (inverse temperature)"),
+ min = -3,
+ max = 3,
+ value = 0,
+ step=.05),
+ h3('Parameter optimality across 40 random runs'),
+ p('When noise parameter is set very low'),
+ img(src = 'optimality_plot.png', style = 'max-width: 100%')
+ ),
+
+ # Show a plot of the generated distribution
+ mainPanel(
+ plotOutput("trialsPlot"),
+ helpText("RW Equation: \n
+ $$p(a_t|s_t) = \\text{logit}^{-1}\\Big(Q_{t-1}(a_{t},s_{t}) + \\epsilon\\big(\\rho r_{t} - Q_{t-1}(a_{t},s_{t})\\big)\\Big)\\cdot (1-\\xi) + \\frac{\\xi}{2}$$"),
+ uiOutput('rw_eq')
+ # ,
+ # tableOutput('runTable')
+ )
+ )
+)
+
+# Define server logic required to draw a histogram
+server <- function(input, output) {
+
+ generateTrials <- function() {
+ p_right <- data.frame(expand.grid(cue=1:2, reward = c(1,5)), pcorrect_if_pressed_r=c(rep(.2,1), rep(.8,1)))
+
+ cue1Indxs <- sample(c(1,3), size = 60, replace = T)
+ cue2Indxs <- sample(c(2,4), size = 60, replace = T)
+
+ manyTrialIndxs <- c(cue1Indxs,cue2Indxs)
+
+ manyTrialIndxsShuffled <- manyTrialIndxs[sample(1:length(manyTrialIndxs),
+ size = length(manyTrialIndxs),
+ replace = F)]
+
+ Trials <- p_right[manyTrialIndxsShuffled,]
+ Trials$crct_if_right <- rbinom(dim(Trials)[1], size = 1, prob = Trials$pcorrect_if_pressed_r)
+ Trials$outcome_r <- Trials$crct_if_right*Trials$reward
+ Trials$outcome_l <- (1-Trials$crct_if_right)*Trials$reward
+ return(Trials)
+ }
+
+ inv_logit <- function(x) exp(x)/(1+exp(x))
+ Phi_approx <- function(x) pnorm(x)
+
+ rw_strategy <- function(trialdf, mu_p){
+ xi <- Phi_approx( mu_p[1])# + sigma[1] * xi_pr[i] )
+ ep <- Phi_approx( mu_p[2])# + sigma[2] * ep_pr[i] )
+ b <- mu_p[3]# + sigma[3] * b_pr; # vectorization
+ rho <- exp( mu_p[4])# + sigma[4] * rho_pr );
+
+ K <- length(unique(trialdf$cue))
+ Tsubj <- dim(trialdf)[1]
+ wv_g <- c(rep(0, K)) # action wegith for go
+ wv_ng <- c(rep(0, K)) # action wegith for nogo
+ qv_g <- c(rep(0, K)) # Q value for go
+ qv_ng <- c(rep(0, K)) # Q value for nogo
+ pGo <- c(rep(0, K)) # prob of go (press)
+
+ trialdf$pressed_r <- NA
+ trialdf$Qgo <- NA
+ trialdf$Qnogo <- NA
+ trialdf$Wgo <- NA
+ trialdf$Wnogo <- NA
+ trialdf$pGo <- NA
+ trialdf$outcome <- NA
+
+ for (t in 1:Tsubj) {
+ wv_g[ trialdf$cue[t] ] <- qv_g[ trialdf$cue[t] ] + b
+ wv_ng[ trialdf$cue[t] ] <- qv_ng[ trialdf$cue[t] ] # qv_ng is always equal to wv_ng (regardless of action)
+ pGo[ trialdf$cue[t] ] = inv_logit( wv_g[ trialdf$cue[t] ] - wv_ng[ trialdf$cue[t] ] )
+ pGo[ trialdf$cue[t] ] = pGo[ trialdf$cue[t] ] * (1 - xi) + xi/2; # noise
+
+ trialdf$pressed_r[t] <- rbinom(n = 1, size = 1, prob = , pGo[ trialdf$cue[t] ]);
+
+ trialdf$Qgo[t] <- qv_g[ trialdf$cue[t] ];
+ trialdf$Qnogo[t] <- qv_ng[ trialdf$cue[t] ];
+ trialdf$Wgo[t] <- wv_g[ trialdf$cue[t] ];
+ trialdf$Wnogo[t] <- wv_ng[ trialdf$cue[t] ];
+ trialdf$pGo[t] <- pGo[ trialdf$cue[t] ];
+
+ # update action values
+ if(trialdf$pressed_r[t] == 1){
+ qv_g[ trialdf$cue[t] ] <- qv_g[ trialdf$cue[t] ] + ep * (rho * trialdf$outcome_r[t] - qv_g[ trialdf$cue[t] ]);
+ trialdf$outcome[t] <- trialdf$outcome_r[t]
+ } else {
+ qv_ng[ trialdf$cue[t] ] <- qv_ng[ trialdf$cue[t] ] + ep * (rho * trialdf$outcome_l[t] - qv_ng[ trialdf$cue[t] ]);
+ trialdf$outcome[t] <- trialdf$outcome_l[t]
+ }
+ } # end of t loop
+ return(trialdf)
+ }
+
+ plot_RW_run <- function(trials, mu_p){
+ single_run <- rw_strategy(trialdf = trials,
+ mu_p = mu_p)
+
+ aplot <- single_run %>%
+ mutate(cue = factor(cue)) %>%
+ group_by(cue) %>%
+ mutate(t = 1:n(), last_outcome = as.numeric( ifelse(lag(pressed_r) == 1 & lag(outcome) == 5, 1,
+ ifelse(lag(pressed_r) == 1 & lag(outcome) == 1, .95,
+ ifelse(lag(pressed_r) == 1 & lag(outcome) == 0, .1,
+ ifelse(lag(pressed_r) == 0 & lag(outcome) == 5, 0,
+ ifelse(lag(pressed_r) == 0 & lag(outcome) == 1, .05,
+ .9)))))),
+ last_press = lag(pressed_r)) %>%
+ ggplot(aes(x = t, y = pGo)) +
+ geom_line(alpha = .1) +
+ geom_line(stat = 'smooth', method = 'gam', formula = y ~ s(x, k = 15, bs = "cr"), alpha = .5) +
+ geom_segment(aes(xend = t, yend = last_outcome), alpha = .1, color = 'black') +
+ geom_point(aes(y = last_outcome, shape = factor(last_press))) +
+ scale_shape_manual(values = c(25,24), name = 'Last press was...', breaks = c(1,0), labels = c('Right', 'left')) +
+ scale_y_continuous(breaks = c(0,.5, 1), labels = c('left', '', 'right'))+
+ geom_point() +
+ facet_wrap(~cue, nrow = 2)+
+ theme(panel.background = element_blank(),
+ # strip.text = element_blank(),
+ strip.background = element_rect(fill = '#eeeeee'))+
+ labs(y = "More likely to press...", x = 'Trial number')
+ return(list(plot = aplot, runData = single_run))
+ }
+
+ someTrials <- generateTrials()
+
+ simulatedTrials <- reactive({
+ plot_RW_run(trials = someTrials,
+ mu_p = c(xi = input$xi, ep = input$ep, b = 0, rho = input$rho))
+ })
+
+ output$trialsPlot <- renderPlot({
+ simulatedTrials()$plot
+ })
+
+ output$rw_eq <- renderUI({
+ xi <- Phi_approx( input$xi )
+ ep <- Phi_approx( input$ep )
+ b <- input$b
+ rho <- exp( input$rho )
+
+ withMathJax(sprintf("RW Equation with transformed values: \n
+ $$p(a_t|s_t) = \\text{logit}^{-1}\\Big(Q_{t-1}(a_{t},s_{t}) + %.02f\\big(%.02f r_{t} - Q_{t-1}(a_{t},s_{t})\\big)\\Big)\\cdot (1-%.02f) + \\frac{%.02f}{2}$$",
+ ep, rho, xi, xi))
+ })
+
+ output$runTable <- renderTable({
+ simulatedTrials()$runData
+ })
+}
+
+# Run the application
+shinyApp(ui = ui, server = server)
+
diff --git a/rsconnect/shinyapps.io/jflournoy/rw_model.dcf b/rsconnect/shinyapps.io/jflournoy/rw_model.dcf
new file mode 100644
index 0000000..419d5a5
--- /dev/null
+++ b/rsconnect/shinyapps.io/jflournoy/rw_model.dcf
@@ -0,0 +1,10 @@
+name: rw_model
+title: rw_model
+account: jflournoy
+server: shinyapps.io
+appId: 188785
+bundleId: 909163
+url: https://jflournoy.shinyapps.io/rw_model/
+when: 1500486751.98579
+asMultiple: FALSE
+asStatic: FALSE
diff --git a/www/optimality_plot.png b/www/optimality_plot.png
new file mode 100644
index 0000000..8a3233c
Binary files /dev/null and b/www/optimality_plot.png differ