fit_exgaussian.r

# iterate over starting parameters and assess density-----------------

do_fit <- function(current_data){

  # mu <- seq(from = 0.2, to = 1.8, by = 0.2)
  # sigma <- seq(from = 0.05, to = 0.65, by = 0.2)
  # tau <- seq(from = 0.05, to = 0.85, by = 0.2)
  
  mu <- 0.5
  sigma <- 0.05
  tau <- 0.10
  
  best_parms <- NULL
  best_fit <- 1e6
  
  for(i in 1:length(mu)){
    for(j in 1:length(sigma)){
      for(k in 1:length(tau)){
        
        current_fit <- optim(par = c(mu[i], 
                                     sigma[j], 
                                     tau[k]), 
                             fn = fit_exg,
                             data = current_data, 
                             min_parms = c(0.01, 0.01, 0.001), 
                             max_parms = c(Inf, Inf, Inf), 
                             control = (parscale = c(1, 0.1, 0.1)))
        
        if(current_fit$value < best_fit){
          best_fit <- current_fit$value
          best_parameters <- round(current_fit$par, 3)
        }
  
      }
    }
  }
  
  return(list(best_parameters = best_parameters, 
              best_fit = best_fit))
    
}

# return the density given a single set of parameter values------------------------
fit_exg <- function(data, 
                    parameters, 
                    min_parms, 
                    max_parms){
  
  
  if((min(parameters - min_parms) < 0) | 
     (min(max_parms - parameters) < 0)){
    return(.Machine$double.xmax)
  }
  
  # unpack starting parameters
  mu <- parameters[1]
  sigma <- parameters[2]
  tau <- parameters[3]
  
  # get the density of the rts with current parameter values
  d <- brms::dexgaussian(data, mu = mu, sigma = sigma, beta = tau)
  
  # return negative log
  d <- -sum(log(d))
  
  return(d)
}