R/family.actuary.R

Defines functions inv.paralogistic paralogistic inv.lomax fisk lomax betaII dagum sinmad dinv.paralogistic dinv.lomax ddagum dparalogistic dfisk dlomax dsinmad dbetaII pinv.paralogistic pinv.lomax pdagum pparalogistic pfisk plomax psinmad qinv.paralogistic qinv.lomax qdagum qparalogistic qfisk qlomax qsinmad rinv.paralogistic rinv.lomax rdagum rparalogistic rfisk rlomax rsinmad dgenbetaII genbetaII genbetaII.Loglikfun4 exponential.mo exponential.mo.control gompertz gompertz.control rgompertz qgompertz pgompertz dgompertz makeham makeham.control rmakeham qmakeham pmakeham dmakeham perks perks.control rperks qperks pperks dperks qbeard dbeard dbeard pmperks dmperks pmbeard dmbeard gumbelII rgumbelII qgumbelII pgumbelII dgumbelII

Documented in betaII dagum dbetaII ddagum dfisk dgenbetaII dgompertz dgumbelII dinv.lomax dinv.paralogistic dlomax dmakeham dparalogistic dperks dsinmad fisk genbetaII genbetaII.Loglikfun4 gompertz gumbelII inv.lomax inv.paralogistic lomax makeham paralogistic pdagum perks pfisk pgompertz pgumbelII pinv.lomax pinv.paralogistic plomax pmakeham pparalogistic pperks psinmad qdagum qfisk qgompertz qgumbelII qinv.lomax qinv.paralogistic qlomax qmakeham qparalogistic qperks qsinmad rdagum rfisk rgompertz rgumbelII rinv.lomax rinv.paralogistic rlomax rmakeham rparalogistic rperks rsinmad sinmad

# These functions are
# Copyright (C) 1998-2024 T.W. Yee, University of Auckland.
# All rights reserved.















dgumbelII <- function(x, scale = 1, shape, log = FALSE) {


  if (!is.logical(log.arg <- log) || length(log) != 1)
    stop("bad input for argument 'log'")
  rm(log)

  LLL <- max(length(x), length(shape), length(scale))
  if (length(x)       != LLL) x       <- rep_len(x,       LLL)
  if (length(shape)   != LLL) shape   <- rep_len(shape,   LLL)
  if (length(scale)   != LLL) scale   <- rep_len(scale,   LLL)


  ans <- x
  index0 <- (x < 0) & is.finite(x) & !is.na(x)

  ans[!index0] <- log(shape[!index0] / scale[!index0]) +
    (shape[!index0] + 1) * log(scale[!index0] / x[!index0]) -
    (x[!index0] / scale[!index0])^(-shape[!index0])
  ans[index0] <- log(0)
  ans[x == Inf] <- log(0)

  if (log.arg) {
  } else {
    ans <- exp(ans)
    ans[index0] <- 0
    ans[x == Inf] <- 0
  }
  ans[shape <= 0 | scale <= 0] <- NaN
  ans
}



pgumbelII <- function(q, scale = 1, shape,
                      lower.tail = TRUE, log.p = FALSE) {

  # 20150121 KaiH
  if (!is.logical(lower.tail) || length(lower.tail ) != 1)
    stop("bad input for argument 'lower.tail'")

  # 20150121 KaiH
  if (!is.logical(log.p) || length(log.p) != 1)
    stop("bad input for argument 'log.p'")

  LLL <- max(length(q), length(shape), length(scale))
  if (length(q)       != LLL) q       <- rep_len(q,       LLL)
  if (length(shape)   != LLL) shape   <- rep_len(shape,   LLL)
  if (length(scale)   != LLL) scale   <- rep_len(scale,   LLL)

  # 20150121 KaiH
  if (lower.tail) {
    if (log.p) {
      ans <- -(q / scale)^(-shape)
      ans[q <= 0 ] <- -Inf
      ans[q == Inf] <- 0
    } else {
      ans <- exp(-(q / scale)^(-shape))
      ans[q <= 0] <- 0
      ans[q == Inf] <- 1
    }
  } else {
    if (log.p) {
      ans <- log(-expm1(-(q / scale)^(-shape)))
      ans[q <= 0] <- 0
      ans[q == Inf] <- -Inf
    } else {
      ans <- -expm1(-(q / scale)^(-shape))
      ans[q <= 0] <- 1
      ans[q == Inf] <- 0
    }
  }
  ans[shape <= 0 | scale <= 0] <- NaN
  ans
}










qgumbelII <- function(p, scale = 1, shape,
                      lower.tail = TRUE, log.p = FALSE) {



  if (!is.logical(lower.tail) || length(lower.tail ) != 1)
    stop("bad input for argument 'lower.tail'")
  if (!is.logical(log.p) || length(log.p) != 1)
    stop("bad input for argument 'log.p'")



  LLL <- max(length(p), length(shape), length(scale))
  if (length(p)       != LLL) p       <- rep_len(p,       LLL)
  if (length(shape)   != LLL) shape   <- rep_len(shape,   LLL)
  if (length(scale)   != LLL) scale   <- rep_len(scale,   LLL)


  if (lower.tail) {
    if (log.p) {
      ln.p <- p
      ans <- scale * (-ln.p)^(-1 / shape)
      ans[ln.p > 0] <- NaN
    } else { # Default
      ans <- scale * (-log(p))^(-1 / shape)
      ans[p < 0] <- NaN
      ans[p == 0] <- 0
      ans[p == 1] <- Inf
      ans[p > 1] <- NaN
    }
  } else {
    if (log.p) {
      ln.p <- p
      ans <- scale * (-log(-expm1(ln.p)))^(-1 / shape)
      ans[ln.p > 0] <- NaN
    } else {
      ans <- scale * (-log1p(-p))^(-1 / shape)
      ans[p < 0] <- NaN
      ans[p == 0] <- Inf
      ans[p == 1] <- 0
      ans[p > 1] <- NaN
    }
  }

  ans[shape <= 0 | scale <= 0] <- NaN
  ans
}


rgumbelII <- function(n, scale = 1, shape) {
  qgumbelII(runif(n), shape = shape, scale = scale)
}









 gumbelII <-
  function(lscale = "loglink", lshape = "loglink",
           iscale = NULL,   ishape = NULL,
           probs.y = c(0.2, 0.5, 0.8),
           perc.out = NULL,  # 50,
           imethod = 1, zero = "shape", nowarning = FALSE) {




  lshape <- as.list(substitute(lshape))
  eshape <- link2list(lshape)
  lshape <- attr(eshape, "function.name")

  lscale <- as.list(substitute(lscale))
  escale <- link2list(lscale)
  lscale <- attr(escale, "function.name")


  if (!is.Numeric(imethod, length.arg = 1,
                  integer.valued = TRUE, positive = TRUE) ||
      imethod > 2)
    stop("argument 'imethod' must be 1 or 2")
  if (!is.Numeric(probs.y, positive  = TRUE) ||
      length(probs.y) < 2 ||
      max(probs.y) >= 1)
    stop("bad input for argument 'probs.y'")
  if (length(perc.out))
    if (!is.Numeric(perc.out, positive  = TRUE) ||
        max(probs.y) >= 100)
    stop("bad input for argument 'perc.out'")


  if (length(ishape))
    if (!is.Numeric(ishape, positive = TRUE))
      stop("argument 'ishape' values must be positive")
  if (length(iscale))
    if (!is.Numeric(iscale, positive = TRUE))
      stop("argument 'iscale' values must be positive")


  new("vglmff",
  blurb = c("Gumbel Type II distribution\n\n",
            "Links:    ",
            namesof("scale", lscale, escale), ", ",
            namesof("shape", lshape, eshape), "\n",
            "Mean:     scale^(1/shape) * gamma(1 - 1 / shape)\n",
            "Variance: scale^(2/shape) * (gamma(1 - 2/shape) - ",
                      "gamma(1 + 1/shape)^2)"),
  constraints = eval(substitute(expression({
    constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
                     M = M, M1 = 2,
                     predictors.names = predictors.names)

  }), list( .zero = zero ))),

  infos = eval(substitute(function(...) {
    list(M1 = 2,
         Q1 = 1,
         dpqrfun = "gumbelII",
         parameters.names = c("scale", "shape"),
         perc.out = .perc.out ,
         zero = .zero )
  }, list( .zero = zero,
           .perc.out = perc.out
         ))),

  initialize = eval(substitute(expression({

    temp5 <-
    w.y.check(w = w, y = y,
              Is.positive.y = TRUE,
              ncol.w.max = Inf,
              ncol.y.max = Inf,
              out.wy = TRUE,
              colsyperw = 1,
              maximize = TRUE)
    w <- temp5$w
    y <- temp5$y

    ncoly <- ncol(y)
    M1 <- 2
    extra$ncoly <- ncoly
    extra$M1 <- M1
    M <- M1 * ncoly


    mynames1 <- param.names("scale", ncoly, skip1 = TRUE)
    mynames2 <- param.names("shape", ncoly, skip1 = TRUE)


    predictors.names <-
        c(namesof(mynames1, .lscale , .escale , tag = FALSE),
          namesof(mynames2, .lshape , .eshape , tag = FALSE))[
          interleave.VGAM(M, M1 = M1)]


    Shape.init <- matrix(if (length( .ishape ))
                             .ishape else 0 + NA,
                         n, ncoly, byrow = TRUE)
    Scale.init <- matrix(if (length( .iscale ))
                             .iscale else 0 + NA,
                         n, ncoly, byrow = TRUE)

    if (!length(etastart)) {
      if (!length( .ishape ) ||
          !length( .iscale )) {
        for (ilocal in 1:ncoly) {

          anyc <- FALSE  # extra$leftcensored|extra$rightcensored
          i11 <- if ( .imethod == 1)
                   anyc else FALSE # can be all data
          probs.y <- .probs.y
          xvec <- log(-log(probs.y))
          fit0 <- lsfit(y  = xvec,
                        x  = log(quantile(y[!i11, ilocal],
                                          probs = probs.y )))


          if (!is.Numeric(Shape.init[, ilocal]))
            Shape.init[, ilocal] <- -fit0$coef["X"]
          if (!is.Numeric(Scale.init[, ilocal]))
            Scale.init[, ilocal] <-
              exp(fit0$coef["Intercept"] / Shape.init[, ilocal])
        }  # ilocal

        etastart <-
          cbind(theta2eta(Scale.init, .lscale , .escale ),
                theta2eta(Shape.init, .lshape , .eshape ))[,
                interleave.VGAM(M, M1 = M1)]
      }
    }
  }), list(
            .lscale = lscale, .lshape = lshape,
            .escale = escale, .eshape = eshape,
            .iscale = iscale, .ishape = ishape,
            .probs.y = probs.y,
            .imethod = imethod ) )),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    Scale <- eta2theta(eta[, c(TRUE, FALSE)], .lscale , .escale )
    Shape <- eta2theta(eta[, c(FALSE, TRUE)], .lshape , .eshape )
    Shape <- as.matrix(Shape)

    if (length( .perc.out ) > 1 && ncol(Shape) > 1)
      stop("argument 'perc.out' should be of length one since ",
           "there are multiple responses")

    if (!length( .perc.out )) {
      return(Scale * gamma(1 - 1 / Shape))
    }

    ans <- if (length( .perc.out ) > 1) {
      qgumbelII(p = matrix( .perc.out / 100, length(Shape),
                           length( .perc.out ), byrow = TRUE),
                shape = Shape, scale = Scale)
    } else {
      qgumbelII(p = .perc.out / 100, Shape, scale = Scale)
    }
    colnames(ans) <- paste0(as.character( .perc.out ), "%")
    ans
  }, list(
           .lscale = lscale, .lshape = lshape,
           .escale = escale, .eshape = eshape,
           .perc.out = perc.out ) )),
  last = eval(substitute(expression({


    M1 <- extra$M1
    misc$link <-
      c(rep_len( .lscale , ncoly),
        rep_len( .lshape , ncoly))[interleave.VGAM(M, M1 = M1)]
    temp.names <- c(mynames1, mynames2)[
        interleave.VGAM(M, M1 = M1)]
    names(misc$link) <- temp.names

    misc$earg <- vector("list", M)
    names(misc$earg) <- temp.names
    for (ii in 1:ncoly) {
      misc$earg[[M1*ii-1]] <- .escale
      misc$earg[[M1*ii  ]] <- .eshape
    }

    misc$M1 <- M1
    misc$imethod <- .imethod
    misc$expected <- TRUE
    misc$multipleResponses <- TRUE
    misc$perc.out <- .perc.out
    misc$true.mu <- FALSE # @fitted is not a true mu


  }), list(
            .lscale = lscale, .lshape = lshape,
            .escale = escale, .eshape = eshape,
            .perc.out = perc.out,
            .imethod = imethod ) )),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    Scale <- eta2theta(eta[, c(TRUE, FALSE)], .lscale , .escale )
    Shape <- eta2theta(eta[, c(FALSE, TRUE)], .lshape , .eshape )
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * dgumbelII(x = y, shape = Shape,
                                  scale = Scale, log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .lscale = lscale, .lshape = lshape,
           .escale = escale, .eshape = eshape
         ) )),
  vfamily = c("gumbelII"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    Scale <- eta2theta(eta[, c(TRUE, FALSE)], .lscale , .escale )
    Shape <- eta2theta(eta[, c(FALSE, TRUE)], .lshape , .eshape )
    okay1 <- all(is.finite(Scale)) && all(0 < Scale) &&
             all(is.finite(Shape)) && all(0 < Shape)
    okay1
  }, list( .lscale = lscale, .lshape = lshape,
           .escale = escale, .eshape = eshape) )),


  simslot = eval(substitute(
  function(object, nsim) {

    pwts <- if (length(pwts <- object@prior.weights) > 0)
              pwts else weights(object, type = "prior")
    if (any(pwts != 1))
      warning("ignoring prior weights")
    eta <- predict(object)
    Scale <- eta2theta(eta[, c(TRUE, FALSE)], .lscale , .escale )
    Shape <- eta2theta(eta[, c(FALSE, TRUE)], .lshape , .eshape )
    rgumbelII(nsim * length(Scale), shape = Shape, scale = Scale)
  }, list( .lscale = lscale, .lshape = lshape,
           .escale = escale, .eshape = eshape
         ) )),



  deriv = eval(substitute(expression({
    M1 <- 2
    Scale <- eta2theta(eta[, c(TRUE, FALSE)], .lscale , .escale )
    Shape <- eta2theta(eta[, c(FALSE, TRUE)], .lshape , .eshape )

    dl.dshape <- 1 / Shape + log(Scale / y) -
                 log(Scale / y) * (Scale / y)^Shape
    dl.dscale <- Shape / Scale -
        (Shape / y) * (Scale / y)^(Shape - 1)


    dscale.deta <- dtheta.deta(Scale, .lscale , .escale )
    dshape.deta <- dtheta.deta(Shape, .lshape , .eshape )

    myderiv <- c(w) * cbind(dl.dscale, dl.dshape) *
                      cbind(dscale.deta, dshape.deta)
    myderiv[, interleave.VGAM(M, M1 = M1)]
  }), list( .lscale = lscale, .lshape = lshape,
            .escale = escale, .eshape = eshape
          ) )),
  weight = eval(substitute(expression({
    EulerM <- -digamma(1.0)


    ned2l.dshape2 <- (1 + trigamma(2) + digamma(2)^2) / Shape^2
    ned2l.dscale2 <-  (Shape / Scale)^2
    ned2l.dshapescale <- digamma(2) / Scale


    wz <- array(c(c(w) * ned2l.dscale2 * dscale.deta^2,
                  c(w) * ned2l.dshape2 * dshape.deta^2,
                  c(w) * ned2l.dshapescale * dscale.deta *
                                             dshape.deta),
                dim = c(n, M / M1, 3))
    wz <- arwz2wz(wz, M = M, M1 = M1)
    wz
  }), list( .lscale = lscale, .lshape = lshape ))))
}







dmbeard <-
  function(x, shape, scale = 1, rho, epsilon, log = FALSE) {


  if (!is.logical(log.arg <- log) || length(log) != 1)
    stop("bad input for argument 'log'")
  rm(log)

  LLL <- max(length(x), length(shape), length(scale),
             length(rho), length(epsilon))
  if (length(x)       != LLL) x       <- rep_len(x,       LLL)
  if (length(shape)   != LLL) shape   <- rep_len(shape,   LLL)
  if (length(scale)   != LLL) scale   <- rep_len(scale,   LLL)
  if (length(rho)     != LLL) rho     <- rep_len(rho,     LLL)
  if (length(epsilon) != LLL) epsilon <- rep_len(epsilon, LLL)


  index0 <- (x < 0)

  ans <- log(epsilon * exp(-x * scale) + shape) +
            (-epsilon * x -
            ((rho * epsilon - 1) / (rho * scale)) *
            (log1p(rho * shape) -
             log(exp(-x * scale) + rho * shape) - scale * x)) -
            log(exp(-x * scale) + shape * rho)

  ans[index0] <- log(0)
  ans[x == Inf] <- log(0)

  if (log.arg) {
  } else {
    ans <- exp(ans)
    ans[index0] <- 0
    ans[x == Inf] <- 0
  }
  ans[shape <= 0 | scale <= 0 | rho <= 0 | epsilon <= 0] <- NaN
  ans
}


pmbeard <- function(q, shape, scale = 1, rho, epsilon) {

  LLL <- max(length(q), length(shape), length(scale),
             length(rho), length(epsilon))
  if (length(q)       != LLL) q       <- rep_len(q,       LLL)
  if (length(shape)   != LLL) shape   <- rep_len(shape,   LLL)
  if (length(scale)   != LLL) scale   <- rep_len(scale,   LLL)
  if (length(rho)     != LLL) rho     <- rep_len(rho,     LLL)
  if (length(epsilon) != LLL) epsilon <- rep_len(epsilon, LLL)


  ans <- -expm1(-epsilon * q -
               ((rho * epsilon - 1) / (rho * scale)) *
               (log1p(rho * shape) -
                log(exp(-scale * q) + rho * shape) - scale * q))
  ans[(q <= 0)] <- 0
  ans[shape <= 0 | scale <= 0 | rho <= 0 | epsilon <= 0] <- NaN
  ans[q == Inf] <- 1
  ans
}









dmperks <-
  function(x, scale = 1, shape, epsilon, log = FALSE) {

  if (!is.logical(log.arg <- log) || length(log) != 1)
    stop("bad input for argument 'log'")
  rm(log)

  LLL <- max(length(x), length(shape), length(scale),
             length(epsilon))
  if (length(x)       != LLL) x       <- rep_len(x,       LLL)
  if (length(shape)   != LLL) shape   <- rep_len(shape,   LLL)
  if (length(scale)   != LLL) scale   <- rep_len(scale,   LLL)
  if (length(epsilon) != LLL) epsilon <- rep_len(epsilon, LLL)


  index0 <- (x < 0)
  ans <- log(epsilon * exp(-x * scale) + shape) +
            (-epsilon * x -
            ((epsilon - 1) / scale) *
            (log1p(shape) -
             log(shape + exp(-x * scale)) -x * scale)) -
            log(exp(-x * scale) + shape)

  ans[index0] <- log(0)
  ans[x == Inf] <- log(0)
  if (log.arg) {
  } else {
    ans <- exp(ans)
    ans[index0] <- 0
    ans[x == Inf] <- 0
  }
  ans[shape <= 0 | scale <= 0 | epsilon <= 0] <- NaN
  ans
}  # dmperks



pmperks <- function(q, scale = 1, shape, epsilon) {

  LLL <- max(length(q), length(shape), length(scale))
  if (length(q)       != LLL) q       <- rep_len(q,       LLL)
  if (length(shape)   != LLL) shape   <- rep_len(shape,   LLL)
  if (length(scale)   != LLL) scale   <- rep_len(scale,   LLL)


  ans <- -expm1(-epsilon * q -
               ((epsilon - 1) / scale) *
               (log1p(shape) -
                log(shape + exp(-q * scale)) - q * scale))

  ans[(q <= 0)] <- 0
  ans[shape <= 0 | scale <= 0] <- NaN
  ans[q == Inf] <- 1
  ans
}














dbeard <- function(x, shape, scale = 1, rho, log = FALSE) {

 warning("does not integrate to unity")

  if (!is.logical(log.arg <- log) || length(log) != 1)
    stop("bad input for argument 'log'")
  rm(log)

  LLL <- max(length(x), length(shape), length(scale),
             length(rho))
  if (length(x)       != LLL) x       <- rep_len(x,       LLL)
  if (length(shape)   != LLL) shape   <- rep_len(shape,   LLL)
  if (length(scale)   != LLL) scale   <- rep_len(scale,   LLL)
  if (length(rho)     != LLL) rho     <- rep_len(rho,     LLL)

  index0 <- (x < 0)
    ans <- log(shape) - x * scale * (rho^(-1 / scale)) +
           log(rho) + log(scale) +
           (rho^(-1 / scale)) * log1p(shape * rho) -
           (1 + rho^(-1 / scale)) *
           log(shape * rho + exp(-x * scale))
    ans[index0] <- log(0)
    ans[x == Inf] <- log(0)


  if (log.arg) {
  } else {
    ans <- exp(ans)
    ans[index0] <- 0
    ans[x == Inf] <- 0
  }
  ans[shape <= 0 | scale <= 0 | rho <= 0] <- NaN
  ans
}






dbeard <- function(x, shape, scale = 1, rho, log = FALSE) {
  alpha <- shape
  beta  <- scale

 warning("does not integrate to unity")

  ret <- ifelse(x <= 0 | beta <= 0, NaN,
                exp(alpha+beta*x)*
                (1+exp(alpha+rho))**(exp(-rho/beta)) /
                (1+exp(alpha+rho+beta*x))**(1+exp(-rho/beta)))
  ret
}



qbeard <- function(x, u = 0.5, alpha = 1, beta = 1,rho = 1) {
  ret <-
    ifelse(x <= 0 | u <= 0 | u >= 1 |
           length(x) != length(u) | beta <= 0,
           NaN,
           (1/beta) * (log((u**(-beta*exp(rho))) *
           (1+exp(alpha+rho+beta*x))-1)-alpha-rho)-x)

  return(ret)
}











dperks <- function(x, scale = 1, shape, log = FALSE) {

  if (!is.logical(log.arg <- log) || length(log) != 1)
    stop("bad input for argument 'log'")
  rm(log)

  LLL <- max(length(x), length(shape), length(scale))
  if (length(x)     != LLL) x     <- rep_len(x,     LLL)
  if (length(shape) != LLL) shape <- rep_len(shape, LLL)
  if (length(scale) != LLL) scale <- rep_len(scale, LLL)

  index0 <- (x < 0)
    ans <- log(shape) - x +
           log1p(shape) / scale -
           (1 + 1 / scale) * log(shape + exp(-x * scale))
    ans[index0] <- log(0)
    ans[x == Inf] <- log(0)

  if (log.arg) {
  } else {
    ans <- exp(ans)
    ans[index0] <- 0
    ans[x == Inf] <- 0
  }
  ans[shape <= 0 | scale <= 0] <- NaN
  ans
}  # dperks



pperks <- function(q, scale = 1, shape,
                   lower.tail = TRUE, log.p = FALSE) {


  if (!is.logical(lower.tail) || length(lower.tail ) != 1)
    stop("bad input for argument 'lower.tail'")

  if (!is.logical(log.p) || length(log.p) != 1)
    stop("bad input for argument 'log.p'")


  LLL <- max(length(q), length(shape), length(scale))
  if (length(q)       != LLL) q       <- rep_len(q,       LLL)
  if (length(shape)   != LLL) shape   <- rep_len(shape,   LLL)
  if (length(scale)   != LLL) scale   <- rep_len(scale,   LLL)

  logS <- -q + (log1p(shape) -
          log(shape + exp(-q * scale))) / scale


  if (lower.tail) {
    if (log.p) {
      ans <- log(-expm1(logS))
      ans[q <= 0 ] <- -Inf
      ans[q == Inf] <- 0
    } else {
      ans <- -expm1(logS)
      ans[q <= 0] <- 0
      ans[q == Inf] <- 1
    }
  } else {
    if (log.p) {
      ans <- logS
      ans[q <= 0] <- 0
      ans[q == Inf] <- -Inf
    } else {
      ans <- exp(logS)
      ans[q <= 0] <- 1
      ans[q == Inf] <- 0
    }
  }

  ans[shape <= 0 | scale <= 0] <- NaN
  ans
}  # pperks


 qperks <-
  function(p, scale = 1, shape, lower.tail = TRUE, log.p = FALSE) {

  if (!is.logical(lower.tail) || length(lower.tail ) != 1)
    stop("bad input for argument 'lower.tail'")

  if (!is.logical(log.p) || length(log.p) != 1)
    stop("bad input for argument 'log.p'")

  LLL <- max(length(p), length(shape), length(scale))
  if (length(p)       != LLL) p       <- rep_len(p,       LLL)
  if (length(shape)   != LLL) shape   <- rep_len(shape,   LLL)
  if (length(scale)   != LLL) scale   <- rep_len(scale,   LLL)


  if (lower.tail) {
    if (log.p) {
      ln.p <- p
      tmp <- scale * log(-expm1(ln.p))
      onemFb <- exp(tmp)
      ans <- (log1p(shape - onemFb) - log(shape) - tmp) / scale
      ans[ln.p > 0] <- NaN
    } else {
      tmp <- scale * log1p(-p)
      onemFb <- exp(tmp)
      ans <- (log1p(shape - onemFb) - log(shape) - tmp) / scale
      ans[p < 0] <- NaN
      ans[p == 0] <- 0
      ans[p == 1] <- Inf
      ans[p > 1] <- NaN
    }
  } else {
    if (log.p) {
      ln.p <- p
      tmp <- scale * ln.p
      onemFb <- exp(tmp)
      ans <- (log1p(shape - onemFb) - log(shape) - tmp) / scale
      ans[ln.p > 0] <- NaN
    } else {
      tmp <- scale * log(p)
      onemFb <- exp(tmp)
      ans <- (log1p(shape - onemFb) - log(shape) - tmp) / scale
      ans[p < 0] <- NaN
      ans[p == 0] <- Inf
      ans[p == 1] <- 0
      ans[p > 1] <- NaN
    }
  }

  ans[shape <= 0 | scale <= 0] <- NaN
  ans
}  # qperks



rperks <- function(n, scale = 1, shape) {
  qperks(runif(n), scale = scale, shape = shape)
}





perks.control <- function(save.weights = TRUE, ...) {
  list(save.weights = save.weights)
}


 perks <-
  function(lscale = "loglink",    lshape = "loglink",
           iscale = NULL,      ishape = NULL,
           gscale = exp(-5:5), gshape = exp(-5:5),
           nsimEIM = 500,
           oim.mean = FALSE,
           zero = NULL, nowarning = FALSE) {



  lshape <- as.list(substitute(lshape))
  eshape <- link2list(lshape)
  lshape <- attr(eshape, "function.name")

  lscale <- as.list(substitute(lscale))
  escale <- link2list(lscale)
  lscale <- attr(escale, "function.name")


  if (!is.Numeric(nsimEIM, length.arg = 1,
                  integer.valued = TRUE))
    stop("bad input for argument 'nsimEIM'")
  if (nsimEIM <= 50)
    warning("argument 'nsimEIM' should be an integer ",
            "greater than 50, say")


  if (length(ishape))
    if (!is.Numeric(ishape, positive = TRUE))
      stop("argument 'ishape' values must be positive")
  if (length(iscale))
    if (!is.Numeric(iscale, positive = TRUE))
      stop("argument 'iscale' values must be positive")




    if (!is.logical(oim.mean) || length(oim.mean) != 1)
      stop("bad input for argument 'oim.mean'")



  new("vglmff",
  blurb = c("Perks' distribution\n\n",
            "Links:    ",
            namesof("scale", lscale, escale), ", ",
            namesof("shape", lshape, eshape), "\n",
  "Median:   qperks(p = 0.5, scale = scale, shape = shape)"),

  constraints = eval(substitute(expression({
    constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
                     M = M, M1 = 2,
                     predictors.names = predictors.names)
  }), list( .zero = zero ))),

  infos = eval(substitute(function(...) {
    list(M1 = 2,
         Q1 = 1,
         dpqrfun = "perks",
         nsimEIM = .nsimEIM ,
         parameters.names = c("scale", "shape"),
         zero = .zero )
  }, list( .zero = zero,
           .nsimEIM = nsimEIM ))),
  initialize = eval(substitute(expression({

    temp5 <-
    w.y.check(w = w, y = y,
              Is.positive.y = TRUE,
              ncol.w.max = Inf,
              ncol.y.max = Inf,
              out.wy = TRUE,
              colsyperw = 1,
              maximize = TRUE)
    w <- temp5$w
    y <- temp5$y



    ncoly <- ncol(y)
    M1 <- 2
    extra$ncoly <- ncoly
    extra$M1 <- M1
    M <- M1 * ncoly


    mynames1 <- param.names("scale", ncoly, skip1 = TRUE)
    mynames2 <- param.names("shape", ncoly, skip1 = TRUE)
    predictors.names <-
        c(namesof(mynames1, .lscale , .escale , tag = FALSE),
          namesof(mynames2, .lshape , .eshape , tag = FALSE))[
          interleave.VGAM(M, M1 = M1)]



    if (!length(etastart)) {

      matH <- matrix(if (length( .ishape )) .ishape else 0 + NA,
                     n, ncoly, byrow = TRUE)
      matC <- matrix(if (length( .iscale )) .iscale else 0 + NA,
                     n, ncoly, byrow = TRUE)

      shape.grid <- .gshape
      scale.grid <- .gscale

      for (spp. in 1:ncoly) {
        yvec <- y[, spp.]
        wvec <- w[, spp.]


        perks.Loglikfun2 <- function(scaleval, shapeval,
                                     y, x, w, extraargs) {
          sum(c(w) * dperks(x = y, shape = shapeval,
                            scale = scaleval, log = TRUE))
        }
        try.this <-
          grid.search2(scale.grid, shape.grid,
                       objfun = perks.Loglikfun2,
                       y = yvec, w = wvec,
                       ret.objfun = TRUE)  # Last value is \ell
        if (!length( .iscale ))
          matC[, spp.] <- try.this["Value1"]
        if (!length( .ishape ))
          matH[, spp.] <- try.this["Value2"]
      }  # spp.

      etastart <-
          cbind(theta2eta(matC, .lscale , .escale ),
                theta2eta(matH, .lshape , .eshape ))[,
                interleave.VGAM(M, M1 = M1)]
    }  # End of !length(etastart)
  }), list( .lscale = lscale, .lshape = lshape,
            .eshape = eshape, .escale = escale,
            .gshape = gshape, .gscale = gscale,
            .ishape = ishape, .iscale = iscale
            ))),

  linkinv = eval(substitute(function(eta, extra = NULL) {
    Scale <- eta2theta(eta[, c(TRUE, FALSE)], .lscale , .escale )
    Shape <- eta2theta(eta[, c(FALSE, TRUE)], .lshape , .eshape )

    qperks(p = 0.5, shape = Shape, scale = Scale)
  }, list( .lscale = lscale, .lshape = lshape,
           .escale = escale, .eshape = eshape ))),
  last = eval(substitute(expression({

    misc$link <-
      c(rep_len( .lscale , ncoly),
        rep_len( .lshape , ncoly))[interleave.VGAM(M, M1 = M1)]
    temp.names <- c(mynames1, mynames2)[interleave.VGAM(M,
                                                        M1 = M1)]
    names(misc$link) <- temp.names

    misc$earg <- vector("list", M)
    names(misc$earg) <- temp.names
    for (ii in 1:ncoly) {
      misc$earg[[M1*ii-1]] <- .escale
      misc$earg[[M1*ii  ]] <- .eshape
    }


    misc$M1 <- M1
    misc$expected <- TRUE
    misc$multipleResponses <- TRUE
    misc$nsimEIM <- .nsimEIM
  }), list( .lscale = lscale, .lshape = lshape,
            .escale = escale, .eshape = eshape,
            .nsimEIM = nsimEIM ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    Scale <- eta2theta(eta[, c(TRUE, FALSE)], .lscale , .escale )
    Shape <- eta2theta(eta[, c(FALSE, TRUE)], .lshape , .eshape )
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * dperks(x = y, shape = Shape,
                               scale = Scale, log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .lscale = lscale, .lshape = lshape,
           .escale = escale, .eshape = eshape ))),
  vfamily = c("perks"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    Scale <- eta2theta(eta[, c(TRUE, FALSE), drop = FALSE],
                       .lscale , .escale )
    Shape <- eta2theta(eta[, c(FALSE, TRUE), drop = FALSE],
                       .lshape , .eshape )
    okay1 <- all(is.finite(Scale)) && all(0 < Scale) &&
             all(is.finite(Shape)) && all(0 < Shape)
    okay1
  }, list( .lscale = lscale, .lshape = lshape,
           .escale = escale, .eshape = eshape) )),





  simslot = eval(substitute(
  function(object, nsim) {

    pwts <- if (length(pwts <- object@prior.weights) > 0)
              pwts else weights(object, type = "prior")
    if (any(pwts != 1))
      warning("ignoring prior weights")
    eta <- predict(object)
    Scale <- eta2theta(eta[, c(TRUE, FALSE)], .lscale , .escale )
    Shape <- eta2theta(eta[, c(FALSE, TRUE)], .lshape , .eshape )
    rperks(nsim * length(Scale), shape = Shape, scale = Scale)
  }, list( .lscale = lscale, .lshape = lshape,
           .escale = escale, .eshape = eshape ))),







  deriv = eval(substitute(expression({
    M1 <- 2
    scale <- eta2theta(eta[, c(TRUE, FALSE), drop = FALSE],
                       .lscale , .escale )
    shape <- eta2theta(eta[, c(FALSE, TRUE), drop = FALSE],
                       .lshape , .eshape )


    temp2 <- exp(y * scale)
    temp3 <- 1 + shape * temp2
    dl.dshape <- 1 / shape + 1 / (scale * (1 + shape)) -
                 (1 + 1 / scale) * temp2 / temp3
    dl.dscale <- y - log1p(shape) / scale^2 +
                 log1p(shape * temp2) / scale^2 -
                 (1 + 1 / scale) * shape * y * temp2 / temp3

    dshape.deta <- dtheta.deta(shape, .lshape , .eshape )
    dscale.deta <- dtheta.deta(scale, .lscale , .escale )

    dthetas.detas <- cbind(dscale.deta, dshape.deta)
    myderiv <- c(w) * cbind(dl.dscale, dl.dshape) * dthetas.detas
    myderiv[, interleave.VGAM(M, M1 = M1)]
  }), list( .lscale = lscale, .lshape = lshape,
            .escale = escale, .eshape = eshape ))),


  weight = eval(substitute(expression({

    NOS <- M / M1
    dThetas.detas <- dthetas.detas[, interleave.VGAM(M, M1 = M1)]

    wz <- matrix(0.0, n, M + M - 1)  # wz is 'tridiagonal'

    ind1 <- iam(NA, NA, M = M1, both = TRUE, diag = TRUE)


    for (spp. in 1:NOS) {
      run.varcov <- 0
      Scale <- scale[, spp.]
      Shape <- shape[, spp.]




      if (FALSE && intercept.only && .oim.mean ) {

 stop("this is wrong")
      temp8 <- (1 + Shape * exp(Scale * y[, spp.]))^2
      nd2l.dadb <- 2 * y[, spp.] * exp(Scale * y[, spp.]) / temp8

      nd2l.dada <- 1 / Shape^2 + 1 / (1 + Shape)^2 -
        2 * exp(2 * Scale * y[, spp.]) / temp8

      nd2l.dbdb <- 2 * Shape * y[, spp.]^2 *
                   exp(Scale * y[, spp.]) / temp8

      ave.oim11 <- weighted.mean(nd2l.dada, w[, spp.])
      ave.oim12 <- weighted.mean(nd2l.dadb, w[, spp.])
      ave.oim22 <- weighted.mean(nd2l.dbdb, w[, spp.])
      run.varcov <- cbind(ave.oim11, ave.oim22, ave.oim12)
    } else {

      for (ii in 1:( .nsimEIM )) {
        ysim <- rperks(n = n, shape = Shape, scale = Scale)
if (ii < 3) {
}

        temp2 <- exp(ysim * Scale)
        temp3 <- 1 + Shape * temp2
        dl.dshape <- 1 / Shape + 1 / (Scale * (1 + Shape)) -
                     (1 + 1 / Scale) * temp2 / temp3
        dl.dscale <- ysim - log1p(Shape) / Scale^2 +
                     log1p(Shape * temp2) / Scale^2 -
                     (1 + 1 / Scale) * Shape * ysim * temp2 / temp3


        temp7 <- cbind(dl.dscale, dl.dshape)
if (ii < 3) {
}
        run.varcov <- run.varcov +
                      temp7[, ind1$row.index] *
                      temp7[, ind1$col.index]
      }
      run.varcov <- cbind(run.varcov / .nsimEIM )

    }



      wz1 <- if (intercept.only)
          matrix(colMeans(run.varcov), n, ncol(run.varcov),
                 byrow = TRUE) else
          run.varcov

      wz1 <- wz1 * dThetas.detas[, M1 * (spp. - 1) + ind1$row] *
                   dThetas.detas[, M1 * (spp. - 1) + ind1$col]


      for (jay in 1:M1)
        for (kay in jay:M1) {
          cptr <- iam((spp. - 1) * M1 + jay,
                      (spp. - 1) * M1 + kay,
                      M = M)
          wz[, cptr] <- wz1[, iam(jay, kay, M = M1)]
        }
    }  # End of for (spp.) loop



    w.wz.merge(w = w, wz = wz, n = n, M = M, ndepy = M / M1)
  }), list( .lscale = lscale,
            .escale = escale,
            .nsimEIM = nsimEIM, .oim.mean = oim.mean ))))
}  # perks()









dmakeham <-
  function(x, scale = 1, shape, epsilon = 0, log = FALSE) {

  if (!is.logical(log.arg <- log) || length(log) != 1)
    stop("bad input for argument 'log'")
  rm(log)

  LLL <- max(length(x), length(shape), length(scale),
             length(epsilon))
  if (length(x)       != LLL) x       <- rep_len(x,       LLL)
  if (length(shape)   != LLL) shape   <- rep_len(shape,   LLL)
  if (length(scale)   != LLL) scale   <- rep_len(scale,   LLL)
  if (length(epsilon) != LLL) epsilon <- rep_len(epsilon, LLL)

  index0 <- (x < 0)
  ans <- log(epsilon * exp(-x * scale) + shape) +
         x * (scale - epsilon) -
         (shape / scale) * expm1(x * scale)
  ans[index0] <- log(0)
  ans[x == Inf] <- log(0)
  if (log.arg) {
  } else {
    ans <- exp(ans)
    ans[index0] <- 0
    ans[x == Inf] <- 0
  }
  ans[shape <= 0 | scale <= 0 | epsilon < 0] <- NaN
  ans
}  # dmakeham



pmakeham <- function(q, scale = 1, shape, epsilon = 0,
                     lower.tail = TRUE, log.p = FALSE) {


  if (!is.logical(lower.tail) || length(lower.tail ) != 1)
    stop("bad input for argument 'lower.tail'")

  if (!is.logical(log.p) || length(log.p) != 1)
    stop("bad input for argument 'log.p'")

  LLL <- max(length(q), length(shape), length(scale),
             length(epsilon))
  if (length(q)       != LLL) q       <- rep_len(q,       LLL)
  if (length(shape)   != LLL) shape   <- rep_len(shape,   LLL)
  if (length(scale)   != LLL) scale   <- rep_len(scale,   LLL)
  if (length(epsilon) != LLL) epsilon <- rep_len(epsilon, LLL)

  if (lower.tail) {
    if (log.p) {
      ans <- log(-expm1(-q * epsilon - (shape / scale) *
                        expm1(scale * q)))
      ans[q <= 0 ] <- -Inf
      ans[q == Inf] <- 0
    } else {
      ans <- -expm1(-q * epsilon - (shape / scale) *
                    expm1(scale * q))
      ans[q <= 0] <- 0
      ans[q == Inf] <- 1
    }
  } else {
    if (log.p) {
      ans <- -q * epsilon - (shape / scale) * expm1(scale * q)
      ans[q <= 0] <- 0
      ans[q == Inf] <- -Inf
    } else {
      ans <- exp(-q * epsilon - (shape / scale) * expm1(scale * q))
      ans[q <= 0] <- 1
      ans[q == Inf] <- 0
    }
  }

  ans[shape <= 0 | scale <= 0 | epsilon < 0] <- NaN
  ans
}  # pmakeham



qmakeham <- function(p, scale = 1, shape, epsilon = 0,
                     lower.tail = TRUE, log.p = FALSE) {

  if (!is.logical(lower.tail) || length(lower.tail ) != 1)
    stop("bad input for argument 'lower.tail'")

  if (!is.logical(log.p) || length(log.p) != 1)
    stop("bad input for argument 'log.p'")

  LLL <- max(length(p), length(shape), length(scale),
             length(epsilon))
  if (length(p)       != LLL) p       <- rep_len(p,       LLL)
  if (length(shape)   != LLL) shape   <- rep_len(shape,   LLL)
  if (length(scale)   != LLL) scale   <- rep_len(scale,   LLL)
  if (length(epsilon) != LLL) epsilon <- rep_len(epsilon, LLL)


  if (lower.tail) {
    if (log.p) {
      ln.p <- p
      ans <- shape / (scale * epsilon) -
        log(-expm1(ln.p)) / epsilon -
        lambertW((shape / epsilon) * exp(shape / epsilon) *
           exp(log(-expm1(ln.p)) * (-scale / epsilon))) / scale
      ans[ln.p == 0] <- Inf
      ans[ln.p > 0] <- NaN
    } else {
      ans <- shape / (scale * epsilon) - log1p(-p) / epsilon -
        lambertW((shape / epsilon) * exp(shape / epsilon) *
                   exp( (-scale / epsilon) * log1p(-p) )) / scale
      ans[p < 0] <- NaN
      ans[p == 0] <- 0
      ans[p == 1] <- Inf
      ans[p > 1] <- NaN
    }
  } else {
    if (log.p) {
      ln.p <- p
      ans <-  shape / (scale * epsilon) - ln.p / epsilon -
        lambertW((shape / epsilon) * exp(shape / epsilon) *
                  exp(ln.p * (-scale / epsilon))) / scale
      ans[ln.p == -Inf] <- Inf
      ans[ln.p > 0] <- NaN
    } else {
      ans <- shape / (scale * epsilon) - log(p) / epsilon -
        lambertW((shape / epsilon) * exp(shape / epsilon) *
                  p^(-scale / epsilon)) / scale
      ans[p < 0] <- NaN
      ans[p == 0] <- Inf
      ans[p == 1] <- 0
      ans[p > 1] <- NaN
    }
  }

  ans[epsilon == 0] <-
    qgompertz(p     =     p[epsilon == 0],
              shape = shape[epsilon == 0],
              scale = scale[epsilon == 0],
              lower.tail = lower.tail,
              log.p = log.p)

  ans[shape <= 0 | scale <= 0 | epsilon < 0] <- NaN
  ans
}  # qmakeham



rmakeham <- function(n, scale = 1, shape, epsilon = 0) {
  qmakeham(runif(n), scale = scale, shape, epsilon = epsilon)
}




makeham.control <- function(save.weights = TRUE, ...) {
  list(save.weights = save.weights)
}



 makeham <-
  function(lscale = "loglink", lshape = "loglink",
           lepsilon = "loglink",
           iscale = NULL,   ishape = NULL,
           iepsilon = NULL,  # 0.3,
           gscale = exp(-5:5),
           gshape = exp(-5:5),
           gepsilon = exp(-4:1),
           nsimEIM = 500,
           oim.mean = TRUE,
           zero = NULL, nowarning = FALSE) {







  lepsil <- lepsilon
  iepsil <- iepsilon


  lshape <- as.list(substitute(lshape))
  eshape <- link2list(lshape)
  lshape <- attr(eshape, "function.name")

  lscale <- as.list(substitute(lscale))
  escale <- link2list(lscale)
  lscale <- attr(escale, "function.name")

  lepsil <- as.list(substitute(lepsil))
  eepsil <- link2list(lepsil)
  lepsil <- attr(eepsil, "function.name")

  if (!is.Numeric(nsimEIM, length.arg = 1,
                  integer.valued = TRUE))
    stop("bad input for argument 'nsimEIM'")
  if (nsimEIM <= 50)
    warning("argument 'nsimEIM' should be an integer ",
            "greater than 50, say")


  if (length(ishape))
    if (!is.Numeric(ishape, positive = TRUE))
      stop("argument 'ishape' values must be positive")
  if (length(iscale))
    if (!is.Numeric(iscale, positive = TRUE))
      stop("argument 'iscale' values must be positive")
  if (length(iepsil))
    if (!is.Numeric(iepsil, positive = TRUE))
      stop("argument 'iepsil' values must be positive")





    if (!is.logical(oim.mean) || length(oim.mean) != 1)
      stop("bad input for argument 'oim.mean'")




  new("vglmff",
  blurb = c("Makeham distribution\n\n",
            "Links:    ",
            namesof("scale",   lscale, escale), ", ",
            namesof("shape",   lshape, eshape), ", ",
            namesof("epsilon", lepsil, eepsil), "\n",
            "Median:   qmakeham(p = 0.5, scale, shape, epsilon)"),

  constraints = eval(substitute(expression({
    constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
                     M = M, M1 = 3,
                     predictors.names = predictors.names)
  }), list( .zero = zero ))),

  infos = eval(substitute(function(...) {
    list(M1 = 3,
         Q1 = 1,
         dpqrfun = "makeham",
         nsimEIM = .nsimEIM,
         parameters.names = c("scale", "shape"),
         zero = .zero )
  }, list( .zero = zero,
           .nsimEIM = nsimEIM ))),
  initialize = eval(substitute(expression({


    temp5 <-
    w.y.check(w = w, y = y,
              Is.positive.y = TRUE,
              ncol.w.max = Inf,
              ncol.y.max = Inf,
              out.wy = TRUE,
              colsyperw = 1,
              maximize = TRUE)
    w <- temp5$w
    y <- temp5$y


    ncoly <- ncol(y)

    M1 <- 3
    extra$ncoly <- ncoly
    extra$M1 <- M1
    M <- M1 * ncoly


    mynames1 <- param.names("scale",   ncoly, skip1 = TRUE)
    mynames2 <- param.names("shape",   ncoly, skip1 = TRUE)
    mynames3 <- param.names("epsilon", ncoly, skip1 = TRUE)
    predictors.names <-
        c(namesof(mynames1, .lscale , .escale , tag = FALSE),
          namesof(mynames2, .lshape , .eshape , tag = FALSE),
          namesof(mynames3, .lepsil , .eepsil , tag = FALSE))[
          interleave.VGAM(M, M1 = M1)]


    if (!length(etastart)) {

      matC <- matrix(if (length( .iscale )) .iscale else 0 + NA,
                     n, ncoly, byrow = TRUE)
      matH <- matrix(if (length( .ishape )) .ishape else 0 + NA,
                     n, ncoly, byrow = TRUE)

      matE <- matrix(if (length( .iepsil )) .iepsil else 0.3,
                     n, ncoly, byrow = TRUE)


      shape.grid <- unique(sort(c( .gshape )))
      scale.grid <- unique(sort(c( .gscale )))



      for (spp. in 1:ncoly) {
        yvec <- y[, spp.]
        wvec <- w[, spp.]


        makeham.Loglikfun2 <- function(scaleval, shapeval,
                                       y, x, w, extraargs) {
          sum(c(w) * dmakeham(x = y, shape = shapeval,
                              epsilon = extraargs$Epsil,
                              scale = scaleval, log = TRUE))
        }
        try.this <-
          grid.search2(scale.grid, shape.grid,
                       objfun = makeham.Loglikfun2,
                       y = yvec, w = wvec,
                       extraargs = list(Epsilon = matE[1, spp.]),
                       ret.objfun = TRUE)  # Last value is \ell
        if (!length( .iscale ))
          matC[, spp.] <- try.this["Value1"]
        if (!length( .ishape ))
          matH[, spp.] <- try.this["Value2"]
      }  # spp.





      epsil.grid <- c( .gepsil )
      for (spp. in 1:ncoly) {
        yvec <- y[, spp.]
        wvec <- w[, spp.]

        makeham.Loglikfun2 <-
          function(epsilval, y, x, w, extraargs) {
          ans <-
          sum(c(w) * dmakeham(x = y, shape = extraargs$Shape,
                              epsilon = epsilval, log = TRUE,
                              scale = extraargs$Scale))
          ans
        }
        Init.epsil <-
            grid.search(epsil.grid,
                        objfun = makeham.Loglikfun2,
                        y = yvec, x = x, w = wvec,
                        extraargs = list(Shape = matH[1, spp.],
                                         Scale = matC[1, spp.]))

        matE[, spp.] <- Init.epsil
      }  # spp.


      etastart <- cbind(theta2eta(matC, .lscale , .escale ),
                        theta2eta(matH, .lshape , .eshape ),
                        theta2eta(matE, .lepsil , .eepsil ))[,
                        interleave.VGAM(M, M1 = M1)]
    }  # End of !length(etastart)
  }),
  list(
       .lshape = lshape, .lscale = lscale, .lepsil = lepsil,
       .eshape = eshape, .escale = escale, .eepsil = eepsil,
       .gshape = gshape, .gscale = gscale, .gepsil = gepsilon,
       .ishape = ishape, .iscale = iscale, .iepsil = iepsil
       ))),

  linkinv = eval(substitute(function(eta, extra = NULL) {
    scale <- eta2theta(eta[, c(TRUE, FALSE, FALSE)],
                       .lscale , .escale )
    shape <- eta2theta(eta[, c(FALSE, TRUE, FALSE)],
                       .lshape , .eshape )
    epsil <- eta2theta(eta[, c(FALSE, FALSE, TRUE)],
                       .lepsil , .eepsil )
    qmakeham(p = 0.5, scale = scale, shape, epsil = epsil)
  },
  list( .lshape = lshape, .lscale = lscale, .lepsil = lepsil,
        .eshape = eshape, .escale = escale, .eepsil = eepsil ))),
  last = eval(substitute(expression({
    M1 <- extra$M1
    misc$link <-
      c(rep_len( .lscale , ncoly),
        rep_len( .lshape , ncoly),
        rep_len( .lepsil , ncoly))[interleave.VGAM(M, M1 = M1)]
    temp.names <- c(mynames1, mynames2, mynames3)[
                    interleave.VGAM(M, M1 = M1)]
    names(misc$link) <- temp.names

    misc$earg <- vector("list", M)
    names(misc$earg) <- temp.names
    for (ii in 1:ncoly) {
      misc$earg[[M1*ii-2]] <- .escale
      misc$earg[[M1*ii-1]] <- .eshape
      misc$earg[[M1*ii  ]] <- .eepsil
    }

    misc$M1 <- M1
    misc$expected <- TRUE
    misc$multipleResponses <- TRUE
    misc$nsimEIM <- .nsimEIM
  }),
  list( .lshape = lshape, .lscale = lscale, .lepsil = lepsil,
        .eshape = eshape, .escale = escale, .eepsil = eepsil,
        .nsimEIM = nsimEIM ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
      scale <- eta2theta(eta[, c(TRUE, FALSE, FALSE)],
                         .lscale , .escale )
      shape <- eta2theta(eta[, c(FALSE, TRUE, FALSE)],
                         .lshape , .eshape )
      epsil <- eta2theta(eta[, c(FALSE, FALSE, TRUE)],
                         .lepsil , .eepsil )
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * dmakeham(y, scale = scale, shape,
                                 epsil = epsil, log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  },
  list( .lshape = lshape, .lscale = lscale, .lepsil = lepsil,
        .eshape = eshape, .escale = escale, .eepsil = eepsil ))),
  vfamily = c("makeham"),
 # ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    Scale <- eta2theta(eta[, c(TRUE, FALSE, FALSE), drop = FALSE],
                       .lscale , .escale )
    shape <- eta2theta(eta[, c(FALSE, TRUE, FALSE), drop = FALSE],
                       .lshape , .eshape )
    epsil <- eta2theta(eta[, c(FALSE, FALSE, TRUE), drop = FALSE],
                       .lepsil , .eepsil )
    okay1 <- all(is.finite(Scale)) && all(0 < Scale) &&
             all(is.finite(shape)) && all(0 < shape) &&
             all(is.finite(epsil)) && all(0 < epsil)
    okay1
  },
  list( .lshape = lshape, .lscale = lscale, .lepsil = lepsil,
        .eshape = eshape, .escale = escale, .eepsil = eepsil ))),




  simslot = eval(substitute(
  function(object, nsim) {

    pwts <- if (length(pwts <- object@prior.weights) > 0)
              pwts else weights(object, type = "prior")
    if (any(pwts != 1))
      warning("ignoring prior weights")
    eta <- predict(object)
    Scale <- eta2theta(eta[, c(TRUE, FALSE, FALSE), drop = FALSE],
                       .lscale , .escale )
    shape <- eta2theta(eta[, c(FALSE, TRUE, FALSE), drop = FALSE],
                       .lshape , .eshape )
    epsil <- eta2theta(eta[, c(FALSE, FALSE, TRUE), drop = FALSE],
                       .lepsil , .eepsil )
    rmakeham(nsim * length(Scale),
             scale = c(Scale), shape = c(shape), epsilon = c(epsil))
  }, list( .lshape = lshape, .lscale = lscale, .lepsil = lepsil,
           .eshape = eshape, .escale = escale, .eepsil = eepsil ))),




  deriv = eval(substitute(expression({
    scale <- eta2theta(eta[, c(TRUE, FALSE, FALSE), drop = FALSE],
                       .lscale , .escale )
    shape <- eta2theta(eta[, c(FALSE, TRUE, FALSE), drop = FALSE],
                       .lshape , .eshape )
    epsil <- eta2theta(eta[, c(FALSE, FALSE, TRUE), drop = FALSE],
                       .lepsil , .eepsil )

    temp2 <- exp(y * scale)
    temp3 <- epsil + shape * temp2
    dl.dshape <- temp2 / temp3 - expm1(y * scale) / scale
    dl.dscale <- shape * y * temp2 / temp3 +
                 shape * expm1(y * scale) / scale^2 -
                 shape * y * temp2 / scale
    dl.depsil <- 1 / temp3 - y

    dshape.deta <- dtheta.deta(shape, .lshape , .eshape )
    dscale.deta <- dtheta.deta(scale, .lscale , .escale )
    depsil.deta <- dtheta.deta(epsil, .lepsil , .eepsil )

    dthetas.detas <- cbind(dscale.deta, dshape.deta, depsil.deta)
    myderiv <- c(w) * cbind(dl.dscale,
                            dl.dshape,
                            dl.depsil) * dthetas.detas
    myderiv[, interleave.VGAM(M, M1 = M1)]
  }),
  list( .lshape = lshape, .lscale = lscale, .lepsil = lepsil,
        .eshape = eshape, .escale = escale, .eepsil = eepsil ))),

  weight = eval(substitute(expression({
    NOS <- M / M1
    dThetas.detas <- dthetas.detas[, interleave.VGAM(M, M1 = M1)]
    wz <- matrix(0.0, n, M + M - 1 + M - 2)  # wz has half-bandwidth 3

    ind1 <- iam(NA, NA, M = M1, both = TRUE, diag = TRUE)  # Use SFS

    for (spp. in 1:NOS) {
      run.varcov <- 0
      Shape <- shape[, spp.]
      Scale <- scale[, spp.]
      Epsil <- epsil[, spp.]

      for (ii in 1:( .nsimEIM )) {
        ysim <- rmakeham(n, scale = Scale, Shape, epsil = Epsil)
        temp2 <- exp(ysim * Scale)
        temp3 <- Epsil + Shape * temp2
        dl.dshape <- temp2 / temp3 - expm1(ysim * Scale) / Scale
        dl.dscale <- Shape * ysim * temp2 / temp3 +
                     Shape * expm1(ysim * Scale) / Scale^2 -
                     Shape * ysim * temp2 / Scale
        dl.depsil <- 1 / temp3 - ysim

        temp7 <- cbind(dl.dscale, dl.dshape, dl.depsil)
        run.varcov <- run.varcov + temp7[, ind1$row.index] *
                                   temp7[, ind1$col.index]
      }
      run.varcov <- cbind(run.varcov / .nsimEIM )


      wz1 <- if (intercept.only)
        matrix(colMeans(run.varcov, na.rm = TRUE),
               n, ncol(run.varcov), byrow = TRUE) else
        run.varcov

      wz1 <- wz1 * dThetas.detas[, M1 * (spp. - 1) + ind1$row] *
                   dThetas.detas[, M1 * (spp. - 1) + ind1$col]
      for (jay in 1:M1)
        for (kay in jay:M1) {  # Now copy wz1 into wz
          cptr <- iam((spp. - 1) * M1 + jay,
                      (spp. - 1) * M1 + kay, M = M)
          wz[, cptr] <- wz1[, iam(jay, kay, M = M1)]
        }
    }  # End of for (spp.) loop
    w.wz.merge(w = w, wz = wz, n = n, M = M, ndepy = M / M1)
  }),
  list( .lshape = lshape, .lscale = lscale, .lepsil = lepsil,
        .eshape = eshape, .escale = escale, .eepsil = eepsil,
        .nsimEIM = nsimEIM, .oim.mean = oim.mean ))))
}  # makeham()










dgompertz <- function(x, scale = 1, shape, log = FALSE) {

  if (!is.logical(log.arg <- log) || length(log) != 1)
    stop("bad input for argument 'log'")
  rm(log)

  LLL <- max(length(x), length(shape), length(scale))
  if (length(x)     != LLL) x     <- rep_len(x,     LLL)
  if (length(shape) != LLL) shape <- rep_len(shape, LLL)
  if (length(scale) != LLL) scale <- rep_len(scale, LLL)


  index0 <- (x < 0)
  index1 <- abs(x * scale) < 0.1 & is.finite(x * scale)
  ans <- log(shape) + x * scale -
         (shape / scale) * (exp(x * scale) - 1)
  ans[index1] <- log(shape[index1]) + x[index1] * scale[index1] -
                 (shape[index1] / scale[index1]) *
                 expm1(x[index1] * scale[index1])
  ans[index0] <- log(0)
  ans[x == Inf] <- log(0)
  if (log.arg) {
  } else {
    ans <- exp(ans)
    ans[index0] <- 0
    ans[x == Inf] <- 0
  }
  ans[shape <= 0 | scale <= 0] <- NaN
  ans
}  # dgompertz



pgompertz <- function(q, scale = 1, shape,
                      lower.tail = TRUE, log.p = FALSE) {


  if (!is.logical(lower.tail) || length(lower.tail ) != 1)
    stop("bad input for argument 'lower.tail'")

  if (!is.logical(log.p) || length(log.p) != 1)
    stop("bad input for argument 'log.p'")

  LLL <- max(length(q), length(shape), length(scale))
  if (length(q)       != LLL) q       <- rep_len(q,       LLL)
  if (length(shape)   != LLL) shape   <- rep_len(shape,   LLL)
  if (length(scale)   != LLL) scale   <- rep_len(scale,   LLL)


  if (lower.tail) {
    if (log.p) {
      ans <- log1p(-exp((-shape / scale) * expm1(scale * q)))
      ans[q <= 0 ] <- -Inf
      ans[q == Inf] <- 0
    } else {
      ans <- -expm1((-shape / scale) * expm1(scale * q))
      ans[q <= 0] <- 0
      ans[q == Inf] <- 1
    }
  } else {
    if (log.p) {
      ans <- (-shape / scale) * expm1(scale * q)
      ans[q <= 0] <- 0
      ans[q == Inf] <- -Inf
    } else {
      ans <- exp((-shape / scale) * expm1(scale * q))
      ans[q <= 0] <- 1
      ans[q == Inf] <- 0
    }
  }
  ans[shape <= 0 | scale <= 0] <- NaN
  ans
}  # pgompertz



qgompertz <- function(p, scale = 1, shape,
                      lower.tail = TRUE, log.p = FALSE) {

  if (!is.logical(lower.tail) || length(lower.tail ) != 1)
    stop("bad input for argument 'lower.tail'")
  if (!is.logical(log.p) || length(log.p) != 1)
    stop("bad input for argument 'log.p'")

  LLL <- max(length(p), length(shape), length(scale))
  if (length(p)       != LLL) p       <- rep_len(p,       LLL)
  if (length(shape)   != LLL) shape   <- rep_len(shape,   LLL)
  if (length(scale)   != LLL) scale   <- rep_len(scale,   LLL)

  if (lower.tail) {
    if (log.p) {
      ln.p <- p
      ans <- log1p((-scale / shape) * log(-expm1(ln.p))) / scale
      ans[ln.p > 0] <- NaN
    } else {
      ans <- log1p((-scale / shape) * log1p(-p)) / scale
      ans[p < 0] <- NaN
      ans[p == 0] <- 0
      ans[p == 1] <- Inf
      ans[p > 1] <- NaN
    }
  } else {
    if (log.p) {
      ln.p <- p
      ans <- log1p((-scale / shape) * ln.p) / scale
      ans[ln.p > 0] <- NaN
    } else {
      ans <- log1p((-scale / shape) * log(p)) / scale
      ans[p < 0] <- NaN
      ans[p == 0] <- Inf
      ans[p == 1] <- 0
      ans[p > 1] <- NaN
    }
  }
  ans[shape <= 0 | scale <= 0] <- NaN
  ans
}  # qgompertz





rgompertz <- function(n, scale = 1, shape) {
  qgompertz(runif(n), scale = scale, shape = shape)
}







gompertz.control <- function(save.weights = TRUE, ...) {
  list(save.weights = save.weights)
}


 gompertz <-
  function(lscale = "loglink", lshape = "loglink",
           iscale = NULL,   ishape = NULL,
           nsimEIM = 500,
           zero = NULL, nowarning = FALSE) {





  lshape <- as.list(substitute(lshape))
  eshape <- link2list(lshape)
  lshape <- attr(eshape, "function.name")

  lscale <- as.list(substitute(lscale))
  escale <- link2list(lscale)
  lscale <- attr(escale, "function.name")



  if (!is.Numeric(nsimEIM, length.arg = 1,
                  integer.valued = TRUE))
    stop("bad input for argument 'nsimEIM'")
  if (nsimEIM <= 50)
    warning("argument 'nsimEIM' should be an integer ",
            "greater than 50, say")


  if (length(ishape))
    if (!is.Numeric(ishape, positive = TRUE))
      stop("argument 'ishape' values must be positive")
  if (length(iscale))
    if (!is.Numeric(iscale, positive = TRUE))
      stop("argument 'iscale' values must be positive")





  new("vglmff",
  blurb = c("Gompertz distribution\n\n",
            "Links:    ",
            namesof("scale", lscale, escale ), ", ",
            namesof("shape", lshape, eshape ), "\n",
            "Median:     scale * log(2 - 1 / shape)"),

  constraints = eval(substitute(expression({
    constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
                     M = M, M1 = 2,
                     predictors.names = predictors.names)
  }), list( .zero = zero ))),

  infos = eval(substitute(function(...) {
    list(M1 = 2,
         Q1 = 1,
         dpqrfun = "gompertz",
         nsimEIM = .nsimEIM,
         parameters.names = c("scale", "shape"),
         zero = .zero )
  }, list( .zero = zero,
           .nsimEIM = nsimEIM ))),
  initialize = eval(substitute(expression({

    temp5 <-
    w.y.check(w = w, y = y,
              Is.positive.y = TRUE,
              ncol.w.max = Inf,
              ncol.y.max = Inf,
              out.wy = TRUE,
              colsyperw = 1,
              maximize = TRUE)
    w <- temp5$w
    y <- temp5$y



    ncoly <- ncol(y)
    M1 <- 2
    extra$ncoly <- ncoly
    extra$M1 <- M1
    M <- M1 * ncoly


    mynames1 <- param.names("scale", ncoly, skip1 = TRUE)
    mynames2 <- param.names("shape", ncoly, skip1 = TRUE)
    predictors.names <-
        c(namesof(mynames1, .lscale , .escale , tag = FALSE),
          namesof(mynames2, .lshape , .eshape , tag = FALSE))[
          interleave.VGAM(M, M1 = M1)]



    if (!length(etastart)) {

      matH <- matrix(if (length( .ishape )) .ishape else 0 + NA,
                     n, ncoly, byrow = TRUE)
      matC <- matrix(if (length( .iscale )) .iscale else 0 + NA,
                     n, ncoly, byrow = TRUE)

      shape.grid <- c(exp(-seq(4, 0.1, len = 07)), 1,
                      exp( seq(0.1, 4, len = 07)))
      scale.grid <- c(exp(-seq(4, 0.1, len = 07)), 1,
                      exp( seq(0.1, 4, len = 07)))

      for (spp. in 1:ncoly) {
        yvec <- y[, spp.]
        wvec <- w[, spp.]


        gompertz.Loglikfun <-
          function(scaleval, y, x, w, extraargs) {
          ans <-
          sum(c(w) * dgompertz(x = y, shape = extraargs$Shape,
                               scale = scaleval, log = TRUE))
          ans
        }

        mymat <- matrix(-1, length(shape.grid), 2)
        for (jlocal in seq_along(shape.grid)) {
          mymat[jlocal, ] <-
            grid.search(scale.grid,
                        objfun = gompertz.Loglikfun,
                        y = yvec, x = x, w = wvec,
                        ret.objfun = TRUE,
               extraargs = list(Shape = shape.grid[jlocal]))
        }
        index.shape <- which(mymat[, 2] == max(mymat[, 2]))[1]

        if (!length( .ishape ))
          matH[, spp.] <- shape.grid[index.shape]
        if (!length( .iscale ))
          matC[, spp.] <- mymat[index.shape, 1]
      }  # spp.

      etastart <- cbind(theta2eta(matC, .lscale , .escale ),
                        theta2eta(matH, .lshape , .eshape ))[,
                        interleave.VGAM(M, M1 = M1)]
    }  # End of !length(etastart)
  }), list( .lshape = lshape, .lscale = lscale,
            .eshape = eshape, .escale = escale,
            .ishape = ishape, .iscale = iscale
          ))),

  linkinv = eval(substitute(function(eta, extra = NULL) {
    scale <- eta2theta(eta[, c(TRUE, FALSE)], .lscale , .escale )
    shape <- eta2theta(eta[, c(FALSE, TRUE)], .lshape , .eshape )
    log1p((scale / shape) * log(2)) / scale
  }, list( .lshape = lshape, .lscale = lscale,
           .eshape = eshape, .escale = escale ))),
  last = eval(substitute(expression({
    M1 <- extra$M1
    misc$link <-
      c(rep_len( .lscale , ncoly),
        rep_len( .lshape , ncoly))[interleave.VGAM(M, M1 = M1)]
    temp.names <- c(mynames1, mynames2)[interleave.VGAM(M,
                                                        M1 = M1)]
    names(misc$link) <- temp.names

    misc$earg <- vector("list", M)
    names(misc$earg) <- temp.names
    for (ii in 1:ncoly) {
      misc$earg[[M1*ii-1]] <- .escale
      misc$earg[[M1*ii  ]] <- .eshape
    }

    misc$M1 <- M1
    misc$expected <- TRUE
    misc$multipleResponses <- TRUE
    misc$nsimEIM <- .nsimEIM
  }), list( .lshape = lshape, .lscale = lscale,
            .eshape = eshape, .escale = escale,
            .nsimEIM = nsimEIM ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    scale <- eta2theta(eta[, c(TRUE, FALSE)], .lscale , .escale )
    shape <- eta2theta(eta[, c(FALSE, TRUE)], .lshape , .eshape )
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * dgompertz(x = y, scale = scale,
                                  shape = shape, log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
    }, list( .lshape = lshape, .lscale = lscale,
             .eshape = eshape, .escale = escale ))),
  vfamily = c("gompertz"),
 # ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    Scale <- eta2theta(eta[, c(TRUE, FALSE), drop = FALSE],
                       .lscale , .escale )
    shape <- eta2theta(eta[, c(FALSE, TRUE), drop = FALSE],
                       .lshape , .eshape )
    okay1 <- all(is.finite(Scale)) && all(0 < Scale) &&
             all(is.finite(shape)) && all(0 < shape)
    okay1
  }, list( .lshape = lshape, .lscale = lscale,
           .eshape = eshape, .escale = escale ))),





  simslot = eval(substitute(
  function(object, nsim) {

    pwts <- if (length(pwts <- object@prior.weights) > 0)
              pwts else weights(object, type = "prior")
    if (any(pwts != 1))
      warning("ignoring prior weights")
    eta <- predict(object)
    Scale <- eta2theta(eta[, c(TRUE, FALSE)], .lscale , .escale )
    Shape <- eta2theta(eta[, c(FALSE, TRUE)], .lshape , .eshape )
    rgompertz(nsim * length(Scale),
              shape = c(Shape), scale = c(Scale))
    }, list( .lshape = lshape, .lscale = lscale,
             .eshape = eshape, .escale = escale ))),



  deriv = eval(substitute(expression({
    M1 <- 2
    scale <- eta2theta(eta[, c(TRUE, FALSE), drop = FALSE],
                       .lscale , .escale )
    shape <- eta2theta(eta[, c(FALSE, TRUE), drop = FALSE],
                       .lshape , .eshape )


    temp2 <- exp(y * scale)
    temp4 <- -expm1(y * scale)
    dl.dshape <- 1 / shape + temp4 / scale
    dl.dscale <- y * (1 - shape * temp2 / scale) -
                 shape * temp4 / scale^2

    dscale.deta <- dtheta.deta(scale, .lscale , .escale )
    dshape.deta <- dtheta.deta(shape, .lshape , .eshape )

    dthetas.detas <- cbind(dscale.deta, dshape.deta)
    myderiv <- c(w) * cbind(dl.dscale, dl.dshape) * dthetas.detas
    myderiv[, interleave.VGAM(M, M1 = M1)]
  }), list( .lshape = lshape, .lscale = lscale,
            .eshape = eshape, .escale = escale ))),


  weight = eval(substitute(expression({

    NOS <- M / M1
    dThetas.detas <- dthetas.detas[, interleave.VGAM(M, M1 = M1)]

    wz <- matrix(0.0, n, M + M - 1)  # wz is 'tridiagonal'

    ind1 <- iam(NA, NA, M = M1, both = TRUE, diag = TRUE)


    for (spp. in 1:NOS) {
      run.varcov <- 0
      Shape <- shape[, spp.]
      Scale <- scale[, spp.]

      for (ii in 1:( .nsimEIM )) {
        ysim <- rgompertz(n = n, shape = Shape, scale = Scale)
if (ii < 3) {
}

        temp2 <- exp(ysim * scale)
        temp4 <- -expm1(ysim * scale)
        dl.dshape <- 1 / shape + temp4 / scale
        dl.dscale <- ysim * (1 - shape * temp2 / scale) -
                     shape * temp4 / scale^2


        temp7 <- cbind(dl.dscale, dl.dshape)
        run.varcov <- run.varcov +
                      temp7[, ind1$row.index] *
                      temp7[, ind1$col.index]
      }
      run.varcov <- cbind(run.varcov / .nsimEIM )

      wz1 <- if (intercept.only)
          matrix(colMeans(run.varcov),
                 nrow = n, ncol = ncol(run.varcov), byrow = TRUE) else
          run.varcov

      wz1 <- wz1 * dThetas.detas[, M1 * (spp. - 1) + ind1$row] *
                   dThetas.detas[, M1 * (spp. - 1) + ind1$col]


      for (jay in 1:M1)
        for (kay in jay:M1) {
          cptr <- iam((spp. - 1) * M1 + jay,
                      (spp. - 1) * M1 + kay,
                      M = M)
          wz[, cptr] <- wz1[, iam(jay, kay, M = M1)]
        }
    }  # End of for (spp.) loop



    w.wz.merge(w = w, wz = wz, n = n, M = M, ndepy = M / M1)
  }), list( .lscale = lscale,
            .escale = escale,
            .nsimEIM = nsimEIM ))))
}  # gompertz()








 dmoe <- function (x, alpha = 1, lambda = 1, log = FALSE) {
  if (!is.logical(log.arg <- log) || length(log) != 1)
    stop("bad input for argument 'log'")
  rm(log)

  LLL <- max(length(x),
             length(alpha),
             length(lambda))
  if (length(x)      != LLL) x      <- rep_len(x,      LLL)
  if (length(alpha)  != LLL) alpha  <- rep_len(alpha,  LLL)
  if (length(lambda) != LLL) lambda <- rep_len(lambda, LLL)

  index0 <- (x < 0)
  if (log.arg) {
    ans <- log(lambda) + (lambda * x) -
           2 * log(expm1(lambda * x) + alpha)
    ans[index0] <- log(0)
  } else {
    ans <- lambda * exp(lambda * x) / (expm1(lambda * x) + alpha)^2
    ans[index0] <- 0
  }
  ans[alpha <= 0 | lambda <= 0] <- NaN
  ans
}



 pmoe <- function (q, alpha = 1, lambda = 1) {
  ret <- ifelse(alpha <= 0 | lambda <= 0, NaN,
                1 - 1 / (expm1(lambda * q) + alpha))
  ret[q < log(2 - alpha) / lambda] <- 0
  ret
}



qmoe <- function (p, alpha = 1, lambda = 1) {
  ifelse(p < 0 | p > 1 | alpha <= 0 | lambda <= 0, NaN,
        log1p(-alpha + 1 / (1 - p)) / lambda)
}



rmoe <- function (n, alpha = 1, lambda = 1) {
  qmoe(p = runif(n), alpha = alpha, lambda = lambda)
}






exponential.mo.control <- function(save.weights = TRUE, ...) {
  list(save.weights = save.weights)
}




 exponential.mo <-
  function(lalpha = "loglink", llambda = "loglink",
           ealpha = list(), elambda = list(),
           ialpha = 1,      ilambda = NULL,
           imethod = 1,
           nsimEIM = 200,
           zero = NULL) {

  stop("fundamentally unable to estimate the parameters as ",
       "the support of the density depends on the parameters")


  lalpha <- as.list(substitute(lalpha))
  ealpha <- link2list(lalpha)
  lalpha <- attr(ealpha, "function.name")

  llambda <- as.list(substitute(llambda))
  elambda <- link2list(llambda)
  llambda <- attr(elambda, "function.name")

  lalpha0 <- lalpha
  ealpha0 <- ealpha
  ialpha0 <- ialpha



  if (!is.Numeric(nsimEIM, length.arg = 1,
                  integer.valued = TRUE))
    stop("bad input for argument 'nsimEIM'")
  if (nsimEIM <= 50)
    warning("argument 'nsimEIM' should be an integer ",
            "greater than 50, say")

  if (length(ialpha0))
    if (!is.Numeric(ialpha0, positive = TRUE))
      stop("argument 'ialpha' values must be positive")
  if (length(ilambda))
    if (!is.Numeric(ilambda, positive = TRUE))
      stop("argument 'ilambda' values must be positive")


  if (!is.Numeric(imethod, length.arg = 1,
                  integer.valued = TRUE, positive = TRUE) ||
     imethod > 2)
    stop("argument 'imethod' must be 1 or 2")



  new("vglmff",
  blurb = c("Marshall-Olkin exponential distribution\n\n",
            "Links:    ",
            namesof("alpha",  lalpha0, ealpha0 ), ", ",
            namesof("lambda", llambda, elambda ), "\n",
            "Median:     log(3 - alpha) / lambda"),

  constraints = eval(substitute(expression({
    constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
                     M = M, M1 = 2,
                     predictors.names = predictors.names)
  }), list( .zero = zero ))),

  infos = eval(substitute(function(...) {
    list(M1 = 2,
         Q1 = 1,
         dpqrfun = "moe",
         nsimEIM = .nsimEIM,
         parameters.names = c("alpha", "lambda"),
         zero = .zero )
  }, list( .zero = zero,
           .nsimEIM = nsimEIM ))),
  initialize = eval(substitute(expression({

    temp5 <-
    w.y.check(w = w, y = y,
              Is.positive.y = TRUE,
              ncol.w.max = Inf,
              ncol.y.max = Inf,
              out.wy = TRUE,
              colsyperw = 1,
              maximize = TRUE)
    w <- temp5$w
    y <- temp5$y



    ncoly <- ncol(y)

    M1 <- 2
    extra$ncoly <- ncoly
    extra$M1 <- M1
    M <- M1 * ncoly


    mynames1 <- param.names("alpha",  ncoly, skip1 = TRUE)
    mynames2 <- param.names("lambda", ncoly, skip1 = TRUE)
    predictors.names <-
        c(namesof(mynames1, .lalpha0 , .ealpha0 , tag = FALSE),
          namesof(mynames2, .llambda , .elambda , tag = FALSE))[
          interleave.VGAM(M, M1 = M1)]



    if (!length(etastart)) {

      matL <- matrix(if (length( .ilambda )) .ilambda else 0,
                     n, ncoly, byrow = TRUE)
      matA <- matrix(if (length( .ialpha0 )) .ialpha0 else 0,
                     n, ncoly, byrow = TRUE)


      for (spp. in 1:ncoly) {
        yvec <- y[, spp.]

        moexpon.Loglikfun <-
          function(lambdaval, y, x, w, extraargs) {
          ans <-
          sum(c(w) * log(dmoe(x = y, alpha = extraargs$alpha,
                              lambda = lambdaval)))
          ans
        }
        Alpha.init <- .ialpha0
        lambda.grid <- seq(0.1, 10.0, len = 21)
        Lambda.init <- grid.search(lambda.grid,
                                   objfun = moexpon.Loglikfun,
                                   y = y, x = x, w = w,
                         extraargs = list(alpha = Alpha.init))

        if (length(mustart)) {
          Lambda.init <- Lambda.init / (1 - Phimat.init)
        }

        if (!length( .ialpha0 ))
          matA[, spp.] <- Alpha0.init
        if (!length( .ilambda ))
          matL[, spp.] <- Lambda.init
      }  # spp.

      etastart <- cbind(theta2eta(matA, .lalpha0, .ealpha0 ),
                        theta2eta(matL, .llambda, .elambda ))[,
                        interleave.VGAM(M, M1 = M1)]
      mustart <- NULL # Since etastart has been computed.
    }  # End of !length(etastart)
  }), list( .lalpha0 = lalpha0, .llambda = llambda,
            .ealpha0 = ealpha0, .elambda = elambda,
            .ialpha0 = ialpha0, .ilambda = ilambda,
            .imethod = imethod
          ))),

  linkinv = eval(substitute(function(eta, extra = NULL) {
    alpha0 <- eta2theta(eta[, c(TRUE, FALSE)], .lalpha0 , .ealpha0 )
    lambda <- eta2theta(eta[, c(FALSE, TRUE)], .llambda , .elambda )
    log(3 - alpha0) / lambda
  }, list( .lalpha0 = lalpha0, .llambda = llambda,
           .ealpha0 = ealpha0, .elambda = elambda ))),
  last = eval(substitute(expression({
    M1 <- extra$M1
    misc$link <-
      c(rep_len( .lalpha0 , ncoly),
        rep_len( .llambda , ncoly))[interleave.VGAM(M, M1 = M1)]
    temp.names <- c(mynames1, mynames2)[interleave.VGAM(M,
                                                        M1 = M1)]
    names(misc$link) <- temp.names

    misc$earg <- vector("list", M)
    names(misc$earg) <- temp.names
    for (ii in 1:ncoly) {
      misc$earg[[M1*ii-1]] <- .ealpha0
      misc$earg[[M1*ii  ]] <- .elambda
    }

    misc$M1 <- M1
    misc$imethod <- .imethod
    misc$expected <- TRUE
    misc$multipleResponses <- TRUE
    misc$nsimEIM <- .nsimEIM
  }), list( .lalpha0 = lalpha0, .llambda = llambda,
            .ealpha0 = ealpha0, .elambda = elambda,
            .nsimEIM = nsimEIM,
            .imethod = imethod ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    alpha0 <- eta2theta(eta[, c(TRUE, FALSE)], .lalpha0 , .ealpha0 )
    lambda <- eta2theta(eta[, c(FALSE, TRUE)], .llambda , .elambda )
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * log(dmoe(x = y, alpha = alpha0,
                                 lambda = lambda))
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
    }, list( .lalpha0 = lalpha0, .llambda = llambda,
             .ealpha0 = ealpha0, .elambda = elambda ))),
  vfamily = c("exponential.mo"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    alpha0 <- eta2theta(eta[, c(TRUE, FALSE), drop = FALSE],
                        .lalpha0 , .ealpha0 )
    lambda <- eta2theta(eta[, c(FALSE, TRUE), drop = FALSE],
                        .llambda , .elambda )
    okay1 <- all(is.finite(alpha0)) && all(0 < alpha0) &&
             all(is.finite(lambda)) && all(0 < lambda)
    okay1
  }, list( .lalpha0 = lalpha0, .llambda = llambda,
           .ealpha0 = ealpha0, .elambda = elambda ))),


  deriv = eval(substitute(expression({
    M1 <- 2
    alpha0 <- eta2theta(eta[, c(TRUE, FALSE), drop = FALSE],
                        .lalpha0 , .ealpha0 )
    lambda <- eta2theta(eta[, c(FALSE, TRUE), drop = FALSE],
                        .llambda , .elambda )

    temp2 <- (expm1(lambda * y) + alpha0)
    dl.dalpha0 <- -2 / temp2
    dl.dlambda <- 1 / lambda + y -
                  2 * y * exp(lambda * y) / temp2

    dalpha0.deta <- dtheta.deta(alpha0, .lalpha0 , .ealpha0 )
    dlambda.deta <- dtheta.deta(lambda, .llambda , .elambda )

    dthetas.detas <- cbind(dalpha0.deta,
                           dlambda.deta)
    myderiv <- c(w) * cbind(dl.dalpha0, dl.dlambda) *
                      dthetas.detas
    myderiv[, interleave.VGAM(M, M1 = M1)]
  }), list( .lalpha0 = lalpha0, .llambda = llambda,
            .ealpha0 = ealpha0, .elambda = elambda ))),


  weight = eval(substitute(expression({

    NOS <- M / M1
    dThetas.detas <- dthetas.detas[, interleave.VGAM(M, M1 = M1)]

    wz <- matrix(0.0, n, M + M - 1)  # wz is 'tridiagonal'

    ind1 <- iam(NA, NA, M = M1, both = TRUE, diag = TRUE)


    for (spp. in 1:NOS) {
      run.varcov <- 0
      Alph <- alpha0[, spp.]
      Lamb <- lambda[, spp.]

      for (ii in 1:( .nsimEIM )) {
        ysim <- rmoe(n = n, alpha = Alph, lambda = Lamb)
if (ii < 3) {
}

        temp2 <- (expm1(lambda * ysim) + alpha0)
        dl.dalpha0 <- -2 / temp2
        dl.dlambda <- 1 / lambda + ysim -
                      2 * ysim * exp(lambda * ysim) / temp2


        temp3 <- cbind(dl.dalpha0, dl.dlambda)
        run.varcov <- run.varcov +
                      temp3[, ind1$row.index] *
                      temp3[, ind1$col.index]
      }
      run.varcov <- cbind(run.varcov / .nsimEIM)

      wz1 <- if (intercept.only)
          matrix(colMeans(run.varcov),
                 n, ncol(run.varcov), byrow = TRUE) else
          run.varcov

      wz1 <- wz1 * dThetas.detas[, M1 * (spp. - 1) + ind1$row] *
                   dThetas.detas[, M1 * (spp. - 1) + ind1$col]


      for (jay in 1:M1)
        for (kay in jay:M1) {
          cptr <- iam((spp. - 1) * M1 + jay,
                      (spp. - 1) * M1 + kay,
                      M = M)
          wz[, cptr] <- wz1[, iam(jay, kay, M = M1)]
        }
    }  # End of for (spp.) loop




    w.wz.merge(w = w, wz = wz, n = n, M = M, ndepy = M / M1)
  }), list( .llambda = llambda,
            .elambda = elambda,
            .nsimEIM = nsimEIM ))))
}  # exponential.mo()









        genbetaII.Loglikfun4 <-
            function(scaleval, shape1.a, shape2.p, shape3.q,
                     y, x, w, extraargs) {
          sum(c(w) * dgenbetaII(x = y, scale = scaleval,
                                shape1.a = shape1.a,
                                shape2.p = shape2.p,
                                shape3.q = shape3.q,
                                log = TRUE))
        }



 genbetaII <-
  function(lscale    = "loglink",
           lshape1.a = "loglink",
           lshape2.p = "loglink",
           lshape3.q = "loglink",
           iscale    = NULL,
           ishape1.a = NULL,
           ishape2.p = NULL,
           ishape3.q = NULL,
           lss       = TRUE,
           gscale    = exp(-5:5),
           gshape1.a = exp(-5:5),
           gshape2.p = exp(-5:5),
           gshape3.q = exp(-5:5),
           zero      = "shape") {





  if (length(lss) != 1 && !is.logical(lss))
    stop("Argument 'lss' not specified correctly")


  if (length(iscale) && !is.Numeric(iscale, positive = TRUE))
    stop("Bad input for argument 'iscale'")

  if (length(ishape1.a) && !is.Numeric(ishape1.a, positive = TRUE))
    stop("Bad input for argument 'ishape1.a'")

  if (length(ishape2.p) && !is.Numeric(ishape2.p, positive = TRUE))
    stop("Bad input for argument 'ishape2.p'")

  if (length(ishape3.q) && !is.Numeric(ishape3.q, positive = TRUE))
    stop("Bad input for argument 'ishape3.q'")



  lscale <- as.list(substitute(lscale))
  escale <- link2list(lscale)
  lscale <- attr(escale, "function.name")

  lshape1.a <- as.list(substitute(lshape1.a))
  eshape1.a <- link2list(lshape1.a)
  lshape1.a <- attr(eshape1.a, "function.name")

  lshape2.p <- as.list(substitute(lshape2.p))
  eshape2.p <- link2list(lshape2.p)
  lshape2.p <- attr(eshape2.p, "function.name")

  lshape3.q <- as.list(substitute(lshape3.q))
  eshape3.q <- link2list(lshape3.q)
  lshape3.q <- attr(eshape3.q, "function.name")


  new("vglmff",
  blurb =
    c("Generalized Beta II distribution \n\n",
      "Links:    ",
      ifelse (lss,
              namesof("scale"   , lscale   , escale),
              namesof("shape1.a", lshape1.a, eshape1.a)), ", ",
      ifelse (lss,
              namesof("shape1.a", lshape1.a, eshape1.a),
              namesof("scale"   , lscale   , escale)), ", ",
      namesof("shape2.p" , lshape2.p, earg = eshape2.p), ", ",
      namesof("shape3.q" , lshape3.q, earg = eshape3.q), "\n",
      "Mean:     scale * gamma(shape2.p + 1/shape1.a) * ",
                "gamma(shape3.q - 1/shape1.a) / ",
                "(gamma(shape2.p) * gamma(shape3.q))"),
  constraints = eval(substitute(expression({
    constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
                     M = M, M1 = 4,
                     predictors.names = predictors.names)
  }), list( .zero = zero ))),
  infos = eval(substitute(function(...) {
    list(M1 = 4,
         Q1 = 1,
         dpqrfun = "genbetaII",
         expected = TRUE,
         zero = .zero ,
         multipleResponses = TRUE,
         parameters.names = if ( .lss )
           c("scale", "shape1.a", "shape2.p", "shape3.q") else
           c("shape1.a", "scale", "shape2.p", "shape3.q"),
         lscale    = .lscale    , lshape1.a = .lshape1.a ,
         escale    = .escale    , eshape1.a = .eshape1.a ,
         lshape2.p = .lshape2.p , lshape3.q = .lshape3.q ,
         eshape2.p = .eshape2.p , eshape3.q = .eshape3.q )
  }, list( .lscale = lscale      , .lshape1.a = lshape1.a,
           .escale = escale      , .eshape1.a = eshape1.a,
           .lshape2.p = lshape2.p, .lshape3.q = lshape3.q,
           .eshape2.p = eshape2.p, .eshape3.q = eshape3.q,
           .lss  = lss ,
           .zero = zero ))),
  initialize = eval(substitute(expression({
    temp5 <- w.y.check(w = w, y = y,
                       Is.positive.y = TRUE,
                       ncol.w.max = Inf,
                       ncol.y.max = Inf,
                       out.wy = TRUE,
                       colsyperw = 1,
                       maximize = TRUE)
    y    <- temp5$y
    w    <- temp5$w
    M1 <- 4  # Number of parameters for one response
    NOS  <- ncoly <- ncol(y)
    M    <- M1*ncol(y)

    scaL.names  <- param.names("scale",     NOS, skip1 = TRUE)
    sha1.names  <- param.names("shape1.a",  NOS, skip1 = TRUE)
    sha2.names  <- param.names("shape2.p",  NOS, skip1 = TRUE)
    sha3.names  <- param.names("shape3.q",  NOS, skip1 = TRUE)

    predictors.names <- c(
      if ( .lss ) {
        c(namesof(scaL.names , .lscale    , earg = .escale    ,
                  tag = FALSE),
          namesof(sha1.names , .lshape1.a , earg = .eshape1.a ,
                  tag = FALSE))
      } else {
        c(namesof(sha1.names , .lshape1.a , earg = .eshape1.a ,
                  tag = FALSE),
          namesof(scaL.names , .lscale    , earg = .escale    ,
                  tag = FALSE))
      },
      namesof(sha2.names , .lshape2.p , .eshape2.p , tag = FALSE),
      namesof(sha3.names , .lshape3.q , .eshape3.q , tag = FALSE))
    predictors.names <- predictors.names[interleave.VGAM(M,
                                                         M1 = M1)]

    if (!length(etastart)) {
      sc.init <-
      aa.init <-
      pp.init <-
      qq.init <- matrix(NA_real_, n, NOS)

      for (spp. in 1:NOS) {  # For each response 'y_spp.'... do:
        yvec <- y[, spp.]
        wvec <- w[, spp.]

          gscale     <- .gscale
          gshape1.a  <- .gshape1.a
          gshape2.p  <- .gshape2.p
          gshape3.q  <- .gshape3.q
        if (length( .iscale    )) gscale    <-
                                      rep_len( .iscale    , NOS)
        if (length( .ishape1.a )) gshape1.a <-
                                      rep_len( .ishape1.a , NOS)
        if (length( .ishape2.p )) gshape2.p <-
                                      rep_len( .ishape2.p , NOS)
        if (length( .ishape3.q )) gshape3.q <-
                                      rep_len( .ishape3.q , NOS)


        try.this <-
          grid.search4(gscale, gshape1.a, gshape2.p, gshape3.q,
                       objfun = genbetaII.Loglikfun4,  # x = x,
                       y = yvec, w = wvec,  # extraargs = NULL,
                       ret.objfun = TRUE)  # Last value is \ell
          sc.init[, spp.] <- try.this["Value1"]
          aa.init[, spp.] <- try.this["Value2"]
          pp.init[, spp.] <- try.this["Value3"]
          qq.init[, spp.] <- try.this["Value4"]
      }  # End of for (spp. ...)


      finite.mean <- 1 < aa.init * qq.init
      COP.use <- 1.15
      while (FALSE && any(!finite.mean)) {
    qq.init[!finite.mean] <- 0.1 + qq.init[!finite.mean] * COP.use
    aa.init[!finite.mean] <- 0.1 + aa.init[!finite.mean] * COP.use
    finite.mean <- 1 < aa.init * qq.init
      }

      etastart <-
        cbind(if ( .lss )
          cbind(theta2eta(sc.init, .lscale    , .escale    ),
                theta2eta(aa.init, .lshape1.a , .eshape1.a )) else
          cbind(theta2eta(aa.init, .lshape1.a , .eshape1.a ),
                theta2eta(sc.init, .lscale    , .escale    )),
          theta2eta(pp.init , .lshape2.p , earg = .eshape2.p ),
          theta2eta(qq.init , .lshape3.q , earg = .eshape3.q ))
      etastart <- etastart[, interleave.VGAM(M, M1 = M1)]
    }  # End of etastart.
  }), list( .lscale    = lscale   , .lshape1.a = lshape1.a,
            .escale    = escale   , .eshape1.a = eshape1.a,
            .iscale    = iscale   , .ishape1.a = ishape1.a,
            .gscale    = gscale   , .gshape1.a = gshape1.a,
            .lshape2.p = lshape2.p, .lshape3.q = lshape3.q,
            .eshape2.p = eshape2.p, .eshape3.q = eshape3.q,
            .ishape2.p = ishape2.p, .ishape3.q = ishape3.q,
            .gshape2.p = gshape2.p, .gshape3.q = gshape3.q,
            .lss = lss ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    M1 <- 4
    NOS <- ncol(eta) / M1
    if ( .lss ) {
      Scale <- eta2theta(eta[, M1*(1:NOS) - 3, drop = FALSE],
                         .lscale ,  earg = .escale )
      aa    <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a )
    } else {
      aa    <- eta2theta(eta[, M1*(1:NOS) - 3, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a )
      Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                         .lscale    ,  earg = .escale )
    }
    parg <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                      .lshape2.p , earg = .eshape2.p )
    qq   <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                      .lshape3.q , earg = .eshape3.q )

    ans <- cbind(Scale * exp(lgamma(parg + 1/aa) +
                             lgamma(qq   - 1/aa) -
                             lgamma(parg) - lgamma(qq)))
    ans
  }, list( .lscale    = lscale   , .lshape1.a = lshape1.a,
           .escale    = escale   , .eshape1.a = eshape1.a,
           .lshape2.p = lshape2.p, .lshape3.q = lshape3.q,
           .eshape2.p = eshape2.p, .eshape3.q = eshape3.q,
           .lss = lss ))),
  last = eval(substitute(expression({
    M1 <- 4

    misc$link <- c(rep_len( if ( .lss ) .lscale else .lshape1.a ,
                           ncoly),
                   rep_len( if ( .lss ) .lshape1.a else .lscale ,
                           ncoly),
                   rep_len( .lshape2.p , ncoly),
                   rep_len( .lshape3.q , ncoly))[
                   interleave.VGAM(M, M1 = M1)]
    temp.names <- if ( .lss ) {
      c(scaL.names, sha1.names, sha2.names, sha3.names)
    } else {
      c(sha1.names, scaL.names, sha2.names, sha3.names)
    }
    names(misc$link) <- temp.names[interleave.VGAM(M, M1 = M1)]

    misc$earg <- vector("list", M)
    names(misc$earg) <- temp.names[interleave.VGAM(M, M1 = M1)]
    for (ii in 1:ncoly) {
      if ( .lss ) {
        misc$earg[[M1*ii-3]] <- .escale
        misc$earg[[M1*ii-2]] <- .eshape1.a
    } else {
        misc$earg[[M1*ii-3]] <- .eshape1.a
        misc$earg[[M1*ii-2]] <- .escale
    }
      misc$earg[[M1*ii-1]] <- .eshape2.p
      misc$earg[[M1*ii  ]] <- .eshape3.q
    }

    misc$expected          <- TRUE
    misc$multipleResponses <- TRUE
  }), list( .lscale    = lscale   , .lshape1.a = lshape1.a,
            .escale    = escale   , .eshape1.a = eshape1.a,
            .lshape2.p = lshape2.p, .lshape3.q = lshape3.q,
            .eshape2.p = eshape2.p, .eshape3.q = eshape3.q,
            .lss = lss ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE,
             eta, extra = NULL, summation = TRUE) {
      M1  <- 4
      NOS <- ncol(eta)/M1
      if ( .lss ) {
        Scale <- eta2theta(eta[, M1*(1:NOS) - 3, drop = FALSE],
                           .lscale    , earg = .escale    )
        aa    <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                           .lshape1.a , earg = .eshape1.a )
      } else {
        aa    <- eta2theta(eta[, M1*(1:NOS) - 3, drop = FALSE],
                           .lshape1.a , earg = .eshape1.a )
        Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                           .lscale    , earg = .escale    )
      }
      parg   <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                          .lshape2.p , earg = .eshape2.p )
      qq     <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                          .lshape3.q , earg = .eshape3.q )
      if (residuals) {
        stop("loglikelihood residuals not implemented yet")
      } else {
        ll.elts <-
          c(w) * dgenbetaII(x = y, scale = Scale, shape1.a = aa,
                            shape2.p = parg, shape3.q = qq,
                            log = TRUE)
        if (summation) {
          sum(ll.elts)
        } else {
          ll.elts
        }
      }
    }, list( .lscale    = lscale   , .lshape1.a = lshape1.a,
             .escale    = escale   , .eshape1.a = eshape1.a,
             .lshape2.p = lshape2.p, .lshape3.q = lshape3.q,
             .eshape2.p = eshape2.p, .eshape3.q = eshape3.q,
             .lss = lss ))),
  vfamily = c("genbetaII"),




  validparams = eval(substitute(function(eta, y, extra = NULL) {
    M1 <- 4
    NOS <- ncol(eta) / M1
    if ( .lss ) {
      Scale <- eta2theta(eta[, M1*(1:NOS) - 3, drop = FALSE],
                         .lscale    , earg = .escale   )
      aa    <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a)
    } else {
      aa    <- eta2theta(eta[, M1*(1:NOS) - 3, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a)
      Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                         .lscale    , earg = .escale )
    }
    parg <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                      .lshape2.p , earg = .eshape2.p)
    qq   <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                      .lshape3.q , earg = .eshape3.q)


    okay1 <- all(is.finite(Scale)) && all(Scale  > 0) &&
             all(is.finite(aa   )) && all(aa     > 0) &&
             all(is.finite(parg )) && all(parg   > 0) &&
             all(is.finite(qq   )) && all(qq     > 0)
    okay.support <- if (okay1)
                    all(-aa * parg < 1 & 1 < aa * qq) else TRUE
    if (!okay.support)
      warning("parameter constraint -a*p < 1 < a*q has  ",
              "been violated; solution may be at the",
              " boundary of the parameter space.")
    okay1 && okay.support
  }, list(  .lscale    = lscale   , .lshape1.a = lshape1.a,
            .escale    = escale   , .eshape1.a = eshape1.a,
            .lshape2.p = lshape2.p, .lshape3.q = lshape3.q,
            .eshape2.p = eshape2.p, .eshape3.q = eshape3.q,
            .lss = lss ))),




  deriv = eval(substitute(expression({
    M1 <- 4
    NOS <- ncol(eta)/M1  # Needed for summary()
    if ( .lss ) {
      Scale <- eta2theta(eta[, M1*(1:NOS) - 3, drop = FALSE],
                         .lscale    , earg = .escale   )
      aa    <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a)
    } else {
      aa    <- eta2theta(eta[, M1*(1:NOS) - 3, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a)
      Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                         .lscale    , earg = .escale )
    }
    parg <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                      .lshape2.p , earg = .eshape2.p)
    qq   <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                      .lshape3.q , earg = .eshape3.q)
    temp1 <- log(y/Scale)
    temp2 <- (y/Scale)^aa
    temp3 <- digamma(parg + qq)
    temp3a <- digamma(parg)
    temp3b <- digamma(qq)
    temp4 <- log1p(temp2)

    dl.dscale <- (aa/Scale) * (-parg + (parg+qq) / (1+1/temp2))
    dl.da <- 1/aa + parg * temp1 - (parg+qq) * temp1 / (1+1/temp2)
    dl.dp <- aa * temp1 + temp3 - temp3a - temp4
    dl.dq <- temp3 - temp3b - temp4

    dscale.deta <- dtheta.deta(Scale, .lscale ,    .escale )
    da.deta     <- dtheta.deta(aa,    .lshape1.a , .eshape1.a )
    dp.deta     <- dtheta.deta(parg,  .lshape2.p , .eshape2.p )
    dq.deta     <- dtheta.deta(qq,    .lshape3.q , .eshape3.q )

    myderiv <- if ( .lss ) {
      c(w) * cbind(dl.dscale * dscale.deta,
                   dl.da * da.deta,
                   dl.dp * dp.deta,
                   dl.dq * dq.deta)
    } else {
      c(w) * cbind(dl.da * da.deta,
                   dl.dscale * dscale.deta,
                   dl.dp * dp.deta,
                   dl.dq * dq.deta)
    }
    myderiv[, interleave.VGAM(M, M1 = M1)]
  }), list(  .lscale    = lscale   , .lshape1.a = lshape1.a,
             .escale    = escale   , .eshape1.a = eshape1.a,
             .lshape2.p = lshape2.p, .lshape3.q = lshape3.q,
             .eshape2.p = eshape2.p, .eshape3.q = eshape3.q,
             .lss = lss ))),
  weight = eval(substitute(expression({
    temp5  <- trigamma(parg + qq)
    temp5a <- trigamma(parg)
    temp5b <- trigamma(qq)

    ned2l.da <- (1 + parg + qq + parg * qq * (temp5a + temp5b +
                (temp3b - temp3a + (parg-qq)/(parg*qq))^2 -
                (parg^2 + qq^2) / (parg*qq)^2)) / (aa^2 * (1+parg+qq))
    ned2l.dscale  <- (aa^2) * parg * qq / ((1+parg+qq) * Scale^2)
    ned2l.dp      <- temp5a - temp5
    ned2l.dq      <- temp5b - temp5
    ned2l.dascale <- (parg - qq - parg * qq *
                     (temp3a -temp3b)) / (Scale*(1 + parg+qq))
    ned2l.dap     <- -(qq   * (temp3a -temp3b) -1) / (aa*(parg+qq))
    ned2l.daq     <- -(parg * (temp3b -temp3a) -1) / (aa*(parg+qq))
    ned2l.dscalep <-  aa * qq   / (Scale*(parg+qq))
    ned2l.dscaleq <- -aa * parg / (Scale*(parg+qq))
    ned2l.dpq     <- -temp5
    wz <- if ( .lss ) {
      array(c(c(w) * ned2l.dscale * dscale.deta^2,
              c(w) * ned2l.da * da.deta^2,
              c(w) * ned2l.dp * dp.deta^2,
              c(w) * ned2l.dq * dq.deta^2,
              c(w) * ned2l.dascale * da.deta * dscale.deta,
              c(w) * ned2l.dap * da.deta * dp.deta,
              c(w) * ned2l.dpq * dp.deta * dq.deta,
              c(w) * ned2l.dscalep * dscale.deta * dp.deta,
              c(w) * ned2l.daq * da.deta * dq.deta,
              c(w) * ned2l.dscaleq * dscale.deta * dq.deta),
            dim = c(n, M/M1, M1*(M1+1)/2))
    } else {
      array(c(c(w) * ned2l.da * da.deta^2,
              c(w) * ned2l.dscale * dscale.deta^2,
              c(w) * ned2l.dp * dp.deta^2,
              c(w) * ned2l.dq * dq.deta^2,
              c(w) * ned2l.dascale * da.deta * dscale.deta,
              c(w) * ned2l.dscalep * dscale.deta * dp.deta,
              c(w) * ned2l.dpq * dp.deta * dq.deta,
              c(w) * ned2l.dap * da.deta * dp.deta,
              c(w) * ned2l.dscaleq * dscale.deta * dq.deta,
              c(w) * ned2l.daq * da.deta * dq.deta),
            dim = c(n, M/M1, M1*(M1+1)/2))
    }
    wz <- arwz2wz(wz, M = M, M1 = M1)
    wz
  }), list( .lscale    = lscale   , .lshape1.a = lshape1.a,
            .escale    = escale   , .eshape1.a = eshape1.a,
            .lshape2.p = lshape2.p, .lshape3.q = lshape3.q,
            .eshape2.p = eshape2.p, .eshape3.q = eshape3.q,
            .lss = lss ))))
}  # genbetaII










dgenbetaII <-
  function(x, scale = 1, shape1.a, shape2.p, shape3.q,
           log = FALSE)  {

  if (!is.logical(log.arg <- log) || length(log) != 1)
    stop("Bad input for argument 'log'")
  rm(log)


    logden <- log(shape1.a) +
      (shape1.a * shape2.p - 1) * log(abs(x)) -
       shape1.a * shape2.p * log(scale) -
       lbeta(shape2.p, shape3.q) -
       (shape2.p + shape3.q) * log1p((abs(x)/scale)^shape1.a)


  if (any(x <= 0) || any(is.infinite(x))) {
    LLL <- max(length(x),        length(scale),
               length(shape1.a), length(shape2.p), length(shape3.q))
    if (length(x)        != LLL) x        <- rep_len(x,        LLL)
    if (length(scale)    != LLL) scale    <- rep_len(scale,    LLL)
    if (length(shape1.a) != LLL) shape1.a <- rep_len(shape1.a, LLL)
    if (length(shape2.p) != LLL) shape2.p <- rep_len(shape2.p, LLL)
    if (length(shape3.q) != LLL) shape3.q <- rep_len(shape3.q, LLL)

    logden[is.infinite(x)] <- log(0)
    logden[x < 0] <- log(0)
    x.eq.0 <- !is.na(x) & (x == 0)
    if (any(x.eq.0)) {
      axp <- shape1.a[x.eq.0] * shape2.p[x.eq.0]
      logden[x.eq.0 & axp <  1] <- log(Inf)
      ind5 <- x.eq.0 & axp == 1
      logden[ind5] <- log(shape1.a[ind5]) -
        shape1.a[ind5] * shape2.p[ind5] * log(scale[ind5]) -
        lbeta(shape2.p[ind5], shape3.q[ind5]) -
        (shape2.p[ind5] + shape3.q[ind5]) *
        log1p((0/scale[ind5])^shape1.a[ind5])
      logden[x.eq.0 & axp >  1] <- log(0)
    }
  }

  if (log.arg) logden else exp(logden)
}  # dgenbetaII





rsinmad <- function(n, scale = 1, shape1.a, shape3.q)
  qsinmad(runif(n), shape1.a = shape1.a, scale = scale,
          shape3.q = shape3.q)


rlomax <- function(n, scale = 1, shape3.q)
  rsinmad(n, scale = scale, shape1.a = 1, shape3.q = shape3.q)


rfisk <- function(n, scale = 1, shape1.a)
  rsinmad(n, scale = scale, shape1.a = shape1.a, shape3.q = 1)


rparalogistic <- function(n, scale = 1, shape1.a)
  rsinmad(n, scale = scale, shape1.a = shape1.a,
          shape3.q = shape1.a)


rdagum <- function(n, scale = 1, shape1.a, shape2.p)
  qdagum(runif(n), scale = scale, shape1.a = shape1.a,
         shape2.p = shape2.p)


rinv.lomax <- function(n, scale = 1, shape2.p)
  rdagum(n, scale = scale, shape1.a = 1, shape2.p = shape2.p)


rinv.paralogistic <- function(n, scale = 1, shape1.a)
  rdagum(n, scale = scale, shape1.a = shape1.a,
         shape2.p = shape1.a)




qsinmad <- function(p, scale = 1, shape1.a, shape3.q,
                    lower.tail = TRUE,
                    log.p = FALSE) {



  if (!is.logical(lower.tail) || length(lower.tail ) != 1)
    stop("bad input for argument 'lower.tail'")

  if (!is.logical(log.p) || length(log.p) != 1)
    stop("bad input for argument 'log.p'")



  LLL <- max(length(p), length(shape1.a), length(scale),
                        length(shape3.q))
  if (length(p)        != LLL) p         <- rep_len(p,         LLL)
  if (length(shape1.a) != LLL) shape1.a  <- rep_len(shape1.a,  LLL)
  if (length(scale)    != LLL) scale     <- rep_len(scale,     LLL)
  if (length(shape3.q) != LLL) shape3.q  <- rep_len(shape3.q,  LLL)


  if (lower.tail) {
    if (log.p) {
      ln.p <- p
      ans <- scale * expm1((-1/shape3.q) *
                           log(-expm1(ln.p)))^(1/shape1.a)
    } else {
      ans <- scale * expm1((-1/shape3.q) * log1p(-p))^(1/shape1.a)
      ans[p == 0] <- 0
      ans[p == 1] <- Inf
    }
  } else {
    if (log.p) {
      ln.p <- p
      ans <- scale * expm1(-ln.p / shape3.q)^(1/shape1.a)
    } else {
      ans <- scale * expm1(-log(p) / shape3.q)^(1/shape1.a)
      ans[p == 0] <- Inf
      ans[p == 1] <- 0
    }
  }

  ans[scale    <= 0 | shape1.a <= 0 | shape3.q <= 0] <- NaN
  ans
}  # qsinmad




qlomax <- function(p, scale = 1, shape3.q,
                   lower.tail = TRUE, log.p = FALSE)
  qsinmad(p, shape1.a = 1, scale = scale, shape3.q = shape3.q,
          lower.tail = lower.tail, log.p = log.p)

qfisk <- function(p, scale = 1, shape1.a,
                  lower.tail = TRUE, log.p = FALSE)
  qsinmad(p, shape1.a = shape1.a, scale = scale, shape3.q = 1,
          lower.tail = lower.tail, log.p = log.p)

qparalogistic <- function(p, scale = 1, shape1.a,
                          lower.tail = TRUE, log.p = FALSE)
  qsinmad(p, shape1.a = shape1.a, scale = scale,
          shape3.q = shape1.a,
          lower.tail = lower.tail, log.p = log.p)





qdagum <- function(p, scale = 1, shape1.a, shape2.p,
                   lower.tail = TRUE, log.p = FALSE) {

  if (!is.logical(lower.tail) || length(lower.tail ) != 1)
    stop("bad input for argument 'lower.tail'")
  if (!is.logical(log.p) || length(log.p) != 1)
    stop("bad input for argument 'log.p'")


  LLL <- max(length(p), length(shape1.a), length(scale),
                        length(shape2.p))
  if (length(p)        != LLL) p         <- rep_len(p,         LLL)
  if (length(shape1.a) != LLL) shape1.a  <- rep_len(shape1.a,  LLL)
  if (length(scale)    != LLL) scale     <- rep_len(scale,     LLL)
  if (length(shape2.p) != LLL) shape2.p  <- rep_len(shape2.p,  LLL)

  if (lower.tail) {
    if (log.p) {
      ln.p <- p
      ans <- scale * (expm1(-ln.p/shape2.p))^(-1/shape1.a)
      ans[ln.p > 0] <- NaN
    } else {
      ans <- scale * (expm1(-log(p)/shape2.p))^(-1/shape1.a)
      ans[p < 0] <- NaN
      ans[p == 0] <- 0
      ans[p == 1] <- Inf
      ans[p > 1] <- NaN
    }
  } else {
    if (log.p) {
      ln.p <- p
      ans <- scale *
        (expm1(-log(-expm1(ln.p))/shape2.p))^(-1/shape1.a)
      ans[ln.p > 0] <- NaN
    } else {
      ans <- scale * (expm1(-log1p(-p)/shape2.p))^(-1/shape1.a)
      ans[p < 0] <- NaN
      ans[p == 0] <- Inf
      ans[p == 1] <- 0
      ans[p > 1] <- NaN
    }
  }

  ans[scale <= 0 | shape1.a <= 0 | shape2.p <= 0] <- NaN
  ans
}  # qdagum




qinv.lomax <- function(p, scale = 1, shape2.p,
                       lower.tail = TRUE, log.p = FALSE)
  qdagum(p, scale = scale, shape1.a = 1, shape2.p = shape2.p,
         lower.tail = lower.tail, log.p = log.p)


qinv.paralogistic <- function(p, scale = 1, shape1.a,
                              lower.tail = TRUE, log.p = FALSE)
  qdagum(p, scale = scale, shape1.a = shape1.a,
  shape2.p = shape1.a,   ##  20150121 Kai; add shape2.p=shape1.a
         lower.tail = lower.tail, log.p = log.p)




psinmad <- function(q, scale = 1, shape1.a, shape3.q,
                    lower.tail = TRUE, log.p = FALSE) {


  if (!is.logical(lower.tail) || length(lower.tail ) != 1)
    stop("bad input for argument 'lower.tail'")

  if (!is.logical(log.p) || length(log.p) != 1)
    stop("bad input for argument 'log.p'")


  LLL <- max(length(q), length(shape1.a), length(scale),
                        length(shape3.q))
  if (length(q)        != LLL) q         <- rep_len(q,         LLL)
  if (length(shape1.a) != LLL) shape1.a  <- rep_len(shape1.a,  LLL)
  if (length(scale)    != LLL) scale     <- rep_len(scale,     LLL)
  if (length(shape3.q) != LLL) shape3.q  <- rep_len(shape3.q,  LLL)

  # 20150121 KaiH
  if (lower.tail) {
    if (log.p) {
      ans <- log1p(-(1 + (q / scale)^shape1.a)^(-shape3.q))
      ans[q <= 0 ] <- -Inf
      ans[q == Inf] <- 0
    } else {
      ans <- exp(log1p(-(1 + (q / scale)^shape1.a)^(-shape3.q)))
      ans[q <= 0] <- 0
      ans[q == Inf] <- 1
    }
  } else {
    if (log.p) {
      ans <- (-shape3.q) * log1p((q / scale)^shape1.a)
      ans[q <= 0] <- 0
      ans[q == Inf] <- -Inf
    } else {
      ans <- (1 + (q / scale)^shape1.a)^(-shape3.q)
      ans[q <= 0] <- 1
      ans[q == Inf] <- 0
    }
  }
  ans[scale <= 0 | shape1.a <= 0 | shape3.q <= 0] <- NaN
  ans
}  # psinmad






plomax <-
  function(q, scale = 1, shape3.q,   # Change the order
           lower.tail = TRUE, log.p = FALSE)
  psinmad(q, shape1.a = 1, scale = scale, shape3.q = shape3.q,
          lower.tail = lower.tail, log.p = log.p)


pfisk <- function(q, scale = 1, shape1.a,
                  lower.tail = TRUE, log.p = FALSE)
  psinmad(q, shape1.a = shape1.a, scale = scale, shape3.q = 1,
          lower.tail = lower.tail, log.p = log.p)


pparalogistic <-
  function(q, scale = 1, shape1.a,   # Change the order
           lower.tail = TRUE, log.p = FALSE)
  psinmad(q, shape1.a = shape1.a, scale = scale,
          shape3.q = shape1.a,  # Add shape3.q = shape1.a
          lower.tail = lower.tail, log.p = log.p)






pdagum <- function(q, scale = 1, shape1.a, shape2.p,
                   lower.tail = TRUE, log.p = FALSE) {



  if (!is.logical(lower.tail) || length(lower.tail ) != 1)
    stop("bad input for argument 'lower.tail'")

  if (!is.logical(log.p) || length(log.p) != 1)
    stop("bad input for argument 'log.p'")


  LLL <- max(length(q), length(shape1.a), length(scale),
                        length(shape2.p))
  if (length(q)        != LLL) q         <- rep_len(q,         LLL)
  if (length(shape1.a) != LLL) shape1.a  <- rep_len(shape1.a,  LLL)
  if (length(scale)    != LLL) scale     <- rep_len(scale,     LLL)
  if (length(shape2.p) != LLL) shape2.p  <- rep_len(shape2.p,  LLL)


  if (lower.tail) {
    if (log.p) {
      ans <- (-shape2.p) * log1p((q/scale)^(-shape1.a))
      ans[q <= 0 ] <- -Inf
      ans[q == Inf] <- 0
    } else {
      ans <- exp( (-shape2.p) * log1p((q/scale)^(-shape1.a)) )
      ans[q <= 0] <- 0
      ans[q == Inf] <- 1
    }
  } else {
    if (log.p) {
      ans <- log1p(-(1 + (q/scale)^(-shape1.a))^(-shape2.p))
      ans[q <= 0] <- 0
      ans[q == Inf] <- -Inf
    } else {
      stop("unfinished")
      ans[q <= 0] <- 1
      ans[q == Inf] <- 0
    }
  }

  ans[shape1.a <= 0 | scale <= 0 | shape2.p <= 0] <- NaN
  ans
}  # pdagum





pinv.lomax <- function(q, scale = 1, shape2.p,
                       lower.tail = TRUE, log.p = FALSE)
  pdagum(q, scale = scale, shape1.a = 1, shape2.p = shape2.p,
         lower.tail = lower.tail, log.p = log.p)


pinv.paralogistic <-
  function(q, scale = 1, shape1.a,
           lower.tail = TRUE, log.p = FALSE)
  pdagum(q, scale = scale, shape1.a = shape1.a,
         shape2.p = shape1.a,
         lower.tail = lower.tail, log.p = log.p)





dbetaII <-
  function(x, scale = 1, shape2.p, shape3.q, log = FALSE)
  dgenbetaII(x = x, scale = scale, shape1.a = 1,
             shape2.p = shape2.p, shape3.q = shape3.q, log = log)




dsinmad <-
  function(x, scale = 1, shape1.a, shape3.q, log = FALSE) {

  if (!is.logical(log.arg <- log) || length(log) != 1)
    stop("bad input for argument 'log'")
  rm(log)

  LLL      <- max(length(x),     length(shape1.a),
                  length(scale), length(shape3.q))
  if (length(x)        != LLL) x         <- rep_len(x,         LLL)
  if (length(shape1.a) != LLL) shape1.a  <- rep_len(shape1.a,  LLL)
  if (length(scale)    != LLL) scale     <- rep_len(scale,     LLL)
  if (length(shape3.q) != LLL) shape3.q  <- rep_len(shape3.q,  LLL)

  Loglik <- rep_len(log(0), LLL)
  xok <- (x > 0) & !is.na(x)  # Avoids log(x) if x<0, and handles NAs
  Loglik[xok] <- log(shape1.a[xok]) + log(shape3.q[xok]) +
                 (shape1.a[xok]-1) * log(x[xok]) -
                shape1.a[xok] * log(scale[xok]) -
      (1 + shape3.q[xok]) *
      log1p((x[xok]/scale[xok])^shape1.a[xok])
  x.eq.0 <- (x == 0) & !is.na(x)
  Loglik[x.eq.0] <- log(shape1.a[x.eq.0]) + log(shape3.q[x.eq.0]) -
                    shape1.a[x.eq.0] * log(scale[x.eq.0])
  Loglik[is.na(x)]  <- NA
  Loglik[is.nan(x)] <- NaN
  Loglik[x == Inf]  <- log(0)

  if (log.arg) Loglik else exp(Loglik)
}



dlomax <- function(x, scale = 1, shape3.q, log = FALSE)
  dsinmad(x, scale = scale, shape1.a = 1, shape3.q = shape3.q,
          log = log)



dfisk <- function(x, scale = 1, shape1.a, log = FALSE)
  dsinmad(x, scale = scale, shape1.a = shape1.a, shape3.q = 1,
          log = log)



dparalogistic <- function(x, scale = 1, shape1.a, log = FALSE)
  dsinmad(x, scale = scale, shape1.a = shape1.a, shape3.q = shape1.a,
          log = log)



ddagum <-
  function(x, scale = 1, shape1.a, shape2.p, log = FALSE) {
  if (!is.logical(log.arg <- log) || length(log) != 1)
    stop("bad input for argument 'log'")
  rm(log)

  LLL <- max(length(x),
             length(shape1.a),
             length(scale),
             length(shape2.p))
  if (length(x)        != LLL) x         <- rep_len(x,         LLL)
  if (length(shape1.a) != LLL) shape1.a  <- rep_len(shape1.a,  LLL)
  if (length(scale)    != LLL) scale     <- rep_len(scale,     LLL)
  if (length(shape2.p) != LLL) shape2.p  <- rep_len(shape2.p,  LLL)

  Loglik <- rep_len(log(0), LLL)
  xok <- (x > 0) & !is.na(x)  # Avoids log(x) if x<0, and handles NAs
  Loglik[xok] <- log(shape1.a[xok]) +
                 log(shape2.p[xok]) +
                 (shape1.a[xok] * shape2.p[xok]-1) * log(    x[xok]) -
                  shape1.a[xok] * shape2.p[xok]    * log(scale[xok]) -
                 (1 + shape2.p[xok]) *
                 log1p((x[xok]/scale[xok])^shape1.a[xok])
  Loglik[shape2.p <= 0] <- NaN

  x.eq.0 <- (x == 0) & !is.na(x)
  Loglik[x.eq.0] <- log(shape1.a[x.eq.0]) +
                    log(shape2.p[x.eq.0]) -
                    shape1.a[x.eq.0] * shape2.p[x.eq.0] *
                    log(scale[x.eq.0])
  Loglik[is.na(x)]  <- NA
  Loglik[is.nan(x)] <- NaN
  Loglik[x == Inf]  <- log(0)

  if (log.arg) Loglik else exp(Loglik)
}



dinv.lomax <- function(x, scale = 1, shape2.p, log = FALSE)
  ddagum(x, scale = scale, shape1.a = 1, shape2.p = shape2.p,
         log = log)


dinv.paralogistic <- function(x, scale = 1, shape1.a, log = FALSE)
  ddagum(x, scale = scale, shape1.a = shape1.a, shape2.p = shape1.a,
         log = log)








 sinmad    <- function(lscale    = "loglink",
                       lshape1.a = "loglink",
                       lshape3.q = "loglink",
                       iscale    = NULL,
                       ishape1.a = NULL,
                       ishape3.q = NULL,
                       imethod   = 1,
                       lss       = TRUE,
                       gscale    = exp(-5:5),
                       gshape1.a = exp(-5:5),
                       gshape3.q = exp(-5:5),
                       probs.y   = c(0.25, 0.50, 0.75),
                       zero      = "shape") {




  if (length(lss) != 1 && !is.logical(lss))
    stop("Argument 'lss' not specified correctly")

  if (!is.Numeric(imethod, length.arg = 1,
                  integer.valued = TRUE,
                  positive = TRUE) || imethod > 2)
    stop("Bad input for argument 'imethod'")

  if (length(iscale) && !is.Numeric(iscale, positive = TRUE))
    stop("Bad input for argument 'iscale'")

  if (length(ishape1.a) && !is.Numeric(ishape1.a, positive = TRUE))
    stop("Bad input for argument 'ishape1.a'")

  if (length(ishape3.q) && !is.Numeric(ishape3.q, positive = TRUE))
    stop("Bad input for argument 'ishape3.q'")

  if (length(probs.y) < 2 || max(probs.y) > 1 ||
        !is.Numeric(probs.y, positive = TRUE))
    stop("Bad input for argument 'probs.y'")

  lscale <- as.list(substitute(lscale))
  escale <- link2list(lscale)
  lscale <- attr(escale, "function.name")

  lshape1.a <- as.list(substitute(lshape1.a))
  eshape1.a <- link2list(lshape1.a)
  lshape1.a <- attr(eshape1.a, "function.name")

  lshape3.q <- as.list(substitute(lshape3.q))
  eshape3.q <- link2list(lshape3.q)
  lshape3.q <- attr(eshape3.q, "function.name")


  new("vglmff",
  blurb =
    c("Singh-Maddala distribution \n\n",
      "Links:    ",
      ifelse (lss,
              namesof("scale"   , lscale   , escale),
              namesof("shape1.a", lshape1.a, eshape1.a)), ", ",
      ifelse (lss,
              namesof("shape1.a", lshape1.a, eshape1.a),
              namesof("scale"   , lscale   , escale)), ", ",
      namesof("shape3.q" , lshape3.q, earg = eshape3.q), "\n",
      "Mean:     scale * gamma(shape2.p + 1/shape1.a) * ",
                "gamma(shape3.q - 1/shape1.a) / ",
                "gamma(shape3.q)"),
  constraints = eval(substitute(expression({
    constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
                     M = M, M1 = 3,
                     predictors.names = predictors.names)
  }), list( .zero = zero ))),
  infos = eval(substitute(function(...) {
    list(M1 = 3,
         Q1 = 1,
         dpqrfun = "sinmad",  # "genbetaII",
         expected = TRUE,
         zero = .zero ,
         multipleResponses = TRUE,
         parameters.names = if ( .lss )
           c("scale", "shape1.a", "shape3.q") else
           c("shape1.a", "scale", "shape3.q"),
         lscale = .lscale ,       lshape1.a = .lshape1.a ,
         escale = .escale ,       eshape1.a = .eshape1.a ,
                                  lshape3.q = .lshape3.q ,
                                  eshape3.q = .eshape3.q )
  }, list( .lscale = lscale      , .lshape1.a = lshape1.a,
           .escale = escale      , .eshape1.a = eshape1.a,
                                   .lshape3.q = lshape3.q,
                                   .eshape3.q = eshape3.q,
           .lss  = lss ,
           .zero = zero ))),

  rqresslot = eval(substitute(
    function(mu, y, w, eta, extra = NULL) {
      M1  <- 3
      NOS <- ncol(eta) / M1
      if ( .lss ) {
        Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                           .lscale    , earg = .escale    )
        aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                           .lshape1.a , earg = .eshape1.a )
      } else {
        aa    <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                           .lshape1.a , earg = .eshape1.a )
        Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                           .lscale    , earg = .escale    )
      }
      parg   <- 1
      qq     <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                          .lshape3.q , earg = .eshape3.q )
    scrambleseed <- runif(1)  # To scramble the seed
      qnorm(psinmad(y, scale = Scale, shape1.a = aa,
                    shape3.q = qq))
    }, list( .lscale    = lscale   , .lshape1.a = lshape1.a,
             .escale    = escale   , .eshape1.a = eshape1.a,
                                     .lshape3.q = lshape3.q,
                                     .eshape3.q = eshape3.q,
             .lss = lss ))),
  initialize = eval(substitute(expression({
    temp5 <- w.y.check(w = w, y = y,
                       Is.positive.y = TRUE,
                       ncol.w.max = Inf,
                       ncol.y.max = Inf,
                       out.wy = TRUE,
                       colsyperw = 1,
                       maximize = TRUE)
    y    <- temp5$y
    w    <- temp5$w
    M1 <- 3  # Number of parameters for one response
    NOS  <- ncoly <- ncol(y)
    M    <- M1*ncol(y)

    scaL.names  <- param.names("scale",     NOS, skip1 = TRUE)
    sha1.names  <- param.names("shape1.a",  NOS, skip1 = TRUE)
    sha3.names  <- param.names("shape3.q",  NOS, skip1 = TRUE)

    predictors.names <- c(
      if ( .lss ) {
        c(namesof(scaL.names , .lscale    , earg = .escale    ,
                  tag = FALSE),
          namesof(sha1.names , .lshape1.a , earg = .eshape1.a ,
                  tag = FALSE))
      } else {
        c(namesof(sha1.names , .lshape1.a , earg = .eshape1.a ,
                  tag = FALSE),
          namesof(scaL.names , .lscale    , earg = .escale    ,
                  tag = FALSE))
      },
      namesof(sha3.names , .lshape3.q , .eshape3.q , tag = FALSE))
    predictors.names <- predictors.names[interleave.VGAM(M,
                                                         M1 = M1)]

    if (!length(etastart)) {
      sc.init <-
      aa.init <-
      qq.init <- matrix(NA_real_, n, NOS)

      for (spp. in 1:NOS) {  # For each response 'y_spp.'... do:
        yvec <- y[, spp.]
        wvec <- w[, spp.]

        if ( .imethod == 1 ) {
          gscale     <- .gscale
          gshape1.a  <- .gshape1.a
          gshape3.q  <- .gshape3.q
          if (length( .iscale    )) gscale    <-
                                        rep_len( .iscale    , NOS)
          if (length( .ishape1.a )) gshape1.a <-
                                        rep_len( .ishape1.a , NOS)
          if (length( .ishape3.q )) gshape3.q <-
                                        rep_len( .ishape3.q , NOS)



        try.this <-
          grid.search4(gscale, gshape1.a, vov3 = 1, gshape3.q,
                       objfun = genbetaII.Loglikfun4,  # x = x,
                       y = yvec, w = wvec,  # extraargs = NULL,
                       ret.objfun = TRUE)  # Last value is \ell
          sc.init[, spp.] <- try.this["Value1"]
          aa.init[, spp.] <- try.this["Value2"]
          qq.init[, spp.] <- try.this["Value4"]
       } else {  # .imethod == 2
          qvec <- .probs.y
          ishape3.q <- if (length( .ishape3.q )) .ishape3.q else 1
          xvec <- log( (1-qvec)^(-1/ ishape3.q ) - 1 )
          fit0 <- lsfit(x = xvec, y = log(quantile(yvec,  qvec)))
          sc.init[, spp.] <- if (length( .iscale )) .iscale else
                             exp(fit0$coef[1])
          aa.init[, spp.] <- if (length( .ishape1.a ))
                             .ishape1.a else 1/fit0$coef[2]
          qq.init[, spp.] <- ishape3.q
        }
      }  # End of for (spp. ...)




      finite.mean <- 1 < aa.init * qq.init
      COP.use <- 1.15
      while (FALSE && any(!finite.mean)) {
    qq.init[!finite.mean] <- 0.1 + qq.init[!finite.mean] * COP.use
    aa.init[!finite.mean] <- 0.1 + aa.init[!finite.mean] * COP.use
    finite.mean <- 1 < aa.init * qq.init
      }

      etastart <-
        cbind(if ( .lss )
          cbind(theta2eta(sc.init, .lscale    , .escale    ),
                theta2eta(aa.init, .lshape1.a , .eshape1.a )) else
          cbind(theta2eta(aa.init, .lshape1.a , .eshape1.a ),
                theta2eta(sc.init, .lscale    , .escale    )),
          theta2eta(qq.init , .lshape3.q , earg = .eshape3.q ))
      etastart <- etastart[, interleave.VGAM(M, M1 = M1)]
    }  # End of etastart.
  }), list( .lscale    = lscale   , .lshape1.a = lshape1.a,
            .escale    = escale   , .eshape1.a = eshape1.a,
            .iscale    = iscale   , .ishape1.a = ishape1.a,
            .gscale    = gscale   , .gshape1.a = gshape1.a,
                                    .lshape3.q = lshape3.q,
                                    .eshape3.q = eshape3.q,
                                    .ishape3.q = ishape3.q,
                                    .gshape3.q = gshape3.q,
            .imethod   = imethod   , .probs.y  = probs.y,
            .lss = lss ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    M1 <- 3
    NOS <- ncol(eta)/M1
    if ( .lss ) {
      Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                         .lscale ,  earg = .escale )
      aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a )
    } else {
      aa    <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a )
      Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lscale    ,  earg = .escale )
    }
    parg <- 1
    qq   <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                      .lshape3.q , earg = .eshape3.q )

    ans <- Scale * exp(lgamma(parg + 1/aa) +
           lgamma(qq   - 1/aa) - lgamma(parg) - lgamma(qq))
    ans[parg + 1/aa <= 0] <- NA
    ans[qq   - 1/aa <= 0] <- NA
    ans[aa          <= 0] <- NA
    ans[Scale       <= 0] <- NA
    ans[qq          <= 0] <- NA
    ans
  }, list( .lscale    = lscale   , .lshape1.a = lshape1.a,
           .escale    = escale   , .eshape1.a = eshape1.a,
                                   .lshape3.q = lshape3.q,
                                   .eshape3.q = eshape3.q,
           .lss = lss ))),
  last = eval(substitute(expression({
    M1 <- 3

    misc$link <- c(rep_len( if ( .lss ) .lscale else .lshape1.a ,
                           ncoly),
                   rep_len( if ( .lss ) .lshape1.a else .lscale ,
                           ncoly),
                   rep_len( .lshape3.q , ncoly))[
                   interleave.VGAM(M, M1 = M1)]
    temp.names <- if ( .lss ) {
      c(scaL.names, sha1.names,             sha3.names)
    } else {
      c(sha1.names, scaL.names,             sha3.names)
    }
    names(misc$link) <- temp.names[interleave.VGAM(M, M1 = M1)]

    misc$earg <- vector("list", M)
    names(misc$earg) <- temp.names[interleave.VGAM(M, M1 = M1)]
    for (ii in 1:ncoly) {
      if ( .lss ) {
        misc$earg[[M1*ii-2]] <- .escale
        misc$earg[[M1*ii-1]] <- .eshape1.a
    } else {
        misc$earg[[M1*ii-2]] <- .eshape1.a
        misc$earg[[M1*ii-1]] <- .escale
    }
      misc$earg[[M1*ii  ]] <- .eshape3.q
    }

    misc$expected          <- TRUE
    misc$multipleResponses <- TRUE
  }), list( .lscale    = lscale   , .lshape1.a = lshape1.a,
            .escale    = escale   , .eshape1.a = eshape1.a,
                                    .lshape3.q = lshape3.q,
                                    .eshape3.q = eshape3.q,
            .lss = lss ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE,
             eta, extra = NULL, summation = TRUE) {
      M1  <- 3
      NOS <- ncol(eta) / M1
      if ( .lss ) {
        Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                           .lscale    , earg = .escale    )
        aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                           .lshape1.a , earg = .eshape1.a )
      } else {
        aa    <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                           .lshape1.a , earg = .eshape1.a )
        Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                           .lscale    , earg = .escale    )
      }
      parg   <- 1
      qq     <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                          .lshape3.q , earg = .eshape3.q )
      if (residuals) {
        stop("loglikelihood residuals not implemented yet")
      } else {
        ll.elts <-
          c(w) * dgenbetaII(y, scale = Scale, shape1.a = aa,
                            shape2.p = parg, shape3.q = qq,
                            log = TRUE)
        if (summation) {
          sum(ll.elts)
        } else {
          ll.elts
        }
      }
    }, list( .lscale    = lscale   , .lshape1.a = lshape1.a,
             .escale    = escale   , .eshape1.a = eshape1.a,
                                     .lshape3.q = lshape3.q,
                                     .eshape3.q = eshape3.q,
             .lss = lss ))),
  vfamily = c("sinmad"),
  simslot = eval(substitute(
  function(object, nsim) {

    pwts <- if (length(pwts <- object@prior.weights) > 0)
              pwts else weights(object, type = "prior")
    if (any(pwts != 1))
      warning("ignoring prior weights")

    eta <- predict(object)
    M1 <- 3
    NOS <- ncol(eta)/M1
    if ( .lss ) {
      Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                         .lscale ,  earg = .escale )
      aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a )
    } else {
      aa    <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a )
      Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lscale    ,  earg = .escale )
    }
    qq   <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                      .lshape3.q , earg = .eshape3.q )

    rsinmad(nsim * length(qq), shape1.a = aa, scale = Scale,
            shape3.q = qq)
  }, list( .lscale    = lscale   , .lshape1.a = lshape1.a,
           .escale    = escale   , .eshape1.a = eshape1.a,
                                   .lshape3.q = lshape3.q,
                                   .eshape3.q = eshape3.q
                      ))),

  validparams = eval(substitute(function(eta, y, extra = NULL) {
    M1 <- 3
    NOS <- ncol(eta) / M1
    if ( .lss ) {
      Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                         .lscale    , earg = .escale   )
      aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a)
    } else {
      aa    <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a)
      Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lscale    , earg = .escale )
    }
    parg <- 1
    qq   <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                      .lshape3.q , earg = .eshape3.q)

    okay1 <- all(is.finite(Scale)) && all(Scale  > 0) &&
             all(is.finite(aa   )) && all(aa     > 0) &&
             all(is.finite(parg )) && all(parg   > 0) &&
             all(is.finite(qq   )) && all(qq     > 0)
    okay.support <- if (okay1)
                    all(-aa * parg < 1 & 1 < aa * qq) else TRUE
    if (!okay.support)
      warning("parameter constraint -a   < 1 < a*q has ",
              "been violated; solution may be at the",
              " boundary of the parameter space.")
    okay1 && okay.support
  }, list(  .lscale    = lscale   , .lshape1.a = lshape1.a,
            .escale    = escale   , .eshape1.a = eshape1.a,
                                    .lshape3.q = lshape3.q,
                                    .eshape3.q = eshape3.q,
            .lss = lss ))),




  deriv = eval(substitute(expression({
    M1 <- 3
    NOS <- ncol(eta)/M1  # Needed for summary()
    if ( .lss ) {
      Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                         .lscale    , earg = .escale   )
      aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a)
    } else {
      aa    <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a)
      Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lscale    , earg = .escale )
    }
    parg <- 1
    qq   <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                      .lshape3.q , earg = .eshape3.q)
    temp1 <- log(y/Scale)
    temp2 <- (y/Scale)^aa
    temp3 <- digamma(parg + qq)
    temp3a <- digamma(parg)
    temp3b <- digamma(qq)
    temp4 <- log1p(temp2)

    dl.dscale <- (aa/Scale) * (-parg + (parg+qq) / (1+1/temp2))
    dl.da <- 1/aa + parg * temp1 - (parg+qq) * temp1 / (1+1/temp2)
    dl.dq <- temp3 - temp3b - temp4

    dscale.deta <- dtheta.deta(Scale, .lscale ,    .escale )
    da.deta     <- dtheta.deta(aa,    .lshape1.a , .eshape1.a )
    dq.deta     <- dtheta.deta(qq,    .lshape3.q , .eshape3.q )

    myderiv <- if ( .lss ) {
      c(w) * cbind(dl.dscale * dscale.deta,
                   dl.da * da.deta,
                   dl.dq * dq.deta)
    } else {
      c(w) * cbind(dl.da * da.deta,
                   dl.dscale * dscale.deta,
                   dl.dq * dq.deta)
    }
    myderiv[, interleave.VGAM(M, M1 = M1)]
  }), list(  .lscale    = lscale   , .lshape1.a = lshape1.a,
             .escale    = escale   , .eshape1.a = eshape1.a,
                                     .lshape3.q = lshape3.q,
                                     .eshape3.q = eshape3.q,
             .lss = lss ))),
  weight = eval(substitute(expression({
    temp5  <- trigamma(parg + qq)
    temp5a <- trigamma(parg)
    temp5b <- trigamma(qq)

    ned2l.da <- (1 + parg + qq + parg * qq * (temp5a + temp5b +
                (temp3b - temp3a + (parg-qq)/(parg*qq))^2 -
          (parg^2 + qq^2) / (parg*qq)^2)) / (aa^2 * (1+parg+qq))
    ned2l.dscale  <- (aa^2) * parg * qq / ((1+parg+qq) * Scale^2)
    ned2l.dq      <- temp5b - temp5
    ned2l.dascale <- (parg - qq - parg * qq *
                     (temp3a -temp3b)) / (Scale*(1 + parg+qq))
    ned2l.daq     <- -(parg * (temp3b -temp3a) -1) / (aa*(parg+qq))
    ned2l.dscaleq <- -aa * parg / (Scale*(parg+qq))
    wz <- if ( .lss ) {
      array(c(c(w) * ned2l.dscale * dscale.deta^2,
              c(w) * ned2l.da * da.deta^2,
              c(w) * ned2l.dq * dq.deta^2,
              c(w) * ned2l.dascale * da.deta * dscale.deta,
              c(w) * ned2l.daq * da.deta * dq.deta,
              c(w) * ned2l.dscaleq * dscale.deta * dq.deta),
            dim = c(n, M/M1, M1*(M1+1)/2))
    } else {
      array(c(c(w) * ned2l.da * da.deta^2,
              c(w) * ned2l.dscale * dscale.deta^2,
              c(w) * ned2l.dq * dq.deta^2,
              c(w) * ned2l.dascale * da.deta * dscale.deta,
              c(w) * ned2l.dscaleq * dscale.deta * dq.deta,
              c(w) * ned2l.daq * da.deta * dq.deta),
            dim = c(n, M/M1, M1*(M1+1)/2))
    }
    wz <- arwz2wz(wz, M = M, M1 = M1)
    wz
  }), list( .lscale    = lscale   , .lshape1.a = lshape1.a,
            .escale    = escale   , .eshape1.a = eshape1.a,
                                    .lshape3.q = lshape3.q,
                                    .eshape3.q = eshape3.q,
            .lss = lss ))))
}  # sinmad






 dagum <-
  function(lscale    = "loglink",
           lshape1.a = "loglink",
           lshape2.p = "loglink",
           iscale    = NULL,
           ishape1.a = NULL,
           ishape2.p = NULL,
           imethod   = 1,
           lss       = TRUE,
           gscale    = exp(-5:5),
           gshape1.a = seq(0.75, 4, by = 0.25),  # exp(-5:5),
           gshape2.p = exp(-5:5),
           probs.y   = c(0.25, 0.50, 0.75),
           zero      = "shape") {





  if (length(lss) != 1 && !is.logical(lss))
    stop("Argument 'lss' not specified correctly")

  if (!is.Numeric(imethod, length.arg = 1,
                  integer.valued = TRUE,
                  positive = TRUE) || imethod > 2)
    stop("Bad input for argument 'imethod'")

  if (length(iscale) && !is.Numeric(iscale, positive = TRUE))
    stop("Bad input for argument 'iscale'")

  if (length(ishape1.a) && !is.Numeric(ishape1.a, positive = TRUE))
    stop("Bad input for argument 'ishape1.a'")

  if (length(ishape2.p) && !is.Numeric(ishape2.p, positive = TRUE))
    stop("Bad input for argument 'ishape2.p'")

  if (length(probs.y) < 2 || max(probs.y) > 1 ||
        !is.Numeric(probs.y, positive = TRUE))
    stop("Bad input for argument 'probs.y'")


  lscale <- as.list(substitute(lscale))
  escale <- link2list(lscale)
  lscale <- attr(escale, "function.name")

  lshape1.a <- as.list(substitute(lshape1.a))
  eshape1.a <- link2list(lshape1.a)
  lshape1.a <- attr(eshape1.a, "function.name")

  lshape2.p <- as.list(substitute(lshape2.p))
  eshape2.p <- link2list(lshape2.p)
  lshape2.p <- attr(eshape2.p, "function.name")


  new("vglmff",
  blurb =
    c("Dagum distribution \n\n",
      "Links:    ",
      ifelse (lss,
              namesof("scale"   , lscale   , escale),
              namesof("shape1.a", lshape1.a, eshape1.a)), ", ",
      ifelse (lss,
              namesof("shape1.a", lshape1.a, eshape1.a),
              namesof("scale"   , lscale   , escale)), ", ",
      namesof("shape2.p" , lshape2.p, earg = eshape2.p), "\n",
      "Mean:     scale * gamma(shape2.p + 1/shape1.a) * ",
                "gamma(1 - 1/shape1.a) / ",
                "gamma(shape2.p)"),
  constraints = eval(substitute(expression({
    constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
                     M = M, M1 = 3,
                     predictors.names = predictors.names)

  }), list( .zero = zero ))),
  infos = eval(substitute(function(...) {
    list(M1 = 3,
         Q1 = 1,
         dpqrfun = "dagum",  # "genbetaII",
         expected = TRUE,
         zero = .zero ,
         multipleResponses = TRUE,
         parameters.names =
           if ( .lss ) c("scale", "shape1.a", "shape2.p") else
                       c("shape1.a", "scale", "shape2.p"),
         lscale    = .lscale    , lshape1.a = .lshape1.a ,
         escale    = .escale    , eshape1.a = .eshape1.a ,
         lshape2.p = .lshape2.p ,
         eshape2.p = .eshape2.p )
  }, list( .lscale = lscale      , .lshape1.a = lshape1.a,
           .escale = escale      , .eshape1.a = eshape1.a,
           .lshape2.p = lshape2.p,
           .eshape2.p = eshape2.p,
           .lss  = lss , .zero = zero ))),

  rqresslot = eval(substitute(
    function(mu, y, w, eta, extra = NULL) {
      M1  <- 3
      NOS <- ncol(eta) / M1
      if ( .lss ) {
        Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                           .lscale    , earg = .escale    )
        aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                           .lshape1.a , earg = .eshape1.a )
      } else {
        aa    <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                           .lshape1.a , earg = .eshape1.a )
        Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                           .lscale    , earg = .escale    )
      }
      parg   <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                          .lshape2.p , earg = .eshape2.p )
    scrambleseed <- runif(1)  # To scramble the seed
      qnorm(pdagum(y, scale = Scale, shape1.a = aa,
                   shape2.p = parg))  # shape3.q = 1
    }, list( .lscale    = lscale   , .lshape1.a = lshape1.a,
             .escale    = escale   , .eshape1.a = eshape1.a,
             .lshape2.p = lshape2.p,
             .eshape2.p = eshape2.p,
             .lss = lss ))),

  initialize = eval(substitute(expression({
    temp5 <- w.y.check(w = w, y = y,
                       Is.positive.y = TRUE,
                       ncol.w.max = Inf,
                       ncol.y.max = Inf,
                       out.wy = TRUE,
                       colsyperw = 1,
                       maximize = TRUE)
    y    <- temp5$y
    w    <- temp5$w
    M1 <- 3  # Number of parameters for one response
    NOS  <- ncoly <- ncol(y)
    M    <- M1*ncol(y)

    scaL.names  <- param.names("scale",     NOS, skip1 = TRUE)
    sha1.names  <- param.names("shape1.a",  NOS, skip1 = TRUE)
    sha2.names  <- param.names("shape2.p",  NOS, skip1 = TRUE)

    predictors.names <- c(
    if ( .lss ) {
    c(namesof(scaL.names , .lscale    , .escale    , tag = FALSE),
      namesof(sha1.names , .lshape1.a , .eshape1.a , tag = FALSE))
    } else {
    c(namesof(sha1.names , .lshape1.a , .eshape1.a , tag = FALSE),
      namesof(scaL.names , .lscale    , .escale    , tag = FALSE))
    },
    namesof(sha2.names , .lshape2.p , .eshape2.p , tag = FALSE))
    predictors.names <- predictors.names[interleave.VGAM(M,
                                                         M1 = M1)]

    if (!length(etastart)) {
      sc.init <-
      aa.init <-
      pp.init <- matrix(NA_real_, n, NOS)

      for (spp. in 1:NOS) {  # For each response 'y_spp.'... do:
        yvec <- y[, spp.]
        wvec <- w[, spp.]

        if ( .imethod == 1 ) {
          gscale     <- .gscale
          gshape1.a  <- .gshape1.a
          gshape2.p  <- .gshape2.p
          if (length( .iscale    )) gscale    <-
                                        rep_len( .iscale    , NOS)
          if (length( .ishape1.a )) gshape1.a <-
                                        rep_len( .ishape1.a , NOS)
          if (length( .ishape2.p )) gshape2.p <-
                                        rep_len( .ishape2.p , NOS)



        try.this <-
          grid.search4(gscale, gshape1.a, gshape2.p, vov4 = 1,
                       objfun = genbetaII.Loglikfun4,  # x = x,
                       y = yvec, w = wvec,  # extraargs = NULL,
                       ret.objfun = TRUE)  # Last value is \ell
          sc.init[, spp.] <- try.this["Value1"]
          aa.init[, spp.] <- try.this["Value2"]
          pp.init[, spp.] <- try.this["Value3"]
       } else {  # .imethod == 2
          qvec <- .probs.y
          ishape2.p <- if (length( .ishape2.p )) .ishape2.p else 1
          xvec <- log( qvec^(-1/ ishape2.p) - 1 )
          fit0 <- lsfit(x = xvec, y = log(quantile(yvec,  qvec)))
          sc.init[, spp.] <- if (length( .iscale )) .iscale else
                             exp(fit0$coef[1])
          aa.init[, spp.] <- if (length( .ishape1.a )) .ishape1.a else
                             -1/fit0$coef[2]
          pp.init[, spp.] <- ishape2.p
      }
    }  # End of for (spp. ...)



      finite.mean <- 1 < aa.init
      COP.use <- 1.15
      while (FALSE && any(!finite.mean)) {
   aa.init[!finite.mean] <- 0.1 + aa.init[!finite.mean] * COP.use
   finite.mean <- 1 < aa.init
      }

      etastart <-
        cbind(if ( .lss )
          cbind(theta2eta(sc.init, .lscale    , .escale    ),
                theta2eta(aa.init, .lshape1.a , .eshape1.a )) else
          cbind(theta2eta(aa.init, .lshape1.a , .eshape1.a ),
                theta2eta(sc.init, .lscale    , .escale    )),
          theta2eta(pp.init , .lshape2.p , earg = .eshape2.p ))
      etastart <- etastart[, interleave.VGAM(M, M1 = M1)]
    }  # End of etastart.
  }), list( .lscale    = lscale   , .lshape1.a = lshape1.a,
            .escale    = escale   , .eshape1.a = eshape1.a,
            .iscale    = iscale   , .ishape1.a = ishape1.a,
            .gscale    = gscale   , .gshape1.a = gshape1.a,
            .lshape2.p = lshape2.p,
            .eshape2.p = eshape2.p,
            .ishape2.p = ishape2.p,
            .gshape2.p = gshape2.p,
            .imethod   = imethod   , .probs.y = probs.y,
            .lss = lss ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    M1 <- 3
    NOS <- ncol(eta)/M1
    if ( .lss ) {
      Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                         .lscale ,  earg = .escale )
      aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a )
    } else {
      aa    <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a )
      Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lscale    ,  earg = .escale )
    }
    parg <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                      .lshape2.p , earg = .eshape2.p )

    qq     <- 1

    ans <- Scale * exp(lgamma(parg + 1/aa) +
                       lgamma(qq   - 1/aa) -
                       lgamma(parg) - lgamma(qq))
    ans[parg + 1/aa <= 0] <- NA
    ans[qq   - 1/aa <= 0] <- NA
    ans[aa          <= 0] <- NA
    ans[Scale       <= 0] <- NA
    ans[parg        <= 0] <- NA
    ans
  }, list( .lscale    = lscale   , .lshape1.a = lshape1.a,
           .escale    = escale   , .eshape1.a = eshape1.a,
           .lshape2.p = lshape2.p,
           .eshape2.p = eshape2.p,
           .lss = lss ))),
  last = eval(substitute(expression({
    M1 <- 3

    misc$link <- c(rep_len( if ( .lss ) .lscale else .lshape1.a ,
                           ncoly),
                   rep_len( if ( .lss ) .lshape1.a else .lscale ,
                           ncoly),
                   rep_len( .lshape2.p , ncoly))[
                   interleave.VGAM(M, M1 = M1)]
    temp.names <- if ( .lss ) {
      c(scaL.names, sha1.names, sha2.names)
    } else {
      c(sha1.names, scaL.names, sha2.names)
    }
    names(misc$link) <- temp.names[interleave.VGAM(M, M1 = M1)]

    misc$earg <- vector("list", M)
    names(misc$earg) <- temp.names[interleave.VGAM(M, M1 = M1)]
    for (ii in 1:ncoly) {
      if ( .lss ) {
        misc$earg[[M1*ii-2]] <- .escale
        misc$earg[[M1*ii-1]] <- .eshape1.a
      } else {
        misc$earg[[M1*ii-2]] <- .eshape1.a
        misc$earg[[M1*ii-1]] <- .escale
      }
      misc$earg[[M1*ii  ]] <- .eshape2.p
    }

    misc$expected          <- TRUE
    misc$multipleResponses <- TRUE
  }), list( .lscale    = lscale   , .lshape1.a = lshape1.a,
            .escale    = escale   , .eshape1.a = eshape1.a,
            .lshape2.p = lshape2.p,
            .eshape2.p = eshape2.p,
            .lss = lss ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE,
             eta, extra = NULL, summation = TRUE) {
      M1  <- 3
      NOS <- ncol(eta) / M1
      if ( .lss ) {
        Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                           .lscale    , earg = .escale    )
        aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                           .lshape1.a , earg = .eshape1.a )
      } else {
        aa    <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                           .lshape1.a , earg = .eshape1.a )
        Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                           .lscale    , earg = .escale    )
      }
      parg   <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                          .lshape2.p , earg = .eshape2.p )
      if (residuals) {
        stop("loglikelihood residuals not implemented yet")
      } else {
        ll.elts <-
          c(w) * dgenbetaII(x = y, scale = Scale, shape1.a = aa,
                            shape2.p = parg, shape3.q = 1,
                            log = TRUE)
        if (summation) {
          sum(ll.elts)
        } else {
          ll.elts
        }
      }
    }, list( .lscale    = lscale   , .lshape1.a = lshape1.a,
             .escale    = escale   , .eshape1.a = eshape1.a,
             .lshape2.p = lshape2.p,
             .eshape2.p = eshape2.p,
             .lss = lss ))),
  vfamily = c("dagum"),

  simslot = eval(substitute(
  function(object, nsim) {

    pwts <- if (length(pwts <- object@prior.weights) > 0)
              pwts else weights(object, type = "prior")
    if (any(pwts != 1))
      warning("ignoring prior weights")
    eta <- predict(object)
    if ( .lss ) {
      Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                         .lscale    , earg = .escale   )
      aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a)
    } else {
      aa    <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a)
      Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lscale    , earg = .escale )
    }
    parg <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                      .lshape2.p , earg = .eshape2.p)
    qq   <- 1
    rdagum(nsim * length(parg), shape1.a = aa, scale = Scale,
           shape2.p = parg)
  }, list(  .lscale    = lscale   , .lshape1.a = lshape1.a,
            .escale    = escale   , .eshape1.a = eshape1.a,
            .lshape2.p = lshape2.p,
            .eshape2.p = eshape2.p,
            .lss = lss ))),



  validparams = eval(substitute(function(eta, y, extra = NULL) {
    M1 <- 3
    NOS <- ncol(eta)/M1  # Needed for summary()
    if ( .lss ) {
      Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                         .lscale    , earg = .escale   )
      aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a)
    } else {
      aa    <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a)
      Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lscale    , earg = .escale )
    }
    parg <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                      .lshape2.p , earg = .eshape2.p)
    qq   <- 1


    okay1 <- all(is.finite(Scale)) && all(Scale  > 0) &&
             all(is.finite(aa   )) && all(aa     > 0) &&
             all(is.finite(parg )) && all(parg   > 0) &&
             all(is.finite(qq   )) && all(qq     > 0)
    okay.support <- if (okay1)
                    all(-aa * parg < 1 & 1 < aa * qq) else TRUE
    if (!okay.support)
      warning("parameter constraint -a*p < 1 < a   has ",
              "been violated; solution may be at",
              " the boundary of the parameter space.")
    okay1 && okay.support
  }, list(  .lscale    = lscale   , .lshape1.a = lshape1.a,
            .escale    = escale   , .eshape1.a = eshape1.a,
            .lshape2.p = lshape2.p,
            .eshape2.p = eshape2.p,
            .lss = lss ))),




  deriv = eval(substitute(expression({
    M1 <- 3
    NOS <- ncol(eta)/M1  # Needed for summary()
    if ( .lss ) {
      Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                         .lscale    , earg = .escale   )
      aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a)
    } else {
      aa    <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a)
      Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lscale    , earg = .escale )
    }
    parg <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                      .lshape2.p , earg = .eshape2.p)
    qq   <- 1
    temp1 <- log(y/Scale)
    temp2 <- (y/Scale)^aa
    temp3 <- digamma(parg + qq)
    temp3a <- digamma(parg)
    temp3b <- digamma(qq)
    temp4 <- log1p(temp2)

    dl.dscale <- (aa/Scale) * (-parg + (parg+qq) / (1+1/temp2))
    dl.da <- 1/aa + parg * temp1 - (parg+qq) * temp1 / (1+1/temp2)
    dl.dp <- aa * temp1 + temp3 - temp3a - temp4

    dscale.deta <- dtheta.deta(Scale, .lscale ,    .escale )
    da.deta     <- dtheta.deta(aa,    .lshape1.a , .eshape1.a )
    dp.deta     <- dtheta.deta(parg,  .lshape2.p , .eshape2.p )

    myderiv <- if ( .lss ) {
      c(w) * cbind(dl.dscale * dscale.deta,
                   dl.da * da.deta,
                   dl.dp * dp.deta)
    } else {
      c(w) * cbind(dl.da * da.deta,
                   dl.dscale * dscale.deta,
                   dl.dp * dp.deta)
    }
    myderiv[, interleave.VGAM(M, M1 = M1)]
  }), list(  .lscale    = lscale   , .lshape1.a = lshape1.a,
             .escale    = escale   , .eshape1.a = eshape1.a,
             .lshape2.p = lshape2.p,
             .eshape2.p = eshape2.p,
             .lss = lss ))),
  weight = eval(substitute(expression({
    temp5  <- trigamma(parg + qq)
    temp5a <- trigamma(parg)
    temp5b <- trigamma(qq)

    ned2l.da <- (1 + parg + qq + parg * qq * (temp5a + temp5b +
                (temp3b - temp3a + (parg-qq)/(parg*qq))^2 -
        (parg^2 + qq^2) / (parg*qq)^2)) / (aa^2 * (1+parg+qq))
    ned2l.dscale  <- (aa^2) * parg * qq / ((1+parg+qq) * Scale^2)
    ned2l.dp      <- temp5a - temp5
    ned2l.dascale <- (parg - qq - parg * qq *
                     (temp3a -temp3b)) / (Scale*(1 + parg+qq))
    ned2l.dap     <- -(qq   * (temp3a -temp3b) -1) / (aa*(parg+qq))
    ned2l.dscalep <-  aa * qq   / (Scale*(parg+qq))
    wz <- if ( .lss ) {
      array(c(c(w) * ned2l.dscale * dscale.deta^2,
              c(w) * ned2l.da * da.deta^2,
              c(w) * ned2l.dp * dp.deta^2,
              c(w) * ned2l.dascale * da.deta * dscale.deta,
              c(w) * ned2l.dap * da.deta * dp.deta,
              c(w) * ned2l.dscalep * dscale.deta * dp.deta),
            dim = c(n, M/M1, M1*(M1+1)/2))
    } else {
      array(c(c(w) * ned2l.da * da.deta^2,
              c(w) * ned2l.dscale * dscale.deta^2,
              c(w) * ned2l.dp * dp.deta^2,
              c(w) * ned2l.dascale * da.deta * dscale.deta,
              c(w) * ned2l.dscalep * dscale.deta * dp.deta,
              c(w) * ned2l.dap * da.deta * dp.deta),
            dim = c(n, M/M1, M1*(M1+1)/2))
    }
    wz <- arwz2wz(wz, M = M, M1 = M1)
    wz
  }), list( .lscale    = lscale   , .lshape1.a = lshape1.a,
            .escale    = escale   , .eshape1.a = eshape1.a,
            .lshape2.p = lshape2.p,
            .eshape2.p = eshape2.p,
            .lss = lss ))))
}  # dagum






 betaII <-
  function(lscale    = "loglink",
           lshape2.p = "loglink",
           lshape3.q = "loglink",
           iscale    = NULL,
           ishape2.p = NULL,
           ishape3.q = NULL,
           imethod   = 1,
           gscale    = exp(-5:5),
           gshape2.p = exp(-5:5),
           gshape3.q = seq(0.75, 4, by = 0.25),  # exp(-5:5),
           probs.y   = c(0.25, 0.50, 0.75),
           zero      = "shape") {




  if (!is.Numeric(imethod, length.arg = 1,
                  integer.valued = TRUE,
                  positive = TRUE) || imethod > 2)
    stop("Bad input for argument 'imethod'")

  if (length(iscale   ) && !is.Numeric(iscale   ,  positive = TRUE))
    stop("Bad input for argument 'iscale'")

  if (length(ishape2.p) && !is.Numeric(ishape2.p, positive = TRUE))
    stop("Bad input for argument 'ishape2.p'")

  if (length(ishape3.q) && !is.Numeric(ishape3.q, positive = TRUE))
    stop("Bad input for argument 'ishape3.q'")

  if (length(probs.y) < 2 || max(probs.y) > 1 ||
        !is.Numeric(probs.y, positive = TRUE))
    stop("Bad input for argument 'probs.y'")


  lscale <- as.list(substitute(lscale))
  escale <- link2list(lscale)
  lscale <- attr(escale, "function.name")

  lshape2.p <- as.list(substitute(lshape2.p))
  eshape2.p <- link2list(lshape2.p)
  lshape2.p <- attr(eshape2.p, "function.name")

  lshape3.q <- as.list(substitute(lshape3.q))
  eshape3.q <- link2list(lshape3.q)
  lshape3.q <- attr(eshape3.q, "function.name")


  new("vglmff",
  blurb =
    c("Beta II distribution \n\n",
      "Links:    ",
      namesof("scale"    , lscale   , earg = escale   ), ", ",
      namesof("shape2.p" , lshape2.p, earg = eshape2.p), ", ",
      namesof("shape3.q" , lshape3.q, earg = eshape3.q), "\n",
      "Mean:     scale * gamma(shape2.p + 1) * ",
                "gamma(shape3.q - 1) / ",
                "(gamma(shape2.p) * gamma(shape3.q))"),
  constraints = eval(substitute(expression({
    constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
                     M = M, M1 = 3,
                     predictors.names = predictors.names)

  }), list( .zero = zero ))),
  infos = eval(substitute(function(...) {
    list(M1 = 3,
         Q1 = 1,
         dpqrfun = "betaII",  # "genbetaII",
         expected = TRUE,
         zero = .zero ,
         multipleResponses = TRUE,
         parameters.names = c("scale", "shape2.p", "shape3.q"),
         lscale    = .lscale    ,
         escale    = .escale    ,
         lshape2.p = .lshape2.p , lshape3.q = .lshape3.q ,
         eshape2.p = .eshape2.p , eshape3.q = .eshape3.q )
  }, list( .lscale = lscale      ,
           .escale = escale      ,
           .lshape2.p = lshape2.p, .lshape3.q = lshape3.q,
           .eshape2.p = eshape2.p, .eshape3.q = eshape3.q,
           .zero = zero ))),
  initialize = eval(substitute(expression({
    temp5 <- w.y.check(w = w, y = y,
                       Is.positive.y = TRUE,
                       ncol.w.max = Inf,
                       ncol.y.max = Inf,
                       out.wy = TRUE,
                       colsyperw = 1,
                       maximize = TRUE)
    y    <- temp5$y
    w    <- temp5$w
    M1 <- 3  # Number of parameters for one response
    NOS  <- ncoly <- ncol(y)
    M    <- M1*ncol(y)

    scaL.names  <- param.names("scale",     NOS, skip1 = TRUE)
    sha2.names  <- param.names("shape2.p",  NOS, skip1 = TRUE)
    sha3.names  <- param.names("shape3.q",  NOS, skip1 = TRUE)

    predictors.names <-
    c(namesof(scaL.names , .lscale    , earg = .escale    , tag = FALSE),
      namesof(sha2.names , .lshape2.p , earg = .eshape2.p , tag = FALSE),
      namesof(sha3.names , .lshape3.q , earg = .eshape3.q , tag = FALSE))
    predictors.names <- predictors.names[interleave.VGAM(M, M1 = M1)]

    if (!length(etastart)) {
      sc.init <-
      pp.init <-
      qq.init <- matrix(NA_real_, n, NOS)

      for (spp. in 1:NOS) {  # For each response 'y_spp.'... do:
        yvec <- y[, spp.]
        wvec <- w[, spp.]

        if ( .imethod == 1 ) {
          gscale     <- .gscale
          gshape2.p  <- .gshape2.p
          gshape3.q  <- .gshape3.q
          if (length( .iscale    )) gscale    <- rep_len( .iscale    , NOS)
          if (length( .ishape2.p )) gshape2.p <- rep_len( .ishape2.p , NOS)
          if (length( .ishape3.q )) gshape3.q <- rep_len( .ishape3.q , NOS)




        try.this <-
          grid.search4(gscale, vov2 = 1, gshape2.p, gshape3.q,
                       objfun = genbetaII.Loglikfun4,  # x = x,
                       y = yvec, w = wvec,  # extraargs = NULL,
                       ret.objfun = TRUE)  # Last value is \ell
          sc.init[, spp.] <- try.this["Value1"]
          pp.init[, spp.] <- try.this["Value3"]
          qq.init[, spp.] <- try.this["Value4"]
      } else {  # .imethod == 2

        sc.init[, spp.] <- if (length( .iscale )) .iscale else {
          qvec <- .probs.y
          ishape3.q <- if (length( .ishape3.q )) .ishape3.q else 1
          xvec <- log( (1-qvec)^(-1/ ishape3.q ) - 1 )
          fit0 <- lsfit(x = xvec, y = log(quantile(yvec, qvec )))
          exp(fit0$coef[1])
        }
        pp.init[, spp.] <- if (length( .ishape2.p )) .ishape2.p else 1.0
        qq.init[, spp.] <- if (length( .ishape3.q )) .ishape3.q else 1.0
       }
     }  # End of for (spp. ...)




      finite.mean <- 1 < qq.init
      COP.use <- 1.15
      while (any(!finite.mean)) {
        qq.init[!finite.mean] <- 0.1 + qq.init[!finite.mean] * COP.use
        finite.mean <- 1 < qq.init
      }

      etastart <-
        cbind(theta2eta(sc.init , .lscale    , earg = .escale    ),
              theta2eta(pp.init , .lshape2.p , earg = .eshape2.p ),
              theta2eta(qq.init , .lshape3.q , earg = .eshape3.q ))
      etastart <- etastart[, interleave.VGAM(M, M1 = M1)]
    }  # End of etastart.
  }), list( .lscale    = lscale   ,
            .escale    = escale   ,
            .iscale    = iscale   ,
            .gscale    = gscale   ,
            .lshape2.p = lshape2.p, .lshape3.q = lshape3.q,
            .eshape2.p = eshape2.p, .eshape3.q = eshape3.q,
            .ishape2.p = ishape2.p, .ishape3.q = ishape3.q,
            .gshape2.p = gshape2.p, .gshape3.q = gshape3.q,
            .imethod   = imethod   , .probs.y = probs.y
                       ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    M1 <- 3
    NOS <- ncol(eta)/M1
    Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                       .lscale ,  earg = .escale )
    aa    <- 1
    parg  <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                      .lshape2.p , earg = .eshape2.p )
    qq    <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                      .lshape3.q , earg = .eshape3.q )

    ans <- Scale * exp(lgamma(parg + 1/aa) +
                       lgamma(qq   - 1/aa) - lgamma(parg) - lgamma(qq))
    ans[parg + 1/aa <= 0] <- NA
    ans[qq   - 1/aa <= 0] <- NA
    ans[Scale       <= 0] <- NA
    ans[parg        <= 0] <- NA
    ans[qq          <= 0] <- NA
    ans
  }, list( .lscale    = lscale   ,
           .escale    = escale   ,
           .lshape2.p = lshape2.p, .lshape3.q = lshape3.q,
           .eshape2.p = eshape2.p, .eshape3.q = eshape3.q
                      ))),
  last = eval(substitute(expression({
    M1 <- 3

    misc$link <-
      c(rep_len( .lscale    , ncoly),
        rep_len( .lshape2.p , ncoly),
        rep_len( .lshape3.q , ncoly))[interleave.VGAM(M, M1 = M1)]
    temp.names <- c(scaL.names,             sha2.names, sha3.names)
    names(misc$link) <- temp.names[interleave.VGAM(M, M1 = M1)]

    misc$earg <- vector("list", M)
    names(misc$earg) <- temp.names[interleave.VGAM(M, M1 = M1)]
    for (ii in 1:ncoly) {
      misc$earg[[M1*ii-2]] <- .escale
      misc$earg[[M1*ii-1]] <- .eshape2.p
      misc$earg[[M1*ii  ]] <- .eshape3.q
    }

    misc$expected          <- TRUE
    misc$multipleResponses <- TRUE
  }), list( .lscale    = lscale   ,
            .escale    = escale   ,
            .lshape2.p = lshape2.p, .lshape3.q = lshape3.q,
            .eshape2.p = eshape2.p, .eshape3.q = eshape3.q
                       ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE,
             eta, extra = NULL, summation = TRUE) {
      M1  <- 3
      NOS <- ncol(eta)/M1
      Scale  <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                          .lscale    , earg = .escale    )
      aa     <- 1
      parg   <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                          .lshape2.p , earg = .eshape2.p )
      qq     <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                          .lshape3.q , earg = .eshape3.q )
      if (residuals) {
        stop("loglikelihood residuals not implemented yet")
      } else {
        ll.elts <-
          c(w) * dgenbetaII(x = y, scale = Scale, shape1.a = aa,
                            shape2.p = parg, shape3.q = qq, log = TRUE)
        if (summation) {
          sum(ll.elts)
        } else {
          ll.elts
        }
      }
    }, list( .lscale    = lscale   ,
             .escale    = escale   ,
             .lshape2.p = lshape2.p, .lshape3.q = lshape3.q,
             .eshape2.p = eshape2.p, .eshape3.q = eshape3.q
                        ))),
  vfamily = c("betaII"),

  validparams = eval(substitute(function(eta, y, extra = NULL) {
    M1 <- 3
    NOS <- ncol(eta)/M1
    Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                       .lscale    , earg = .escale    )
    aa    <- 1
    parg  <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                       .lshape2.p , earg = .eshape2.p )
    qq    <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                       .lshape3.q , earg = .eshape3.q )


    okay1 <- all(is.finite(Scale)) && all(Scale  > 0) &&
             all(is.finite(aa   )) && all(aa     > 0) &&
             all(is.finite(parg )) && all(parg   > 0) &&
             all(is.finite(qq   )) && all(qq     > 0)
    okay.support <- if (okay1)
                    all(-aa * parg < 1 & 1 < aa * qq) else TRUE
    if (!okay.support)
      warning("parameter constraint -p < 1 <   q has",
              " been violated; solution may be at the ",
              " boundary of the parameter space.")
    okay1 && okay.support
  }, list(  .lscale    = lscale   ,
            .escale    = escale   ,
            .lshape2.p = lshape2.p, .lshape3.q = lshape3.q,
            .eshape2.p = eshape2.p, .eshape3.q = eshape3.q ))),




  deriv = eval(substitute(expression({
    M1 <- 3
    NOS <- ncol(eta)/M1  # Needed for summary()
    Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE],
                       .lscale    , earg = .escale    )
    aa    <- 1
    parg  <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                       .lshape2.p , earg = .eshape2.p )
    qq    <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                       .lshape3.q , earg = .eshape3.q )
    temp1 <- log(y/Scale)
    temp2 <- (y/Scale)^aa
    temp3 <- digamma(parg + qq)
    temp3a <- digamma(parg)
    temp3b <- digamma(qq)
    temp4 <- log1p(temp2)

    dl.dscale <- (aa/Scale) * (-parg + (parg+qq) / (1+1/temp2))
    dl.dp <- aa * temp1 + temp3 - temp3a - temp4
    dl.dq <- temp3 - temp3b - temp4

    dscale.deta <- dtheta.deta(Scale, .lscale ,    .escale )
    dp.deta     <- dtheta.deta(parg,  .lshape2.p , .eshape2.p )
    dq.deta     <- dtheta.deta(qq,    .lshape3.q , .eshape3.q )

    myderiv <-
      c(w) * cbind(dl.dscale * dscale.deta,
                   dl.dp * dp.deta,
                   dl.dq * dq.deta)
    myderiv[, interleave.VGAM(M, M1 = M1)]
  }), list(  .lscale    = lscale   ,
             .escale    = escale   ,
             .lshape2.p = lshape2.p, .lshape3.q = lshape3.q,
             .eshape2.p = eshape2.p, .eshape3.q = eshape3.q
                        ))),
  weight = eval(substitute(expression({
    temp5  <- trigamma(parg + qq)
    temp5a <- trigamma(parg)
    temp5b <- trigamma(qq)

    ned2l.dscale  <- (aa^2) * parg * qq / ((1+parg+qq) * Scale^2)
    ned2l.dp      <- temp5a - temp5
    ned2l.dq      <- temp5b - temp5
    ned2l.dscalep <-  aa * qq   / (Scale*(parg+qq))
    ned2l.dscaleq <- -aa * parg / (Scale*(parg+qq))
    ned2l.dpq     <- -temp5
    wz <-
      array(c(c(w) * ned2l.dscale * dscale.deta^2,
              c(w) * ned2l.dp * dp.deta^2,
              c(w) * ned2l.dq * dq.deta^2,
              c(w) * ned2l.dscalep * dscale.deta * dp.deta,
              c(w) * ned2l.dpq * dp.deta * dq.deta,  # Switched!!
              c(w) * ned2l.dscaleq * dscale.deta * dq.deta),
            dim = c(n, M/M1, M1*(M1+1)/2))
    wz <- arwz2wz(wz, M = M, M1 = M1)
    wz
  }), list( .lscale    = lscale   ,
            .escale    = escale   ,
            .lshape2.p = lshape2.p, .lshape3.q = lshape3.q,
            .eshape2.p = eshape2.p, .eshape3.q = eshape3.q
                       ))))
}  # betaII






 lomax <-
  function(lscale    = "loglink",
           lshape3.q = "loglink",
           iscale    = NULL,
           ishape3.q = NULL,
           imethod   = 1,
           gscale    = exp(-5:5),
           gshape3.q = seq(0.75, 4, by = 0.25),  # exp(-5:5),
           probs.y   = c(0.25, 0.50, 0.75),
           zero      = "shape") {




  if (!is.Numeric(imethod, length.arg = 1,
                  integer.valued = TRUE,
                  positive = TRUE) || imethod > 2)
    stop("Bad input for argument 'imethod'")

  if (length(iscale) && !is.Numeric(iscale, positive = TRUE))
    stop("Bad input for argument 'iscale'")

  if (length(ishape3.q) && !is.Numeric(ishape3.q, positive = TRUE))
    stop("Bad input for argument 'ishape3.q'")

  if (length(probs.y) < 2 || max(probs.y) > 1 ||
        !is.Numeric(probs.y, positive = TRUE))
    stop("Bad input for argument 'probs.y'")


  lscale <- as.list(substitute(lscale))
  escale <- link2list(lscale)
  lscale <- attr(escale, "function.name")

  lshape3.q <- as.list(substitute(lshape3.q))
  eshape3.q <- link2list(lshape3.q)
  lshape3.q <- attr(eshape3.q, "function.name")


  new("vglmff",
  blurb =
    c("Lomax distribution \n\n",
      "Links:    ",
      namesof("scale"    , lscale   , earg = escale   ), ", ",
      namesof("shape3.q" , lshape3.q, earg = eshape3.q), "\n",
      "Mean:     scale / (shape3.q - 1)"),
  constraints = eval(substitute(expression({
    constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
                     M = M, M1 = 2,
                     predictors.names = predictors.names)

  }), list( .zero = zero ))),
  infos = eval(substitute(function(...) {
    list(M1 = 2,
         Q1 = 1,
         dpqrfun = "lomax",  # "genbetaII",
         expected = TRUE,
         zero = .zero ,
         multipleResponses = TRUE,
         parameters.names = c("scale", "shape3.q"),
         lscale = .lscale ,
         escale = .escale ,
                                  lshape3.q = .lshape3.q ,
                                  eshape3.q = .eshape3.q )
  }, list( .lscale = lscale      ,
           .escale = escale      ,
                                   .lshape3.q = lshape3.q,
                                   .eshape3.q = eshape3.q,
           .zero = zero ))),

  rqresslot = eval(substitute(
    function(mu, y, w, eta, extra = NULL) {
      M1  <- 2
      NOS <- ncol(eta) / M1
      Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lscale    , earg = .escale    )
      aa    <- 1
      parg  <- 1
      qq    <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lshape3.q , earg = .eshape3.q )
    scrambleseed <- runif(1)  # To scramble the seed
      qnorm(plomax(y, scale = Scale, shape3.q = qq))
    }, list( .lscale    = lscale   ,
             .escale    = escale   ,
                                     .lshape3.q = lshape3.q,
                                     .eshape3.q = eshape3.q
           ))),
  initialize = eval(substitute(expression({
    temp5 <- w.y.check(w = w, y = y,
                       Is.positive.y = TRUE,
                       ncol.w.max = Inf,
                       ncol.y.max = Inf,
                       out.wy = TRUE,
                       colsyperw = 1,
                       maximize = TRUE)
    y    <- temp5$y
    w    <- temp5$w
    M1 <- 2  # Number of parameters for one response
    NOS  <- ncoly <- ncol(y)
    M    <- M1*ncol(y)

    scaL.names  <- param.names("scale",     NOS, skip1 = TRUE)
    sha3.names  <- param.names("shape3.q",  NOS, skip1 = TRUE)

    predictors.names <-
      c(namesof(scaL.names , .lscale    , .escale    , tag = FALSE),
        namesof(sha3.names , .lshape3.q , .eshape3.q , tag = FALSE))
    predictors.names <- predictors.names[interleave.VGAM(M,
                                                         M1 = M1)]

    if (!length(etastart)) {
      sc.init <-
      qq.init <- matrix(NA_real_, n, NOS)

      for (spp. in 1:NOS) {  # For each response 'y_spp.'... do:
        yvec <- y[, spp.]
        wvec <- w[, spp.]

        if ( .imethod == 1) {
          gscale     <- .gscale
          gshape3.q  <- .gshape3.q
          if (length( .iscale    )) gscale    <-
                                        rep_len( .iscale    , NOS)
          if (length( .ishape3.q )) gshape3.q <-
                                        rep_len( .ishape3.q , NOS)



        try.this <-
          grid.search4(gscale, vov2 = 1, vov3 = 1, gshape3.q,
                       objfun = genbetaII.Loglikfun4,  # x = x,
                       y = yvec, w = wvec,  # extraargs = NULL,
                       ret.objfun = TRUE)  # Last value is \ell
          sc.init[, spp.] <- try.this["Value1"]
          qq.init[, spp.] <- try.this["Value4"]
       } else {  # .imethod == 2
          qvec <- .probs.y
          iscale <- if (length( .iscale )) .iscale else 1
          xvec <- log1p( quantile(yvec / iscale, probs = qvec) )
          fit0 <- lsfit(xvec, y = -log1p(-qvec), intercept = FALSE)
          sc.init[, spp.] <- iscale
           qq.init[, spp.] <- if (length( .ishape3.q ))
                                  .ishape3.q else
                             fit0$coef
        }
      }  # End of for (spp. ...)

      finite.mean <- 1 < qq.init
      COP.use <- 1.15
      while (FALSE && any(!finite.mean)) {
        qq.init[!finite.mean] <- 0.1 + qq.init[!finite.mean] * COP.use
        finite.mean <- 1 < qq.init
      }

      etastart <-
        cbind(theta2eta(sc.init, .lscale    , earg = .escale    ),
              theta2eta(qq.init, .lshape3.q , earg = .eshape3.q ))
      etastart <- etastart[, interleave.VGAM(M, M1 = M1)]
    }  # End of etastart.
  }), list( .lscale    = lscale   ,
            .escale    = escale   ,
            .iscale    = iscale   ,
            .gscale    = gscale   ,
                                    .lshape3.q = lshape3.q,
                                    .eshape3.q = eshape3.q,
                                    .ishape3.q = ishape3.q,
                                    .gshape3.q = gshape3.q,
            .imethod   = imethod   , .probs.y  = probs.y
                       ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    M1 <- 2
    NOS <- ncol(eta)/M1
    Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                       .lscale    , earg = .escale    )
    aa    <- 1
    parg  <- 1
    qq    <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                       .lshape3.q , earg = .eshape3.q )

    ans <- Scale * exp(lgamma(parg + 1/aa) +
                       lgamma(qq   - 1/aa) -
                       lgamma(parg) - lgamma(qq))
    ans[parg + 1/aa <= 0] <- NA
    ans[qq   - 1/aa <= 0] <- NA
    ans[Scale       <= 0] <- NA
    ans[qq          <= 0] <- NA
    ans
  }, list( .lscale    = lscale   ,
           .escale    = escale   ,
                                   .lshape3.q = lshape3.q,
                                   .eshape3.q = eshape3.q
                      ))),
  last = eval(substitute(expression({
    M1 <- 2

    misc$link <-
      c(rep_len( .lscale    , ncoly),
        rep_len( .lshape3.q , ncoly))[interleave.VGAM(M, M1 = M1)]
    temp.names <-
      c(scaL.names,                         sha3.names)
    names(misc$link) <- temp.names[interleave.VGAM(M, M1 = M1)]

    misc$earg <- vector("list", M)
    names(misc$earg) <- temp.names[interleave.VGAM(M, M1 = M1)]
    for (ii in 1:ncoly) {
      misc$earg[[M1*ii-1]] <- .escale
      misc$earg[[M1*ii  ]] <- .eshape3.q
    }

    misc$expected          <- TRUE
    misc$multipleResponses <- TRUE
  }), list( .lscale    = lscale   ,
            .escale    = escale   ,
                                    .lshape3.q = lshape3.q,
                                    .eshape3.q = eshape3.q
                       ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE,
             eta, extra = NULL, summation = TRUE) {
      M1  <- 2
      NOS <- ncol(eta) / M1
      Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lscale    , earg = .escale    )
      aa    <- 1
      parg  <- 1
      qq    <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lshape3.q , earg = .eshape3.q )
      if (residuals) {
        stop("loglikelihood residuals not implemented yet")
      } else {
        ll.elts <-
          c(w) * dgenbetaII(x = y, scale = Scale, shape1.a = aa,
                            shape2.p = parg, shape3.q = qq,
                            log = TRUE)
        if (summation) {
          sum(ll.elts)
        } else {
          ll.elts
        }
      }
    }, list( .lscale    = lscale   ,
             .escale    = escale   ,
                                     .lshape3.q = lshape3.q,
                                     .eshape3.q = eshape3.q
                        ))),
  vfamily = c("lomax"),
  simslot = eval(substitute(
  function(object, nsim) {

    pwts <- if (length(pwts <- object@prior.weights) > 0)
              pwts else weights(object, type = "prior")
    if (any(pwts != 1))
      warning("ignoring prior weights")

    eta <- predict(object)
    M1 <- 2
    NOS <- ncol(eta)/M1
    Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                       .lscale ,  earg = .escale )
    qq    <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                       .lshape3.q , earg = .eshape3.q )

    rlomax(nsim * length(qq), scale = Scale, shape3.q = qq)
  }, list( .lscale    = lscale   ,
           .escale    = escale   ,
                                   .lshape3.q = lshape3.q,
                                   .eshape3.q = eshape3.q
                      ))),

  validparams = eval(substitute(function(eta, y, extra = NULL) {
    M1 <- 2
    NOS <- ncol(eta) / M1
    Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                       .lscale    , earg = .escale   )
    aa    <- 1
    parg  <- 1
    qq    <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                       .lshape3.q , earg = .eshape3.q)
    okay1 <- all(is.finite(Scale)) && all(Scale  > 0) &&
             all(is.finite(qq   )) && all(qq     > 0)
    okay.support <- if (okay1)
                    all(-aa * parg < 1 & 1 < aa * qq) else TRUE
    if (!okay.support)
      warning("parameter constraint 1 < q has",
              " been violated; solution may be at",
              " the boundary of the parameter space.")
    okay1 && okay.support
  }, list(  .lscale    = lscale   ,
            .escale    = escale   ,
                                    .lshape3.q = lshape3.q,
                                    .eshape3.q = eshape3.q ))),




  deriv = eval(substitute(expression({
    M1 <- 2
    NOS <- ncol(eta)/M1  # Needed for summary()
    Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                       .lscale    , earg = .escale   )
    aa    <- 1
    parg  <- 1
    qq    <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                       .lshape3.q , earg = .eshape3.q)
    temp1 <- log(y/Scale)
    temp2 <- (y/Scale)^aa
    temp3 <- digamma(parg + qq)
    temp3a <- digamma(parg)
    temp3b <- digamma(qq)
    temp4 <- log1p(temp2)

    dl.dscale <- (aa/Scale) * (-parg + (parg+qq) / (1+1/temp2))
    dl.dq <- temp3 - temp3b - temp4

    dscale.deta <- dtheta.deta(Scale, .lscale ,    earg = .escale )
    dq.deta     <- dtheta.deta(qq,    .lshape3.q , earg = .eshape3.q )

    myderiv <-
      c(w) * cbind(dl.dscale * dscale.deta,
                   dl.dq * dq.deta)
    myderiv[, interleave.VGAM(M, M1 = M1)]
  }), list(  .lscale    = lscale   ,
             .escale    = escale   ,
                                     .lshape3.q = lshape3.q,
                                     .eshape3.q = eshape3.q
                        ))),
  weight = eval(substitute(expression({
    temp5  <- trigamma(parg + qq)
    temp5a <- trigamma(parg)
    temp5b <- trigamma(qq)

    ned2l.dscale  <- (aa^2) * parg * qq / ((1+parg+qq) * Scale^2)
    ned2l.dq      <- temp5b - temp5
    ned2l.dscaleq <- -aa * parg / (Scale*(parg+qq))
    wz <-
      array(c(c(w) * ned2l.dscale * dscale.deta^2,
              c(w) * ned2l.dq * dq.deta^2,
              c(w) * ned2l.dscaleq * dscale.deta * dq.deta),
            dim = c(n, M/M1, M1*(M1+1)/2))
    wz <- arwz2wz(wz, M = M, M1 = M1)
    wz
  }), list( .lscale    = lscale   ,
            .escale    = escale   ,
                                    .lshape3.q = lshape3.q,
                                    .eshape3.q = eshape3.q
                       ))))
}  # lomax







 fisk <-
  function(lscale    = "loglink",
           lshape1.a = "loglink",
           iscale    = NULL,
           ishape1.a = NULL,
           imethod   = 1,
           lss       = TRUE,
           gscale    = exp(-5:5),
           gshape1.a = seq(0.75, 4, by = 0.25),  # exp(-5:5),
           probs.y   = c(0.25, 0.50, 0.75),
           zero      = "shape") {





  if (length(lss) != 1 && !is.logical(lss))
    stop("Argument 'lss' not specified correctly")

  if (!is.Numeric(imethod, length.arg = 1,
                  integer.valued = TRUE,
                  positive = TRUE) || imethod > 2)
    stop("Bad input for argument 'imethod'")

  if (length(iscale) && !is.Numeric(iscale, positive = TRUE))
    stop("Bad input for argument 'iscale'")

  if (length(ishape1.a) && !is.Numeric(ishape1.a, positive = TRUE))
    stop("Bad input for argument 'ishape1.a'")

  if (length(probs.y) < 2 || max(probs.y) > 1 ||
        !is.Numeric(probs.y, positive = TRUE))
    stop("Bad input for argument 'probs.y'")


  lscale <- as.list(substitute(lscale))
  escale <- link2list(lscale)
  lscale <- attr(escale, "function.name")

  lshape1.a <- as.list(substitute(lshape1.a))
  eshape1.a <- link2list(lshape1.a)
  lshape1.a <- attr(eshape1.a, "function.name")


  new("vglmff",
  blurb =
    c("Fisk distribution \n\n",
      "Links:    ",
      ifelse (lss,
              namesof("scale"   , lscale   , escale),
              namesof("shape1.a", lshape1.a, eshape1.a)), ", ",
      ifelse (lss,
              namesof("shape1.a", lshape1.a, eshape1.a),
              namesof("scale"   , lscale   , escale)), "\n",
      "Mean:     scale * gamma(1 + 1/shape1.a) * ",
                "gamma(1 - 1/shape1.a)"),
  constraints = eval(substitute(expression({
    constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
                     M = M, M1 = 2,
                     predictors.names = predictors.names)

  }), list( .zero = zero ))),
  infos = eval(substitute(function(...) {
    list(M1 = 2,
         Q1 = 1,
         dpqrfun = "fisk",  # "genbetaII",
         expected = TRUE,
         zero = .zero ,
         multipleResponses = TRUE,
         parameters.names = if ( .lss )
           c("scale", "shape1.a") else
           c("shape1.a", "scale"),
         lscale = .lscale ,       lshape1.a = .lshape1.a ,
         escale = .escale ,       eshape1.a = .eshape1.a )
  }, list( .lscale = lscale      , .lshape1.a = lshape1.a,
           .escale = escale      , .eshape1.a = eshape1.a,
           .lss  = lss ,
           .zero = zero ))),

  rqresslot = eval(substitute(
    function(mu, y, w, eta, extra = NULL) {
      M1  <- 2
      NOS <- ncol(eta) / M1
    if ( .lss ) {
      Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lscale ,  earg = .escale )
      aa    <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lshape1.a , earg = .eshape1.a )
    } else {
      aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a )
      Scale <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lscale    ,  earg = .escale )
    }
    parg <- 1
    qq   <- 1
    scrambleseed <- runif(1)  # To scramble the seed
    qnorm(pfisk(y, scale = Scale, shape1.a = aa))
    }, list( .lscale    = lscale   , .lshape1.a = lshape1.a,
             .escale    = escale   , .eshape1.a = eshape1.a,
             .lss = lss ))),
  initialize = eval(substitute(expression({
    temp5 <- w.y.check(w = w, y = y,
                       Is.positive.y = TRUE,
                       ncol.w.max = Inf,
                       ncol.y.max = Inf,
                       out.wy = TRUE,
                       colsyperw = 1,
                       maximize = TRUE)
    y    <- temp5$y
    w    <- temp5$w
    M1 <- 2  # Number of parameters for one response
    NOS  <- ncoly <- ncol(y)
    M    <- M1*ncol(y)

    scaL.names  <- param.names("scale",     NOS, skip1 = TRUE)
    sha1.names  <- param.names("shape1.a",  NOS, skip1 = TRUE)

    predictors.names <- if ( .lss ) {
      c(namesof(scaL.names , .lscale    , .escale    , tag = FALSE),
        namesof(sha1.names , .lshape1.a , .eshape1.a , tag = FALSE))
    } else {
      c(namesof(sha1.names , .lshape1.a , .eshape1.a , tag = FALSE),
        namesof(scaL.names , .lscale    , .escale    , tag = FALSE))
    }
    predictors.names <- predictors.names[interleave.VGAM(M,
                                                         M1 = M1)]

    if (!length(etastart)) {
      sc.init <-
      aa.init <- matrix(NA_real_, n, NOS)

      for (spp. in 1:NOS) {  # For each response 'y_spp.'... do:
        yvec <- y[, spp.]
        wvec <- w[, spp.]

        if ( .imethod == 1 ) {
          gscale     <- .gscale
          gshape1.a  <- .gshape1.a
          if (length( .iscale    )) gscale    <-
                                        rep_len( .iscale    , NOS)
          if (length( .ishape1.a )) gshape1.a <-
                                        rep_len( .ishape1.a , NOS)




        try.this <-
          grid.search4(gscale, gshape1.a, vov3 = 1, vov4 = 1,
                       objfun = genbetaII.Loglikfun4,  # x = x,
                       y = yvec, w = wvec,  # extraargs = NULL,
                       ret.objfun = TRUE)  # Last value is \ell
          sc.init[, spp.] <- try.this["Value1"]
          aa.init[, spp.] <- try.this["Value2"]
       } else {  # .imethod == 2
          qvec <- .probs.y
          iscale <- if (length( .iscale )) .iscale else 1
          xvec <- log( quantile(yvec / iscale, probs = qvec) )
           fit0 <- lsfit(x = xvec, y = logitlink(qvec),
                         intercept = FALSE)
          sc.init[, spp.] <- iscale
           aa.init[, spp.] <- if (length( .ishape1.a ))
                                  .ishape1.a else
                             fit0$coef
        }
      }  # End of for (spp. ...)


      finite.mean <- 1 < aa.init
      COP.use <- 1.15
      while (FALSE && any(!finite.mean)) {
          aa.init[!finite.mean] <- 0.1 +
              aa.init[!finite.mean] * COP.use
        finite.mean <- 1 < aa.init
      }

      etastart <-
        if ( .lss )
          cbind(theta2eta(sc.init, .lscale    , .escale    ),
                theta2eta(aa.init, .lshape1.a , .eshape1.a )) else
          cbind(theta2eta(aa.init, .lshape1.a , .eshape1.a ),
                theta2eta(sc.init, .lscale    , .escale    ))
      etastart <- etastart[, interleave.VGAM(M, M1 = M1)]
    }  # End of etastart.
  }), list( .lscale    = lscale   , .lshape1.a = lshape1.a,
            .escale    = escale   , .eshape1.a = eshape1.a,
            .iscale    = iscale   , .ishape1.a = ishape1.a,
            .gscale    = gscale   , .gshape1.a = gshape1.a,
            .imethod   = imethod   , .probs.y  = probs.y,
            .lss = lss ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    M1 <- 2
    NOS <- ncol(eta)/M1
    if ( .lss ) {
      Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lscale ,  earg = .escale )
      aa    <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lshape1.a , earg = .eshape1.a )
    } else {
      aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a )
      Scale <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lscale    ,  earg = .escale )
    }
    parg <- 1
    qq   <- 1

    ans <- Scale * exp(lgamma(parg + 1/aa) +
                       lgamma(qq   - 1/aa) -
                       lgamma(parg) - lgamma(qq))
    ans[parg + 1/aa <= 0] <- NA
    ans[qq   - 1/aa <= 0] <- NA
    ans[Scale       <= 0] <- NA
    ans[aa          <= 0] <- NA
    ans
  }, list( .lscale    = lscale   , .lshape1.a = lshape1.a,
           .escale    = escale   , .eshape1.a = eshape1.a,
           .lss = lss ))),
  last = eval(substitute(expression({
    M1 <- 2

    misc$link <- c(rep_len( if ( .lss ) .lscale else .lshape1.a ,
                           ncoly),
                   rep_len( if ( .lss ) .lshape1.a else .lscale ,
                           ncoly))[
                   interleave.VGAM(M, M1 = M1)]
    temp.names <- if ( .lss ) {
      c(scaL.names, sha1.names)
    } else {
      c(sha1.names, scaL.names)
    }
    names(misc$link) <- temp.names[interleave.VGAM(M, M1 = M1)]

    misc$earg <- vector("list", M)
    names(misc$earg) <- temp.names[interleave.VGAM(M, M1 = M1)]
    for (ii in 1:ncoly)
      if ( .lss ) {
        misc$earg[[M1*ii-1]] <- .escale
        misc$earg[[M1*ii  ]] <- .eshape1.a
      } else {
        misc$earg[[M1*ii-1]] <- .eshape1.a
        misc$earg[[M1*ii  ]] <- .escale
      }
  }), list( .lscale    = lscale   , .lshape1.a = lshape1.a,
            .escale    = escale   , .eshape1.a = eshape1.a,
            .lss = lss ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE,
             eta, extra = NULL, summation = TRUE) {
      M1  <- 2
      NOS <- ncol(eta)/M1
      if ( .lss ) {
        Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                           .lscale    , earg = .escale    )
        aa    <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                           .lshape1.a , earg = .eshape1.a )
      } else {
        aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                           .lshape1.a , earg = .eshape1.a )
        Scale <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                           .lscale    , earg = .escale    )
      }
      if (residuals) {
        stop("loglikelihood residuals not implemented yet")
      } else {
        ll.elts <-
            c(w) * dfisk(x = y, scale = Scale, shape1.a = aa,
                         log = TRUE)
        if (summation) {
          sum(ll.elts)
        } else {
          ll.elts
        }
      }
    }, list( .lscale    = lscale   , .lshape1.a = lshape1.a,
             .escale    = escale   , .eshape1.a = eshape1.a,
             .lss = lss ))),
  vfamily = c("fisk"),
  simslot = eval(substitute(
  function(object, nsim) {

    pwts <- if (length(pwts <- object@prior.weights) > 0)
              pwts else weights(object, type = "prior")
    if (any(pwts != 1))
      warning("ignoring prior weights")

    eta <- predict(object)
    M1 <- 2
    NOS <- ncol(eta)/M1
    if ( .lss ) {
      Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lscale ,  earg = .escale )
      aa    <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lshape1.a , earg = .eshape1.a )
    } else {
      aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a )
      Scale <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lscale    ,  earg = .escale )
    }
    rfisk(nsim * length(aa), shape1.a = aa, scale = Scale)
  }, list( .lscale    = lscale   , .lshape1.a = lshape1.a,
           .escale    = escale   , .eshape1.a = eshape1.a,
           .lss = lss ))),

  validparams = eval(substitute(function(eta, y, extra = NULL) {
    M1 <- 2
    NOS <- ncol(eta) / M1
    if ( .lss ) {
      Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lscale    , earg = .escale   )
      aa    <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lshape1.a , earg = .eshape1.a)
    } else {
      aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a)
      Scale <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lscale    , earg = .escale )
    }
    parg <- 1
    qq   <- 1

    okay1 <- all(is.finite(Scale)) && all(Scale  > 0) &&
             all(is.finite(aa   )) && all(aa     > 0) &&
             all(is.finite(parg )) && all(parg   > 0) &&
             all(is.finite(qq   )) && all(qq     > 0)
    okay.support <- if (okay1)
                    all(-aa * parg < 1 & 1 < aa * qq) else TRUE
    if (!okay.support)
      warning("parameter constraint -a   < 1 < a   has",
              " been violated; solution may be at",
              " the boundary of the parameter space.")
    okay1 && okay.support
  }, list(  .lscale    = lscale   , .lshape1.a = lshape1.a,
            .escale    = escale   , .eshape1.a = eshape1.a,
            .lss = lss ))),




  deriv = eval(substitute(expression({
    M1 <- 2
    NOS <- ncol(eta)/M1  # Needed for summary()
    if ( .lss ) {
      Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lscale    , earg = .escale   )
      aa    <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lshape1.a , earg = .eshape1.a)
    } else {
      aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a)
      Scale <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lscale    , earg = .escale )
    }
    parg <- 1
    qq   <- 1
    temp1 <- log(y/Scale)
    temp2 <- (y/Scale)^aa

    dl.dscale <- (aa/Scale) * (-parg + (parg+qq) / (1+1/temp2))
    dl.da <- 1/aa + parg * temp1 - (parg+qq) * temp1 / (1+1/temp2)

    dscale.deta <- dtheta.deta(Scale, .lscale ,    earg = .escale )
    da.deta     <- dtheta.deta(aa,    .lshape1.a , earg = .eshape1.a )

    myderiv <- if ( .lss ) {
      c(w) * cbind(dl.dscale * dscale.deta,
                   dl.da * da.deta)
    } else {
      c(w) * cbind(dl.da * da.deta,
                   dl.dscale * dscale.deta)
    }
    myderiv[, interleave.VGAM(M, M1 = M1)]
  }), list(  .lscale    = lscale   , .lshape1.a = lshape1.a,
             .escale    = escale   , .eshape1.a = eshape1.a,
             .lss = lss ))),
  weight = eval(substitute(expression({
    ned2l.da <- (1 + (-6 + pi^2) / 9) / aa^2
    ned2l.dscale  <- ((aa / Scale)^2) / 3


    wz <- matrix(0, n, M)  # Diagonal
    if ( .lss ) {
      wz[, c(TRUE, FALSE)] <- ned2l.dscale * dscale.deta^2
      wz[, c(FALSE, TRUE)] <- ned2l.da * da.deta^2
    } else {
      wz[, c(TRUE, FALSE)] <- ned2l.da * da.deta^2
      wz[, c(FALSE, TRUE)] <- ned2l.dscale * dscale.deta^2
    }
    c(w) * wz
  }), list( .lscale    = lscale   , .lshape1.a = lshape1.a,
            .escale    = escale   , .eshape1.a = eshape1.a,
            .lss = lss ))))
}  # fisk






 inv.lomax <- function(lscale    = "loglink",
                       lshape2.p = "loglink",
                       iscale    = NULL,
                       ishape2.p = NULL,
                       imethod   = 1,
                       gscale    = exp(-5:5),
                       gshape2.p = exp(-5:5),
                       probs.y   = c(0.25, 0.50, 0.75),
                       zero      = "shape2.p") {






  if (!is.Numeric(imethod, length.arg = 1,
                  integer.valued = TRUE,
                  positive = TRUE) || imethod > 2)
    stop("Bad input for argument 'imethod'")

  if (length(iscale) && !is.Numeric(iscale, positive = TRUE))
    stop("Bad input for argument 'iscale'")

  if (length(ishape2.p) && !is.Numeric(ishape2.p, positive = TRUE))
    stop("Bad input for argument 'ishape2.p'")

  if (length(probs.y) < 2 || max(probs.y) > 1 ||
        !is.Numeric(probs.y, positive = TRUE))
    stop("Bad input for argument 'probs.y'")


  lscale <- as.list(substitute(lscale))
  escale <- link2list(lscale)
  lscale <- attr(escale, "function.name")

  lshape2.p <- as.list(substitute(lshape2.p))
  eshape2.p <- link2list(lshape2.p)
  lshape2.p <- attr(eshape2.p, "function.name")


  new("vglmff",
  blurb =
    c("Inverse Lomax distribution \n\n",
      "Links:    ",
      namesof("scale"    , lscale   , earg = escale),    ", ",
      namesof("shape2.p" , lshape2.p, earg = eshape2.p), "\n",
      "Mean:     does not exist"),
  constraints = eval(substitute(expression({
    constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
                     M = M, M1 = 2,
                     predictors.names = predictors.names)

  }), list( .zero = zero ))),
  infos = eval(substitute(function(...) {
    list(M1 = 2,
         Q1 = 1,
         dpqrfun = "inv.lomax",  # "genbetaII",
         expected = TRUE,
         zero = .zero ,
         multipleResponses = TRUE,
         parameters.names = c("scale", "shape2.p"),
         lscale    = .lscale    ,
         escale    = .escale    ,
         lshape2.p = .lshape2.p ,
         eshape2.p = .eshape2.p )
  }, list( .lscale = lscale      ,
           .escale = escale      ,
           .lshape2.p = lshape2.p,
           .eshape2.p = eshape2.p,
           .zero = zero ))),

  rqresslot = eval(substitute(
    function(mu, y, w, eta, extra = NULL) {
    M1 <- 2
    NOS <- ncol(eta) / M1
    Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                       .lscale    , earg = .escale    )
    parg  <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                       .lshape2.p , earg = .eshape2.p )

    scrambleseed <- runif(1)  # To scramble the seed
    qnorm(pinv.lomax(y, scale = Scale, shape2.p = parg))
    }, list( .lscale  = lscale , .lshape2.p = lshape2.p,
             .escale  = escale , .eshape2.p = eshape2.p ))),
  initialize = eval(substitute(expression({
    temp5 <- w.y.check(w = w, y = y,
                       Is.positive.y = TRUE,
                       ncol.w.max = Inf,
                       ncol.y.max = Inf,
                       out.wy = TRUE,
                       colsyperw = 1,
                       maximize = TRUE)
    y    <- temp5$y
    w    <- temp5$w
    M1 <- 2  # Number of parameters for one response
    NOS  <- ncoly <- ncol(y)
    M    <- M1*ncol(y)

    scaL.names  <- param.names("scale",     NOS, skip1 = TRUE)
    sha2.names  <- param.names("shape2.p",  NOS, skip1 = TRUE)

    predictors.names <-
    c(namesof(scaL.names , .lscale    , .escale    , tag = FALSE),
      namesof(sha2.names , .lshape2.p , .eshape2.p , tag = FALSE))
    predictors.names <- predictors.names[interleave.VGAM(M,
                                                         M1 = M1)]

    if (!length(etastart)) {
      sc.init <-
      pp.init <- matrix(NA_real_, n, NOS)

      for (spp. in 1:NOS) {  # For each response 'y_spp.'... do:
        yvec <- y[, spp.]
        wvec <- w[, spp.]

        if ( .imethod == 1 ) {
          gscale     <- .gscale
          gshape2.p  <- .gshape2.p
          if (length( .iscale    )) gscale    <-
                                        rep_len( .iscale    , NOS)
          if (length( .ishape2.p )) gshape2.p <-
                                        rep_len( .ishape2.p , NOS)





        try.this <-
          grid.search4(gscale, vov2 = 1, gshape2.p, vov4 = 1,
                       objfun = genbetaII.Loglikfun4,  # x = x,
                       y = yvec, w = wvec,  # extraargs = NULL,
                       ret.objfun = TRUE)  # Last value is \ell
          sc.init[, spp.] <- try.this["Value1"]
          pp.init[, spp.] <- try.this["Value3"]
       } else {  # .imethod == 2
          qvec <- .probs.y
          ishape2.p <- if (length( .ishape2.p )) .ishape2.p else 1
          xvec <- log( qvec^(-1/ ishape2.p) - 1 )
          fit0 <- lsfit(x = xvec, y = log(quantile(yvec,  qvec)))
          sc.init[, spp.] <- if (length( .iscale )) .iscale else
                             exp(fit0$coef[1])
          pp.init[, spp.] <- ishape2.p
      }
    }  # End of for (spp. ...)


      etastart <-
        cbind(theta2eta(sc.init, .lscale    , earg = .escale    ),
              theta2eta(pp.init, .lshape2.p , earg = .eshape2.p ))
      etastart <- etastart[, interleave.VGAM(M, M1 = M1)]
    }  # End of etastart.
  }), list( .lscale    = lscale   ,
            .escale    = escale   ,
            .iscale    = iscale   ,
            .gscale    = gscale   ,
            .lshape2.p = lshape2.p,
            .eshape2.p = eshape2.p,
            .ishape2.p = ishape2.p,
            .gshape2.p = gshape2.p,
            .imethod   = imethod   , .probs.y = probs.y
                       ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {

    M1 <- 2
    NOS <- ncol(eta) / M1
    Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                       .lscale    , earg = .escale    )
    parg  <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                       .lshape2.p , earg = .eshape2.p )

    qinv.lomax(p = 0.5, scale = Scale, shape2.p = parg)
  }, list( .lscale    = lscale   ,
           .escale    = escale   ,
           .lshape2.p = lshape2.p,
           .eshape2.p = eshape2.p
                      ))),
  last = eval(substitute(expression({
    M1 <- 2

    misc$link <-
      c(rep_len( .lscale    , ncoly),
        rep_len( .lshape2.p , ncoly))[interleave.VGAM(M, M1 = M1)]
    temp.names <- c(scaL.names, sha2.names)
    names(misc$link) <- temp.names[interleave.VGAM(M, M1 = M1)]

    misc$earg <- vector("list", M)
    names(misc$earg) <- temp.names[interleave.VGAM(M, M1 = M1)]
    for (ii in 1:ncoly) {
      misc$earg[[M1*ii-1]] <- .escale
      misc$earg[[M1*ii  ]] <- .eshape2.p
    }

    misc$expected          <- TRUE
    misc$multipleResponses <- TRUE
  }), list( .lscale    = lscale   ,
            .escale    = escale   ,
            .lshape2.p = lshape2.p,
            .eshape2.p = eshape2.p
                       ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE,
             eta, extra = NULL, summation = TRUE) {
      M1  <- 2
      NOS <- ncol(eta)/M1
      Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lscale    , earg = .escale    )
      parg  <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lshape2.p , earg = .eshape2.p )
      if (residuals) {
        stop("loglikelihood residuals not implemented yet")
      } else {
        ll.elts <-
          c(w) * dgenbetaII(x = y, scale = Scale, shape1.a = 1,
                            shape2.p = parg, shape3.q = 1,
                            log = TRUE)
        if (summation) {
          sum(ll.elts)
        } else {
          ll.elts
        }
      }
    }, list( .lscale    = lscale   ,
             .escale    = escale   ,
             .lshape2.p = lshape2.p,
             .eshape2.p = eshape2.p
                        ))),
  vfamily = c("inv.lomax"),

  simslot = eval(substitute(
  function(object, nsim) {

    pwts <- if (length(pwts <- object@prior.weights) > 0)
              pwts else weights(object, type = "prior")
    if (any(pwts != 1))
      warning("ignoring prior weights")
    eta <- predict(object)
    Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                       .lscale    , earg = .escale    )
    parg  <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                       .lshape2.p , earg = .eshape2.p )
    aa   <- 1
    qq   <- 1
    rinv.lomax(nsim * length(Scale), scale = Scale, shape2.p = parg)
  }, list(  .lscale    = lscale   ,
            .escale    = escale   ,
            .lshape2.p = lshape2.p,
            .eshape2.p = eshape2.p
                       ))),



  validparams = eval(substitute(function(eta, y, extra = NULL) {
    M1 <- 2
    NOS <- ncol(eta) / M1
    Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                       .lscale    , earg = .escale    )
    parg  <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                       .lshape2.p , earg = .eshape2.p )
    qq   <- 1
    aa   <- 1
    okay1 <- all(is.finite(Scale)) && all(Scale  > 0) &&
             all(is.finite(aa   )) && all(aa     > 0) &&
             all(is.finite(parg )) && all(parg   > 0) &&
             all(is.finite(qq   )) && all(qq     > 0)
      okay.support <- if (okay1)
                      all(-aa * parg < 1              ) else TRUE
    if (!okay.support)
      warning("parameter constraint -a*p < 1 has",
              " been violated; solution may be at",
              " the boundary of the parameter space.")
    okay1 && okay.support
  }, list(  .lscale    = lscale   ,
            .escale    = escale   ,
            .lshape2.p = lshape2.p,
            .eshape2.p = eshape2.p
         ))),




  deriv = eval(substitute(expression({
    M1 <- 2
    NOS <- ncol(eta)/M1  # Needed for summary()
    Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                       .lscale    , earg = .escale    )
    parg  <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                       .lshape2.p , earg = .eshape2.p )
    qq   <- 1
    aa   <- 1
    temp1 <- log(y/Scale)
    temp2 <- (y/Scale)^aa
    temp3 <- digamma(parg + qq)
    temp3a <- digamma(parg)
    temp3b <- digamma(qq)
    temp4 <- log1p(temp2)

    dl.dscale <- (aa/Scale) * (-parg + (parg+qq) / (1+1/temp2))
    dl.dp <- aa * temp1 + temp3 - temp3a - temp4

    dscale.deta <- dtheta.deta(Scale, .lscale ,    earg = .escale )
    dp.deta     <- dtheta.deta(parg,  .lshape2.p , earg = .eshape2.p )

    myderiv <-
      c(w) * cbind(dl.dscale * dscale.deta,
                   dl.dp * dp.deta)
    myderiv[, interleave.VGAM(M, M1 = M1)]
  }), list(  .lscale    = lscale   ,
             .escale    = escale   ,
             .lshape2.p = lshape2.p,
             .eshape2.p = eshape2.p
                        ))),
  weight = eval(substitute(expression({
    temp5  <- trigamma(parg + qq)
    temp5a <- trigamma(parg)
    temp5b <- trigamma(qq)

    ned2l.dscale  <- (aa^2) * parg * qq / ((1+parg+qq) * Scale^2)
    ned2l.dp      <- temp5a - temp5
    ned2l.dscalep <-  aa * qq   / (Scale*(parg+qq))
    wz <-
      array(c(c(w) * ned2l.dscale * dscale.deta^2,
              c(w) * ned2l.dp * dp.deta^2,
              c(w) * ned2l.dscalep * dscale.deta * dp.deta),
            dim = c(n, M/M1, M1*(M1+1)/2))
    wz <- arwz2wz(wz, M = M, M1 = M1)
    wz
  }), list( .lscale    = lscale   ,
            .escale    = escale   ,
            .lshape2.p = lshape2.p,
            .eshape2.p = eshape2.p
                       ))))
}  # inv.lomax






 paralogistic <-
  function(lscale    = "loglink",
           lshape1.a = "loglink",
           iscale    = NULL,
           ishape1.a = NULL,
           imethod   = 1,
           lss       = TRUE,
           gscale    = exp(-5:5),
           gshape1.a = seq(0.75, 4, by = 0.25),  # exp(-5:5),
           probs.y   = c(0.25, 0.50, 0.75),
           zero      = "shape") {





  if (length(lss) != 1 && !is.logical(lss))
    stop("Argument 'lss' not specified correctly")

  if (!is.Numeric(imethod, length.arg = 1,
                  integer.valued = TRUE,
                  positive = TRUE) || imethod > 2)
    stop("Bad input for argument 'imethod'")

  if (length(iscale) && !is.Numeric(iscale, positive = TRUE))
    stop("Bad input for argument 'iscale'")

  if (length(ishape1.a) && !is.Numeric(ishape1.a, positive = TRUE))
    stop("Bad input for argument 'ishape1.a'")

  if (length(probs.y) < 2 || max(probs.y) > 1 ||
        !is.Numeric(probs.y, positive = TRUE))
    stop("Bad input for argument 'probs.y'")


  lscale <- as.list(substitute(lscale))
  escale <- link2list(lscale)
  lscale <- attr(escale, "function.name")

  lshape1.a <- as.list(substitute(lshape1.a))
  eshape1.a <- link2list(lshape1.a)
  lshape1.a <- attr(eshape1.a, "function.name")


  new("vglmff",
  blurb =
    c("Paralogistic distribution \n\n",
      "Links:    ",
      ifelse (lss,
              namesof("scale"   , lscale   , escale),
              namesof("shape1.a", lshape1.a, eshape1.a)), ", ",
      ifelse (lss,
              namesof("shape1.a", lshape1.a, eshape1.a),
              namesof("scale"   , lscale   , escale)), "\n",
      "Mean:     scale * gamma(1 + 1/shape1.a) * ",
                "gamma(shape1.a - 1/shape1.a) / ",
                "gamma(shape1.a)"),
  constraints = eval(substitute(expression({
    constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
                     M = M, M1 = 2,
                     predictors.names = predictors.names)

  }), list( .zero = zero ))),
  infos = eval(substitute(function(...) {
    list(M1 = 2,
         Q1 = 1,
         dpqrfun = "paralogistic",  # "genbetaII",
         expected = TRUE,
         zero = .zero ,
         multipleResponses = TRUE,
         parameters.names = if ( .lss )
           c("scale", "shape1.a") else
           c("shape1.a", "scale"),
         lscale = .lscale ,       lshape1.a = .lshape1.a ,
         escale = .escale ,       eshape1.a = .eshape1.a )
  }, list( .lscale = lscale      , .lshape1.a = lshape1.a,
           .escale = escale      , .eshape1.a = eshape1.a,
           .lss  = lss ,
           .zero = zero ))),

  rqresslot = eval(substitute(
    function(mu, y, w, eta, extra = NULL) {
      M1  <- 2
      NOS <- ncol(eta) / M1
    if ( .lss ) {
      Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lscale ,  earg = .escale )
      aa    <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lshape1.a , earg = .eshape1.a )
    } else {
      aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a )
      Scale <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lscale    ,  earg = .escale )
    }
    scrambleseed <- runif(1)  # To scramble the seed
    qnorm(pparalogistic(y, scale = Scale, shape1.a = aa))
    }, list( .lscale    = lscale   , .lshape1.a = lshape1.a,
             .escale    = escale   , .eshape1.a = eshape1.a,
             .lss = lss))),
  initialize = eval(substitute(expression({
    temp5 <- w.y.check(w = w, y = y,
                       Is.positive.y = TRUE,
                       ncol.w.max = Inf,
                       ncol.y.max = Inf,
                       out.wy = TRUE,
                       colsyperw = 1,
                       maximize = TRUE)
    y    <- temp5$y
    w    <- temp5$w
    M1 <- 2  # Number of parameters for one response
    NOS  <- ncoly <- ncol(y)
    M    <- M1*ncol(y)

    scaL.names  <- param.names("scale",     NOS, skip1 = TRUE)
    sha1.names  <- param.names("shape1.a",  NOS, skip1 = TRUE)

    predictors.names <-
      if ( .lss ) {
        c(namesof(scaL.names , .lscale    , earg = .escale    ,
                  tag = FALSE),
          namesof(sha1.names , .lshape1.a , earg = .eshape1.a ,
                  tag = FALSE))
      } else {
        c(namesof(sha1.names , .lshape1.a , earg = .eshape1.a ,
                  tag = FALSE),
          namesof(scaL.names , .lscale    , earg = .escale    ,
                  tag = FALSE))
      }
    predictors.names <- predictors.names[interleave.VGAM(M,
                                                         M1 = M1)]

    if (!length(etastart)) {
      sc.init <-
      aa.init <- matrix(NA_real_, n, NOS)

      for (spp. in 1:NOS) {  # For each response 'y_spp.'... do:
        yvec <- y[, spp.]
        wvec <- w[, spp.]

        if ( .imethod == 1 ) {
          gscale     <- .gscale
          gshape1.a  <- .gshape1.a
          if (length( .iscale    )) gscale    <-
                                        rep_len( .iscale    , NOS)
          if (length( .ishape1.a )) gshape1.a <-
                                        rep_len( .ishape1.a , NOS)


        paralogistic.Loglikfun2 <-
            function(scaleval, shape1.a,
                     y, x, w, extraargs) {
          sum(c(w) * dgenbetaII(x = y, scale = scaleval,
                                shape1.a = shape1.a,
                                shape2.p = 1,
                                shape3.q = shape1.a, log = TRUE))
        }

        try.this <-
          grid.search2(gscale, gshape1.a, # vov3 = 1, vov4 = gshape1.a,
                       objfun = paralogistic.Loglikfun2,  # x = x,
                       y = yvec, w = wvec,  # extraargs = NULL,
                       ret.objfun = TRUE)  # Last value is \ell
          sc.init[, spp.] <- try.this["Value1"]
          aa.init[, spp.] <- try.this["Value2"]
       } else {  # .imethod == 2
          qvec <- .probs.y
          ishape3.q <- if (length( .ishape1.a )) .ishape1.a else 1
          xvec <- log( (1-qvec)^(-1/ ishape3.q) - 1 )
          fit0 <- lsfit(x = xvec, y = log(quantile(yvec,  qvec)))
          sc.init[, spp.] <- if (length( .iscale )) .iscale else
                             exp(fit0$coef[1])
           aa.init[, spp.] <- if (length( .ishape1.a ))
                                  .ishape1.a else
                             1/fit0$coef[2]
        }
      }  # End of for (spp. ...)

      finite.mean <- 1 < aa.init * aa.init
      COP.use <- 1.15
      while (FALSE && any(!finite.mean)) {
    aa.init[!finite.mean] <- 0.1 + aa.init[!finite.mean] * COP.use
    finite.mean <- 1 < aa.init * aa.init
      }

      etastart <- if ( .lss )
          cbind(theta2eta(sc.init, .lscale    , .escale    ),
                theta2eta(aa.init, .lshape1.a , .eshape1.a )) else
          cbind(theta2eta(aa.init, .lshape1.a , .eshape1.a ),
                theta2eta(sc.init, .lscale    , .escale    ))
      etastart <- etastart[, interleave.VGAM(M, M1 = M1)]
    }  # End of etastart.
  }), list( .lscale    = lscale   , .lshape1.a = lshape1.a,
            .escale    = escale   , .eshape1.a = eshape1.a,
            .iscale    = iscale   , .ishape1.a = ishape1.a,
            .gscale    = gscale   , .gshape1.a = gshape1.a,
            .imethod   = imethod  , .probs.y  = probs.y,
            .lss = lss ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    M1 <- 2
    NOS <- ncol(eta) / M1
    if ( .lss ) {
      Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lscale ,  earg = .escale )
      aa    <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lshape1.a , earg = .eshape1.a )
    } else {
      aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a )
      Scale <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lscale    ,  earg = .escale )
    }
    parg <- 1
    qq   <- aa

    ans <- Scale * exp(lgamma(parg + 1/aa) +
                       lgamma(qq   - 1/aa) -
                       lgamma(parg) - lgamma(qq))
    ans[parg + 1/aa <= 0] <- NA
    ans[qq   - 1/aa <= 0] <- NA
    ans[aa          <= 0] <- NA
    ans[Scale       <= 0] <- NA
    ans
  }, list( .lscale    = lscale   , .lshape1.a = lshape1.a,
           .escale    = escale   , .eshape1.a = eshape1.a,
           .lss = lss ))),
  last = eval(substitute(expression({
    M1 <- 2

    misc$link <- c(rep_len(if ( .lss ) .lscale else .lshape1.a ,
                           ncoly),
                   rep_len(if ( .lss ) .lshape1.a else .lscale ,
                           ncoly))[
                   interleave.VGAM(M, M1 = M1)]
    temp.names <- if ( .lss ) {
      c(scaL.names, sha1.names)
    } else {
      c(sha1.names, scaL.names)
    }
    names(misc$link) <- temp.names[interleave.VGAM(M, M1 = M1)]

    misc$earg <- vector("list", M)
    names(misc$earg) <- temp.names[interleave.VGAM(M, M1 = M1)]
    for (ii in 1:ncoly)
      if ( .lss ) {
        misc$earg[[M1*ii-1]] <- .escale
        misc$earg[[M1*ii  ]] <- .eshape1.a
      } else {
        misc$earg[[M1*ii-1]] <- .eshape1.a
        misc$earg[[M1*ii  ]] <- .escale
      }

    misc$expected          <- TRUE
    misc$multipleResponses <- TRUE
  }), list( .lscale    = lscale   , .lshape1.a = lshape1.a,
            .escale    = escale   , .eshape1.a = eshape1.a,
            .lss = lss ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE,
             eta, extra = NULL, summation = TRUE) {
      M1  <- 2
      NOS <- ncol(eta)/M1
      if ( .lss ) {
        Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                           .lscale    , earg = .escale    )
        aa    <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                           .lshape1.a , earg = .eshape1.a )
      } else {
        aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                           .lshape1.a , earg = .eshape1.a )
        Scale <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                           .lscale    , earg = .escale    )
      }
      parg   <- 1
      qq     <- aa
      if (residuals) {
        stop("loglikelihood residuals not implemented yet")
      } else {
        ll.elts <-
          c(w) * dgenbetaII(x = y, scale = Scale, shape1.a = aa,
                            shape2.p = parg, shape3.q = aa,
                            log = TRUE)
        if (summation) {
          sum(ll.elts)
        } else {
          ll.elts
        }
      }
    }, list( .lscale    = lscale   , .lshape1.a = lshape1.a,
             .escale    = escale   , .eshape1.a = eshape1.a,
             .lss = lss ))),
  vfamily = c("paralogistic"),
  simslot = eval(substitute(
  function(object, nsim) {

    pwts <- if (length(pwts <- object@prior.weights) > 0)
              pwts else weights(object, type = "prior")
    if (any(pwts != 1))
      warning("ignoring prior weights")

    eta <- predict(object)
    M1 <- 2
    NOS <- ncol(eta)/M1
    if ( .lss ) {
      Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lscale ,  earg = .escale )
      aa    <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lshape1.a , earg = .eshape1.a )
    } else {
      aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a )
      Scale <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lscale    ,  earg = .escale )
    }
    rparalogistic(nsim * length(Scale), shape1.a = aa, scale = Scale)
  }, list( .lscale    = lscale   , .lshape1.a = lshape1.a,
           .escale    = escale   , .eshape1.a = eshape1.a,
           .lss = lss ))),

  validparams = eval(substitute(function(eta, y, extra = NULL) {
    M1 <- 2
    NOS <- ncol(eta) / M1
    if ( .lss ) {
      Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lscale    , earg = .escale   )
      aa    <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lshape1.a , earg = .eshape1.a)
    } else {
      aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a)
      Scale <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lscale    , earg = .escale )
    }
    parg <- 1
    qq   <- aa
    okay1 <- all(is.finite(Scale)) && all(Scale  > 0) &&
             all(is.finite(aa   )) && all(aa     > 0) &&
             all(is.finite(parg )) && all(parg   > 0) &&
             all(is.finite(qq   )) && all(qq     > 0)
    okay.support <- if (okay1)
                    all(-aa * parg < 1 & 1 < aa * qq) else TRUE
    if (!okay.support)
      warning("parameter constraint -a   < 1 < a*a has",
              " been violated; solution may be at",
              " the boundary of the parameter space.")
    okay1 && okay.support
  }, list(  .lscale    = lscale   , .lshape1.a = lshape1.a,
            .escale    = escale   , .eshape1.a = eshape1.a,
            .lss = lss ))),




  deriv = eval(substitute(expression({
    M1 <- 2
    NOS <- ncol(eta)/M1  # Needed for summary()
    if ( .lss ) {
      Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lscale    , earg = .escale   )
      aa    <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lshape1.a , earg = .eshape1.a)
    } else {
      aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a)
      Scale <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lscale    , earg = .escale )
    }
    parg <- 1
    qq   <- aa

    temp1 <- log(y/Scale)
    temp2 <- (y/Scale)^aa
    temp3 <- digamma(parg + qq)
    temp3a <- digamma(parg)
    temp3b <- digamma(qq)
    temp4 <- log1p(temp2)

    dl.dscale <- (aa/Scale) * (-parg + (parg+qq) / (1+1/temp2))
    dl.da <- 1/aa + parg * temp1 - (parg+qq) * temp1 / (1+1/temp2)

    dscale.deta <- dtheta.deta(Scale, .lscale ,    .escale )
    da.deta     <- dtheta.deta(aa,    .lshape1.a , .eshape1.a )

    myderiv <- if ( .lss ) {
      c(w) * cbind(dl.dscale * dscale.deta,
                   dl.da * da.deta)
    } else {
      c(w) * cbind(dl.da * da.deta,
                   dl.dscale * dscale.deta)
    }
    myderiv[, interleave.VGAM(M, M1 = M1)]
  }), list(  .lscale    = lscale   , .lshape1.a = lshape1.a,
             .escale    = escale   , .eshape1.a = eshape1.a,
             .lss = lss ))),
  weight = eval(substitute(expression({
    temp5  <- trigamma(parg + qq)
    temp5a <- trigamma(parg)
    temp5b <- trigamma(qq)

    ned2l.da <- (1 + parg + qq + parg * qq * (temp5a + temp5b +
                (temp3b - temp3a + (parg-qq)/(parg*qq))^2 -
         (parg^2 + qq^2) / (parg*qq)^2)) / (aa^2 * (1+parg+qq))
    ned2l.dscale  <- (aa^2) * parg * qq / ((1+parg+qq) * Scale^2)
    ned2l.dascale <- (parg - qq - parg * qq *
                     (temp3a -temp3b)) / (Scale*(1 + parg+qq))
    wz <- if ( .lss ) {
      array(c(c(w) * ned2l.dscale * dscale.deta^2,
              c(w) * ned2l.da * da.deta^2,
              c(w) * ned2l.dascale * da.deta * dscale.deta),
            dim = c(n, M/M1, M1*(M1+1)/2))
    } else {
      array(c(c(w) * ned2l.da * da.deta^2,
              c(w) * ned2l.dscale * dscale.deta^2,
              c(w) * ned2l.dascale * da.deta * dscale.deta),
            dim = c(n, M/M1, M1*(M1+1)/2))
    }
    wz <- arwz2wz(wz, M = M, M1 = M1)
    wz
  }), list( .lscale    = lscale   , .lshape1.a = lshape1.a,
            .escale    = escale   , .eshape1.a = eshape1.a,
            .lss = lss))))
}  # paralogistic






 inv.paralogistic <-
   function(lscale    = "loglink",
            lshape1.a = "loglink",
            iscale    = NULL,
            ishape1.a = NULL,
            imethod   = 1,
            lss       = TRUE,
            gscale    = exp(-5:5),
            gshape1.a = seq(0.75, 4, by = 0.25),  # exp(-5:5),
            probs.y   = c(0.25, 0.50, 0.75),
            zero      = "shape") {




  if (length(lss) != 1 && !is.logical(lss))
    stop("Argument 'lss' not specified correctly")

  if (!is.Numeric(imethod, length.arg = 1,
                  integer.valued = TRUE,
                  positive = TRUE) || imethod > 2)
    stop("Bad input for argument 'imethod'")

  if (length(iscale) && !is.Numeric(iscale, positive = TRUE))
    stop("Bad input for argument 'iscale'")

  if (length(ishape1.a) && !is.Numeric(ishape1.a, positive = TRUE))
    stop("Bad input for argument 'ishape1.a'")

  if (length(probs.y) < 2 || max(probs.y) > 1 ||
        !is.Numeric(probs.y, positive = TRUE))
    stop("Bad input for argument 'probs.y'")


  lscale <- as.list(substitute(lscale))
  escale <- link2list(lscale)
  lscale <- attr(escale, "function.name")

  lshape1.a <- as.list(substitute(lshape1.a))
  eshape1.a <- link2list(lshape1.a)
  lshape1.a <- attr(eshape1.a, "function.name")


  new("vglmff",
  blurb =
    c("Inverse paralogistic distribution \n\n",
      "Links:    ",
      ifelse (lss,
              namesof("scale"   , lscale   , escale),
              namesof("shape1.a", lshape1.a, eshape1.a)), ", ",
      ifelse (lss,
              namesof("shape1.a", lshape1.a, eshape1.a),
              namesof("scale"   , lscale   , escale)), "\n",
      "Mean:     scale * gamma(shape1.a + 1/shape1.a) * ",
                "gamma(1 - 1/shape1.a) / ",
                "gamma(shape1.a)"),
  constraints = eval(substitute(expression({
    constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
                     M = M, M1 = 2,
                     predictors.names = predictors.names)

  }), list( .zero = zero ))),
  infos = eval(substitute(function(...) {
    list(M1 = 2,
         Q1 = 1,
         dpqrfun = "inv.paralogistic",  # "genbetaII",
         expected = TRUE,
         zero = .zero ,
         multipleResponses = TRUE,
         parameters.names = if ( .lss )
           c("scale", "shape1.a") else
           c("shape1.a", "scale"),
         lscale    = .lscale    , lshape1.a = .lshape1.a ,
         escale    = .escale    , eshape1.a = .eshape1.a )
  }, list( .lscale = lscale      , .lshape1.a = lshape1.a,
           .escale = escale      , .eshape1.a = eshape1.a,
           .lss  = lss ,
           .zero = zero ))),

  rqresslot = eval(substitute(
    function(mu, y, w, eta, extra = NULL) {
      M1  <- 2
      NOS <- ncol(eta) / M1
    if ( .lss ) {
      Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lscale ,  earg = .escale )
      aa    <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lshape1.a , earg = .eshape1.a )
    } else {
      aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a )
      Scale <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lscale    ,  earg = .escale )
    }
    parg <- aa
    qq   <- 1
    scrambleseed <- runif(1)  # To scramble the seed
    qnorm(pinv.paralogistic(y, scale = Scale, shape1.a = aa))
    }, list( .lscale    = lscale   , .lshape1.a = lshape1.a,
             .escale    = escale   , .eshape1.a = eshape1.a,
             .lss = lss ))),
  initialize = eval(substitute(expression({
    temp5 <- w.y.check(w = w, y = y,
                       Is.positive.y = TRUE,
                       ncol.w.max = Inf,
                       ncol.y.max = Inf,
                       out.wy = TRUE,
                       colsyperw = 1,
                       maximize = TRUE)
    y    <- temp5$y
    w    <- temp5$w
    M1 <- 2  # Number of parameters for one response
    NOS  <- ncoly <- ncol(y)
    M    <- M1*ncol(y)

    scaL.names  <- param.names("scale",     NOS, skip1 = TRUE)
    sha1.names  <- param.names("shape1.a",  NOS, skip1 = TRUE)

    predictors.names <- if ( .lss ) {
    c(namesof(scaL.names , .lscale    , .escale    , tag = FALSE),
      namesof(sha1.names , .lshape1.a , .eshape1.a , tag = FALSE))
    } else {
    c(namesof(sha1.names , .lshape1.a , .eshape1.a , tag = FALSE),
      namesof(scaL.names , .lscale    , .escale    , tag = FALSE))
   }
    predictors.names <- predictors.names[interleave.VGAM(M, M1 = M1)]

    if (!length(etastart)) {
      sc.init <-
      aa.init <- matrix(NA_real_, n, NOS)

      for (spp. in 1:NOS) {  # For each response 'y_spp.'... do:
        yvec <- y[, spp.]
        wvec <- w[, spp.]

        if ( .imethod == 1 ) {
          gscale     <- .gscale
          gshape1.a  <- .gshape1.a
          if (length( .iscale    )) gscale    <-
                                        rep_len( .iscale    , NOS)
          if (length( .ishape1.a )) gshape1.a <-
                                        rep_len( .ishape1.a , NOS)


        inv.paralogistic.Loglikfun2 <-
            function(scaleval, shape1.a,
                     y, x, w, extraargs) {
          sum(c(w) * dgenbetaII(x = y, scale = scaleval,
                                shape1.a = shape1.a,
                                shape2.p = shape1.a,
                                shape3.q = 1, log = TRUE))
        }


        try.this <-
          grid.search2(gscale, gshape1.a, # vov3 = 1, vov4 = gshape1.a,
                       objfun = inv.paralogistic.Loglikfun2,  # x = x,
                       y = yvec, w = wvec,  # extraargs = NULL,
                       ret.objfun = TRUE)  # Last value is \ell
          sc.init[, spp.] <- try.this["Value1"]
          aa.init[, spp.] <- try.this["Value2"]
        } else {  # .imethod == 2
          qvec <- .probs.y
          ishape2.p <- if (length( .ishape1.a )) .ishape1.a else 1
          xvec <- log( qvec^(-1/ ishape2.p) - 1 )
          fit0 <- lsfit(x = xvec, y = log(quantile(yvec,  qvec)))
          sc.init[, spp.] <- if (length( .iscale )) .iscale else
                             exp(fit0$coef[1])
            aa.init[, spp.] <- if (length( .ishape1.a ))
                                   .ishape1.a else
                             -1/fit0$coef[2]
        }
      }  # End of for (spp. ...)



      finite.mean <- 1 < aa.init
      COP.use <- 1.15
      while (FALSE && any(!finite.mean)) {
    aa.init[!finite.mean] <- 0.1 + aa.init[!finite.mean] * COP.use
    finite.mean <- 1 < aa.init
      }

      etastart <- if ( .lss )
          cbind(theta2eta(sc.init, .lscale    , .escale    ),
                theta2eta(aa.init, .lshape1.a , .eshape1.a )) else
          cbind(theta2eta(aa.init, .lshape1.a , .eshape1.a ),
                theta2eta(sc.init, .lscale    , .escale    ))
      etastart <- etastart[, interleave.VGAM(M, M1 = M1)]
    }  # End of etastart.
  }), list( .lscale    = lscale   , .lshape1.a = lshape1.a,
            .escale    = escale   , .eshape1.a = eshape1.a,
            .iscale    = iscale   , .ishape1.a = ishape1.a,
            .gscale    = gscale   , .gshape1.a = gshape1.a,
            .imethod   = imethod  , .probs.y   = probs.y,
            .lss = lss ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    M1 <- 2
    NOS <- ncol(eta) / M1
    if ( .lss ) {
      Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lscale ,  earg = .escale )
      aa    <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lshape1.a , earg = .eshape1.a )
    } else {
      aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a )
      Scale <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lscale    ,  earg = .escale )
    }
    parg <- aa
    qq   <- 1

    ans <- Scale * exp(lgamma(parg + 1/aa) +
                       lgamma(qq   - 1/aa) -
                       lgamma(parg) - lgamma(qq))
    ans[parg + 1/aa <= 0] <- NA
    ans[qq   - 1/aa <= 0] <- NA
    ans[aa          <= 0] <- NA
    ans[Scale       <= 0] <- NA
    ans
  }, list( .lscale    = lscale   , .lshape1.a = lshape1.a,
           .escale    = escale   , .eshape1.a = eshape1.a,
           .lss = lss ))),
  last = eval(substitute(expression({
    M1 <- 2

    misc$link <- c(rep_len(if ( .lss ) .lscale else .lshape1.a ,
                           ncoly),
                   rep_len(if ( .lss ) .lshape1.a else .lscale ,
                           ncoly))[
                   interleave.VGAM(M, M1 = M1)]
    temp.names <- if ( .lss ) {
      c(scaL.names, sha1.names)
    } else {
      c(sha1.names, scaL.names)
    }
    names(misc$link) <- temp.names[interleave.VGAM(M, M1 = M1)]

    misc$earg <- vector("list", M)
    names(misc$earg) <- temp.names[interleave.VGAM(M, M1 = M1)]
    for (ii in 1:ncoly)
      if ( .lss ) {
        misc$earg[[M1*ii-1]] <- .escale
        misc$earg[[M1*ii  ]] <- .eshape1.a
      } else {
        misc$earg[[M1*ii-1]] <- .eshape1.a
        misc$earg[[M1*ii  ]] <- .escale
      }

    misc$expected          <- TRUE
    misc$multipleResponses <- TRUE
  }), list( .lscale    = lscale   , .lshape1.a = lshape1.a,
            .escale    = escale   , .eshape1.a = eshape1.a,
            .lss = lss ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE,
             eta, extra = NULL, summation = TRUE) {
      M1  <- 2
      NOS <- ncol(eta)/M1
      if ( .lss ) {
        Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                           .lscale    , earg = .escale    )
        aa    <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                           .lshape1.a , earg = .eshape1.a )
      } else {
        aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                           .lshape1.a , earg = .eshape1.a )
        Scale <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                           .lscale    , earg = .escale    )
      }
      parg   <- aa
      if (residuals) {
        stop("loglikelihood residuals not implemented yet")
      } else {
        ll.elts <-
          c(w) * dgenbetaII(x = y, scale = Scale, shape1.a = aa,
                            shape2.p = aa, shape3.q = 1,
                            log = TRUE)
        if (summation) {
          sum(ll.elts)
        } else {
          ll.elts
        }
      }
    }, list( .lscale    = lscale   , .lshape1.a = lshape1.a,
             .escale    = escale   , .eshape1.a = eshape1.a,
             .lss = lss ))),
  vfamily = c("inv.paralogistic"),

  simslot = eval(substitute(
  function(object, nsim) {

    pwts <- if (length(pwts <- object@prior.weights) > 0)
              pwts else weights(object, type = "prior")
    if (any(pwts != 1))
      warning("ignoring prior weights")
    eta <- predict(object)
    if ( .lss ) {
      Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lscale    , earg = .escale   )
      aa    <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lshape1.a , earg = .eshape1.a)
    } else {
      aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a)
      Scale <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lscale    , earg = .escale )
    }
    parg <- aa
    qq   <- 1
    rinv.paralogistic(nsim * length(Scale), shape1.a = aa, scale = Scale)
  }, list(  .lscale    = lscale   , .lshape1.a = lshape1.a,
            .escale    = escale   , .eshape1.a = eshape1.a,
            .lss = lss ))),



  validparams = eval(substitute(function(eta, y, extra = NULL) {
    M1 <- 2
    NOS <- ncol(eta) / M1
    if ( .lss ) {
      Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lscale    , earg = .escale   )
      aa    <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lshape1.a , earg = .eshape1.a)
    } else {
      aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a)
      Scale <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lscale    , earg = .escale )
    }
    parg <- aa
    qq   <- 1
    okay1 <- all(is.finite(Scale)) && all(Scale  > 0) &&
             all(is.finite(aa   )) && all(aa     > 0) &&
             all(is.finite(parg )) && all(parg   > 0) &&
             all(is.finite(qq   )) && all(qq     > 0)
    okay.support <- if (okay1)
                    all(-aa * parg < 1 & 1 < aa * qq) else TRUE
    if (!okay.support)
      warning("parameter constraint -a*a < 1 < a   has",
              " been violated; solution may be at the ",
              " boundary of the parameter space.")
    okay1 && okay.support
  }, list(  .lscale    = lscale   , .lshape1.a = lshape1.a,
            .escale    = escale   , .eshape1.a = eshape1.a,
            .lss = lss ))),




  deriv = eval(substitute(expression({
    M1 <- 2
    NOS <- ncol(eta)/M1  # Needed for summary()
    if ( .lss ) {
      Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lscale    , earg = .escale   )
      aa    <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lshape1.a , earg = .eshape1.a)
    } else {
      aa    <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE],
                         .lshape1.a , earg = .eshape1.a)
      Scale <- eta2theta(eta[, M1*(1:NOS)    , drop = FALSE],
                         .lscale    , earg = .escale )
    }
    parg <- aa
    qq   <- 1
    temp1 <- log(y/Scale)
    temp2 <- (y/Scale)^aa
    temp3 <- digamma(parg + qq)
    temp3a <- digamma(parg)
    temp3b <- digamma(qq)
    temp4 <- log1p(temp2)

    dl.dscale <- (aa/Scale) * (-parg + (parg+qq) / (1+1/temp2))
    dl.da <- 1/aa + parg * temp1 - (parg+qq) * temp1 / (1+1/temp2)

    dscale.deta <- dtheta.deta(Scale, .lscale ,    earg = .escale )
    da.deta     <- dtheta.deta(aa,    .lshape1.a , earg = .eshape1.a )

    myderiv <- if ( .lss ) {
      c(w) * cbind(dl.dscale * dscale.deta,
                   dl.da * da.deta)
    } else {
      c(w) * cbind(dl.da * da.deta,
                   dl.dscale * dscale.deta)
    }
    myderiv[, interleave.VGAM(M, M1 = M1)]
  }), list(  .lscale    = lscale   , .lshape1.a = lshape1.a,
             .escale    = escale   , .eshape1.a = eshape1.a,
             .lss = lss ))),
  weight = eval(substitute(expression({
    temp5  <- trigamma(parg + qq)
    temp5a <- trigamma(parg)
    temp5b <- trigamma(qq)

    ned2l.da <- (1 + parg + qq + parg * qq * (temp5a + temp5b +
                (temp3b - temp3a + (parg-qq)/(parg*qq))^2 -
                (parg^2 + qq^2) / (parg*qq)^2)) / (aa^2 * (1+parg+qq))
    ned2l.dscale  <- (aa^2) * parg * qq / ((1+parg+qq) * Scale^2)
    ned2l.dascale <- (parg - qq - parg * qq *
                     (temp3a -temp3b)) / (Scale*(1 + parg+qq))
    wz <- if ( .lss ) {
      array(c(c(w) * ned2l.dscale * dscale.deta^2,
              c(w) * ned2l.da * da.deta^2,
              c(w) * ned2l.dascale * da.deta * dscale.deta),
            dim = c(n, M/M1, M1*(M1+1)/2))
    } else {
      array(c(c(w) * ned2l.da * da.deta^2,
              c(w) * ned2l.dscale * dscale.deta^2,
              c(w) * ned2l.dascale * da.deta * dscale.deta),
            dim = c(n, M/M1, M1*(M1+1)/2))
    }
    wz <- arwz2wz(wz, M = M, M1 = M1)
    wz
  }), list( .lscale    = lscale   , .lshape1.a = lshape1.a,
            .escale    = escale   , .eshape1.a = eshape1.a,
            .lss = lss ))))
}  # inv.paralogistic

Try the VGAM package in your browser

Any scripts or data that you put into this service are public.

VGAM documentation built on Sept. 18, 2024, 9:09 a.m.