YEAB (version 1.0.6)

gaussian_fit: Gaussian + ramp fit with LM algorithm

Description

Gaussian + ramp fit with LM algorithm

Usage

gaussian_fit(
  responses,
  time_vec,
  par = list(a = 0.1, d = 0.1, t0 = 18, b = 10, c = 1),
  max.iter = 500
)

Value

a numeric vector of coefficients, the same as the par argument

Arguments

responses

numeric, vector of response or response rate

time_vec

numeric, time bins

par

a list of parameters for the gaussian + linear; see Buhusi, C. V., Perera, D., & Meck, W. H. (2005) for an explanation

max.iter

numeric, max number of iterations

Details

Ver Buhusi, C. V., Perera, D., & Meck, W. H. (2005). Memory for timing visual and auditory signals in albino and pigmented rats.

This algorithm uses the nonlinear least squares nls.lm (Levenberg–Marquardt) from the minpack.lm package

Examples

# Function to create synthetic data
g_plus_lin <- function(par, time) {
  par$a * exp(-0.5 * ((time - par$t0) / par$b)**2) + par$c * (time - par$t0) + par$d
}
# real params
pars <- list(a = 20, t0 = 20, b = 10, c = 0.2, d = 1)
# time vector for simulation
ti <- seq(0, 60, 0.1)
# time vector for sampling with 2 sec of resolution
ti_data <- seq(0, 60, 2)
# r(t) real
y_curva <- g_plus_lin(par = pars, ti)
# r(t) sampled with noise
y_data <- g_plus_lin(par = pars, ti_data) + rnorm(length(ti_data), 0, sd = 2)
# param estimation
par_est <- gaussian_fit(responses = y_data, t = ti_data, par = pars, max.iter = 10500)
par_est
# fitted curve
y_hat <- g_plus_lin(par_est |> as.list(), ti)
# plot results

plot(ti,
  y_curva,
  type = "l",
  col = "blue",
  lwd = 2,
  ylim = c(0, max(y_curva, y_data)),
  xlab = "Time in trial",
  ylab = "R(t)",
)
points(ti_data, y_data, pch = 21, bg = "red", cex = 1.2)
lines(
  ti,
  y_hat,
  col = "green2",
  lwd = 2
)
legend(
  "topright",
  legend = c("real", "real + noise", "ajuste nls.lm"),
  lty = c(1, 0, 1),
  pch = c(NA, 21),
  pt.bg = c(NA, "red"),
  col = c("blue", 1, "green2"),
  pt.cex = 0.9,
  cex = 0.6
)