Setting parameter values

Let’s simulate 4 individuals. We have to define their behavioural rules and parameters. Let’s assume that behavioural rules for all 4 individuals are the same, where they employ the Rescorla-Wagner learning with the softmax choice policy, while assume that their parameters are different from each other.

# individual parameters
true_alpha <- c(0.3, 0.6, 0.3, 0.6) # learning rate (step size)
true_beta <- c(5, 5, 15, 15) # inverse temperature of the softmax function

Running the model

Now we have defined both task and individuals. It’s time to run the model.²

# Simulation
# The initial status of the four agents
# First, I prefer to define empty lists or matrices where data will be stored later
Q <- array(dim = c(num_options, horizon, group_size))
Q[,1,] <- Q_initial # Q_initial = 1.25
netChoiceProb <- array(dim = c(num_options, horizon, group_size))
netChoiceProb[,1,] <- 1/num_options
choices <- matrix(nrow=group_size, ncol=horizon)
payoffs <- matrix(nrow=group_size, ncol=horizon)
riskyChoiceProb <- matrix(nrow=group_size, ncol=horizon)

# looping from t = 1 to t = horizon
for (t in 1:horizon) {
  # each individual chooses one option based on his/her choice probability that is netChoiceProb[,t,]
  choices[,t] <- apply(netChoiceProb[,t,] # apply() allows us to apply the 1st element to the 'function()' repeatedly 
                       , MARGIN = 2
                       , function(netChoiceProbInstance){ 
                                sample(1:num_options
                                       , 1
                                       , prob=netChoiceProbInstance
                                       , replace=FALSE
                                ) 
                         }
                        )
  
  # each subject earns payoffs whose amount depends upon choices[,t]
  payoffs[,t] <- mapply(rnorm, 1, mean_list[choices[,t]], sd_list[choices[,t]]) # payoff drawn from `rnorm(mean, sd)`
  
  if (t < horizon) {
    # Rescorla-Wagner updating rule: Q_next <- Q + alpha * (payoff - Q)
    # The function `rescorla_wagner_updating()` is defined in the first chunk above
    Q <- rescorla_wagner_updating(t
                                 , numoptions
                                 , group_size
                                 , horizon
                                 , true_alpha
                                 , choices[,t]
                                 , payoffs[,t]
                                 , Q
                                 )
  
    
    
    # Softmax choice base solely on Q values
    # Calculate "exp( beta*Q_k )" for each option k
    Q_exp = ( Q[,t+1,] * rep(true_beta, each = num_options) ) %>% exp() # %>% is called a pipe operator (you can also use |> )
    # Then, calculating the softmax probability,
    # that is, exp( beta*Q_k )/(exp( beta*Q_1 ) + exp( beta*Q_2 )) for each slot
    softmaxMatrix = Q_exp %>% 
      apply(MARGIN=1, divideVector, denominator = apply(Q_exp, MARGIN=2, sum)) %>% 
      t()
    
    # Finally, we want to store values of softmax
    for (i in 1:group_size) {
      netChoiceProb[, t+1,i] <- softmaxMatrix[,i]
    }
  }
}
# Hooooray! A simulation run is done by now. 
# Next, we need to record this result
# First, let's summarize this run: 
for(i in 1:group_size) {
   riskyChoiceProb[i,] <- netChoiceProb[2,,i]
}

speudo_data <- data.frame(riskyChoiceProb = riskyChoiceProb) %>% 
  cbind(data.frame(choices = choices)) %>% 
  cbind(data.frame(payoffs = payoffs)) %>% 
  tidyr::pivot_longer(everything()
                      , names_to = c(".value", "trial")
                      , names_pattern = "([A-Za-z]+).(\\d+)"
                      )
speudo_data$indiv <- rep(1:group_size, each = horizon)
speudo_data$trial <- speudo_data$trial %>% as.character() %>% as.numeric() # trial was assumed to be a factor value

speudo_data %>% ggplot() + 
  geom_line(aes(trial, riskyChoiceProb, group = indiv), colour='grey')+ 
  geom_point(aes(trial, choices-1, colour = as.factor(choices))) + 
  geom_hline(yintercept = 0.5, linetype = 'dashed', colour = 'blue') + 
  scale_x_continuous(name="Trial", limits=c(1, horizon))  + 
  scale_color_viridis_d(name="Choice")  + 
  facet_wrap(~indiv) + 
  labs(y = 'Probability of choosing the risky arm', title = 'Simulated behaivour') +
  theme_test() 

knitr::kable(head(speudo_data, 12), "simple")

Calculating log-likelihood under a set of parameters

The next step is fitting the learning model with this pseudo data. Parameter estimation basically consists of the following two steps.

1. Setting a parameter value
2. Calculating log-likelihood

That’s it. Let me walk through them one by one.

alpha_test <- c(0.5, 0.5, 0.5, 0.5) # learning rate values 
beta_test <- c(5, 5, 5, 5) # inverse temperature
alpha_test

## [1] 0.5 0.5 0.5 0.5

beta_test

## [1] 5 5 5 5

# picking up one individual
this_individual <- 1
data_this_individual <- speudo_data %>% dplyr::filter(indiv == this_individual) # picking up the focal agent
this_choices <- data_this_individual$choices # this agent's choices observed
this_payoffs <- data_this_individual$payoffs # this agent's payoffs observed
this_alpha <- alpha_test[this_individual] #  this agent's alpha assumed
this_beta <- beta_test[this_individual] # this agent's alpha assumed

# initial setting
Q <- matrix(rep(NA, num_options * horizon), nrow = num_options)
logChoiceProb_estimated <- matrix(rep(NA, num_options * horizon), nrow = num_options)
Q[,1] <- Q_initial # Note: this initial value can also be a subject of estimation
logChoiceProb_estimated[,1] <- log(1/num_options) # random choice at the 1st trial
logLik <- rep(0, horizon)

# calculating log-likelihood
for (t in 1:horizon) {
  # recording the log likelihood of the choice made at this trial
  logLik[t] <- logChoiceProb_estimated[this_choices[t], t]
  
  # calculating the choice probability for the following trial
  if (t < horizon) {
    # copying Q values to the next trial
    Q[,t+1] <- Q[,t] 
    
    # Rescorla-Wagner Updating for the chosen option 
    Q[this_choices[t],t+1] <- 
      Q[this_choices[t],t] + 
      this_alpha * (this_payoffs[t] - Q[this_choices[t],t])
    
    # Softmax transformation (log based)
    logChoiceProb_estimated[,t+1] <- 
      this_beta * Q[, t+1] - log( exp(this_beta * Q[, t+1]) %>% sum() )
  }
}

# the sum of the negative log-likelihood
sum(-logLik)

## [1] 37.57145

# calculating log-likelihood
logLik_random <- rep(0, horizon)
for (t in 1:horizon) {
  # a probability to choose an option is 1/num_options 
  logLik_random[t] <- log(1/num_options) 
}

# the sum of the negative log-likelihood
sum(-logLik_random)

## [1] 69.31472

The brute force method of parameter estimation

Theoretically, you can always explore for the better combination by thoroughly exploring the entire (or reasonably wide) ranges of parameters, an approach called a brute force method. Let’s explore alpha for 0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1 and beta for 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20.³

# data storage
first_individual_logliks <- 
  data.frame(
    alpha = rep(seq(0, 1, 0.05), each = 21)
    , beta = rep(seq(0, 20, 1), 21)
    , logLik = rep(NA, 21^2)
  )

# The brute force method of parameter estimation
for (this_alpha in seq(0, 1, 0.05)) {
  for (this_beta in seq(0, 20, 1)) {
    # initial setting
    Q <-matrix(rep(NA, num_options * horizon), nrow = num_options)
    logChoiceProb_estimated <- 
      matrix(rep(NA, num_options * horizon), nrow = num_options)
    Q[,1] <- Q_initial 
    logChoiceProb_estimated[,1] <- log(1/num_options) 
    logLik <- rep(0, horizon)
    
    # calculating log-likelihood
    for (t in 1:horizon) {
      # recording the log likelihood of the choice made at this trial
      logLik[t] <- logChoiceProb_estimated[this_choices[t], t]
      
      # calculating the choice probability for the following trial
      if (t < horizon) {
        # copying Q values to the next trial
        Q[,t+1] <- Q[,t] 
        
        # Rescorla-Wagner Updating for the chosen option 
        Q[this_choices[t],t+1] <- 
          Q[this_choices[t],t] + 
          this_alpha * (this_payoffs[t] - Q[this_choices[t],t])
        
        # Softmax transformation (log based)
        logChoiceProb_estimated[,t+1] <- 
          this_beta * Q[, t+1] - log( exp(this_beta * Q[, t+1]) %>% sum() )
      }
    }
    
    # the sum of log-likelihood
    first_individual_logliks$logLik[which(first_individual_logliks$alpha == this_alpha & first_individual_logliks$beta == this_beta)] <- sum(logLik)
  }
}

# which combination of the parameter gave the maximum likelihood?
maximum_likelihood_position <- which(first_individual_logliks$logLik == (first_individual_logliks$logLik %>% max()))

best_fit_values <- first_individual_logliks[maximum_likelihood_position, ]

first_individual_logliks %>% ggplot() +
  geom_tile(aes(alpha, beta, fill = logLik)) +
  geom_point(data=data.frame(alpha=true_alpha[1], beta=true_beta[1]), mapping=aes(alpha,beta), size = 3, colour='red')+ # true value
  geom_text(data=data.frame(alpha=true_alpha[1], beta=true_beta[1]), mapping=aes(alpha,beta), hjust=-0.2, vjust=0.3, label="True value", colour = 'red') +
  geom_point(data=best_fit_values, mapping=aes(alpha,beta), size = 3, colour='blue')+ # best fitting value
  geom_text(data=best_fit_values, mapping=aes(alpha,beta), hjust=-0.2, vjust=0.3, colour='blue', label="MLE") +
  scale_fill_viridis_c('Log likelihood')+
  labs(title = 'The brute force parameter estimation for the agent 1\n(The red and blues points are the true and estimated values, respectively)')+
  theme_classic()

That’s it. How hard could it be? This maximum log likelihood method has been widely used, but often criticized for its unreliability. Thought the brute force method above gave a rough estimation which was more or less correct, it suggests a wide range of likely \(\beta\) values (i.e. the yellow-ish area ranges widely along the \(\beta\)-axis in the low \(\alpha\) area.).

This is why we need more sophisticated method to reliably infer the experimental data (e.g. a hierarchical Bayesian method).