R/single-lba.r

In rtdists: Response Time Distributions

#' Single accumulator of linear ballistic accumulator (LBA)
#' 
#' Density, distribution function, and random generation for a single accumulator of the LBA model with the following parameters: \code{A} (upper value of starting point), \code{b} (response threshold), \code{t0} (non-decision time), and driftrate (\code{v}). All functions are available with different distributions underlying the drift rate: Normal (\code{norm}), Gamma (\code{gamma}), Frechet (\code{frechet}), and log normal (\code{lnorm}).
#' 
#' @param rt a vector of RTs.
#' @param n desired number of observations (scalar integer).
#' @param A start point interval or evidence in accumulator before beginning of decision process. Start point varies from trial to trial in the interval [0, \code{A}] (uniform distribution). Average amount of evidence before evidence accumulation across trials is \code{A}/2.
#' @param b response threshold. (\code{b} - \code{A}/2) is a measure of "response caution". 
#' @param t0 non-decision time or response time constant (in seconds). Lower bound for the duration of all non-decisional processes (encoding and response execution).
#' @param st0 variability of non-decision time, such that \code{t0} is uniformly distributed between \code{t0} and \code{t0} + \code{st0}. Only available in random number generation functions \code{rlba_}.
#' 
#' @param mean_v,sd_v mean and standard deviation of normal distribution for drift rate (\code{norm}). See \code{\link{Normal}}
#' @param shape_v,rate_v,scale_v shape, rate, and scale of gamma (\code{gamma}) and scale and shape of Frechet (\code{frechet}) distributions for drift rate. See \code{\link{GammaDist}} or \code{\link[evd]{frechet}}. For Gamma, scale = 1/shape and shape = 1/scale.
#' @param meanlog_v,sdlog_v mean and standard deviation of lognormal distribution on the log scale for drift rate (\code{lnorm}). See \code{\link{Lognormal}}.
#' 
#' @param posdrift logical. Should driftrates be forced to be positive? Default is \code{TRUE}. (Uses truncated normal for random generation).
#' @param robust logical. Should robust normal distributions be used for \code{norm} and \code{lnorm}? Can be helpful in rare cases but is approximately three times slower than the non-robust versions. Default is \code{FALSE}.
#' 
#' 
#' @details These functions are mainly for internal purposes. We do not recommend to use them. Use the high-level functions described in \code{/link{LBA}} instead.
#' 
#' @return All functions starting with a \code{d} return the density (PDF), all functions starting with \code{p} return the distribution function (CDF), and all functions starting with \code{r} return random response times and responses (in a \code{matrix}).
#' 
#' @note Density (i.e., \code{dlba_}), distribution (i.e., \code{plba_}), and random derivative (i.e., \code{rlba_}) functions are vectorized for all parameters (i.e., in case parameters are not of the same length as \code{rt}, parameters are recycled). Furthermore, the random derivative functions also accept a matrix of length \code{n} in which each column corresponds to a accumulator specific value (see \code{\link{rLBA}} for a more user-friendly way).
#' 
#' @references 
#' 
#' Brown, S. D., & Heathcote, A. (2008). The simplest complete model of choice response time: Linear ballistic accumulation. \emph{Cognitive Psychology}, 57(3), 153-178. doi:10.1016/j.cogpsych.2007.12.002
#' 
#' Donkin, C., Averell, L., Brown, S., & Heathcote, A. (2009). Getting more from accuracy and response time data: Methods for fitting the linear ballistic accumulator. \emph{Behavior Research Methods}, 41(4), 1095-1110. doi:10.3758/BRM.41.4.1095
#' 
#' Heathcote, A., & Love, J. (2012). Linear deterministic accumulator models of simple choice. \emph{Frontiers in Psychology}, 3, 292. doi:10.3389/fpsyg.2012.00292
#' 
#' @importFrom evd rfrechet dfrechet pfrechet
#' @importFrom msm rtnorm
#' @importFrom gsl gamma_inc
#' @importFrom stats dgamma dlnorm dnorm pgamma plnorm pnorm rgamma rlnorm rnorm runif
#' 
#' 
#' @name single-LBA
#' 
#' @example examples/examples.slba.R
#' 
NULL

# protected normal desity and cdf
pnormP <- function(x,mean=0,sd=1,lower.tail=TRUE) ifelse(abs(x)<7,pnorm(x, mean=mean, sd=sd,lower.tail=lower.tail),ifelse(x<0,0,1))
dnormP <- function(x,mean=0,sd=1) ifelse(abs(x)<7,dnorm(x,mean=mean,sd=sd),0)

make_r <- function(drifts, n,b,A,n_v,t0,st0=0) {
  drifts <- drifts[1:n,]
  drifts[drifts<0] <- 0
  if (is.null(dim(A))) starts <- matrix(runif(min=0,max=A,n=n*n_v),ncol=n_v,byrow=TRUE)
  else starts <- apply(A, c(1,2), function(x) runif(min=0, max = x, 1))
  if (is.null(dim(b))) ttf <- t((b-t(starts)))/drifts
  else ttf <- (b-starts)/drifts
  rt <- apply(ttf+t0+runif(min=0,max=st0,n=n),1,min)
  resp <- apply(ttf+t0,1,which.min)
  bad <- !is.finite(rt)
  if (any(bad)) {
    warning(paste(sum(bad),"infinite RTs removed and less than", n, "rts returned"))
    resp <- resp[!bad]
    rt <- rt[!bad]
  }
  cbind(rt=rt,response=resp)
}

rem_t0 <- function(rt, t0) pmax(rt - t0, 0)

check_single_arg <- function(...) {
  mc <- match.call()
  vars <- all.vars(mc)
  arguments <- list(...)
  for(i in seq_along(arguments)) {
    if (length(arguments[[i]]) != 1) stop(paste(vars[i], "needs to be of length 1!"))
    if (!is.numeric(arguments[[i]]) | !is.finite(arguments[[i]])) stop(paste(vars[i], "needs to be numeric and finite!"))
  }
}

check_vector <- function(...) {
  mc <- match.call()
  vars <- all.vars(mc)
  dots <- list(...)
  for (i in seq_along(dots)) {
    if ((vars[i] == "rt") && (any(dots[[i]] < 0))) stop("rt needs to contain only positive values.") 
    if (!is.vector(dots[[i]], "numeric")) stop(paste(vars[[i]], "needs to be a numeric vector!"))
    if (length(dots[[i]]) < 1) stop(paste(vars[[i]], "needs to have a length >= 1."))
  }
}

error_message_b_smaller_A <- "b cannot be smaller than A!"

####### Normal:

#' @rdname single-LBA
#' @export dlba_norm
dlba_norm <- function(rt,A,b, t0, mean_v, sd_v, posdrift=TRUE, robust = FALSE) {
  #check_single_arg(A=A, b=b, t0=t0, mean_v=mean_v, sd_v=sd_v)
  check_vector(rt, A, b, t0, mean_v, sd_v)
  # bring all arguments to length of rt
  nn <- length(rt)
  A <- rep(A, length.out = nn)
  b <- rep(b, length.out = nn)
  t0 <- rep(t0, length.out = nn)
  mean_v <- rep(mean_v, length.out = nn)
  sd_v <- rep(sd_v, length.out = nn)
  if (any(b < A)) stop(error_message_b_smaller_A)
  dlba_norm_core(rt = rt, A = A, b = b, t0 = t0, mean_v = mean_v, sd_v = sd_v, posdrift = posdrift, robust = robust, nn = nn)
}

## this functions expects all arguments to have the samel length (which is nn)
dlba_norm_core <- function(rt,A,b, t0, mean_v, sd_v, posdrift=TRUE, robust = FALSE, nn) {
  if (robust) { # robust == TRUE uses robust versions of the normal distributions
    pnorm1 <- pnormP
    dnorm1 <- dnormP
  } else {
    pnorm1 <- pnorm
    dnorm1 <- dnorm
  }
  rt <- rem_t0(rt, t0) # rmove t0 from rt
  if (posdrift) denom <- pmax(pnorm1(mean_v/sd_v),1e-10) else denom <- rep(1, nn)
  
  if (any(A<1e-10, na.rm = TRUE)) {
    # for A<1e-10 save results in out_A
    A_small <- A<1e-10
    out_A <- pmax(0, ((b[A_small]/rt[A_small]^2)*dnorm1(b[A_small]/rt[A_small],mean_v[A_small],sd=sd_v[A_small]))/denom[A_small], na.rm = TRUE) 
    # calculate other results into out_o
    zs <- rt[!A_small]*sd_v[!A_small]
    zu <- rt[!A_small]*mean_v[!A_small]
    chiminuszu <- b[!A_small]-zu
    chizu <- chiminuszu/zs
    chizumax <- (chiminuszu-A[!A_small])/zs
    out_o <- pmax(0,(mean_v[!A_small]*(pnorm1(chizu)-pnorm1(chizumax)) + sd_v[!A_small]*(dnorm1(chizumax)-dnorm1(chizu)))/(A[!A_small]*denom[!A_small]))
    # combine out_A and out_o
    out <- numeric(nn)
    out[!A_small] <- out_o
    out[A_small] <- out_A
    return(out) 
  } else {
    zs <- rt*sd_v
    zu <- rt*mean_v
    chiminuszu <- b-zu
    chizu <- chiminuszu/zs
    chizumax <- (chiminuszu-A)/zs
    return(pmax(0,(mean_v*(pnorm1(chizu)-pnorm1(chizumax)) + sd_v*(dnorm1(chizumax)-dnorm1(chizu)))/(A*denom), na.rm=TRUE)) 
  }
}

#' @rdname single-LBA
#' @export plba_norm
plba_norm <- function(rt,A,b,t0,mean_v, sd_v,posdrift=TRUE, robust = FALSE) {
  check_vector(rt, A, b, t0, mean_v, sd_v)
  nn <- length(rt)
  A <- rep(A, length.out = nn)
  b <- rep(b, length.out = nn)
  t0 <- rep(t0, length.out = nn)
  mean_v <- rep(mean_v, length.out = nn)
  sd_v <- rep(sd_v, length.out = nn)
  if (any(b < A)) stop(error_message_b_smaller_A)
  plba_norm_core(rt = rt, A = A, b = b, t0 = t0, mean_v = mean_v, sd_v = sd_v, posdrift = posdrift, robust = robust, nn = nn)
}

plba_norm_core <- function(rt,A,b,t0,mean_v, sd_v,posdrift=TRUE, robust = FALSE, nn) {
  if (robust) { # robust == TRUE uses robust versions of the normal distributions
    pnorm1 <- pnormP
    dnorm1 <- dnormP
  } else {
    pnorm1 <- pnorm
    dnorm1 <- dnorm  
  }

  rt <- rem_t0(rt, t0)
  if (posdrift) denom <- pmax(pnorm1(mean_v/sd_v),1e-10) else denom <- 1
  
  if (any(A<1e-10, na.rm = TRUE)) {
    # for A<1e-10 save results in out_A
    A_small <- A<1e-10
    out_A <- pmin(1, pmax(0, (pnorm1(b[A_small]/rt[A_small],mean=mean_v[A_small],sd=sd_v[A_small],lower.tail=FALSE))/denom[A_small], na.rm=TRUE))

    # calculate other results into out_o
    zs <- rt[!A_small]*sd_v[!A_small]
    zu <- rt[!A_small]*mean_v[!A_small]
    chiminuszu <- b[!A_small]-zu
    xx <- chiminuszu-A[!A_small]
    chizu <- chiminuszu/zs
    chizumax <- xx/zs
    tmp1 <- zs*(dnorm1(chizumax)-dnorm1(chizu))
    tmp2 <- xx*pnorm1(chizumax)-chiminuszu*pnorm1(chizu)
    out_o <- pmin(pmax(0,(1+(tmp1+tmp2)/A[!A_small])/denom[!A_small]), 1)
    
    # combine out_A and out_o
    out <- numeric(nn)
    out[!A_small] <- out_o
    out[A_small] <- out_A
    return(out)
  } else {
    zs <- rt*sd_v
    zu <- rt*mean_v
    chiminuszu <- b-zu
    xx <- chiminuszu-A
    chizu <- chiminuszu/zs
    chizumax <- xx/zs
    tmp1 <- zs*(dnorm1(chizumax)-dnorm1(chizu))
    tmp2 <- xx*pnorm1(chizumax)-chiminuszu*pnorm1(chizu)
    return(pmin(pmax(0,(1+(tmp1+tmp2)/A)/denom, na.rm=TRUE), 1))
  }
}

#' @rdname single-LBA
#' @export rlba_norm
rlba_norm <- function(n,A,b,t0,mean_v, sd_v, st0=0,posdrift=TRUE) {
  #check_single_arg(n, A, b, t0, st0)
  check_single_arg(n)
  if (any(b < A)) stop(error_message_b_smaller_A)
  n_v <- max(length(mean_v), length(sd_v))
  if (posdrift) drifts <- matrix(rtnorm(n=n*n_v, mean=mean_v, sd=sd_v, lower=0),ncol=n_v,byrow=TRUE)  
  else drifts <- matrix(rnorm(n=n*n_v, mean=mean_v, sd=sd_v),ncol=n_v,byrow=TRUE)
  make_r(drifts=drifts, n=n, b=b,A=A, n_v=n_v, t0=t0, st0=st0)
}


####### Gamma:

#' @rdname single-LBA
#' @export dlba_gamma
dlba_gamma <- function(rt,A,b,t0,shape_v,rate_v, scale_v) {
  
  if (!missing(rate_v) && !missing(scale_v)) stop("specify 'rate_v' or 'scale_v', but not both")
  if (missing(rate_v)) rate_v <- 1/scale_v
  
  check_vector(rt, A, b=b, t0, shape_v, rate_v)
  
  nn <- length(rt)
  A <- rep(A, length.out = nn)
  b <- rep(b, length.out = nn)
  t0 <- rep(t0, length.out = nn)
  shape_v <- rep(shape_v, length.out = nn)
  rate_v <- rep(rate_v, length.out = nn)
  if (any(b < A)) stop(error_message_b_smaller_A)
  
  dlba_gamma_core(rt=rt,A=A,b=b,t0=t0, shape_v=shape_v, rate_v=rate_v, nn=nn)
  
}


dlba_gamma_core <- function(rt,A,b,t0,shape_v, rate_v, nn) {
  rt <- rem_t0(rt, t0)
  min <- (b-A)/rt
  max <- b/rt
  
  
  if (any(A<1e-10, na.rm = TRUE)) {
    # for A<1e-10 save results in out_A
    A_small<- A<1e-10
    
    out_A<-  pmax(0, (b[A_small]/rt[A_small]^2)*(dgamma(b[A_small]/rt[A_small],shape=shape_v[A_small],rate=rate_v[A_small])), na.rm = TRUE)
    
    min <- (b[!A_small]-A[!A_small])/rt[!A_small]
    max <- b[!A_small]/rt[!A_small]
    
    Gmax <- pgamma(max, shape_v[!A_small], rate=rate_v[!A_small])
    Gmin <- pgamma(min, shape_v[!A_small], rate=rate_v[!A_small])
    Gmax2 <- pgamma(max, (shape_v[!A_small]+1), rate=rate_v[!A_small])
    Gmin2 <- pgamma(min, (shape_v[!A_small]+1), rate=rate_v[!A_small])
    zgamma <- ( ((Gmax2-Gmin2)*gamma(shape_v[!A_small]+1))/((Gmax-Gmin)*rate_v[!A_small]*gamma(shape_v[!A_small])) )
    
    diffG <- function(rt,point,shape_v, rate_v) {
      (-point/(rt^2))*dgamma(point/rt,shape_v[!A_small],rate = rate_v)
    } #NB:point refers to the constants b OR b-A.
    u <- (Gmax2-Gmin2)
    v <- (Gmax-Gmin)
    udash <- (diffG(rt[!A_small], b[!A_small], shape_v[!A_small]+1, rate_v[!A_small])- diffG(rt[!A_small], (b[!A_small]-A[!A_small]), shape_v[!A_small]+1, rate_v[!A_small]))
    vdash <- (diffG(rt[!A_small], b[!A_small], shape_v[!A_small], rate_v[!A_small])- diffG(rt[!A_small], (b[!A_small]-A[!A_small]), shape_v[!A_small], rate_v[!A_small]))
    const <- gamma(shape_v[!A_small]+1)/(rate_v[!A_small]*gamma(shape_v[!A_small]))
    diffzgamma <- ((udash*v - vdash*u)/(v^2))*const #quotient rule
    term1 <- (Gmax - Gmin)*(zgamma + (rt*diffzgamma))
    term2 <- diffG(rt,b,shape_v[!A_small],rate_v[!A_small])*((zgamma*rt)-b)
    term3 <- diffG(rt,(b-A),shape_v[!A_small],rate_v[!A_small])*(b-A-(zgamma*rt))
    out_o <- ((term1+term2+term3)/A)
    out_o[!is.finite(out_o)] <- 0 # Set NaN or -Inf or Inf to pdf=0
    
    out<- numeric(nn)
    out[!A_small] <- out_o
    out[A_small] <- out_A
    
    return(pmax(0, out))
  } else{
    
    Gmax <- pgamma(max, shape_v, rate=rate_v)
    Gmin <- pgamma(min, shape_v, rate=rate_v)
    Gmax2 <- pgamma(max, (shape_v+1), rate=rate_v)
    Gmin2 <- pgamma(min, (shape_v+1), rate=rate_v)
    zgamma <- ( ((Gmax2-Gmin2)*gamma(shape_v+1))/((Gmax-Gmin)*rate_v*gamma(shape_v)) )
    
    diffG <- function(rt,point,shape_v, rate_v) {
      (-point/(rt^2))*dgamma(point/rt,shape_v,rate = rate_v)
    } #NB:point refers to the constants b OR b-A.
    u <- (Gmax2-Gmin2)
    v <- (Gmax-Gmin)
    udash <- (diffG(rt, b, shape_v+1, rate_v)- diffG(rt, (b-A), shape_v+1, rate_v))
    vdash <- (diffG(rt, b, shape_v, rate_v)- diffG(rt, (b-A), shape_v, rate_v))
    const <- gamma(shape_v+1)/(rate_v*gamma(shape_v))
    diffzgamma <- ((udash*v - vdash*u)/(v^2))*const #quotient rule
    term1 <- (Gmax - Gmin)*(zgamma + (rt*diffzgamma))
    term2 <- diffG(rt,b,shape_v,rate_v)*((zgamma*rt)-b)
    term3 <- diffG(rt,(b-A),shape_v,rate_v)*(b-A-(zgamma*rt))
    out.value <- ((term1+term2+term3)/A)
    out.value[!is.finite(out.value)] <- 0 # Set NaN or -Inf or Inf to pdf=0
    return(pmax(0, out.value))
  }
}


#' @rdname single-LBA
#' @export plba_gamma  
plba_gamma <- function(rt,A,b,t0,shape_v, rate_v, scale_v) {
  if (!missing(rate_v) && !missing(scale_v)) stop("specify 'rate_v' or 'scale_v', but not both")
  if (missing(rate_v)) rate_v <- 1/scale_v
  
  check_vector(rt, A, b=b, t0, shape_v, rate_v)
  
  nn <- length(rt)
  A <- rep(A, length.out = nn)
  b <- rep(b, length.out = nn)
  t0 <- rep(t0, length.out = nn)
  shape_v <- rep(shape_v, length.out = nn)
  rate_v <- rep(rate_v, length.out = nn)
  if (any(b < A)) stop(error_message_b_smaller_A)
  
  plba_gamma_core(rt=rt,A=A,b=b,t0=t0,shape_v=shape_v, rate_v=rate_v, nn=nn)
}

plba_gamma_core <- function(rt,A,b,t0,shape_v, rate_v, nn) {
  
  rt <- rem_t0(rt, t0)
  if (any(A<1e-10, na.rm = TRUE)) {
    # for A<1e-10 save results in out_A
    A_small<- A<1e-10
    
    out_A <- pmin(1, pmax(0, (pgamma(b[A_small]/rt[A_small],shape=shape_v[A_small],rate=rate_v[A_small],lower.tail=FALSE)), na.rm=TRUE))
    
    min <- (b[!A_small]-A[!A_small])/rt[!A_small]
    max <- b[!A_small]/rt[!A_small]
    Gmax <- pgamma(max, shape_v[!A_small], rate=rate_v[!A_small])
    Gmin <- pgamma(min, shape_v[!A_small], rate=rate_v[!A_small])
    Gmax2 <- pgamma(max, (shape_v[!A_small]+1), rate=rate_v[!A_small])
    Gmin2 <- pgamma(min, (shape_v[!A_small]+1), rate=rate_v[!A_small])
    zgamma <- ((Gmax2-Gmin2)*gamma(shape_v[!A_small]+1))/((Gmax-Gmin)*rate_v[!A_small]*gamma(shape_v[!A_small]))
    
    term1 <- ((rt*zgamma) - b)/A
    term2 <- (b-A-(rt*zgamma))/A
    pmax <- pgamma(max, shape_v[!A_small], rate = rate_v[!A_small])
    pmin <- pgamma(min, shape_v[!A_small], rate = rate_v[!A_small])
    out_o <- (1 + pmax*term1 + pmin*term2)
    out_o[rt==Inf] <- 1 # term1=Inf and term2=-Inf cancel in this case
    out_o[!is.finite(out_o)] <- 0 # Set NaN or -Inf to CDF=0
    
    
    out<- numeric(nn)
    out[!A_small] <- out_o
    out[A_small] <- out_A
    return(pmin(pmax(0, out), 1))
    
  } else {
    min <- (b-A)/rt
    max <- b/rt
    Gmax <- pgamma(max, shape_v, rate=rate_v)
    Gmin <- pgamma(min, shape_v, rate=rate_v)
    Gmax2 <- pgamma(max, (shape_v+1), rate=rate_v)
    Gmin2 <- pgamma(min, (shape_v+1), rate=rate_v)
    zgamma <- ((Gmax2-Gmin2)*gamma(shape_v+1))/((Gmax-Gmin)*rate_v*gamma(shape_v)) 
    
    term1 <- ((rt*zgamma) - b)/A
    term2 <- (b-A-(rt*zgamma))/A
    pmax <- pgamma(max, shape_v, rate = rate_v)
    pmin <- pgamma(min, shape_v, rate = rate_v)
    out.value <- (1 + pmax*term1 + pmin*term2)
    out.value[rt==Inf] <- 1 # term1=Inf and term2=-Inf cancel in this case
    out.value[!is.finite(out.value)] <- 0 # Set NaN or -Inf to CDF=0
    return(pmin(pmax(0, out.value), 1))
  }
}

#' @rdname single-LBA
#' @export rlba_gamma
rlba_gamma <- function(n,A,b,t0,shape_v, rate_v, scale_v, st0=0) {
  check_single_arg(n)
  if (any(b < A)) stop(error_message_b_smaller_A)
  if (!missing(rate_v) && !missing(scale_v)) stop("specify 'rate_v' or 'scale_v', but not both")
  if (missing(rate_v)) rate_v <- 1/scale_v
  n_v <- max(length(shape_v), length(rate_v))  
  drifts <- matrix(rgamma(n=n*n_v,shape = shape_v,rate = rate_v),ncol=n_v,byrow=TRUE)
  make_r(drifts=drifts, n=n, b=b,A=A, n_v = n_v, t0=t0, st0=st0)
}


####### Frechet:

#' @rdname single-LBA
#' @export dlba_frechet
dlba_frechet <- function(rt,A,b,t0,shape_v, scale_v) {

  check_vector(rt, A, b, t0, shape_v, scale_v)
  nn <- length(rt)
  A <- rep(A, length.out = nn)
  b <- rep(b, length.out = nn)
  t0 <- rep(t0, length.out = nn)
  shape_v <- rep(shape_v, length.out = nn)
  scale_v <- rep(scale_v, length.out = nn)
  if (any(b < A)) stop(error_message_b_smaller_A)
  
  dlba_frechet_core(rt=rt,A=A,b=b,t0=t0,shape_v=shape_v, scale_v=scale_v, nn=nn)
}

dlba_frechet_core <- function(rt,A,b,t0,shape_v, scale_v, nn) {
  
  rt <- rem_t0(rt, t0)
  
  ps <- cbind(b, b-A, scale_v,shape_v)
  ps_below_zero <- apply(ps, 1, function(x) any(x <= 0))
  
  # rt <- pmax(rt,0) #not needed, see rem_t0
  t_old <- rt
  
  out <- numeric(nn)
  
  if (sum(!ps_below_zero) > 0) {
    rt <- rt[!ps_below_zero]
    A <- A[!ps_below_zero]
    b <- b[!ps_below_zero]
    t0 <- t0[!ps_below_zero]
    shape_v <- shape_v[!ps_below_zero]
    scale_v <- scale_v[!ps_below_zero]
  
    min <- (b-A)/rt
    max <- b/rt
    Gmax <- pfrechet(max, loc=0, scale=scale_v, shape=shape_v)
    Gmin <- pfrechet(min, loc=0, scale=scale_v, shape=shape_v)
    D <- Gmax - Gmin
    gam <- gamma_inc(1-(1/shape_v), (1/scale_v*max)^(-shape_v))-gamma_inc(1-(1/shape_v), (1/scale_v*min)^(-shape_v))
    zfrechet <- gam/(1/scale_v*D)
    diffG1 <- ((-b/(rt^2))*dfrechet(b/rt, loc=0, scale=scale_v, shape=shape_v))
    diffG2 <- ((-(b-A)/(rt^2))*dfrechet((b-A)/rt, loc=0, scale=scale_v, shape=shape_v))    
    diffD <- diffG1 - diffG2    
    diffgam <- (-shape_v*(((1/scale_v*b)^(-shape_v+1))/(rt^(-shape_v+2)))*exp(-(1/scale_v*b/rt)^(-shape_v))) - (-shape_v*(((1/scale_v*(b-A))^(-shape_v+1))/(rt^(-shape_v+2)))*exp(-(1/scale_v*(b-A)/rt)^(-shape_v)))
    diffzfrechet <- ((1/scale_v)^(-1))*(((-D^(-2))*diffD)*gam + (diffgam*(D^(-1))))
    term1 <- (Gmax - Gmin)*(zfrechet + (rt*diffzfrechet))
    term2 <- diffG1*((zfrechet*rt)-b)
    term3 <- diffG2*(b-A-(zfrechet*rt))
    out.value <- ((term1+term2+term3)/A)
    out.value[!is.finite(out.value)] <- 0 # Set NaN or -Inf or Inf to pdf=0
    out[!ps_below_zero] <- out.value
  }
  return(pmax(0, out))
}

#' @rdname single-LBA
#' @export plba_frechet
plba_frechet <- function(rt,A,b,t0,shape_v, scale_v) {
  check_vector(rt, A, b, t0, shape_v, scale_v)
  nn <- length(rt)
  A <- rep(A, length.out = nn)
  b <- rep(b, length.out = nn)
  t0 <- rep(t0, length.out = nn)
  shape_v <- rep(shape_v, length.out = nn)
  scale_v <- rep(scale_v, length.out = nn)
  if (any(b < A)) stop(error_message_b_smaller_A)
  
  plba_frechet_core(rt=rt,A=A,b=b,t0=t0,shape_v=shape_v, scale_v=scale_v, nn=nn)
}

plba_frechet_core <- function(rt,A,b,t0,shape_v, scale_v, nn) {  
  rt <- rem_t0(rt, t0)
  
  ps <- cbind(b, b-A, scale_v,shape_v)
  ps_below_zero <- apply(ps, 1, function(x) any(x <= 0))
  
  # rt <- pmax(rt,0) #not needed, see rem_t0
  t_old <- rt
  
  out <- numeric(nn)
  
  if (sum(!ps_below_zero) > 0) {
    
    rt <- rt[!ps_below_zero]
    A <- A[!ps_below_zero]
    b <- b[!ps_below_zero]
    t0 <- t0[!ps_below_zero]
    shape_v <- shape_v[!ps_below_zero]
    scale_v <- scale_v[!ps_below_zero]
    
    # rt <- pmax(rt,0) #not needed, see rem_t0
    min <- (b-A)/rt
    max <- b/rt
    pmax <- pfrechet(max, loc=0, scale=scale_v, shape=shape_v)
    pmin <- pfrechet(min, loc=0, scale=scale_v, shape=shape_v)
    zfrechet <- (gamma_inc(1-(1/shape_v),(1/scale_v*max)^(-shape_v))-gamma_inc(1-(1/shape_v),(1/scale_v*min)^(-shape_v)))/(1/scale_v*(pmax-pmin))    
    term1 <- ((rt*zfrechet) - b)/A
    term2 <- (b-A-(rt*zfrechet))/A 
    out.value <- (1 + pmax*term1 + pmin*term2)
    out.value[rt==Inf] <- 1 # term1=Inf and term2=-Inf cancel in this case
    out.value[!is.finite(out.value)] <- 0 # Set NaN or -Inf to CDF=0
    out[!ps_below_zero] <- out.value
  }
  
  return(pmin(pmax(0, out), 1))
}


#' @rdname single-LBA
#' @export rlba_frechet
rlba_frechet <- function(n,A,b,t0,shape_v, scale_v,st0=0){
  check_single_arg(n)
  if (any(b < A)) stop(error_message_b_smaller_A)
  n_v <- max(length(shape_v), length(scale_v))
  drifts <- matrix(rfrechet(n=n*n_v, loc=0, scale=scale_v, shape=shape_v),ncol=n_v,byrow=TRUE)
  
  make_r(drifts=drifts, n=n, b=b,A=A, n_v=n_v, t0=t0, st0=st0)
}


####### Log-Normal:

#' @rdname single-LBA
#' @export dlba_lnorm
dlba_lnorm <- function(rt,A,b,t0,meanlog_v, sdlog_v, robust = FALSE) {
  check_vector(rt, A, b, t0, meanlog_v, sdlog_v)
  nn <- length(rt)
  A <- rep(A, length.out = nn)
  b <- rep(b, length.out = nn)
  t0 <- rep(t0, length.out = nn)
  meanlog_v <- rep(meanlog_v, length.out = nn)
  sdlog_v <- rep(sdlog_v, length.out = nn)
  if (any(b < A)) stop(error_message_b_smaller_A)
  
  dlba_lnorm_core(rt=rt,A=A,b=b,t0=t0,meanlog_v=meanlog_v, sdlog_v=sdlog_v, robust = robust, nn=nn)
}

dlba_lnorm_core <- function(rt,A,b,t0,meanlog_v, sdlog_v, robust=FALSE, nn) {
  if (robust) { # robust == TRUE uses robust versions of the normal distributions
    pnorm1 <- pnormP
    dnorm1 <- dnormP
  } else {
    pnorm1 <- pnorm
    dnorm1 <- dnorm  
  }
  
  rt <- rem_t0(rt, t0)
  
  if (any(A<1e-10, na.rm = TRUE)) {
    # for A<1e-10 save results in out_A
    A_small<- A<1e-10
    out_A <- pmax(0, (b[A_small]/rt[A_small]^2)*(dlnorm(b[A_small]/rt[A_small],meanlog_v[A_small],sdlog=sdlog_v[A_small])), na.rm = TRUE)
    
    # calculate other results into out_o
    #Should there also be a check that A_small has at least one FALSE to bother running this code? -Angus; Probably not necessary, Henrik.
    min <- (b-A)/rt
    max <- b/rt
    
    zlognorm <- (exp(meanlog_v[!A_small]+(sdlog_v[!A_small]^2)/2)*(pnorm1((log(max)-meanlog_v[!A_small]-(sdlog_v[!A_small]^2))/sdlog_v[!A_small])-pnorm1((log(max)-meanlog_v[!A_small]-(sdlog_v[!A_small]^2))/sdlog_v[!A_small]))) / (pnorm1((log(max)-meanlog_v[!A_small])/sdlog_v[!A_small])-pnorm1((log(max)-meanlog_v[!A_small])/sdlog_v[!A_small]))
    Gmax <- plnorm(max,meanlog=meanlog_v[!A_small],sdlog=sdlog_v[!A_small])
    Gmin <- plnorm(max,meanlog=meanlog_v[!A_small],sdlog=sdlog_v[!A_small])
    
    u <- (pnorm1((log(max)-meanlog_v[!A_small]-(sdlog_v[!A_small])^2)/sdlog_v[!A_small])-pnorm1((log(max)-meanlog_v[!A_small]-(sdlog_v[!A_small])^2)/sdlog_v[!A_small]))
    v <- (pnorm1((log(max)-meanlog_v[!A_small])/sdlog_v[!A_small])-pnorm1((log(max)-meanlog_v[!A_small])/sdlog_v[!A_small]))
    
    udash <- (((-1/(sdlog_v[!A_small]*rt[!A_small]))*dnorm1((log(b[!A_small]/rt[!A_small])-meanlog_v[!A_small]-(sdlog_v[!A_small])^2)/sdlog_v[!A_small])) - ((-1/(sdlog_v[!A_small]*rt[!A_small]))*dnorm1((log((b[!A_small]-A[!A_small])/rt[!A_small])-meanlog_v[!A_small]-(sdlog_v[!A_small])^2)/sdlog_v[!A_small])))
    vdash <- (((-1/(sdlog_v[!A_small]*rt[!A_small]))*dnorm1((log(b[!A_small]/rt[!A_small])-meanlog_v[!A_small])/sdlog_v[!A_small])) - ((-1/(sdlog_v[!A_small]*rt[!A_small]))*dnorm1((log((b[!A_small]-A[!A_small])/rt[!A_small])-meanlog_v[!A_small])/sdlog_v[!A_small])))
    const <- exp(meanlog_v[!A_small]+((sdlog_v[!A_small])^2)/2)
    
    diffzlognorm <- ((udash*v - vdash*u)/(v^2))*const #quotient rule
    term1 <- (Gmax - Gmin)*(zlognorm + (rt[!A_small]*diffzlognorm))
    term2 <- ((-b[!A_small]/(rt[!A_small]^2))*dlnorm(b[!A_small]/rt[!A_small],meanlog=meanlog_v[!A_small],sdlog=sdlog_v[!A_small]))*((zlognorm*rt[!A_small])-b[!A_small])
    term3 <- (b[!A_small]-A[!A_small]-(zlognorm*rt[!A_small]))*((-(b[!A_small]-A[!A_small])/(rt[!A_small]^2))*dlnorm((b[!A_small]-A[!A_small])/rt[!A_small],meanlog=meanlog_v[!A_small],sdlog=sdlog_v[!A_small]))
    out_o<- ((term1+term2+term3)/A[!A_small])
    out_o[!is.finite(out_o)] <- 0 # Set NaN or -Inf or Inf to pdf=0
    
    out<- numeric(nn)
    out[!A_small] <- out_o
    out[A_small] <- out_A
    
    return(pmax(0, out))
    
  } else{
    
    min <- (b-A)/rt
    max <- b/rt
    
    zlognorm <- (exp(meanlog_v+(sdlog_v^2)/2)*(pnorm1((log(max)-meanlog_v-(sdlog_v^2))/sdlog_v)-pnorm1((log(min)-meanlog_v-(sdlog_v^2))/sdlog_v))) / (pnorm1((log(max)-meanlog_v)/sdlog_v)-pnorm1((log(min)-meanlog_v)/sdlog_v))
    Gmax <- plnorm(max,meanlog=meanlog_v,sdlog=sdlog_v) 
    Gmin <- plnorm(min,meanlog=meanlog_v,sdlog=sdlog_v)
    
    u <- (pnorm1((log(max)-meanlog_v-(sdlog_v)^2)/sdlog_v)-pnorm1((log(min)-meanlog_v-(sdlog_v)^2)/sdlog_v))
    v <- (pnorm1((log(max)-meanlog_v)/sdlog_v)-pnorm1((log(min)-meanlog_v)/sdlog_v))
    
    udash <- (((-1/(sdlog_v*rt))*dnorm1((log(b/rt)-meanlog_v-(sdlog_v)^2)/sdlog_v)) - ((-1/(sdlog_v*rt))*dnorm1((log((b-A)/rt)-meanlog_v-(sdlog_v)^2)/sdlog_v)))
    vdash <- (((-1/(sdlog_v*rt))*dnorm1((log(b/rt)-meanlog_v)/sdlog_v)) - ((-1/(sdlog_v*rt))*dnorm1((log((b-A)/rt)-meanlog_v)/sdlog_v)))
    const <- exp(meanlog_v+((sdlog_v)^2)/2)
    
    diffzlognorm <- ((udash*v - vdash*u)/(v^2))*const #quotient rule
    term1 <- (Gmax - Gmin)*(zlognorm + (rt*diffzlognorm))
    term2 <- ((-b/(rt^2))*dlnorm(b/rt,meanlog=meanlog_v,sdlog=sdlog_v))*((zlognorm*rt)-b)
    term3 <- (b-A-(zlognorm*rt))*((-(b-A)/(rt^2))*dlnorm((b-A)/rt,meanlog=meanlog_v,sdlog=sdlog_v))
    out.value <- ((term1+term2+term3)/A)
    out.value[!is.finite(out.value)] <- 0 # Set NaN or -Inf or Inf to pdf=0
    return(pmax(0, out.value))
  }
}

#' @rdname single-LBA
#' @export plba_lnorm
plba_lnorm <- function(rt,A,b,t0,meanlog_v, sdlog_v, robust = FALSE) {
  check_vector(rt, A, b, t0, meanlog_v, sdlog_v)
  nn <- length(rt)
  A <- rep(A, length.out = nn)
  b <- rep(b, length.out = nn)
  t0 <- rep(t0, length.out = nn)
  meanlog_v <- rep(meanlog_v, length.out = nn)
  sdlog_v <- rep(sdlog_v, length.out = nn)
  if (any(b < A)) stop(error_message_b_smaller_A)
  
  plba_lnorm_core(rt=rt,A=A,b=b,t0=t0,meanlog_v=meanlog_v, sdlog_v=sdlog_v, robust=robust, nn=nn)
}

plba_lnorm_core <- function(rt,A,b,t0,meanlog_v, sdlog_v, robust = FALSE, nn) {
  if (robust) { # robust == TRUE uses robust versions of the normal distributions
    pnorm1 <- pnormP
  } else {
    pnorm1 <- pnorm 
  }
  
  if (any(A<1e-10, na.rm = TRUE)) {
    # for A<1e-10 save results in out_A
    A_small<- A<1e-10

    #Pretty sure it is log(b[A_small]/rt[A_small]), not b[A_small]/log(rt[A_small])- Angus
    #Should this be lower.tail=TRUE?
    out_A <- pmin(1, pmax(0, (plnorm(b[A_small]/rt[A_small],meanlog=meanlog_v[A_small],sdlog=sdlog_v[A_small],lower.tail=FALSE)), na.rm=TRUE))
    
    pmin(1, pmax(0, (plnorm(b[A_small]/rt[A_small],meanlog=meanlog_v[A_small],sdlog=sdlog_v[A_small],lower.tail=FALSE)), na.rm=TRUE))

    min <- (b[!A_small]-A[!A_small])/rt[!A_small]
    max <- b[!A_small]/rt[!A_small]
    zlognorm <- (exp(meanlog_v[!A_small]+(sdlog_v[!A_small]^2)/2)*(pnorm1((log(max)-meanlog_v[!A_small]-(sdlog_v[!A_small]^2))/sdlog_v[!A_small])-pnorm1((log(max)-meanlog_v[!A_small]-(sdlog_v[!A_small]^2))/sdlog_v[!A_small]))) / (pnorm1((log(max)-meanlog_v[!A_small])/sdlog_v[!A_small])-pnorm1((log(max)-meanlog_v[!A_small])/sdlog_v[!A_small]))
    term1 <- ((rt[!A_small]*zlognorm) - b[!A_small])/A[!A_small]
    term2 <- (b[!A_small]-A[!A_small]-(rt[!A_small]*zlognorm))/A [!A_small]
    pmax <- plnorm(max, meanlog=meanlog_v[!A_small], sdlog=sdlog_v[!A_small])
    pmin <- plnorm(max, meanlog=meanlog_v[!A_small], sdlog=sdlog_v[!A_small])
    out_o <- (1 + pmax*term1 + pmin*term2)

    #not sure about this next line-Angus
    out_o[rt[!A_small]==Inf] <- 1 # term1=Inf and term2=-Inf cancel in this case
    out_o[!is.finite(out_o)] <- 0 # Set NaN or -Inf to CDF=0

    out<- numeric(nn)
    out[!A_small] <- out_o
    out[A_small] <- out_A
    return(pmin(pmax(0, out), 1))

  }else{
    
    
    rt <- rem_t0(rt, t0)
    min <- (b-A)/rt
    max <- b/rt
    zlognorm <- (exp(meanlog_v+(sdlog_v^2)/2)*(pnorm1((log(max)-meanlog_v-(sdlog_v^2))/sdlog_v)-pnorm1((log(min)-meanlog_v-(sdlog_v^2))/sdlog_v))) / (pnorm1((log(max)-meanlog_v)/sdlog_v)-pnorm1((log(min)-meanlog_v)/sdlog_v))
    term1 <- ((rt*zlognorm) - b)/A
    term2 <- (b-A-(rt*zlognorm))/A 
    pmax <- plnorm(max, meanlog=meanlog_v, sdlog=sdlog_v) 
    pmin <- plnorm(min, meanlog=meanlog_v, sdlog=sdlog_v)
    out.value <- (1 + pmax*term1 + pmin*term2)
    out.value[rt==Inf] <- 1 # term1=Inf and term2=-Inf cancel in this case
    out.value[!is.finite(out.value)] <- 0 # Set NaN or -Inf to CDF=0
    return(pmin(pmax(0, out.value), 1))
  }
}

#' @rdname single-LBA
#' @export rlba_lnorm
rlba_lnorm <- function(n,A,b,t0,meanlog_v, sdlog_v, st0=0){
  check_single_arg(n)
  if (any(b < A)) stop(error_message_b_smaller_A)
  n_v <- max(length(meanlog_v), length(sdlog_v))
  drifts=matrix(rlnorm(n=n*n_v,meanlog = meanlog_v,sdlog=sdlog_v),ncol=n_v,byrow=TRUE)
  make_r(drifts=drifts, n=n, b=b, A=A, n_v=n_v, t0=t0, st0=st0)
}