#' The Ex-Wald family
#'
#' @author Freddy Hernandez, \email{fhernanb@unal.edu.co}
#'
#' @description
#' The function \code{ExWALD()} defines the Ex-wALD distribution, three-parameter
#' continuous distribution for a \code{gamlss.family} object to be used in 
#' GAMLSS fitting using the function \code{gamlss()}.
#'
#' @param mu.link defines the mu.link, with "log" link as the default for the mu parameter.
#' @param sigma.link defines the sigma.link, with "log" link as the default for the sigma parameter.
#' @param nu.link defines the nu.link, with "log" link as the default for the nu parameter.
#'
#' @references
#' Schwarz, W. (2001). The ex-Wald distribution as a descriptive model 
#' of response times. Behavior Research Methods, 
#' Instruments, & Computers, 33, 457-469.
#' 
#' 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.
#'
#' @seealso \link{dExWALD}.
#'
#' @details
#' The Ex-Wald distribution with parameters \eqn{\mu}, \eqn{\sigma} 
#' and \eqn{\nu} has density given by
#'
#' \eqn{f(x |\mu, \sigma, \nu) = \frac{1}{\nu} \exp(\frac{-x}{\nu} + \sigma(\mu-k)) F_W(x|k, \sigma) \, \text{for} \, k \geq 0}
#'
#' \eqn{f(x |\mu, \sigma, \nu) = \frac{1}{\nu} \exp\left( \frac{-(\sigma-\mu)^2}{2x} \right) Re \left( w(k^\prime \sqrt{x/2} + \frac{\sigma i}{\sqrt{2x}}) \right)  \, \text{for} \, k < 0}
#'
#' where \eqn{k=\sqrt{\mu^2-\frac{2}{\nu}}}, 
#' \eqn{k^\prime=\sqrt{\frac{2}{\nu}-\mu^2}} and
#' \eqn{F_W} corresponds to the cumulative function of 
#' the Wald distribution. 
#' 
#' More details about those expressions
#' can be found on page 680 from Heathcote (2004).
#'
#' @returns Returns a gamlss.family object which can be used to fit a 
#' Ex-WALD distribution in the \code{gamlss()} function.
#'
#' @example examples/examples_ExWALD.R
#'
#' @importFrom gamlss.dist checklink
#' @importFrom gamlss rqres.plot
#' @export
ExWALD <- function(mu.link="log", 
                   sigma.link="log", 
                   nu.link="log") {
  
  mstats <- checklink("mu.link", "Ex-Wald",
                      substitute(mu.link), c("log"))
  dstats <- checklink("sigma.link", "Ex-Wald",
                      substitute(sigma.link), c("log"))
  vstats <- checklink("nu.link", "Ex-Wald",
                      substitute(nu.link), c("log"))
  
  structure(list(family     = c("ExWALD", "Ex-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 derivates
                 
                 dldm = function(y, mu, sigma, nu) {
                   dm <- gamlss::numeric.deriv(dExWALD(y, mu, sigma, nu, log=TRUE),
                                               theta="mu",
                                               delta=0.001)
                   dldm <- as.vector(attr(dm, "gradient"))
                   dldm
                 },
                 
                 dldd = function(y, mu, sigma, nu) {
                   dd <- gamlss::numeric.deriv(dExWALD(y, mu, sigma, nu, log=TRUE),
                                               theta="sigma",
                                               delta=0.001)
                   dldd <- as.vector(attr(dd, "gradient"))
                   dldd
                 },
                 
                 dldv = function(y, mu, sigma, nu) {
                   dv <- gamlss::numeric.deriv(dExWALD(y, mu, sigma, nu, log=TRUE),
                                               theta="nu",
                                               delta=0.001)
                   dldv <- as.vector(attr(dv, "gradient"))
                   dldv
                 },

                 # Second derivates
                 
                 d2ldm2 = function(y, mu, sigma, nu) {
                   dm <- gamlss::numeric.deriv(dExWALD(y, mu, sigma, nu, log=TRUE),
                                                 theta="mu",
                                                 delta=0.001)
                   dldm   <- as.vector(attr(dm, "gradient"))
                   d2ldm2 <- - dldm*dldm
                   d2ldm2 <- ifelse(d2ldm2 < -1e-15, d2ldm2, -1e-15)
                   d2ldm2
                 },
                 
                 d2ldmdd = function(y, mu, sigma, nu) {
                   dm <- gamlss::numeric.deriv(dExWALD(y, mu, sigma, nu, log=TRUE),
                                                 theta="mu",
                                                 delta=0.001)
                   dldm <- as.vector(attr(dm, "gradient"))
                   dd   <- gamlss::numeric.deriv(dExWALD(y, mu, sigma, nu, log=TRUE),
                                                 theta="sigma",
                                                 delta=0.001)
                   dldd    <- as.vector(attr(dd, "gradient"))
                   d2ldmdd <- - dldm*dldd
                   d2ldmdd <- ifelse(d2ldmdd < -1e-15, d2ldmdd, -1e-15)
                   d2ldmdd
                 },
                 
                 d2ldmdv = function(y, mu, sigma, nu) {
                   dm <- gamlss::numeric.deriv(dExWALD(y, mu, sigma, nu, log=TRUE),
                                                 theta="mu",
                                                 delta=0.001)
                   dldm <- as.vector(attr(dm, "gradient"))
                   dv   <- gamlss::numeric.deriv(dExWALD(y, mu, sigma, nu, log=TRUE),
                                                 theta="nu",
                                                 delta=0.001)
                   dldv    <- as.vector(attr(dv, "gradient"))
                   d2ldmdv <- - dldm*dldv
                   d2ldmdv <- ifelse(d2ldmdv < -1e-15, d2ldmdv, -1e-15)
                   d2ldmdv
                 },
                 
                 d2ldd2  = function(y, mu, sigma, nu) {
                   dd  <- gamlss::numeric.deriv(dExWALD(y, mu, sigma, nu, log=TRUE),
                                                 theta="sigma",
                                                 delta=0.001)
                   dldd   <- as.vector(attr(dd, "gradient"))
                   d2ldd2 <- - dldd*dldd
                   d2ldd2 <- ifelse(d2ldd2 < -1e-15, d2ldd2, -1e-15)
                   d2ldd2
                 },
                 
                 d2ldddv = function(y, mu, sigma, nu) {
                   dd <- gamlss::numeric.deriv(dExWALD(y, mu, sigma, nu, log=TRUE),
                                                 theta="sigma",
                                                 delta=0.001)
                   dldd <- as.vector(attr(dd, "gradient"))
                   dv   <- gamlss::numeric.deriv(dExWALD(y, mu, sigma, nu, log=TRUE),
                                                 theta="nu",
                                                 delta=0.001)
                   dldv    <- as.vector(attr(dv, "gradient"))
                   d2ldddv <- - dldd*dldv
                   d2ldddv <- ifelse(d2ldddv < -1e-15, d2ldddv, -1e-15)
                   d2ldddv
                 },
                 
                 d2ldv2  = function(y, mu, sigma, nu) {
                   dv   <- gamlss::numeric.deriv(dExWALD(y, mu, sigma, nu, log=TRUE),
                                                 theta="nu",
                                                 delta=0.001)
                   dldv   <- as.vector(attr(dv, "gradient"))
                   d2ldv2 <- - dldv*dldv
                   d2ldv2 <- ifelse(d2ldv2 < -1e-15, d2ldv2, -1e-15)
                   d2ldv2
                 },
                 
                 G.dev.incr = function(y, mu, sigma, nu, ...) -2*dExWALD(y, mu, sigma, nu, log=TRUE),
                 rqres      = expression(rqres(pfun="pExWALD", type="Continuous", y=y, mu=mu, sigma=sigma, nu=nu)),
                 
                 mu.initial    = expression(mu    <- rep(fitexw(y)[1], length(y))),
                 sigma.initial = expression(sigma <- rep(fitexw(y)[2], length(y))),
                 nu.initial    = expression(nu    <- rep(fitexw(y)[3], length(y))),
                 
                 mu.valid    = function(mu)    all(mu > 0),
                 sigma.valid = function(sigma) all(sigma > 0),
                 nu.valid    = function(nu)    all(nu > 0),
                 
                 y.valid = function(y) all(y > 0),
                 
                 mean     = function(mu, sigma, nu) {nu + sigma/mu}, 
                 variance = function(mu, sigma, nu) {nu^2 + sigma/mu^3}
  ),
  class = c("gamlss.family", "family"))
}
#' Auxiliar function for the Ex-Wald distribution
#' @description This function generates start values.
#' @param x a value for x.
#' @param p a value for p.
#' @return returns a vector with starting values.
#' @keywords internal
#' @importFrom stats sd var
#' @export
exwstpt <- function(x, p=0.5) {
  t <- p*sd(x)
  m <- sqrt((mean(x)-t)/(var(x)-t^2))
  a <- m*(mean(x)-t)
  res <- c(m, a, t)
  return(res)
}
#' Auxiliar pwald function for the Ex-Wald distribution
#' @description This function emulates the pWALD function.
#' @param w a value for w.
#' @param m a value for m.
#' @param a a value for a.
#' @param s a value for s.
#' @return returns a vector with cumulative probabilities.
#' @keywords internal
#' @export
pwald <- function(w, m, a, s=0) {
  w <- w - s; 
  sqrtw <- sqrt(w); 
  k1 <- (m*w-a)/sqrtw; 
  k2 <- (m*w+a)/sqrtw
  p1 <- exp(2*a*m); 
  p2 <- pnorm(-k2); 
  bad <- (p1==Inf) | (p2==0); 
  p <- p1*p2
  p[bad] <- (exp(-(k1[bad]^2)/2 - 0.94/(k2[bad]^2))/(k2[bad]*((2*pi)^.5)))
  p + pnorm(k1)
}
#' Auxiliar function to generate random values for Ex-Wald distribution
#' @description This function generates another form to generate values.
#' @param x a value for x.
#' @param y a value for y.
#' @return returns a vector with values.
#' @keywords internal
#' @export
rew <- function(x, y) {
  uv <- uandv(y, x)
  exp(y^2-x^2)*(cos(2*x*y)*(1-uv[,1])+sin(2*x*y)*uv[,2])
}
#' Auxiliar den_exw function for the Ex-Wald distribution
#' @description This function emulates the dExWALD function.
#' @param r value for r.
#' @param m value for m.
#' @param a value for a.
#' @param t value for t.
#' @return returns a vector with probabilities.
#' @keywords internal
#' @export
den_exw <- function(r, m, a, t) {
  k <- (m^2 - (2/t))
  if (k < 0) {
    res <- exp(m*a - (a^2)/(2*r) - r*(m^2)/2)*
      rew(sqrt(-r*k/2), a/sqrt(2*r))/t
  } else {
    k <- sqrt(k)
    res <- pwald(r,k,a)*exp(a*(m-k) - (r/t))/t
  }
  return(res)
}
#' Auxiliar function to obtain minus loglik for Ex-Wald 
#' @description This function calculates minus loglik.
#' @param p a vector with the parameters.
#' @param x a vector with the random sample.
#' @return returns the minus loglik.
#' @keywords internal
#' @export
negllexw <- function(p, x) {
  -sum(log(den_exw(x, p[1], p[2], p[3])))
}
#' Auxiliar function for the Ex-Wald distribution
#' @description This function generates starting values.
#' @param rt a vector with the random sample.
#' @param p a value for p.
#' @param start an optional start vector.
#' @param scaleit logical value to scale.
#' @return returns a vector with cumulative probabilities.
#' @keywords internal
#' @importFrom stats nlminb
#' @export
fitexw <- function(rt, p=0.5,
                   start=exwstpt(rt, p),
                   scaleit=TRUE){
  if (scaleit) 
    fit <- nlminb(start=start,
                  lower=c(1e-8, 1e-8, 1),
                  objective=negllexw,
                  x=rt,
                  scale=(1/start),
                  control=list(eval.max=400, iter.max=300))
  else 
    fit <- nlminb(start=start,
                  lower=c(1e-8, 1e-8, 1),
                  objective=negllexw,
                  x=rt,
                  control=list(eval.max=400, iter.max=300, 
                               scale.tol=1e-8))

  return(fit$par)
}