R/family.others.R

Defines functions eCDF.vglm qtplot.extlogF1 fix.crossing.vglm is.crossing.vglm extlogF1 extlogF1.control dextlogF rzipfmb qzipfmb pzipfmb dzipfmb AR1.gammas AR1EIM pzoibetabinom pzoibetabinom.ab rzoibetabinom rzoibetabinom.ab dzoibetabinom dzoibetabinom.ab log1pexp log1mexp qzoabeta pzoabeta rzoabeta dzoabeta tpnff3 tpnff rtpn qtpn pos ptpn dtpn rweibull3 qweibull3 pweibull3 dweibull3 explogff explogff.control rexplog qexplog pexplog dexplog expgeometric expgeometric.control rexpgeom qexpgeom pexpgeom dexpgeom genrayleigh genrayleigh.control rgenray qgenray pgenray dgenray exppoisson rexppois pexppois qexppois dexppois

Documented in AR1EIM AR1.gammas dexpgeom dexpgeom dexplog dexppois dextlogF dgenray dweibull3 dzipfmb dzoabeta dzoibetabinom dzoibetabinom.ab eCDF.vglm expgeometric expgeometric explogff exppoisson extlogF1 extlogF1.control fix.crossing.vglm genrayleigh is.crossing.vglm log1mexp log1pexp pexpgeom pexpgeom pexplog pexppois pgenray pweibull3 pzipfmb pzoabeta pzoibetabinom pzoibetabinom.ab qexpgeom qexpgeom qexplog qexppois qgenray qweibull3 qzipfmb qzoabeta rexpgeom rexpgeom rexplog rexppois rgenray rweibull3 rzipfmb rzoabeta rzoibetabinom rzoibetabinom.ab

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


















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


  N <- max(length(x), length(shape), length(rate))
  if (length(x)      != N) x      <- rep_len(x,     N)
  if (length(shape)  != N) shape  <- rep_len(shape, N)
  if (length(rate)   != N) rate   <- rep_len(rate,  N)
  logdensity <- rep_len(log(0), N)

  xok <- (0 < x)

  logdensity[xok] <- log(shape[xok]) + log(rate[xok]) -
                     log1p(-exp(-shape[xok])) - shape[xok] -
                     rate[xok] * x[xok] + shape[xok] *
                     exp(-rate[xok] * x[xok])

  logdensity[shape <= 0] <- NaN
  logdensity[rate <= 0] <- NaN
  if (log.arg) logdensity else exp(logdensity)
}






qexppois<- function(p, rate = 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'")

  if (lower.tail) {
    if (log.p) {
      ln.p <- p
      ans <- -log(log(exp(ln.p) *
                      (-expm1(shape)) + exp(shape)) / shape) / rate
      ans[ln.p > 0] <- NaN
    } else {
      ans <- -log(log(p * (-expm1(shape)) + exp(shape)) / shape) / rate
      ans[p < 0] <- NaN
      ans[p > 1] <- NaN
    }
  } else {
    if (log.p) {
      ln.p <- p
      ans <- -log(log(expm1(ln.p) *
                      expm1(shape) + exp(shape)) / shape) / rate
      ans[ln.p > 0] <- NaN
    } else {
      ans <- -log(log(p * expm1(shape) + 1) / shape) / rate
      ans[p < 0] <- NaN
      ans[p > 1] <- NaN
    }
  }
  ans[(shape <= 0) | (rate <= 0)] <- NaN
  ans
}














pexppois<- function(q, rate = 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'")

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




rexppois <- function(n, rate = 1, shape) {
  ans <- -log(log(runif(n) * (-expm1(shape)) +
         exp(shape)) / shape) / rate
  ans[(shape <= 0) | (rate <= 0)] <- NaN
  ans
}









 exppoisson <- function(lrate = "loglink", lshape = "loglink",
                        irate = 2.0, ishape = 1.1,
                        zero = NULL) {

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

  lratee <- as.list(substitute(lrate))
  eratee <- link2list(lratee)
  lratee <- attr(eratee, "function.name")


  iratee <- irate




  if (length(ishape) &&
      !is.Numeric(ishape, positive = TRUE))
    stop("bad input for argument 'ishape'")
  if (length(iratee) &&
      !is.Numeric(iratee, positive = TRUE))
    stop("bad input for argument 'irate'")

  ishape[abs(ishape - 1) < 0.01] = 1.1


  new("vglmff",
  blurb = c("Exponential Poisson distribution \n \n",
            "Links:    ",
            namesof("rate",  lratee, earg = eratee), ", ",
            namesof("shape", lshape, earg = eshape), "\n",
            "Mean:     shape/(expm1(shape) * rate)) * ",
                      "genhypergeo(c(1, 1), c(2, 2), shape)"),


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

  infos = eval(substitute(function(...) {
    list(M1 = 2,
         Q1 = 1,
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = c("rate", "shape"),
         lrate  = .lratee ,
         lshape = .lshape ,
         zero = .zero )
  }, list( .zero = zero, .lratee = lratee, .lshape = lshape ))),


  initialize = eval(substitute(expression({

    temp5 <-
    w.y.check(w = w, y = y,
              out.wy = TRUE,
              maximize = TRUE)
    w <- temp5$w
    y <- temp5$y


    predictors.names <- c(
      namesof("rate",  .lratee , earg = .eratee , short = TRUE),
      namesof("shape", .lshape , earg = .eshape , short = TRUE))

    if (!length(etastart)) {
      ratee.init <- if (length( .iratee ))
              rep_len( .iratee , n) else
              stop("Need to input a value into argument 'iratee'")
      shape.init <- if (length( .ishape ))
                      rep_len( .ishape , n) else
                      (1/ratee.init - mean(y)) / ((y *
                      exp(-ratee.init * y))/n)


      ratee.init <- rep_len(weighted.mean(ratee.init, w = w), n)

      etastart <-
        cbind(theta2eta(ratee.init, .lratee , earg = .eratee ),
              theta2eta(shape.init, .lshape , earg = .eshape ))


    }
  }), list( .lshape = lshape, .lratee = lratee,
            .ishape = ishape, .iratee = iratee,
            .eshape = eshape, .eratee = eratee))),

  linkinv = eval(substitute(function(eta, extra = NULL) {
    ratee <- eta2theta(eta[, 1], .lratee , earg = .eratee )
    shape <- eta2theta(eta[, 2], .lshape , earg = .eshape )





    qexppois(p = 0.5, rate = ratee, shape = shape)
  }, list( .lshape = lshape, .lratee = lratee,
           .eshape = eshape, .eratee = eratee))),

  last = eval(substitute(expression({
    misc$link <-    c( rate = .lratee , shape = .lshape )
    misc$earg <- list( rate = .eratee , shape = .eshape )

    misc$expected <- TRUE
    misc$multipleResponses <- FALSE
  }), list( .lshape = lshape, .lratee = lratee,
            .eshape = eshape, .eratee = eratee))),

  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    ratee <- eta2theta(eta[, 1], .lratee , earg = .eratee )
    shape <- eta2theta(eta[, 2], .lshape , earg = .eshape )

    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * dexppois(x = y, shape = shape, rate = ratee,
                                 log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .lratee = lratee , .lshape = lshape ,
           .eshape = eshape , .eratee = eratee ))),

  vfamily = c("exppoisson"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    ratee <- eta2theta(eta[, 1], .lratee , earg = .eratee )
    shape <- eta2theta(eta[, 2], .lshape , earg = .eshape )
    okay1 <- all(is.finite(ratee)) && all(0 < ratee) &&
             all(is.finite(shape)) && all(0 < shape)
    okay1
  }, list( .lratee = lratee , .lshape = lshape ,
           .eshape = eshape , .eratee = eratee ))),

  deriv = eval(substitute(expression({
    ratee <- eta2theta(eta[, 1], .lratee , earg = .eratee )
    shape <- eta2theta(eta[, 2], .lshape , earg = .eshape )

    dl.dratee <- 1/ratee - y - y * shape * exp(-ratee * y)
    dl.dshape <- 1/shape - 1/expm1(shape) - 1 + exp(-ratee * y)
    dratee.deta <- dtheta.deta(ratee, .lratee , earg = .eratee )
    dshape.deta <- dtheta.deta(shape, .lshape , earg = .eshape )
    c(w) * cbind(dl.dratee * dratee.deta,
                 dl.dshape * dshape.deta)
  }), list( .lshape = lshape, .lratee = lratee,
            .eshape = eshape, .eratee = eratee ))),

  weight = eval(substitute(expression({

    temp1 <- -expm1(-shape)

    ned2l.dshape2 <- (1 + exp(2 * shape) - shape^2 * exp(shape) - 2 *
                      exp(shape)) / (shape * temp1)^2


    ned2l.dratee2 <- 1 / ratee^2 - (shape^2 * exp(-shape) / (4 *
                    ratee^2 * temp1)) *
                    genhypergeo(c(2, 2, 2), c(3, 3, 3), shape)

    ned2l.drateeshape <- (shape * exp(-shape) / (4 * ratee * temp1)) *
                           genhypergeo(c(2, 2), c(3, 3), shape)

    wz <- matrix(0, n, dimm(M))
    wz[, iam(1, 1, M)] <- dratee.deta^2 * ned2l.dratee2
    wz[, iam(1, 2, M)] <- dratee.deta * dshape.deta * ned2l.drateeshape
    wz[, iam(2, 2, M)] <- dshape.deta^2 * ned2l.dshape2
    c(w) * wz
  }), list( .zero = zero ))))
}










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



  N <- max(length(x), length(shape), length(scale))
  if (length(x)      != N) x      <- rep_len(x,     N)
  if (length(shape)  != N) shape  <- rep_len(shape, N)
  if (length(scale)  != N) scale  <- rep_len(scale, N)
  logdensity <- rep_len(log(0), N)

  if (any(xok <- (x > 0))) {
    temp1 <- x[xok] / scale[xok]
    logdensity[xok] <- log(2) + log(shape[xok]) + log(x[xok]) -
                       2 * log(scale[xok]) - temp1^2  +
                       (shape[xok] - 1) * log1p(-exp(-temp1^2))
  }
  logdensity[(shape <= 0) | (scale <= 0)] <- NaN
  logdensity[is.infinite(x)] <- log(0)  # 20141209 KaiH
  if (log.arg) {
    logdensity
  } else {
     exp(logdensity)
  }
}






pgenray <- 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'")

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









qgenray <- 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'")

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






rgenray <- function(n, scale = 1, shape) {
  ans <- qgenray(runif(n), shape = shape, scale = scale)
  ans[(shape <= 0) | (scale <= 0)] <- NaN
  ans
}




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


 genrayleigh <-
  function(lscale = "loglink", lshape = "loglink",
           iscale = NULL,   ishape = NULL,
           tol12 = 1.0e-05,
           nsimEIM = 300, zero = 2) {

  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 (length(ishape) &&
      !is.Numeric(ishape, positive = TRUE))
    stop("bad input for argument 'ishape'")
  if (length(iscale) &&
      !is.Numeric(iscale, positive = TRUE))
    stop("bad input for argument 'iscale'")

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



  new("vglmff",
  blurb = c("Generalized Rayleigh distribution\n",
            "Links:    ",
            namesof("scale", lscale, earg = escale), ", ",
            namesof("shape", lshape, earg = eshape), "\n"),
  constraints = eval(substitute(expression({
    constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
                                predictors.names = predictors.names,
                                M1 = 2)
  }), list( .zero = zero ))),

  infos = eval(substitute(function(...) {
    list(M1 = 2,
         Q1 = 1,
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = c("scale", "shape"),
         nsimEIM = .nsimEIM ,
         lscale = .lscale ,
         lshape = .lshape ,
         zero = .zero )
  }, list( .zero = zero, .lscale = lscale, .lshape = lshape,
           .nsimEIM = nsimEIM ))),


  initialize = eval(substitute(expression({

    temp5 <-
    w.y.check(w = w, y = y,
              out.wy = TRUE,
              maximize = TRUE)
    w <- temp5$w
    y <- temp5$y





    predictors.names <- c(
      namesof("scale", .lscale , earg = .escale , short = TRUE),
      namesof("shape", .lshape , earg = .eshape , short = TRUE))

    if (!length(etastart)) {
      genrayleigh.Loglikfun <- function(scale, y, x, w, extraargs) {
        temp1 <- y / scale
        shape <- -1 / weighted.mean(log1p(-exp(-temp1^2)), w = w)

        ans <- sum(c(w) * (log(2) + log(shape) + log(y) -
                           2 * log(scale) - temp1^2  +
                           (shape - 1) * log1p(-exp(-temp1^2))))
        ans
      }
      scale.grid <- seq(0.2 * stats::sd(c(y)),
                        5.0 * stats::sd(c(y)), len = 29)
      scale.init <- if (length( .iscale )) .iscale else
                    grid.search(scale.grid, objfun = genrayleigh.Loglikfun,
                                y = y, x = x, w = w)
      scale.init <- rep_len(scale.init, length(y))

      shape.init <- if (length( .ishape )) .ishape else
                    -1 / weighted.mean(log1p(-exp(-(y/scale.init)^2)),
                     w = w)
      shape.init <- rep_len(shape.init, length(y))

      etastart <- cbind(theta2eta(scale.init, .lscale , earg = .escale ),
                        theta2eta(shape.init, .lshape , earg = .eshape ))

        }
    }), list( .lscale = lscale, .lshape = lshape,
              .iscale = iscale, .ishape = ishape,
              .escale = escale, .eshape = eshape))),

  linkinv = eval(substitute(function(eta, extra = NULL) {
    Scale <- eta2theta(eta[, 1], .lscale , earg = .escale )
    shape <- eta2theta(eta[, 2], .lshape , earg = .eshape )
    qgenray(p = 0.5, shape = shape, scale = Scale)
  }, list( .lshape = lshape, .lscale = lscale,
           .eshape = eshape, .escale = escale ))),

  last = eval(substitute(expression({
    misc$link <-    c(scale = .lscale, shape = .lshape )

    misc$earg <- list(scale = .escale, shape = .eshape )

    misc$expected <- TRUE
    misc$nsimEIM <- .nsimEIM
    misc$multipleResponses <- FALSE
  }), 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[, 1], .lscale , earg = .escale )
    shape <- eta2theta(eta[, 2], .lshape , earg = .eshape )

    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * dgenray(x = y, shape = shape,
                                scale = Scale, log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .lshape = lshape , .lscale = lscale ,
           .eshape = eshape , .escale = escale ))),

  vfamily = c("genrayleigh"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    Scale <- eta2theta(eta[, 1], .lscale , earg = .escale )
    shape <- eta2theta(eta[, 2], .lshape , earg = .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 ))),

  deriv = eval(substitute(expression({
    Scale <- eta2theta(eta[, 1], .lscale , earg = .escale )
    shape <- eta2theta(eta[, 2], .lshape , earg = .eshape )
    dscale.deta <- dtheta.deta(Scale, .lscale , earg = .escale )
    dshape.deta <- dtheta.deta(shape, .lshape , earg = .eshape )
    dthetas.detas <- cbind(dscale.deta, dshape.deta)

    temp1 <- y / Scale
    temp2 <- exp(-temp1^2)
    temp3 <- temp1^2 / Scale
    AAA   <- 2 * temp1^2 / Scale  # 2 * y^2 / Scale^3
    BBB   <- -expm1(-temp1^2)     # denominator
    dl.dshape <- 1/shape + log1p(-temp2)
    dl.dscale <- -2 / Scale + AAA * (1 - (shape - 1) * temp2 / BBB)

    dl.dshape[!is.finite(dl.dshape)] =
      max(dl.dshape[is.finite(dl.dshape)])

    answer <- c(w) * cbind(dl.dscale, dl.dshape) * dthetas.detas
    answer
  }), list( .lshape = lshape , .lscale = lscale,
            .eshape = eshape,  .escale = escale ))),

  weight = eval(substitute(expression({


    run.varcov <- 0
    ind1 <- iam(NA, NA, M = M, both = TRUE, diag = TRUE)
    for (ii in 1:( .nsimEIM )) {
      ysim <- rgenray(n = n, shape = shape, scale = Scale)

      temp1 <- ysim / Scale
      temp2 <- exp(-temp1^2)  # May be 1 if ysim is very close to 0.
      temp3 <- temp1^2 / Scale
      AAA   <- 2 * temp1^2 / Scale  # 2 * y^2 / Scale^3
      BBB   <- -expm1(-temp1^2)     # denominator
      dl.dshape <- 1/shape + log1p(-temp2)
      dl.dscale <- -2 / Scale + AAA * (1 - (shape - 1) * temp2 / BBB)

      dl.dshape[!is.finite(dl.dshape)] <- max(
      dl.dshape[is.finite(dl.dshape)])

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

    wz <- if (intercept.only)
        matrix(colMeans(run.varcov, na.rm = FALSE),
               n, ncol(run.varcov), byrow = TRUE) else run.varcov
    wz <- wz * dthetas.detas[, ind1$row] * dthetas.detas[, ind1$col]
    c(w) * wz
  }), list( .lshape = lshape , .lscale = lscale,
            .eshape = eshape,  .escale = escale,
            .tol12 = tol12, .nsimEIM = nsimEIM ))))
}










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


  N <- max(length(x), length(scale), length(shape))
  if (length(x)      != N) x      <- rep_len(x,     N)
  if (length(scale)  != N) scale  <- rep_len(scale, N)
  if (length(shape)  != N) shape  <- rep_len(shape, N)
  logdensity <- rep_len(log(0), N)

  if (any(xok <- (x > 0))) {
    temp1 <- -x[xok] / scale[xok]
    logdensity[xok] <- -log(scale[xok]) + log1p(-shape[xok]) +
                       temp1 - 2 * log1p(-shape[xok] * exp(temp1))
  }

  logdensity[(scale <= 0) | (shape <= 0) | (shape >= 1)] <- NaN
  if (log.arg) {
    logdensity
  } else {
     exp(logdensity)
  }
}


pexpgeom <- function(q, scale = 1, shape) {
  temp1 <- -q / scale
  ans <- -expm1(temp1) / (1 - shape * exp(temp1))
  ans[q <= 0] <- 0
  ans[(scale <= 0) | (shape <= 0) | (shape >= 1)] <- NaN
  ans
}


qexpgeom <- function(p, scale = 1, shape) {
  ans <- (-scale) * log((p - 1) / (p * shape - 1))
  ans[(scale <= 0) | (shape <= 0) | (shape >= 1)] <- NaN
  ans[p < 0] <- NaN
  ans[p > 1] <- NaN
  ans[p == 0] <- 0
  ans[p == 1] <- Inf
  ans
}


rexpgeom <- function(n, scale = 1, shape) {
  ans <- qexpgeom(runif(n), shape = shape, scale = scale)
  ans[(scale <= 0) | (shape <= 0) | (shape >= 1)] <- NaN
  ans
}






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


 expgeometric <- function(lscale = "loglink", lshape = "logitlink",
                          iscale = NULL,   ishape = NULL,
                          tol12 = 1.0e-05, zero = 1,
                          nsimEIM = 400) {


  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 (length(ishape))
    if (!is.Numeric(ishape, positive = TRUE) || any(ishape >= 1))
      stop("bad input for argument 'ishape'")

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



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


  new("vglmff",
  blurb = c("Exponential geometric distribution\n\n",
            "Links:    ",
            namesof("scale", lscale, earg = escale), ", ",
            namesof("shape", lshape, earg = eshape), "\n",
            "Mean:     ", "(shape - 1) * log(1 - ",
            "shape) / (shape / scale)"),

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

  infos = eval(substitute(function(...) {
    list(M1 = 2,
         Q1 = 1,
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = c("scale", "shape"),
         nsimEIM = .nsimEIM ,
         lscale = .lscale ,
         lshape = .lshape ,
         zero = .zero )
  }, list( .zero = zero, .lscale = lscale, .lshape = lshape,
           .nsimEIM = nsimEIM ))),

  initialize = eval(substitute(expression({

    temp5 <-
    w.y.check(w = w, y = y,
              out.wy = TRUE,
              maximize = TRUE)
    w <- temp5$w
    y <- temp5$y



    predictors.names <- c(
      namesof("scale", .lscale , earg = .escale , short = TRUE),
      namesof("shape", .lshape , earg = .eshape , short = TRUE))

    if (!length(etastart)) {

      scale.init <- if (is.Numeric( .iscale , positive = TRUE)) {
                      rep_len( .iscale , n)
                    } else {
                      stats::sd(c(y))  # The papers scale parameter beta
                    }

      shape.init <- if (is.Numeric( .ishape , positive = TRUE)) {
                      rep_len( .ishape , n)
                    } else {
                      rep_len(2 - exp(median(y)/scale.init), n)
                    }
      shape.init[shape.init >= 0.95] <- 0.95
      shape.init[shape.init <= 0.05] <- 0.05


      etastart <-
        cbind(theta2eta(scale.init, .lscale , earg = .escale ),
              theta2eta(shape.init, .lshape , earg = .eshape ))

    }
   }), list( .lscale = lscale, .lshape = lshape,
             .iscale = iscale, .ishape = ishape,
             .escale = escale, .eshape = eshape))),

  linkinv = eval(substitute(function(eta, extra = NULL) {
    Scale <- eta2theta(eta[, 1], .lscale , earg = .escale )
    shape <- eta2theta(eta[, 2], .lshape , earg = .eshape )

    (shape - 1) * log1p(-shape) / (shape / Scale)

  }, list( .lscale = lscale, .lshape = lshape,
           .escale = escale, .eshape = eshape ))),

  last = eval(substitute(expression({
    misc$link <-    c(scale = .lscale , shape = .lshape )

    misc$earg <- list(scale = .escale , shape = .eshape )

    misc$expected <- TRUE
    misc$nsimEIM <- .nsimEIM
    misc$multipleResponses <- FALSE
  }), 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[, 1], .lscale , earg = .escale )
    shape <- eta2theta(eta[, 2], .lshape , earg = .eshape )

    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * dexpgeom(x = y, scale = Scale, shape = shape,
                                 log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .lscale = lscale, .lshape = lshape,
           .escale = escale, .eshape = eshape))),

  vfamily = c("expgeometric"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    Scale <- eta2theta(eta[, 1], .lscale , earg = .escale )
    shape <- eta2theta(eta[, 2], .lshape , earg = .eshape )
    okay1 <- all(is.finite(Scale)) && all(0 < Scale) &&
             all(is.finite(shape)) && all(0 < shape & shape < 1)
    okay1
  }, list( .lscale = lscale, .lshape = lshape,
           .escale = escale, .eshape = eshape))),

  deriv = eval(substitute(expression({
    Scale <- eta2theta(eta[, 1], .lscale , earg = .escale )
    shape <- eta2theta(eta[, 2], .lshape , earg = .eshape )

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

    dscale.deta <- dtheta.deta(Scale, .lscale , earg = .escale )
    dshape.deta <- dtheta.deta(shape, .lshape , earg = .eshape )
    dthetas.detas <- cbind(dscale.deta, dshape.deta)

    answer <- c(w) * cbind(dl.dscale, dl.dshape) * dthetas.detas
    answer
  }), list( .lscale = lscale , .lshape = lshape,
            .escale = escale,  .eshape = eshape ))),

  weight = eval(substitute(expression({









      run.varcov <- 0
      ind1 <- iam(NA, NA, M = M, both = TRUE, diag = TRUE)

      if (length( .nsimEIM )) {
        for (ii in 1:( .nsimEIM )) {
          ysim <- rexpgeom(n, scale=Scale, shape=shape)

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

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

      run.varcov <- run.varcov / .nsimEIM

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

      wz <- wz * dthetas.detas[, ind1$row] *
                 dthetas.detas[, ind1$col]
    }

    c(w) * wz
  }), list( .nsimEIM = nsimEIM ))))
}











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


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

  logdensity <- rep_len(log(0), N)
  if (any(xok <- (x > 0))) {
    temp1 <- -x[xok] / scale[xok]
    logdensity[xok] <- -log(-log(shape[xok])) - log(scale[xok]) +
                       log1p(-shape[xok]) + temp1 -
                       log1p(-(1-shape[xok]) * exp(temp1))
  }

  logdensity[(scale <= 0) | (shape <= 0) | (shape >= 1)] <- NaN
  if (log.arg) {
    logdensity
  } else {
     exp(logdensity)
  }
}


pexplog <- function(q, scale = 1, shape) {
  ans <- 1 - log1p(-(1-shape) * exp(-q / scale)) / log(shape)
  ans[q <= 0] <- 0
  ans[(scale <= 0) | (shape <= 0) | (shape >= 1)] <- NaN
  ans
}



qexplog <- function(p, scale = 1, shape) {


  ans <- -scale * (log1p(-shape^(1.0 - p)) - log1p(-shape))

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



rexplog <- function(n, scale = 1, shape) {
  ans <- qexplog(runif(n), scale = scale, shape = shape)
  ans[(scale <= 0) | (shape <= 0) | (shape >= 1)] <- NaN
  ans
}









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


 explogff <- function(lscale = "loglink", lshape = "logitlink",
                      iscale = NULL,   ishape = NULL,
                      tol12 = 1.0e-05, zero = 1,
                      nsimEIM = 400) {

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

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


  if (length(ishape))
    if (!is.Numeric(ishape, positive = TRUE) ||
        any(ishape >= 1))
      stop("bad input for argument 'ishape'")

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



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


  new("vglmff",
  blurb = c("Exponential logarithmic distribution\n\n",
            "Links:    ",
            namesof("scale", lscale, earg = escale), ", ",
            namesof("shape", lshape, earg = eshape), "\n",
            "Mean:     ", "(-polylog(2, 1 - p) * scale) / log(shape)"),

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

  infos = eval(substitute(function(...) {
    list(M1 = 2,
         Q1 = 1,
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = c("scale", "shape"),
         nsimEIM = .nsimEIM ,
         lscale = .lscale ,
         lshape = .lshape ,
         zero = .zero )
  }, list( .zero = zero, .lscale = lscale, .lshape = lshape,
           .nsimEIM = nsimEIM ))),


  initialize = eval(substitute(expression({

    temp5 <-
    w.y.check(w = w, y = y,
              out.wy = TRUE,
              maximize = TRUE)
    w <- temp5$w
    y <- temp5$y



    predictors.names <- c(
      namesof("scale", .lscale , earg = .escale , short = TRUE),
      namesof("shape", .lshape , earg = .eshape , short = TRUE))

    if (!length(etastart)) {

      scale.init <- if (is.Numeric( .iscale , positive = TRUE)) {
                     rep_len( .iscale , n)
                   } else {
                     stats::sd(c(y))
                   }

      shape.init <- if (is.Numeric( .ishape , positive = TRUE)) {
                     rep_len( .ishape , n)
                   } else {
                      rep_len((exp(median(y)/scale.init) - 1)^2, n)
                   }
      shape.init[shape.init >= 0.95] <- 0.95
      shape.init[shape.init <= 0.05] <- 0.05


      etastart <-
        cbind(theta2eta(scale.init, .lscale , earg = .escale ),
              theta2eta(shape.init, .lshape , earg = .eshape ))

    }
   }), list( .lscale = lscale, .lshape = lshape,
             .iscale = iscale, .ishape = ishape,
             .escale = escale, .eshape = eshape))),

  linkinv = eval(substitute(function(eta, extra = NULL) {
    scale <- eta2theta(eta[, 1], .lscale , earg = .escale )
    shape <- eta2theta(eta[, 2], .lshape , earg = .eshape )



    qexplog(p = 0.5, shape = shape, scale = scale)

  }, list( .lscale = lscale, .lshape = lshape,
           .escale = escale, .eshape = eshape ))),

  last = eval(substitute(expression({
    misc$link <-    c(scale = .lscale , shape = .lshape )

    misc$earg <- list(scale = .escale , shape = .eshape )

    misc$expected <- TRUE
    misc$nsimEIM <- .nsimEIM
    misc$multipleResponses <- FALSE
  }), 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[, 1], .lscale , earg = .escale )
    shape <- eta2theta(eta[, 2], .lshape , earg = .eshape )


    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * dexplog(x = y, scale = Scale,
                                shape = shape, log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .lscale = lscale, .lshape = lshape,
           .escale = escale, .eshape = eshape))),

  vfamily = c("explogff"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    Scale <- eta2theta(eta[, 1], .lscale , earg = .escale )
    shape <- eta2theta(eta[, 2], .lshape , earg = .eshape )
    okay1 <- all(is.finite(Scale)) && all(0 < Scale) &&
             all(is.finite(shape)) && all(0 < shape & shape < 1)
    okay1
  }, list( .lscale = lscale, .lshape = lshape,
           .escale = escale, .eshape = eshape))),

  deriv = eval(substitute(expression({
    Scale <- eta2theta(eta[, 1], .lscale , earg = .escale )
    shape <- eta2theta(eta[, 2], .lshape , earg = .eshape )

     temp2 <- exp(-y / Scale)
     temp3 <- y / Scale^2
     temp4 <- 1 - shape
     dl.dscale <- (-1 / Scale) + temp3 + (temp4 * temp3 *
                  temp2) / (1 - temp4 * temp2)
     dl.dshape <- -1 / (shape * log(shape)) - 1 / temp4 -
                  temp2 / (1 - temp4 * temp2)

    dscale.deta <- dtheta.deta(Scale, .lscale , earg = .escale )
    dshape.deta <- dtheta.deta(shape, .lshape , earg = .eshape )
    dthetas.detas <- cbind(dscale.deta, dshape.deta)

    answer <- c(w) * cbind(dl.dscale, dl.dshape) * dthetas.detas
    answer
  }), list( .lscale = lscale , .lshape = lshape,
            .escale = escale,  .eshape = eshape ))),

  weight = eval(substitute(expression({



        run.varcov <- 0
        ind1 <- iam(NA, NA, M = M, both = TRUE, diag = TRUE)

        if (length( .nsimEIM )) {
          for (ii in 1:( .nsimEIM )) {
            ysim <- rexplog(n, scale = Scale, shape = shape)

            temp2 <- exp(-ysim / Scale)
            temp3 <- ysim / Scale^2
            temp4 <- 1 - shape
            dl.dscale <- (-1 / Scale) + temp3 + (temp4 * temp3 *
                         temp2) / (1 - temp4 * temp2)
            dl.dshape <- -1 / (shape * log(shape)) - 1 / temp4 -
                         temp2 / (1 - temp4 * temp2)

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

          run.varcov <- run.varcov / .nsimEIM

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

          wz <- wz * dthetas.detas[, ind1$row] *
                    dthetas.detas[, ind1$col]
        }

    c(w) * wz
  }), list( .nsimEIM = nsimEIM ))))
}












dweibull3 <- function(x, location = 0, scale = 1, shape,
                      log = FALSE) {

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


  dweibull(x = x - location, shape = shape,
           scale = scale, log = log.arg)
}

pweibull3 <- function(q, location = 0, scale = 1, shape) {
  pweibull(q = q - location, scale = scale, shape = shape)
}


qweibull3 <- function(p, location = 0, scale = 1, shape) {
  location + qweibull(p = p, shape = shape, scale = scale)
}


rweibull3 <- function(n, location = 0, scale = 1, shape) {
  location + rweibull(n = n, shape = shape, scale = scale)
}









   ### Two-piece normal (TPN) family


dtpn <- function(x, location = 0, scale = 1, skewpar = 0.5,
                 log.arg = FALSE) {


  if (any(skewpar <= 0 |
          skewpar >= 1 |
          scale   <= 0 ,
           na.rm = TRUE))
    stop("some parameters out of bound")

  N <- max(length(x), length(location), length(scale),
           length(skewpar))
  if (length(x)        != N) x         <- rep_len(x,        N)
  if (length(scale)    != N) scale     <- rep_len(scale,    N)
  if (length(location) != N) location  <- rep_len(location, N)
  if (length(skewpar)  != N) skewpar   <- rep_len(skewpar,  N)

  zedd <- (x - location) / scale

  log.s1 <-  -zedd^2 / (8 * skewpar^2)
  log.s2 <-  -zedd^2 / (8 * (1 - skewpar)^2)

  logdensity <- log.s1
  logdensity[zedd > 0] <- log.s2[zedd > 0]

  logdensity <- logdensity -log(scale) - log(sqrt(2 * pi))

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

ptpn <- function(q, location = 0, scale = 1, skewpar = 0.5) {

  if (any(skewpar <= 0 |
          skewpar >= 1 |
          scale   <= 0 ,
          na.rm = TRUE))
    stop("some parameters out of bound")


 zedd <- (q - location) / scale

  s1 <- 2 * skewpar * pnorm(zedd, sd = 2 * skewpar)  #/ scale
  s2 <- skewpar + (1 - skewpar) *
        pgamma(zedd^2 / (8 * (1-skewpar)^2), 0.5)

ans <- rep_len(0.0, length(zedd))
ans[zedd <= 0] <- s1[zedd <= 0]
ans[zedd > 0] <- s2[zedd > 0]

ans
}



pos <- function(x) ifelse(x > 0, x, 0.0)


qtpn <- function(p, location = 0, scale = 1, skewpar = 0.5) {

  pp = p
  if (any(pp      <= 0 |
          pp      >= 1 |
          skewpar <= 0 |
          skewpar >= 1 |
          scale   <= 0 ,
             na.rm = TRUE))
    stop("some parameters out of bound")
    # Recycle the vectors to equal lengths
  LLL <- max(length(pp), length(location), length(scale),
             length(skewpar))
  if (length(pp)       != LLL) pp       <- rep_len(pp,       LLL)
  if (length(location) != LLL) location <- rep_len(location, LLL)
  if (length(scale)    != LLL) scale    <- rep_len(scale,    LLL)
  if (length(skewpar)  != LLL) skewpar  <- rep_len(skewpar,  LLL)

  qtpn <- rep_len(NA_real_, length(LLL))
  qtpn <- qnorm(pp / (2 * skewpar), sd = 2 * skewpar)
  qtpn[pp > skewpar] <- sqrt(8 * ( 1 - skewpar)^2 *
                        qgamma(pos( pp - skewpar) / (
                        1 - skewpar),.5))[pp > skewpar]

   qtpn * scale + location

}





rtpn <- function(n, location = 0, scale = 1, skewpar = 0.5) {


  qtpn(p = runif(n), location = location,
       scale = scale, skewpar = skewpar)
}





tpnff <- function(llocation = "identitylink", lscale = "loglink",
                  pp = 0.5, method.init = 1,  zero = 2) {
  if (!is.Numeric(method.init, length.arg = 1,
                  integer.valued = TRUE, positive = TRUE) ||
      method.init > 4)
     stop("argument 'imethod' must be 1 or 2 or 3 or 4")

  if (!is.Numeric(pp, length.arg = 1, positive = TRUE))
    stop("bad input for argument 'pp'")


  llocat <- as.list(substitute(llocation))
  elocat <- link2list(llocat)
  llocat <- attr(elocat, "function.name")

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



  new("vglmff",
  blurb = c("Two-piece normal distribution \n\n",
            "Links: ",
            namesof("location",  llocat,  earg = elocat), ", ",
            namesof("scale",     lscale,  earg = escale), "\n\n",
            "Mean: "),
  constraints = eval(substitute(expression({
    constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
                                predictors.names = predictors.names,
                                M1 = 2)
  }), list( .zero = zero ))),

  infos = eval(substitute(function(...) {
    list(M1 = 2,
         Q1 = 1,
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = c("location", "scale"),
         llocation = .llocat ,
         lscale    = .lscale ,
         zero = .zero )
  }, list( .zero = zero,
           .llocat = llocat,
           .lscale = lscale ))),

  initialize = eval(substitute(expression({

    temp5 <-
    w.y.check(w = w, y = y,
              out.wy = TRUE,
              maximize = TRUE)
    w <- temp5$w
    y <- temp5$y



    predictors.names <-
       c(namesof("location", .llocat , earg = .elocat , tag = FALSE),
         namesof("scale",    .lscale , earg = .escale , tag = FALSE))




    if (!length(etastart)) {
        junk <- lm.wfit(x = x, y = c(y), w = c(w))
        scale.y.est <-
          sqrt( sum(c(w) * junk$resid^2) / junk$df.residual )
        location.init <- if ( .llocat == "loglink")
          pmax(1/1024, y) else {

        if ( .method.init == 3) {
          rep_len(weighted.mean(y, w), n)
        } else if ( .method.init == 2) {
          rep_len(median(rep(y, w)), n)
        } else if ( .method.init == 1) {
          junk$fitted
        } else {
          y
        }
      }
      etastart <- cbind(
           theta2eta(location.init,  .llocat , earg = .elocat ),
           theta2eta(scale.y.est,    .lscale , earg = .escale ))
    }
  }), list( .llocat = llocat, .lscale = lscale,
            .elocat = elocat, .escale = escale,
            .method.init = method.init ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    eta2theta(eta[, 1], .llocat , earg = .elocat )
  }, list( .llocat = llocat,
           .elocat = elocat, .escale = escale ))),
  last = eval(substitute(expression({
    misc$link     <-    c("location" = .llocat , "scale" = .lscale )

    misc$earg     <- list("location" = .elocat , "scale" = .escale )

    misc$expected <- TRUE
    misc$pp       <- .pp
    misc$method.init <- .method.init
    misc$multipleResponses <- FALSE
  }), list( .llocat = llocat, .lscale = lscale,
            .elocat = elocat, .escale = escale,
            .pp     = pp,        .method.init = method.init ))),

  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    location <- eta2theta(eta[, 1], .llocat , earg = .elocat )
    myscale  <- eta2theta(eta[, 2], .lscale , earg = .escale )
    ppay     <- .pp
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * dtpn(y, skewpar = ppay, location = location,
                             scale = myscale, log.arg = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .llocat = llocat, .lscale = lscale,
           .elocat = elocat, .escale = escale,
           .pp      = pp ))),
  vfamily = c("tpnff"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    mylocat <- eta2theta(eta[, 1], .llocat ,  earg = .elocat )
    myscale <- eta2theta(eta[, 2], .lscale ,  earg = .escale )
    okay1 <- all(is.finite(mylocat)) &&
             all(is.finite(myscale)) && all(0 < myscale)
    okay1
  }, list( .llocat = llocat, .lscale = lscale,
           .elocat = elocat, .escale = escale,
           .pp      = pp ))),
  deriv = eval(substitute(expression({
    mylocat <- eta2theta(eta[, 1], .llocat ,  earg = .elocat )
    myscale <- eta2theta(eta[, 2], .lscale ,  earg = .escale )
    mypp    <- .pp

    zedd <- (y - mylocat) / myscale
    cond2 <-    (zedd > 0)

    dl.dlocat        <-  zedd / (4 * mypp^2)  # cond1
    dl.dlocat[cond2] <- (zedd / (4 * (1 - mypp)^2))[cond2]
    dl.dlocat        <- dl.dlocat / myscale

    dl.dscale        <-  zedd^2 / (4 * mypp^2)
    dl.dscale[cond2] <- (zedd^2 / (4 * (1 - mypp)^2))[cond2]
    dl.dscale        <- (-1 + dl.dscale) / myscale


    dlocat.deta <- dtheta.deta(mylocat, .llocat, earg = .elocat)
    dscale.deta <- dtheta.deta(myscale, .lscale, earg = .escale)

    ans <- c(w) * cbind(dl.dlocat * dlocat.deta,
                        dl.dscale * dscale.deta)
    ans
  }), list( .llocat = llocat, .lscale = lscale,
            .elocat = elocat, .escale = escale,
            .pp      = pp ))),
  weight = eval(substitute(expression({
    wz   <- matrix(0, n, M)  # diag matrix; y is one-col too
    temp10 <- mypp * (1 - mypp)
    ned2l.dlocat2        <- 1 / ((4 * temp10) * myscale^2)
    ned2l.dscale2        <- 2 /  myscale^2


    wz[, iam(1, 1, M)] <- ned2l.dlocat2 * dlocat.deta^2
    wz[, iam(2, 2, M)] <- ned2l.dscale2 * dscale.deta^2
  # wz[, iam(3, 3, M)] <- ned2l.dskewpar2 * dskewpa.deta^2
  # wz[, iam(1, 3, M)] <- ned2l.dlocatdskewpar * dskewpar.deta * dlocat.deta
    c(w) * wz
  }))))
}



  ########################################################################


tpnff3 <- function(llocation = "identitylink",
                    lscale   = "loglink",
                    lskewpar = "identitylink",
                    method.init = 1,  zero = 2)
{
  if (!is.Numeric(method.init, length.arg = 1,
                  integer.valued = TRUE, positive = TRUE) ||
      method.init > 4)
    stop("argument 'imethod' must be 1 or 2 or 3 or 4")



  llocat <- as.list(substitute(llocation))
  elocat <- link2list(llocat)
  llocat <- attr(elocat, "function.name")

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

  lskewp <- as.list(substitute(lskewpar))
  eskewp <- link2list(lskewp)
  lskewp <- attr(eskewp, "function.name")






  new("vglmff",
  blurb = c("Two-piece normal distribution \n\n",
            "Links: ",
            namesof("location", llocat, earg = elocat), ", ",
            namesof("scale",    lscale, earg = escale),  ", ",
            namesof("skewpar",  lskewp, earg = eskewp),  "\n\n",
            "Mean: "),
  constraints = eval(substitute(expression({
    constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
                                predictors.names = predictors.names,
                                M1 = 2)
  }), list( .zero = zero ))),

  infos = eval(substitute(function(...) {
    list(M1 = 2,
         Q1 = 1,
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = c("location", "scale", "skewpar"),
         llocation = .llocat ,
         lscale    = .lscale ,
         lskewpar  = .lskewp ,
         zero = .zero )
  }, list( .zero = zero,
           .llocat = llocat,
           .lscale = lscale,
           .lskewp = lskewp ))),


  initialize = eval(substitute(expression({

    temp5 <-
    w.y.check(w = w, y = y,
              out.wy = TRUE,
              maximize = TRUE)
    w <- temp5$w
    y <- temp5$y


    predictors.names <-
       c(namesof("location", .llocat, earg = .elocat, tag = FALSE),
         namesof("scale",    .lscale, earg = .escale, tag = FALSE),
         namesof("skewpar",  .lskewp, earg = .eskewp, tag = FALSE))

    if (!length(etastart)) {
      junk = lm.wfit(x = x, y = c(y), w = c(w))
      scale.y.est <- sqrt(sum(c(w) * junk$resid^2) / junk$df.residual)
      location.init <- if ( .llocat == "loglink") pmax(1/1024, y) else {
        if ( .method.init == 3) {
          rep_len(weighted.mean(y, w), n)
        } else if ( .method.init == 2) {
          rep_len(median(rep(y, w)), n)
        } else if ( .method.init == 1) {
          junk$fitted
        } else {
          y
        }
      }
      skew.l.in <- sum((y < location.init)) / length(y)
      etastart <- cbind(
           theta2eta(location.init, .llocat,   earg = .elocat),
           theta2eta(scale.y.est,   .lscale,   earg = .escale),
           theta2eta(skew.l.in,     .lskewp, earg = .escale))
    }
  }), list( .llocat = llocat, .lscale = lscale, .lskewp = lskewp,
            .elocat = elocat, .escale = escale, .eskewp = eskewp,

            .method.init=method.init ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    eta2theta(eta[, 1], .llocat, earg = .elocat)
  }, list( .llocat = llocat,
           .elocat = elocat, .escale = escale ))),
  last = eval(substitute(expression({
    misc$link     <-     c("location" = .llocat,
                           "scale"    = .lscale,
                           "skewpar"  = .lskewp)

    misc$earg     <-  list("location" = .elocat,
                           "scale"    = .escale,
                           "skewpar"  = .eskewp)

    misc$expected <- TRUE
         misc$method.init <- .method.init
  }), list( .llocat = llocat, .lscale = lscale, .lskewp = lskewp,
            .elocat = elocat, .escale = escale, .eskewp = eskewp,
                    .method.init = method.init ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {

   locat    <- eta2theta(eta[, 1], .llocat , earg = .elocat )
   myscale  <- eta2theta(eta[, 2], .lscale , earg = .escale )
   myskew   <- eta2theta(eta[, 3], .lskewp , earg = .eskewp )

    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * dtpn(y, location = locat,  scale = myscale,
                             skewpar = myskew, log.arg = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
   }
  }, list( .llocat = llocat, .lscale = lscale, .lskewp = lskewp,
           .elocat = elocat, .escale = escale, .eskewp = eskewp
           ))),
  vfamily = c("tpnff3"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    mylocat <- eta2theta(eta[, 1], .llocat , earg = .elocat )
    myscale <- eta2theta(eta[, 2], .lscale , earg = .escale )
    myskew  <- eta2theta(eta[, 3], .lskewp , earg = .eskewp )
    okay1 <- all(is.finite(mylocat)) &&
             all(is.finite(myscale)) && all(0 < myscale) &&
             all(is.finite(myskew ))
    okay1
  }, list( .llocat = llocat, .lscale = lscale, .lskewp = lskewp,
           .elocat = elocat, .escale = escale, .eskewp = eskewp
         ))),
  deriv = eval(substitute(expression({
    mylocat <- eta2theta(eta[, 1], .llocat,   earg = .elocat)
    myscale <- eta2theta(eta[, 2], .lscale,   earg = .escale)
    myskew  <- eta2theta(eta[, 3], .lskewp, earg = .eskewp)


    zedd <- (y - mylocat) / myscale
    cond2 <-    (zedd > 0)

    dl.dlocat        <-  zedd / (4 * myskew^2)  # cond1
    dl.dlocat[cond2] <- (zedd / (4 * (1 - myskew)^2))[cond2]
    dl.dlocat        <- dl.dlocat / myscale

    dl.dscale        <-  zedd^2 / (4 * myskew^2)
    dl.dscale[cond2] <- (zedd^2 / (4 * (1 - myskew)^2))[cond2]
    dl.dscale        <- (-1 + dl.dscale) / myscale

    dl.dskewpar      <-     zedd^2 /  (4 * myskew^3)
    dl.dskewpar[cond2] <- (-zedd^2 /  (4 * (1 - myskew)^3))[cond2]



    dlocat.deta <- dtheta.deta(mylocat, .llocat, earg = .elocat)
    dscale.deta <- dtheta.deta(myscale, .lscale, earg = .escale)
    dskewpar.deta <- dtheta.deta(myskew, .lskewp, earg = .eskewp)
    ans <-
    c(w) * cbind(dl.dlocat * dlocat.deta,
              dl.dscale * dscale.deta,
              dl.dskewpar * dskewpar.deta
              )
    ans
  }), list( .llocat = llocat, .lscale = lscale, .lskewp = lskewp,
            .elocat = elocat, .escale = escale, .eskewp = eskewp
            ))),
  weight = eval(substitute(expression({
    wz <- matrix(NA_real_, n, dimm(M))  # diag matrix; y is one-col too

    temp10 <- myskew * (1 - myskew)

    ned2l.dlocat2        <- 1 / ((4 * temp10) * myscale^2)
    ned2l.dscale2        <- 2 /  myscale^2
    ned2l.dskewpar2      <- 3 / temp10
    ned2l.dlocatdskewpar <- (-2 * sqrt(2)) / (temp10 * sqrt(pi) *
                             myscale)

    wz[, iam(1, 1,M)] <- ned2l.dlocat2 * dlocat.deta^2
    wz[, iam(2, 2,M)] <- ned2l.dscale2 * dscale.deta^2
    wz[, iam(3, 3,M)] <- ned2l.dskewpar2 * dskewpar.deta^2
    wz[, iam(1, 3,M)] <- ned2l.dlocatdskewpar * dskewpar.deta *
                         dlocat.deta

    ans
    c(w) * wz
  }))))
}





dzoabeta <- function(x, shape1, shape2, pobs0 = 0,
                     pobs1 = 0, log = FALSE, tol = .Machine$double.eps) {
  log.arg <- log
  rm(log)
  LLL <- max(length(x), length(shape1),
             length(shape2), length(pobs0), length(pobs1))
  if (LLL != length(x))      x      <- rep_len(x,      LLL)
  if (LLL != length(shape1)) shape1 <- rep_len(shape1, LLL)
  if (LLL != length(shape2)) shape2 <- rep_len(shape2, LLL)
  if (LLL != length(pobs0))  pobs0  <- rep_len(pobs0,  LLL)
  if (LLL != length(pobs1))  pobs1  <- rep_len(pobs1,  LLL)
  ans <- rep_len(NA_real_, LLL)

  k1 <- (pobs0 < -tol | pobs1 < -tol |
    (pobs0 + pobs1) > (1 + tol))
  k4 <- is.na(pobs0) | is.na(pobs1)
  ans[!k4 & !k1] <- dbeta(x[!k4 & !k1],
                          shape1[!k4 & !k1],
                          shape2[!k4 & !k1], log = TRUE) +
                    log1p(-(pobs0[!k4 & !k1] + pobs1[!k4 & !k1]))
  k2 <- x == 0 & pobs0 > 0 & !is.na(x)
  k3 <- x == 1 & pobs1 > 0 & !is.na(x)
  ans[k2 & !k4 & !k1] <- log(pobs0[k2 & !k4 & !k1])
  ans[k3 & !k4 & !k1] <- log(pobs1[k3 & !k4 & !k1])
  if (!log.arg) ans <- exp(ans)
  if (any(k1 & !k4)) {
    ans[k1 & !k4] <- NaN
    warning("NaNs produced")
  }
  ans
}



rzoabeta <- function(n, shape1, shape2, pobs0 = 0, pobs1 = 0,
                     tol = .Machine$double.eps) {
  use.n <- if ((length.n <- length(n)) > 1) {
    length.n
  } else {
    if (!is.Numeric(n, integer.valued = TRUE, length.arg = 1,
                    positive = TRUE)) {
      stop("bad input for argument 'n'")
    } else {
      n
    }
  }
  shape1 <- rep_len(shape1, use.n)
  shape2 <- rep_len(shape2, use.n)
  pobs0  <- rep_len(pobs0,  use.n)
  pobs1  <- rep_len(pobs1,  use.n)
  random.number <- runif(use.n)
  ans <- rep_len(NA_real_, use.n)
  k5 <- (pobs0 < -tol | pobs1 < -tol |
           (pobs0 + pobs1) > (1 + tol))
  k4 <- is.na(pobs0) | is.na(pobs1)
  ans[!k4] <- qzoabeta(random.number[!k4], shape1 = shape1,
                       shape2 = shape2, pobs0 = pobs0,
                       pobs1 = pobs1)
  if (any(k5 & !k4)) {
    ans[k5 & !k4] <- NaN
    warning("NaNs produced")
  }
  ans
}



pzoabeta <- function(q, shape1, shape2, pobs0 = 0, pobs1 = 0,
                     lower.tail = TRUE, log.p = FALSE,
                     tol = .Machine$double.eps) {
  LLL <- max(length(q), length(shape1),
             length(shape2), length(pobs0), length(pobs1))
  if (LLL != length(q))      q      <- rep_len(q,      LLL)
  if (LLL != length(shape1)) shape1 <- rep_len(shape1, LLL)
  if (LLL != length(shape2)) shape2 <- rep_len(shape2, LLL)
  if (LLL != length(pobs0))  pobs0  <- rep_len(pobs0,  LLL)
  if (LLL != length(pobs1))  pobs1  <- rep_len(pobs1,  LLL)
  k3 <- (pobs0 < -tol | pobs1 < -tol |
        (pobs0 + pobs1) > (1 + tol))
  k4 <- is.na(pobs0) | is.na(pobs1)
  ans <- rep_len(NA_real_, LLL)
  ans[!k3 & !k4] <- pbeta(q[!k3 & !k4],
                          shape1[!k3 & !k4],
                          shape2[!k3 & !k4], log.p = TRUE) +
    log1p(-(pobs0[!k3 & !k4] + pobs1[!k3 & !k4]))
  ans <- exp(ans)
  k1 <- q >= 0 & !is.na(q)
  k2 <- q >= 1 & !is.na(q)
  ans[k1 & !k3 & !k4] <- ans[k1 & !k3 & !k4] +
    pobs0[k1 & !k3 & !k4]
  ans[k2 & !k3 & !k4] <- ans[k2 & !k3 & !k4] +
    pobs1[k2 & !k3 & !k4]
  if (!lower.tail & log.p) {
    ans <- log1p(-ans)
  } else {
    if (!lower.tail)
      ans <- 1 - ans
    if (log.p)
      ans <- log(ans)
  }
  if (any(k3 & !k4)) {
    ans[k3 & !k4] <- NaN
    warning("NaNs produced")
  }
  ans
}



qzoabeta <- function(p, shape1, shape2, pobs0 = 0, pobs1 = 0,
                     lower.tail = TRUE, log.p = FALSE,
                     tol = .Machine$double.eps) {
  LLL <- max(length(p), length(shape1),
             length(shape2), length(pobs0), length(pobs1))
  if (LLL != length(p))      p      <- rep_len(p,      LLL)
  if (LLL != length(shape1)) shape1 <- rep_len(shape1, LLL)
  if (LLL != length(shape2)) shape2 <- rep_len(shape2, LLL)
  if (LLL != length(pobs0))  pobs0  <- rep_len(pobs0,  LLL)
  if (LLL != length(pobs1))  pobs1  <- rep_len(pobs1,  LLL)
  k0 <- (pobs0 < -tol | pobs1 < -tol |
        (pobs0 + pobs1) > (1 + tol))
  k4 <- is.na(pobs0) | is.na(pobs1)
  ans <- rep_len(NA_real_, LLL)
  if (!lower.tail & log.p) {
    p <- -expm1(p)
  } else{
    if (!lower.tail)
      p <- 1 - p
    if (log.p) {
      p <- exp(p)
    }
  }
  k1 <- p >= 0 & p <= pobs0 & !is.na(p)
  k2 <- p > pobs0 & p < (1 - pobs1) & !is.na(p)
  k3 <- p >= (1 - pobs1) & p <= 1 & !is.na(p)
  ans[k1 & !k0 & !k4] <- 0
  ans[k2 & !k0 & !k4] <-
    qbeta((p[k2 & !k0 & !k4] -
           pobs0[k2 & !k0 & !k4]) / (1 - pobs0[k2 & !k0 & !k4] -
           pobs1[k2 & !k0 & !k4]),
           shape1 = shape1[k2 & !k0 & !k4],
           shape2 = shape2[k2 & !k0 & !k4])
  ans[k3 & !k0 & !k4] <- 1
  if (any(k0 & !k4)) {
    ans[k3 & !k4] <- NaN
    warning("NaNs produced")
  }
  ans
}





log1mexp <- function(x) {
  if (any(x < 0 & !is.na(x)))
    stop("Inputs need to be non-negative!")
  ifelse(x <= log(2), log(-expm1(-x)), log1p(-exp(-x)))
}


log1pexp <- function(x){

  ifelse(x <= -37, exp(x),
         ifelse(x <= 18, log1p(exp(x)),
                ifelse(x <= 33, x + exp(-x), x)))
}






dzoibetabinom.ab <- function(x, size, shape1, shape2, pstr0 = 0,
                             pstrsize = 0, log = FALSE) {
  log.arg <- log
  rm(log)
  LLL <- max(length(x), length(size), length(shape1),
             length(shape2), length(pstr0), length(pstrsize))
  if (LLL != length(x))        x        <- rep_len(x,        LLL)
  if (LLL != length(size))     size     <- rep_len(size,     LLL)
  if (LLL != length(shape1))   shape1   <- rep_len(shape1,   LLL)
  if (LLL != length(shape2))   shape2   <- rep_len(shape2,   LLL)
  if (LLL != length(pstr0))    pstr0    <- rep_len(pstr0,    LLL)
  if (LLL != length(pstrsize)) pstrsize <- rep_len(pstrsize, LLL)
  ans <- rep_len(NA_real_, LLL)
  k1 <- pstr0 < 0 | pstrsize < 0 |
           (pstr0 + pstrsize) > 1
  k <- is.na(size) | is.na(shape1) | is.na(shape2) |
    is.na(pstr0) | is.na(pstrsize) | is.na(x)
  if (sum(!k & !k1) > 0) {
    ans[!k & !k1] <-
      dbetabinom.ab(x[!k & !k1], size[!k & !k1], shape1[!k & !k1],
                    shape2[!k & !k1], log = TRUE) +
      log1p(-(pstr0[!k & !k1]+pstrsize[!k & !k1]))
    if (!log.arg) ans <- exp(ans)
  }
  k2 <- x == 0 & pstr0 > 0
  k3 <- x == size & pstrsize > 0
  if (sum(k2 & !k & !k1) > 0)
    ans[k2 & !k & !k1] <- pstr0[k2 & !k & !k1] +
      ans[k2 & !k & !k1]
  if (sum(k3 & !k & !k1) > 0)
    ans[k3 & !k & !k1] <- pstrsize[k3 & !k & !k1] +
      ans[k3 & !k & !k1]
  if (any(k1 & !k)) {
    ans[k1 & !k] <- NaN
    warning("NaNs produced")
  }
  ans
}



dzoibetabinom <- function(x, size, prob, rho = 0, pstr0 = 0,
                          pstrsize = 0, log = FALSE) {
  dzoibetabinom.ab(x, size, shape1 = prob * (1 - rho) / rho,
                   shape2 = (1 - prob) * (1 - rho) / rho,
                   pstr0 = pstr0, pstrsize = pstrsize, log = log)
}



rzoibetabinom.ab <- function(n, size, shape1, shape2,
                             pstr0 = 0, pstrsize = 0) {
  use.n <- if ((length.n <- length(n)) > 1) {
    length.n
  } else {
    if (!is.Numeric(n, integer.valued = TRUE, length.arg = 1,
                    positive = TRUE)) {
      stop("bad input for argument 'n'")
    } else {
      n
    }
  }
  size     <- rep_len(size,     use.n)
  shape1   <- rep_len(shape1,   use.n)
  shape2   <- rep_len(shape2,   use.n)
  pstr0    <- rep_len(pstr0,    use.n)
  pstrsize <- rep_len(pstrsize, use.n)
  ans      <- rep_len(NA_real_, use.n)
  k <- is.na(size) | is.na(shape1) | is.na(shape2) |
       is.na(pstr0) | is.na(pstrsize)
  k1 <- pstr0 < 0 | pstrsize < 0 |
    (pstr0 + pstrsize) > 1
  random.number <- runif(use.n)
  k2 <- random.number[!k] < pstr0[!k]
  k3 <- pstr0[!k] <= random.number[!k] &
    random.number[!k] <= (1 - pstrsize[!k])
  k4 <- (1 - pstrsize[!k]) < random.number[!k]
  if (sum(k2 & !k1 & !k) > 0)
    ans[k2 & !k1 & !k] <- 0
  if (sum(k3 & !k1 & !k) > 0)
    ans[k3 & !k1 & !k] <- rbetabinom.ab(sum(k3 & !k1 & !k),
                                        size =  size[k3 & !k1 & !k],
                                        shape1 = shape1[k3 & !k1 & !k],
                                        shape2 = shape2[k3 & !k1 & !k])
  if (sum(k4 & !k1 & !k) > 0)
    ans[k4 & !k1 & !k] <- size[k4 & !k1 & !k]
  ans
}



rzoibetabinom <- function(n, size, prob, rho = 0, pstr0 = 0,
                          pstrsize = 0) {
  rzoibetabinom.ab(n, size, shape1 = prob * (1 - rho) / rho,
                   shape2 = (1 - prob) * (1 - rho) / rho,
                   pstr0 = pstr0,
                   pstrsize = pstrsize)
}



pzoibetabinom.ab <- function(q, size, shape1, shape2, pstr0 = 0,
                             pstrsize = 0, lower.tail = TRUE,
                             log.p = FALSE) {
  LLL <- max(length(q), length(size), length(shape1),
             length(shape2), length(pstr0), length(pstrsize))
  if (LLL != length(q))        q        <- rep_len(q,        LLL)
  if (LLL != length(size))     size     <- rep_len(size,     LLL)
  if (LLL != length(shape1))   shape1   <- rep_len(shape1,   LLL)
  if (LLL != length(shape2))   shape2   <- rep_len(shape2,   LLL)
  if (LLL != length(pstr0))    pstr0    <- rep_len(pstr0,    LLL)
  if (LLL != length(pstrsize)) pstrsize <- rep_len(pstrsize, LLL)
  ans <- rep_len(NA_real_, LLL)
  k <- is.na(size) | is.na(shape1) | is.na(shape2) |
    is.na(pstr0) | is.na(pstrsize) | is.na(q)
  k1 <- pstr0 < 0 | pstrsize < 0 |
    (pstr0 + pstrsize) > 1
  if (sum(!k1 & !k) > 0)
    ans[!k & !k1] <-
      pbetabinom.ab(q[!k & !k1], size[!k & !k1],
                    shape1[!k & !k1], shape2[!k & !k1], log.p = TRUE) +
      log1p(-(pstr0[!k & !k1] + pstrsize[!k & !k1]))
  ans <- exp(ans)
  k2 <- q >= 0
  k3 <- q >= size
  if (sum(k2 & !k1 & !k) > 0)
    ans[k2 & !k & !k1] <- ans[k2 & !k & !k1] +
      pstr0[k2 & !k & !k1]
  if (sum(k3 & !k1 & !k) > 0)
    ans[k3 & !k & !k1] <- ans[k3 & !k & !k1] +
      pstrsize[k3 & !k & !k1]
  if (!lower.tail & log.p) {
    ans <- log1p(-ans)
  } else {
    if (!lower.tail)
      ans <- 1 - ans
    if (log.p)
      ans <- log(ans)
  }
  if (any(!k & k1)) {
    ans[!k & k1] <- NaN
    warning("NaNs produced")
  }
  ans
}


pzoibetabinom <- function(q, size, prob, rho,
                          pstr0 = 0, pstrsize = 0,
                        lower.tail = TRUE, log.p = FALSE) {
  pzoibetabinom.ab(q, size, shape1 = prob * (1 - rho) / rho,
                 shape2 = (1 - prob) * (1 - rho) / rho,
                 pstr0 = pstr0, pstrsize = pstrsize,
                 lower.tail = lower.tail, log.p = log.p)
}














  AR1EIM<- function(x = NULL,
                    var.arg   = NULL,
                    p.drift   = NULL,
                    WNsd      = NULL,
                    ARcoeff1  = NULL,
                    eps.porat = 1e-2) {

    if (!is.matrix(x))
      stop("Argument 'x' must be a matrix.")

    yy   <- x
    M    <- 3
    nn   <- nrow(x)
    nn0  <- numeric(0)
    NOS  <- ncol(x)

    if (!is.matrix(WNsd))
      WNsd <- matrix(WNsd, nrow = nn, ncol = NOS, byrow = TRUE)

    if (!is.matrix(ARcoeff1))
      ARcoeff1 <- matrix(ARcoeff1, nrow = nn, ncol = NOS, byrow = TRUE)

    if (!is.Numeric(eps.porat, length.arg = 1) || eps.porat < 0 ||
        eps.porat > 1e-2)
      stop("Bad input for argument 'eps.porat'.")

    sdTSR   <- colMeans(WNsd)
    sdTSv   <- colMeans(WNsd)
    drift.v <- rep(p.drift, NOS)[1:NOS]
    Aux11   <- (NOS > 1)
    the1v   <- colMeans(ARcoeff1)
    JFin    <- array(0.0, dim = c(nn, NOS, M + (M - 1) + (M - 2) ))

    for (spp in 1:NOS) {

      x <- yy[, spp]
      the1    <- the1v[spp]
      drift.p <- drift.v[spp]
      sdTS    <- sdTSv[spp]

      r <- numeric(nn)
      r <- AR1.gammas(x = x, lags = nn - 1)
      r[nn] <- r[1]

      s0 <- numeric(nn)
      s1 <- numeric(nn)
      s1 <- if (var.arg) (the1^(0:(nn - 1))) / (1 - the1^2) else
             2 * (the1^(0:(nn - 1))) * sdTS / (1 - the1^2)

      s2    <- numeric(nn)
      help1 <- c(0:(nn - 1))
      s2    <- help1 * (the1^(help1 - 1)) * (sdTS^2) / (1 - the1^2) +
                   2 * (sdTS^2) * (the1^(help1 + 1)) / (1 - the1^2)^2
      sMat <- cbind(s0, s1, s2)

      J  <- array(NA_real_,
                  dim = c(length(the1) + 2, length(the1) + 2, nn))
      Jp <- array(NA_real_,
                  dim = c(length(the1) + 2, length(the1) + 2, nn))

      alpha    <- numeric(nn)
      alpha[1] <- 1
      delta    <- r[1]
      eta      <- matrix(NA_real_, nrow = nn, ncol = M)
      eta[1, ] <- cbind(s0[1], s1[1], s2[1])

      psi      <- matrix(0, nrow = nn, ncol = length(the1) + 2)
      psi[1, ] <- cbind(s0[1], s1[1], s2[1]) / r[1]

      u0  <- rep(1/(1 - sign(the1v[1]) * min(0.975, abs(the1v[1]))), nn )
      u1  <- rep(drift.p/(1 - the1)^2, nn)
      uMat <- cbind(u0, rep(0, nn), u1)

      aux1 <- matrix(sMat[1, ],
                     nrow = 2 + length(the1),
                     ncol = 2 + length(the1), byrow = TRUE)
      diag(aux1) <- sMat[1, ]
      J[, , 1]   <- Jp[, , 1] <- aux1 * t(aux1) / (2 * r[1]^2)
      J[1, 1, 1] <- Jp[1, 1, 1] <- 1 / sdTS^2
      JFin[1, spp, 1:M] <- Jp[, , 1][row(Jp[, , 1]) == col(Jp[, , 1])]
      Neps.porat <- 1.819*eps.porat*(1e-10)

      dk   <- matrix(NA_real_, nrow = 1, ncol = length(the1) + 2)
      eR   <- matrix(NA_real_, nrow = 1, ncol = length(the1) + 2)
      cAux2 <- d55 <- numeric(nn); d55[1] <- 0.1

      for (jay in 1:(nn - 1)) {

        cAux <- as.numeric(alpha[1:jay] %*%
                    r[2:(jay + 1)][length(r[2:(jay + 1)]):1])/delta

        dk <- alpha[1:jay] %*%
         sMat[2:(jay + 1), , drop = FALSE][length(sMat[2:(jay + 1)]):1, ]

        delta        <- delta * (1 - cAux^2)
        d55[jay + 1] <- cAux^2

        if ((d55[jay + 1] < eps.porat*1e-2) || (jay > 1e1)) {
          nn0 <- jay
          break
        }

        eta[jay + 1, ] <- dk
        tAux <- numeric(jay + 1)
        tAux <- alpha[1:(jay + 1)] -
                      cAux * alpha[1:(jay + 1)][(jay + 1):1]
        alpha[1:(jay + 1)] <- tAux[1:(jay + 1)]

        eR <- alpha[1:(jay + 1)][(jay + 1):1] %*%
          eta[1:(jay + 1), , drop = FALSE]

        tAux <- eta[1:(jay + 1), ] -
          cAux * eta[1:(jay + 1), ][(jay + 1):1, ]

        eta[1:(jay + 1), ] <- tAux

        AuxE <- matrix(eR, nrow = jay + 1, ncol = M, byrow = TRUE)
        Aux3 <- matrix(alpha[1:(jay + 1)][(jay + 1):1],
                       nrow = jay + 1, ncol = M, byrow = FALSE)
        Aux4 <- matrix(alpha[1:(jay + 1)],
                       nrow = jay + 1, ncol = M, byrow = FALSE)
                  tAux <- psi[1:(jay + 1), ] -
                      cAux * psi[1:(jay + 1), ][(jay + 1):1, ] +
                         AuxE * (Aux3 - cAux * Aux4) / delta

        if (any(dim(psi[1:(jay + 1), ])) != any(dim(tAux)) )
          stop("Invalids 'psi' and 'tAux'.")

        psi[1:(jay + 1), ] <- tAux
        fk <- alpha[1:(jay + 1)] %*% eta[1:(jay + 1), ]
        gk <- alpha[1:(jay + 1)][(jay + 1):1] %*% uMat[1:(jay + 1), ]

        Auxf <- matrix(fk, nrow = M, ncol = M, byrow = FALSE)
        Auxg <- matrix(gk, nrow = M, ncol = M, byrow = FALSE)
        J[, , jay + 1] <-
          J[, , jay] + t(eta[1:(jay + 1), ]) %*% psi[1:(jay + 1), ] /
                      delta - 0.5 * Auxf * t(Auxf) / delta^2 +
                          Auxg * t(Auxg) / delta

        Jp[, , jay + 1] <- J[, , jay + 1] - J[, , jay]
        JFin[jay + 1, spp , 1:M ] <-
          Jp[, , jay + 1][col(Jp[, , jay + 1]) == row(Jp[, , jay + 1])]

        helpC <- numeric(0)
        for (kk in 1:(M - 1))  {
          TF1 <- ( col(Jp[, , jay + 1]) >= row(Jp[, , jay + 1]) )
          TF2 <- (abs(col(Jp[, , jay + 1]) - row(Jp[, , jay + 1])) == kk )
          helpC <- c(helpC, Jp[, , jay + 1][TF1 & TF2])
        }
        rm(TF1, TF2)

        JFin[jay + 1, spp , -(1:M) ] <- helpC
      }

      if (length(nn0))
        for (kk in nn0:(nn - 1)) {
          J[, , kk + 1] <- J[, , nn0] + (kk - nn0 + 1) * Jp[, , nn0]
          Jp[, , kk + 1] <- J[, , kk + 1] - J[, , kk]

          JFin[kk + 1, spp , 1:M ] <-
            Jp[, , kk + 1][col(Jp[, , kk + 1]) == row(Jp[, , kk + 1])]

          helpC <- numeric(0)
          for (ll in 1:(M - 1))  {
           TF1 <- ( col(Jp[, , kk + 1]) >= row(Jp[, , kk + 1]) )
           TF2 <- (abs(col(Jp[, , kk + 1]) - row(Jp[, , kk + 1])) == ll)
            helpC <- c(helpC, Jp[, , kk + 1][TF1 & TF2])
          }
          rm(TF1, TF2)
          JFin[kk + 1, spp , -(1:M) ] <- helpC
        }
      JFin[which(JFin <= Neps.porat)] <-
        abs( JFin[which(JFin <= Neps.porat)])
    }

    JFin

  } # End






AR1.gammas <- function(x, y = NULL, lags = 1) {
  xx  <- matrix(x, ncol = 1)
  nx  <- nrow(xx)

  if (lags < 0 || !(is.Numeric(lags, integer.valued = TRUE)))
    stop("'lags' must be a positive integer.")

  if (length(y)) {
    yy <- matrix(y, ncol = 1)
    ny <- nrow(yy)
    if (nx != ny)
      stop("Number of rows differs.") else
        n <- nx
  } else {
    yy <- xx
    n  <- nrow(xx)
  }

  myD <- numeric(lags + 1)
  myD[1] <- if (length(y)) cov(xx, yy) else cov(xx, xx)  # i.e. var(xx)
  if (lags > 0)
    for (ii in 1:lags)
      myD[ii + 1]  <- cov(xx[-(1:ii), 1], yy[1:(n - ii) , 1])

  myD
}





dzipfmb <- function(x, shape, start = 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(start), length(shape))
  if (length(x)     != LLL) x     <- rep_len(x,     LLL)
  if (length(start) != LLL) start <- rep_len(start, LLL)
  if (length(shape) != LLL) shape <- rep_len(shape, LLL)
  
  bad0 <- !is.finite(shape) | !is.finite(start) |
          shape <= 0 | 1 <= shape |
          start != round(start) | start < 1
  bad <- bad0 | !is.finite(x) | x < start | x != round(x)

  logpdf <- x + shape + start
  
  if (any(!bad)) {
    logpdf[!bad] <- lbeta(    x[!bad] - shape[!bad], shape[!bad] + 1) - 
                    lbeta(start[!bad] - shape[!bad], shape[!bad])
  }
  
  logpdf[!bad0 & is.infinite(x)] <- log(0)
  logpdf[!bad0 & x < start     ] <- log(0)
  logpdf[!bad0 & x != round(x) ] <- log(0)
  logpdf[ bad0] <- NaN

  if (log.arg) logpdf else exp(logpdf)
} # dzipfmb()



pzipfmb <-
  function(q, shape, start = 1, 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'")
  
  LLL <- max(length(shape), length(q), length(start))
  if (length(q)     != LLL) q     <- rep_len(q,     LLL)
  if (length(shape) != LLL) shape <- rep_len(shape, LLL)
  if (length(start) != LLL) start <- rep_len(start, LLL)
  
  q.use <- pmax(start, floor(q + 1))
  
  bad0 <- !is.finite(shape) | !is.finite(start) |
          shape <= 0 | 1 <= shape |
          start != round(start) | start < 1
  bad <- bad0 | !is.finite(q.use)

  ans <- q + shape + start

  log.S.short <- lgamma(start[!bad]) + lgamma(q.use[!bad] - shape[!bad]) -
                 lgamma(q.use[!bad]) - lgamma(start[!bad] - shape[!bad])
  
  ans[!bad] <-
  if (lower.tail) {
    if (log.p) { # lower.tail = T, log.p = T
      log1p(-exp(log.S.short))
    } else {  # lower.tail = T, log.p = F
      -expm1(log.S.short)
    }
  } else { 
    if (log.p) {  # lower.tail = F, log.p = T
      log.S.short
    } else {      # lower.tail = F, log.p = F
      exp(log.S.short)
    }
  }

  ans[!bad0 & is.infinite(q.use) & start < q.use] <-
    if (lower.tail) {if (log.p) log(1) else 1} else
    {if (log.p) log(0) else 0}

  ans[bad0] <- NaN
  ans
}  # pzipfmb()






qzipfmb <- function(p, shape, start = 1) {

  LLL <- max(length(p), length(shape), length(start))
  if (length(p)     != LLL) p     <- rep_len(p,     LLL)
  if (length(shape) != LLL) shape <- rep_len(shape, LLL)
  if (length(start) != LLL) start <- rep_len(start, LLL)
  ans <- p + shape + start

  bad0 <- !is.finite(shape) | !is.finite(start) |
          shape <= 0 | 1 <= shape |
          start != round(start) | start < 1
  bad <- bad0 | !is.finite(p) | p <= 0 | 1 <= p

  lo <- rep_len(start, LLL) - 0.5
  approx.ans <- lo  # True at lhs
  hi <- 2 * lo + 10.5
  dont.iterate <- bad
  done <- dont.iterate | p <= pzipfmb(hi, shape, start = start)
  iter <- 0
  max.iter <- round(log2(.Machine$double.xmax)) - 2
  max.iter <- round(log2(1e300)) - 2
  while (!all(done) && iter < max.iter) {
    lo[!done] <- hi[!done]
    hi[!done] <- 2 * hi[!done] + 10.5  # Bug fixed
    done[!done] <- (p[!done] <= pzipfmb(hi[!done],
                    shape = shape[!done], start[!done]))
    iter <- iter + 1
  }

  foo <- function(q, shape, start, p)
    pzipfmb(q, shape = shape, start) - p

  lhs <- dont.iterate |
         p <= dzipfmb(start, shape = shape, start)
  
  approx.ans[!lhs] <-
    bisection.basic(foo, lo[!lhs], hi[!lhs], tol = 1/16,
                    shape = shape[!lhs],
                    start = start[!lhs], p = p[!lhs])
  faa <- floor(approx.ans[!lhs])
  tmp <-
    ifelse(pzipfmb(faa, shape = shape[!lhs], start[!lhs]) < p[!lhs] &
           p[!lhs] <= pzipfmb(faa+1, shape = shape[!lhs], start[!lhs]),
           faa+1, faa)
  ans[!lhs] <- tmp

  vecTF <- !bad0 & !is.na(p) &
           p <= dzipfmb(start, shape = shape, start)
  ans[vecTF] <- start[vecTF]

  ans[!bad0 & !is.na(p) & p == 0] <- start[!bad0 & !is.na(p) & p == 0]
  ans[!bad0 & !is.na(p) & p == 1] <- Inf
  ans[!bad0 & !is.na(p) & p <  0] <- NaN
  ans[!bad0 & !is.na(p) & p >  1] <- NaN
  ans[ bad0] <- NaN
  ans
}  # qzipfmb



rzipfmb <- function(n, shape, start = 1) {
  qzipfmb(runif(n), shape, start = start)
}  # rzipfmb












dextlogF <- function(x, lambda, tau,
                   location = 0, scale = 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(lambda), length(tau),
             length(location), length(scale))
  if (length(x)        != LLL) x        <- rep_len(x,        LLL)
  if (length(lambda)   != LLL) lambda   <- rep_len(lambda,   LLL)
  if (length(tau)      != LLL) tau      <- rep_len(tau,      LLL)
  if (length(location) != LLL) location <- rep_len(location, LLL)
  if (length(scale)    != LLL) scale    <- rep_len(scale,    LLL)
  
  bad0 <- !is.finite(lambda) | !is.finite(tau) | 
          !is.finite(location) | !is.finite(scale) | 
          lambda <= 0 | tau < 0 | 1 < tau | scale <= 0
  bad <- bad0 | !is.finite(x)

  logpdf <- x + lambda + tau + location + scale
  
  if (any(!bad)) {
    zedd <- (x[!bad] - location[!bad]) / scale[!bad]
    logpdf[!bad] <- (1 - tau[!bad]) * zedd -
      lambda[!bad] * log1p(exp(zedd / lambda[!bad])) -
      lbeta((1 - tau[!bad]) * lambda[!bad],
                 tau[!bad]  * lambda[!bad]) -
      log(lambda[!bad] * scale[!bad])
  }
  
  logpdf[!bad0 & is.infinite(x)] <- log(0)
  logpdf[ bad0] <- NaN

  if (log.arg) logpdf else exp(logpdf)
}  # dextlogf()









extlogF1.control <-
  function(stepsize = 0.5,
           maxit = 100,
           ...) {
  list(stepsize = stepsize,
       maxit    = maxit)
}



 extlogF1 <-
  function(tau = c(0.25, 0.5, 0.75),  # NULL,  # \in (0, 1)
           parallel = TRUE ~ 0,  # FALSE,
           seppar = 0,
           tol0 = -0.001,  # Negative means relative, + means absolute
           llocation = "identitylink",
           ilocation = NULL,
           lambda.arg = NULL,  # 0.1,  # NULL means an adaptive value
           scale.arg = 1,  # Best to leave this alone
           ishrinkage = 0.95,
           digt = 4,
           idf.mu = 3,
           imethod = 1) {


  apply.parint.locat <- FALSE


  if (!is.Numeric(seppar, length.arg = 1) || seppar < 0)
    stop("bad input for argument 'seppar'")
  if (!is.Numeric(tol0, length.arg = 1))  # || tol0 < 0
    stop("bad input for argument 'tol0'")

  if (!is.Numeric(tau, positive = TRUE))
    stop("bad input for argument 'tau'")
  if (any(1 <= tau))
    stop("argument 'tau' must have values in (0, 1)")
  if (length(tau) > 1 && any(diff(tau) <= 0))
    stop("argument 'tau' must be an increasing sequence")

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


  llocation <- llocation
  llocat <- as.list(substitute(llocation))
  elocat <- link2list(llocat)
  llocat <- attr(elocat, "function.name")

  ilocat <- ilocation


  if (!is.Numeric(ishrinkage, length.arg = 1) ||
     ishrinkage < 0 ||
     ishrinkage > 1)
    stop("bad input for argument 'ishrinkage'")
  if (!is.Numeric(scale.arg, positive = TRUE))
    stop("bad input for argument 'scale.arg'")




  new("vglmff",
  blurb = c("One-parameter extended log-F distribution\n\n",
            "Links:      ",
            namesof("location", llocat, earg = elocat),
            "\n", "\n",
            "Quantiles:  location(tau)"),

  constraints = eval(substitute(expression({

    onemat <- matrix(1, M, 1)
    constraints.orig <- constraints


    cm1.locat <- diag(M)
    cmk.locat <- onemat
    con.locat <- cm.VGAM(cmk.locat,
                         x = x, bool = .parallel ,
                         constraints = constraints.orig,
                         apply.int   = .apply.parint.locat ,
                         cm.default           = cm1.locat,
                         cm.intercept.default = cm1.locat)
    
    constraints <- con.locat

  }), list( .parallel = parallel, .seppar = seppar, .tol0 = tol0,
            .apply.parint.locat = apply.parint.locat ))),



  infos = eval(substitute(function(...) {
    list(M1 = NA,  # 1,
         Q1 = NA,  # 1,
         tau        = .tau ,
         scale.arg  = .scale.arg ,
         lambda.arg = .lambda.arg ,
         seppar     = .seppar ,
         tol0.arg   = as.vector( .tol0 ),  # Could be negative
         digt       = .digt ,
         expected   = TRUE,
         multipleResponses = TRUE,
         parameters.names  = "location")
  }, list( .lambda.arg = lambda.arg, .scale.arg = scale.arg,
           .digt = digt,
           .tau = tau, .seppar = seppar, .tol0 = tol0 ))),
  initialize = eval(substitute(expression({

    temp5 <-
    w.y.check(w = w, y = y,
              ncol.w.max = if (length( .tau ) > 1) 1 else Inf,
              ncol.y.max = if (length( .tau ) > 1) 1 else Inf,
              out.wy = TRUE,
              maximize = TRUE)
    w <- temp5$w
    y <- temp5$y



    extra$ncoly <- ncoly <- ncol(y)
    if ((ncoly > 1) && (length( .tau ) > 1 ||
        length( .scale.arg ) > 1))
      stop("response must be a vector if 'tau' or 'scale.arg' ",
           "has a length greater than one")



    if ((midspread <- diff(quantile(y, probs = c(0.25, 0.75)))) == 0)
      stop("could not work out an adaptive 'lambda'") else                
    lambda <- if (is.null( .lambda.arg )) {
    muxfactor <- 25
    extra$midspread <- midspread
         as.vector((midspread / .scale.arg) / muxfactor)
    } else {
      as.vector( .lambda.arg )
    }
    extra$lambda.arg <- lambda
    tol0 <- as.vector( .tol0 )  # A negative value means relative
    if (tol0 < 0)
      tol0 <- as.vector(abs(tol0) * midspread)
    extra$tol0 <- tol0  # A positive value means absolute

    extra$M <- M <- max(length( .scale.arg ),
                        ncoly,
                        length( .tau ))  # Recycle
    extra$scale.arg <- rep_len( .scale.arg , M)
    extra$tau       <- rep_len( .tau       , M)
    extra$n <- n


    extra$tau.names <- tau.names <-
      paste("(tau = ", round(extra$tau, digits = .digt), ")", sep = "")
    extra$Y.names <- Y.names <- if (ncoly > 1) dimnames(y)[[2]] else "y"
    if (is.null(Y.names) || any(Y.names == ""))
      extra$Y.names <- Y.names <- paste("y", 1:ncoly, sep = "")
    extra$y.names <- y.names <-
      if (ncoly > 1) paste(Y.names, tau.names, sep = "") else tau.names

    extra$individual <- FALSE

    mynames1 <- param.names("location", M, skip1 = TRUE)
    predictors.names <-
        c(namesof(mynames1, .llocat , earg = .elocat , tag = FALSE))

    seppar <- as.vector( .seppar )
    if (seppar > 0 && M > 1) {
    }  # seppar

    locat.init <- matrix(0, n, M)
    if (!length(etastart)) {

      for (jay in 1:M) {
        y.use <- if (ncoly > 1) y[, jay] else y
        locat.init[, jay] <- if ( .imethod == 1) {
           quantile(y.use, probs = extra$tau[jay])
        } else if ( .imethod == 2) {
          weighted.mean(y.use, w[, min(jay, ncol(w))])
        } else if ( .imethod == 3) {
          median(y.use)
        } else if ( .imethod == 4) {
          Fit5 <- vsmooth.spline(x = x[, min(ncol(x), 2)],
                                 y = y.use, w = w, df = .idf.mu )
          c(predict(Fit5, x = x[, min(ncol(x), 2)])$y)
        } else {
          use.this <- weighted.mean(y.use, w[, min(jay, ncol(w))])
          (1 - .ishrinkage ) * y.use + .ishrinkage * use.this
        }

        if (length( .ilocat )) {
          locat.init <- matrix( .ilocat , n, M, byrow = TRUE)
        }

        if ( .llocat == "loglink") locat.init <- abs(locat.init)
        etastart <-
          cbind(theta2eta(locat.init, .llocat , earg = .elocat ))
      }
    }
    }), list( .imethod = imethod, .seppar = seppar, .tol0 = tol0,
              .idf.mu = idf.mu, .lambda.arg = lambda.arg,
              .ishrinkage = ishrinkage, .digt = digt,
              .elocat = elocat, .scale.arg = scale.arg,
              .llocat = llocat, .tau = tau,
              .ilocat = ilocat ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    locat <- eta2theta(eta, .llocat , earg = .elocat )
    if (length(locat) > extra$n)
      dimnames(locat) <- list(dimnames(eta)[[1]], extra$y.names)
    locat
  }, list( .elocat = elocat,
           .llocat = llocat, .scale.arg = scale.arg,
           .tau = tau, .lambda.arg = lambda.arg ))),
  last = eval(substitute(expression({
    misc$link <- setNames(rep_len( .llocat , M), mynames1)

    misc$earg <- vector("list", M)
    names(misc$earg) <- names(misc$link)
    for (ii in 1:M) {
      misc$earg[[ii]] <- ( .elocat )
    }

    extra$eCDF <- numeric(M)
    locat <- as.matrix(locat)
    for (ii in 1:M) {
      y.use <- if (ncoly > 1) y[, ii] else y
      extra$eCDF[ii] <-
        100 * weighted.mean(y.use <= locat[, ii], w[, min(ii, ncol(w))])
    }
    names(extra$eCDF) <- y.names


    extra$scale.arg <- ( .scale.arg )
    }), list( .elocat = elocat,
              .llocat = llocat, .scale.arg = scale.arg, .tau = tau ))),

  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL, summation = TRUE) {
    ymat <- matrix(y, extra$n, extra$M)
    locat <- eta2theta(eta, .llocat , earg = .elocat )
    taumat <- matrix(extra$tau,       extra$n, extra$M, byrow = TRUE)
    Scale  <- matrix(extra$scale.arg, extra$n, extra$M, byrow = TRUE)
    lambda  <- extra$lambda.arg

    Ans <- if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * dextlogF(x = c(ymat), location = c(locat),
                                 scale = c(Scale), tau = c(taumat),
                                 lambda = lambda, log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }

    n <- nrow(eta)
    M <- NCOL(eta)
    seppar <- as.vector( .seppar )
    tol0 <- extra$tol0  # Absolute
    if (seppar > 0 && M > 1) {
      negdiffmat <- as.matrix(locat[, -ncol(locat)] - locat[, -1])
      negdiffmat <- matrix(pmax(0, tol0 + negdiffmat)^3, n, M-1)
      Adjustment <- seppar * sum(negdiffmat)
      Ans <- Ans - Adjustment
    }  # seppar

    Ans 
  }, list( .elocat = elocat, .scale.arg = scale.arg, .tau = tau,
           .llocat = llocat, .seppar = seppar ))),
  vfamily = c("extlogF1"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    locat <- eta2theta(eta, .llocat , earg = .elocat )
    okay1 <- all(is.finite(locat))
    okay1
  }, list( .elocat = elocat, .scale.arg = scale.arg, .tau = tau,
           .llocat = llocat ))),


  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)
    extra  <- object@extra
    locat  <- eta2theta(eta, .llocat , .elocat )
    Scale  <- matrix(extra$scale.arg, extra$n, extra$M, byrow = TRUE)
    taumat <- matrix(extra$tau,       extra$n, extra$M, byrow = TRUE)
    lambda <- extra$lambda.arg
    relogF(nsim * length(Scale), location = c(locat),
           lambda = lambda, scale = c(Scale), tau = c(taumat))
  }, list( .elocat = elocat, .scale.arg = scale.arg, .tau = tau,
           .llocat = llocat ))),

  deriv = eval(substitute(expression({
    seppar <- as.vector( .seppar )
    tol0 <- extra$tol0  # Absolute
    locat <- eta2theta(eta, .llocat , earg = .elocat )

    ymat <- matrix(y, n, M)
    Scale <- matrix(extra$scale.arg, n, M, byrow = TRUE)

    lambda <- extra$lambda.arg

    taumat <- matrix(extra$tau, n, M, byrow = TRUE)
    zedd <- (ymat - locat) / Scale

    dl.dlocat <- (taumat - 1 + plogis(zedd / lambda)) / Scale
    dlocat.deta <- dtheta.deta(locat, .llocat , earg = .elocat )

    sep.penalize <- seppar > 0 && M > 1
    if (sep.penalize) {
      dl.dlocat[, 1] <- dl.dlocat[, 1] - 3 * seppar *
                        pmax(0, tol0 + locat[, 1] - locat[, 2])^2
      dl.dlocat[, M] <- dl.dlocat[, M] + 3 * seppar *
                        pmax(0, tol0 + locat[, M-1] - locat[, M])^2
      if (M > 2)
      for (jay in 2:(M-1))
        dl.dlocat[, jay] <- dl.dlocat[, jay] - 3 * seppar * (
          pmax(0, tol0 + locat[, jay  ] - locat[, jay+1])^2 -
          pmax(0, tol0 + locat[, jay-1] - locat[, jay  ])^2)
    }  # seppar

    c(w) * cbind(dl.dlocat * dlocat.deta)
  }), list( .elocat = elocat, .seppar = seppar,
            .llocat = llocat, .scale.arg = scale.arg, .tau = tau ))),

  weight = eval(substitute(expression({
    ned2l.dlocat2 <- taumat * (1 - taumat) / ((1 + lambda) * Scale^2)

    if (sep.penalize) {
      ned2l.dlocat2[, 1] <- ned2l.dlocat2[, 1] + 6 * seppar *
                            pmax(0, tol0 + locat[,   1] - locat[, 2])
      ned2l.dlocat2[, M] <- ned2l.dlocat2[, M] + 6 * seppar *
                            pmax(0, tol0 + locat[, M-1] - locat[, M])
      if (M > 2)
        for (jay in 2:(M-1))
          ned2l.dlocat2[, jay] <-
          ned2l.dlocat2[, jay] + 6 * seppar * (
            pmax(0, tol0 + locat[, jay  ] - locat[, jay+1]) +
            pmax(0, tol0 + locat[, jay-1] - locat[, jay  ]))
      rhs.mat <- -6 * seppar * pmax(0, tol0 + locat[, -M] - locat[, -1])
      rhs.mat <- matrix(c(rhs.mat), n, M-1)  # Right dimension
    }  # seppar

    rhs.mat <- if (sep.penalize)
      rhs.mat * dlocat.deta[, -M] * dlocat.deta[, -1] else NULL
    wz <- cbind(ned2l.dlocat2 * dlocat.deta^2,
                rhs.mat)
    c(w) * wz
  }), list( .elocat = elocat, .scale.arg = scale.arg,
            .llocat = llocat ))))
}  # extlogF1()








 setClass("extlogF1",                 contains = "vglmff")



setMethod("showvglmS4VGAM",
          signature(VGAMff = "extlogF1"),
  function(object,
           VGAMff,
           ...) {
  cat("\nQuantiles:\n")
  print(eCDF(as(object, "vglm")))
  invisible(object)
  })




setMethod("showsummaryvglmS4VGAM",  signature(VGAMff = "extlogF1"),
  function(object,
           VGAMff,
           ...) {
  cat("Quantiles:\n")
  print(eCDF(as(object, "vglm"), all = TRUE))
  cat("\n")

  invisible(object)
})








is.crossing.vglm <- function(object, ...) {
  if (!is(object, "vglm"))
    stop("the object is not a VGLM")
  if (any(object@family@vfamily == "lms.bcn"))
    return(FALSE)
  if (!any(object@family@vfamily == "extlogF1"))
    stop("the object does not have a 'extlogF1' family function")
  locat <- fitted(object)
  if (ncol(locat) == 1) {
    TRUE
  } else {
    diffmat <- locat[, -1] - locat[, -ncol(locat)]
    crossq <- any(diffmat < 0)
    crossq
  }
}  # is.crossing.vglm



if (!isGeneric("is.crossing"))
  setGeneric("is.crossing", function(object, ...)
             standardGeneric("is.crossing"),
             package = "VGAM")


setMethod("is.crossing",  "vglm", function(object, ...)
    is.crossing.vglm(object, ...))







 fix.crossing.vglm <-
     function(object, maxit = 100, trace = FALSE,  # etastart = NULL,
              ...) {
  if (!is.crossing(object) ||
      any(object@family@vfamily == "lms.bcn"))
    return(object)

  if (!any(object@family@vfamily == "extlogF1"))
    stop("the object does not have a 'extlogF1' family function")

  object.save <- object
  M <- npred(object.save)
  if (M > 1 &
      !all(is.parallel(object.save)) &
      nrow(coef(object.save, matrix = TRUE)) > 1) {
    if (!has.intercept(object.save))
      stop("the object must have an intercept term")
    for (jay in 1:M) {  # M is an upperbound
      Hlist <- constraints(object, type = "term")
      locat <- fitted(object) 
      diffmat <- locat[,           -1, drop = FALSE] -
                 locat[, -ncol(locat), drop = FALSE]
      crossq <- any(diffmat < 0)
      if (crossq) {
        min.index <- which.min(apply(diffmat, 2, min))
        Hk <- Hlist[[2]]  # Next covariate past the intercept
        use.min.index <- 1 + which(cumsum(colSums(Hk)) == min.index)
        if (ncol(Hk) == 1) break
        Hk[, use.min.index - 1] <- Hk[, use.min.index - 1] +
                                   Hk[, use.min.index]
        Hk <- Hk[, -use.min.index, drop = FALSE]
        for (kay in 2:length(Hlist)) {  # Omit the intercept
          Hlist[[kay]] <- Hk
        }  # kay
      } else {  # crossq
        break
      }
      if (is.numeric(maxit))
        object@control$maxit <- maxit  # Overwrite this possibly
      if (is.logical(trace))
        object@control$trace <- trace  # Overwrite this possibly
      object <-
        vglm(formula = as.formula(formula(object)),
             family = extlogF1(tau = object.save@extra$tau,
                             lambda.arg = object.save@extra$lambda.arg,
                             scale.arg = object.save@extra$scale.arg,
                             llocation = linkfun(object.save)[1],
                             parallel = FALSE),
             data = get(object.save@misc$dataname),
             control = object@control,  # maxit, trace, etc.
             constraints = Hlist)  # New constraints; new object too
    }  # jay
  }  # M > 1 and Hk not all parallel

  object  # A new object without any crossing quantiles
}  # fix.crossing



if (!isGeneric("fix.crossing"))
  setGeneric("fix.crossing", function(object, ...)
             standardGeneric("fix.crossing"),
             package = "VGAM")


setMethod("fix.crossing",  "vglm", function(object, ...)
    fix.crossing.vglm(object, ...))







if (FALSE)
qtplot.extlogF1 <- function(lambda,
                          tau = c(0.25, 0.5, 0.75),
                          location = 0, scale = 1,
                          eta = NULL) {

    
  lp <- length(tau)
  lambda <- rep(lambda, length = lp)
  answer <- matrix(NA_real_, nrow(eta), lp,
                   dimnames = list(dimnames(eta)[[1]],
                                   as.character(tau)))
  for (ii in 1:lp) {
    answer[, ii] <- qtlogF(tau      = tau[ii],
                           lambda   = lambda[ii],
                           location = location,  # eta[, ii],
                           scale    = scale)
  }
  answer
}








 eCDF.vglm <-
     function(object, all = FALSE,
              ...) {

  okayfuns <- c("lms.bcn", "extlogF1")
  if (!any(object@family@vfamily %in% okayfuns))
    stop("the object does not have an empirical CDF defined on it")

  
  if (any(object@family@vfamily == "extlogF1")) {
    Ans <- object@extra$eCDF / 100
    if (all) {
      Ans <- cbind(Ans, object@extra$tau)
      rownames(Ans) <- rep(" ", NROW(Ans))  # NULL  #
      colnames(Ans) <- c("ecdf", "tau")
    }
  }  # "extlogF1"


  if (any(object@family@vfamily == "lms.bcn")) {
    Ans <- numeric(length(object@extra$percentiles))
    locat <- as.matrix(fitted(object))
    M <- npred(object)
    y <- depvar(object)
    w <- weights(object, type = "prior")
    for (ii in 1:M) {
      y.use <- if (ncol(y) > 1) y[, ii] else y
      Ans[ii] <- 1 *  # Was 100, but now 1
        weighted.mean(y.use <= locat[, ii], w[, min(ii, ncol(w))])
    }
    if (all) {
      Ans <- cbind(Ans, object@extra$percentiles / 100)  # tau, really
      rownames(Ans) <- rep(" ", NROW(Ans))  # NULL  #
      colnames(Ans) <- c("ecdf", "tau")
    } else {
      digt <- 4  # extlogF1() default
      tau.names <-
        paste("(tau = ", round(object@extra$percentiles / 100,
                               digits = digt), ")",
            sep = "")
      Y.names <- if (ncol(y) > 1) dimnames(y)[[2]] else "y"
      if (is.null(Y.names) || any(Y.names == ""))
        Y.names <- paste("y", 1:ncol(y), sep = "")
      y.names <- if (ncol(y) > 1)
        paste(Y.names, tau.names, sep = "") else tau.names
      names(Ans) <- y.names
    }
  }  # "lms.bcn"

  Ans
}  # eCDF.vglm



if (!isGeneric("eCDF"))
  setGeneric("eCDF", function(object, ...)
             standardGeneric("eCDF"),
             package = "VGAM")


setMethod("eCDF",  "vglm", function(object, ...)
    eCDF.vglm(object, ...))




         

Try the VGAM package in your browser

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

VGAM documentation built on Jan. 16, 2021, 5:21 p.m.