rtdists_lba

library(rtdists)

set.seed(42)

# declare correct & error drift rate
correct_drift <- 2.4
error_drift <- 1 - correct_drift

# generate RTs with the diffusion model
rt1 <- rLBA(500, A = 0.5, b = 1, t0 = 0.5, 
            mean_v = c(correct_drift, 1 - correct_drift), 
            sd_v = c(1, 1))

# check accuracy
prop.table(table(rt1$response))


# declare the function to fit the LBA model (from rtdists package)
objective_fun <- function(par, rt, response, distribution = "norm") {
  # simple parameters
  spar <- par[!grepl("[12]$", names(par))]  
  
  # distribution parameters:
  dist_par_names <- unique(sub("[12]$", "", grep("[12]$" ,names(par), value = TRUE)))
  dist_par <- vector("list", length = length(dist_par_names))
  names(dist_par) <- dist_par_names
  for (i in dist_par_names) dist_par[[i]] <- as.list(unname(par[grep(i, names(par))]))
  dist_par$sd_v <- c(1, dist_par$sd_v) 
  
  # get summed log-likelihood:
  d <- do.call(dLBA, args = c(rt=list(rt), response=list(response), spar, dist_par, 
                              distribution=distribution, silent=TRUE))
  if (any(d == 0)) return(1e6)
  else return(-sum(log(d)))
}


#--- 
## search of parameter space

# variables to keep track of best objective fit & the best parameters
best_fit <- Inf
best_parms <- numeric(6)


# search 500 random starting points
for(i in 1:500){
  
  print(i)
  
  # initial parameters
  init_par <- c(runif(3, 0, 0.5), runif(3, 0.5, 2))
  
  # ensure b is larger than A
  init_par[2] <- init_par[1] + init_par[2]
  names(init_par) <- c("A", "b", "t0", "mean_v1", "mean_v2", "sd_v2")
  
  # run the model fit
  fit <- nlminb(objective_fun, start = init_par, rt=rt1$rt, response=rt1$response, lower = 0)
  
  # if this led to a better fit than the best so far, store parameters
  if(fit$objective < best_fit){
    best_fit <- fit$objective
    best_parms <- init_par
  }
  
}

# display the best-fitting parmeters
best_parms