ReactionTime:

# ---------------------------------------------------------------------------------------
# This function generates and fit the Shifted Wald distribution
# using GAMLSS framewok as function of
#
# m > 0 = mu (rate of evidence accrual)
# a > 0 = sigma (response threshold)
# s > 0 = nu (the shift)
# We use the parameterization given by Heathcote, A. (2004). Fitting Wald and
# ex-Wald distributions to response time data: An example using functions for
# the S-PLUS package. Behavior Research Methods, Instruments, &
# Computers,36,678-694.
#
# Create by Raydonal Ospina, 2016
# Modified by:
# Raydonal  2016
# ---------------------------------------------------------------------------------------

#----------------------------------------------------------------------------------------
shiftWALD <-function (mu.link = "log", sigma.link = "log", nu.link ="identity")
{
    mstats <- checklink("mu.link", "Shift-Wald", substitute(mu.link), c("1/mu^2", "inverse", "log", "identity", "own"))
    dstats <- checklink("sigma.link", "Shift-Wald", substitute(sigma.link),  c("1/mu^2", "inverse", "log", "identity", "own"))
    vstats <- checklink("nu.link", "Shift-Wald", substitute(nu.link), c("inverse", "log", "identity", "own"))

    structure(
          list(family = c("shiftWALD", "Shift-Wald"),
           parameters = list(mu=TRUE, sigma=TRUE, nu=TRUE),

           nopar = 3,

           type = "Continuous",

           mu.link = as.character(substitute(mu.link)),

           sigma.link = as.character(substitute(sigma.link)),

           nu.link = as.character(substitute(nu.link)),

           mu.linkfun = mstats$linkfun,

           sigma.linkfun = dstats$linkfun,

           nu.linkfun = vstats$linkfun,

           mu.linkinv = mstats$linkinv,

           sigma.linkinv = dstats$linkinv,

           nu.linkinv = vstats$linkinv,

           mu.dr = mstats$mu.eta,

           sigma.dr = dstats$mu.eta,

           nu.dr = vstats$mu.eta,

              #first derivate of log-density respect to mu
                dldm = function(y, mu, sigma, nu) {
                                dldm = sigma-mu*(y-nu)
                                dldm
                                }, #OK

              #second derivate of log-density respect to mu
              d2ldm2 = function(y,nu) {
                                dldm2 = nu-y #sigma/mu
                                dldm2
                                }, # -y, #OK

              #first derivate of log-density respect to sigma
                dldd = function(y, mu, sigma, nu) {
                                dldd = mu + (1/sigma) - (sigma/(y-nu))
                                dldd
                                }, # OK

              #second derivate of log-density respect to sigma
              d2ldd2 = function(y, sigma, nu) {
                                d2ldd2 =   (-1/sigma^2) - (1/(y-nu))
                                d2ldd2
                                },   # (-1/sigma^2 - 1/y), # OK,


              #first derivate log-density respect to nu
              dldv = function(y, mu, sigma, nu) {
                              dldv = 3/(2*(y-nu)) - (sigma^2/ (2*(y-nu)^2)) + (mu^2/2)
                              dldv
                              },

              #second derivate log-density respect to nu
              d2ldv2 = function(y, sigma, nu) {
              d2ldv2 = (3/2 - (sigma^2/(y-nu)))/ (y-nu)^2
              d2ldv2
              },


              #partial derivate of log-density respect to mu and sigma
              d2ldmdd = function(y,mu,sigma) {
              d2ldmdd = rep(1,length(y))
              d2ldmdd
              },

              #partial derivate of log-density respect to mu and nu
              d2ldmdv = function(mu) {
              d2ldmdv = mu
              d2ldmdv
              },

              #partial derivate of log-density respect to sigma and nu
              d2ldddv = function(y, sigma, nu) {
              d2ldddv = -sigma/(y-nu)^2
              d2ldddv
              },

         G.dev.incr  = function(y, mu, sigma, nu, ...) -2*dshiftWALD(y, mu, sigma, nu, log=TRUE),

          rqres = expression(rqres(pfun="pshiftWALD", type="Continuous", y=y, mu=mu, sigma=sigma, nu=nu)),

#         # Initial values based on moments
        mu.initial = expression( mu <- rep(sqrt(mean(y)/var(y)), length(y))),
        sigma.initial = expression(sigma <- rep(sqrt(mean(y)/var(y))*mean(y), length(y))),
        nu.initial = expression(nu <- rep(0.9*min(y), length(y))),

               y.valid = function(y)  all(y > 0),
              mu.valid = function(mu) all(mu > 0),
           sigma.valid = function(sigma)  all(sigma > 0),
              nu.valid = function(nu) all(nu <= min(y) )

            ),
            class = c("gamlss.family","family"))
}
#----------------------------------------------------------------------------------------

#----------------------------------------------------------------------------------------
# Shifted Wald density
dshiftWALD <-function(x, mu = 1, sigma = 1, nu=0, log=FALSE) # OK
{        if (any(mu < 0))  stop(paste("mu must be positive", "\n", ""))
         if (any(sigma < 0))  stop(paste("sigma must be positive", "\n", ""))
         if (any(nu > min(x)))  stop(paste("nu must be less than the minimum observed x ", "\n", ""))
         if (any(x < 0))  stop(paste("x must be positive", "\n", ""))
         x = x-nu
         log.lik <- log(sigma) -(sigma-mu*x)^2/(2*x) -0.5*log(2*pi) -((3/2)*log(x))
         if(log==FALSE) fy  <- exp(log.lik) else fy <- log.lik
         fy
}



#----------------------------------------------------------------------------------------
# Shifted Wald cumulative density
pshiftWALD <- function(q, mu = 1, sigma = 1, nu=0, lower.tail = TRUE, log.p = FALSE) # OK
  {    #  browser()
    if (any(mu < 0))  stop(paste("mu must be positive", "\n", ""))
    if (any(sigma < 0))  stop(paste("sigma must be positive", "\n", ""))

    if (any(q < 0))  stop(paste("y must be positive", "\n", ""))
    if (any(nu > min(q)))  stop(paste("nu must be less than the minimum observed x ", "\n", ""))
      lq <- length(q)
   sigma <- rep(sigma, length = lq)
      mu <- rep(mu, length = lq)

    q = q-nu
    cdf1 <- pnorm(((mu*q) - sigma)/sqrt(q))

   lcdf2 <- ( 2*mu*sigma)+pnorm(- ( ((mu*q) + sigma)/sqrt(q)),log.p=TRUE)
        cdf <- cdf1+ exp(lcdf2)

    if(lower.tail==TRUE) cdf  <- cdf else  cdf <- 1-cdf
    if(log.p==FALSE) cdf  <- cdf else  cdf <- log(cdf)
    cdf
   }
#----------------------------------------------------------------------------------------

#
#
# #----------------------------------------------------------------------------------------
# # Shifted Wald Quantile function - Precisa ser corrigida - Ainda não tenho certeza que este totalmente correta
#
# qshiftWALD <- function(p, mu=1, sigma=1, nu=0,  lower.tail = TRUE, log.p = FALSE)
#  {
#     #---functions--------------------------------------------
#        h1 <- function(q)
#        {
#      pshiftWALD(q , mu = mu[i], sigma = sigma[i], nu=nu[i])-p[i]
#        }
#        h <- function(q)
#        {
#      pshiftWALD(q , mu = mu[i], sigma = sigma[i], nu=nu[i])
#        }
#      #-------------------------------------------------------
#     if (any(mu <= 0))  stop(paste("mu must be positive", "\n", ""))
#     if (any(sigma <= 0))  stop(paste("sigma must be positive", "\n", ""))
#     #if (any(nu < 0))  stop(paste("nu must be positive", "\n", ""))
#     if (log.p==TRUE) p <- exp(p) else p <- p
#     if (lower.tail==TRUE) p <- p else p <- 1-p
#     if (any(p < 0)|any(p > 1))  stop(paste("p must be between 0 and 1", "\n", ""))
#         lp <-  max(length(p),length(mu),length(sigma))
#           p <- rep(p, length = lp)
#       sigma <- rep(sigma, length = lp)
#          mu <- rep(mu, length = lp)
#           q <- rep(0,lp)
#          for (i in  seq(along=p))
#          {
#          if (h(mu[i])<p[i])
#           {
#            interval <- c(mu[i], mu[i]+sigma[i])
#            j <-2
#            while (h(interval[2]) < p[i])
#               {interval[2]<- mu[i]+j*sigma[i]
#               j<-j+1
#               }
#            }
#           else
#            {
#            interval <-  interval <- c(.Machine$double.xmin, mu[i])
#            }
#         q[i] <- uniroot(h1, interval)$root
#          }
#     q
#    }
# #----------------------------------------------------------------------------------------




qshiftWALD <- function(p, mu=1, sigma=1, nu=0,  lower.tail = TRUE, log.p = FALSE)
{


  if(length(n)>1) n <- length(n);
  if(length(mu)>1 && length(mu)!=n) mu <- rep(mu,length=n)
  if(length(sigma)>1 && length(sigma)!=n) lambda <- rep(sigma,length=n)
  y2 <- rchisq(n,1); y2onm <- y2/mu; u <- runif(n)
  r1 <- (2*sigma + y2onm - sqrt(y2onm*(4*sigma+y2onm)))/(2*mu)
  r2 <- (sigma/mu)^2/r1; ifelse(u < sigma/(sigma+mu*r1), nu+r1, nu+r2)




  #---functions--------------------------------------------
  h1 <- function(q)
  {
    pshiftWALD(q , mu = mu[i], sigma = sigma[i], nu=nu[i])-p[i]
  }
  h <- function(q)
  {
    pshiftWALD(q , mu = mu[i], sigma = sigma[i], nu=nu[i])
  }
  #-------------------------------------------------------
  if (any(mu <= 0))  stop(paste("mu must be positive", "\n", ""))
  if (any(sigma <= 0))  stop(paste("sigma must be positive", "\n", ""))
  #if (any(nu < 0))  stop(paste("nu must be positive", "\n", ""))
  if (log.p==TRUE) p <- exp(p) else p <- p
  if (lower.tail==TRUE) p <- p else p <- 1-p
  if (any(p < 0)|any(p > 1))  stop(paste("p must be between 0 and 1", "\n", ""))
  lp <-  max(length(p),length(mu),length(sigma))
  p <- rep(p, length = lp)
  sigma <- rep(sigma, length = lp)
  mu <- rep(mu, length = lp)
  q <- rep(0,lp)
  for (i in  seq(along=p))
  {
    if (h(mu[i])<p[i])
    {
      interval <- c(mu[i], mu[i]+sigma[i])
      j <-2
      while (h(interval[2]) < p[i])
      {interval[2]<- mu[i]+j*sigma[i]
       j<-j+1
      }
    }
    else
    {
      interval <-  interval <- c(.Machine$double.xmin, mu[i])
    }
    q[i] <- uniroot(h1, interval)$root
  }
  q
}




#----------------------------------------------------------------------------------------
# Shifted Wald random function adapted from pp. 79-80, Dagpunar, J. (1988).
# Principles of Random Variate Generation. Clarendon Press, Oxford.
rwald2 <- function(n,mu,sigma,nu=0) {
  if(length(n)>1) n <- length(n);
  if(length(mu)>1 && length(mu)!=n) mu <- rep(mu,length=n)
  if(length(sigma)>1 && length(sigma)!=n) lambda <- rep(sigma,length=n)
  y2 <- rchisq(n,1); y2onm <- y2/mu; u <- runif(n)
  r1 <- (2*sigma + y2onm - sqrt(y2onm*(4*sigma+y2onm)))/(2*mu)
  r2 <- (sigma/mu)^2/r1; ifelse(u < sigma/(sigma+mu*r1), nu+r1, nu+r2)
}


#----------------------------------------------------------------------------------------
rshiftWALD <- function(n, mu=1, sigma=1, nu=0, ...) # OK
  {
  if (any(mu <= 0))  stop(paste("mu must be positive", "\n", ""))
  if (any(sigma <= 0))  stop(paste("sigma must be positive", "\n", ""))
  if (any(nu > min(x)))  stop(paste("nu must be less than the minimum observed x ", "\n", ""))
  if (any(n <= 0))  stop(paste("n must be a positive integer", "\n", ""))
    n <- ceiling(n)
    p <- runif(n)
    r <- qshiftWALD(p,mu=mu,sigma=sigma, nu=nu, ...)
    r
  }
#----------------------------------------------------------------------------------------