R/family.aunivariate.R

Defines functions diffzeta rdiffzeta qdiffzeta pdiffzeta ddiffzeta diffzeta rdiffzeta qdiffzeta pdiffzeta ddiffzeta zipf rzipf qzipf pzipf dzipf gharmonic gharmonic2 zetaff rzeta qzeta pzeta dzeta topple rtopple qtopple ptopple dtopple zoabetaR betaprime betaR betaff truncgeometric qbenf qbenf pbenf rbenf dbenf logF dlogF nefghs dnefghs slash slash.control pslash dslash dpoislindley lindley rlind plind dlind yulesimon yulesimon.control ryules qyules pyules dyules skellam skellam.control rskellam dskellam riceff riceff.control qrice price marcumQ rrice drice kumar pkumar qkumar rkumar dkumar rhzeta qhzeta phzeta dhzeta hzeta hzeta.control

Documented in betaff betaprime betaR dbenf ddiffzeta dhzeta diffzeta dkumar dlind dlogF drice dskellam dslash dtopple dyules dzeta dzipf gharmonic gharmonic2 hzeta kumar lindley logF marcumQ pbenf pdiffzeta phzeta pkumar plind price pslash ptopple pyules pzeta pzipf qbenf qdiffzeta qhzeta qkumar qrice qtopple qyules qzeta qzipf rbenf rdiffzeta rhzeta riceff rkumar rlind rrice rskellam rtopple ryules rzeta rzipf skellam slash topple truncgeometric yulesimon zetaff zipf zoabetaR

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














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



hzeta <-
  function(lshape = "logloglink", ishape = NULL, nsimEIM = 100) {

  stopifnot(ishape > 0)
  stopifnot(nsimEIM > 10,
            length(nsimEIM) == 1,
            nsimEIM == round(nsimEIM))



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


  new("vglmff",
  blurb = c("Haight's Zeta distribution f(y) = (2y-1)^(-shape)",
            " - (2y+1)^(-shape),\n",
            "    shape>0, y = 1, 2,....\n\n",
            "Link:    ",
            namesof("shape", lshape, earg = eshape), "\n\n",
            "Mean:     (1-2^(-shape)) * zeta(shape) if shape>1",
            "\n",
        "Variance: (1-2^(1-shape)) * zeta(shape-1) - mean^2 if",
        " shape>2"),
  infos = eval(substitute(function(...) {
    list(M1 = 1,
         Q1 = 1,
         dpqrfun = "hzeta",
         expected = FALSE,
         multipleResponses = FALSE,
         parameters.names = c("shape"),
         lshape = .lshape ,
         nsimEIM = .nsimEIM )
  }, list( .nsimEIM = nsimEIM, .lshape = lshape ))),

  rqresslot = eval(substitute(
    function(mu, y, w, eta, extra = NULL) {
    Shape <- eta2theta(eta, .lshape , earg = .eshape )
    scrambleseed <- runif(1)  # To scramble the seed
    qnorm(runif(length(y),
                phzeta(y - 1, Shape),
                phzeta(y    , Shape)))
  }, list( .lshape = lshape, .eshape = eshape ))),


  initialize = eval(substitute(expression({

    w.y.check(w = w, y = y,
              Is.integer.y = TRUE,
              Is.positive.y = TRUE)


    predictors.names <-
      namesof("shape", .lshape , earg = .eshape , tag = FALSE)

    if (!length(etastart)) {
      a.init <- if (length( .ishape)) .ishape else {
        if ((meany <- weighted.mean(y, w)) < 1.5) 3.0 else
        if (meany < 2.5) 1.4 else 1.1
      }
      a.init <- rep_len(a.init, n)
      etastart <- theta2eta(a.init, .lshape , earg = .eshape )
    }
  }), list( .lshape = lshape, .eshape = eshape,
            .ishape = ishape ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    mu <- (1-2^(-shape)) * zeta(shape)
    mu[shape <= 1] <- Inf
    mu
  }, list( .lshape = lshape, .eshape = eshape ))),
  last = eval(substitute(expression({
    misc$link <-    c(shape = .lshape)

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

    misc$nsimEIM <- .nsimEIM

  }), list( .lshape = lshape, .eshape = eshape,
            .nsimEIM = nsimEIM ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * dhzeta(x = y, shape = shape, log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .lshape = lshape, .eshape = eshape ))),
  vfamily = c("hzeta"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    okay1 <- all(is.finite(shape)) && all(0 < shape)
    okay1
  }, list( .lshape = lshape, .eshape = eshape ))),




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

    pwts <- if (length(pwts <- object@prior.weights) > 0)
              pwts else weights(object, type = "prior")
    if (any(pwts != 1))
      warning("ignoring prior weights")
    eta <- predict(object)
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    rhzeta(nsim * length(shape), shape = shape)
  }, list( .lshape = lshape, .eshape = eshape ))),



  deriv = eval(substitute(expression({
    shape <- eta2theta(eta, .lshape , earg = .eshape )

    dshape.deta <- dtheta.deta(shape, .lshape , earg = .eshape )

    d3 <- deriv3(~ log((2*y-1)^(-shape) - (2*y+1)^(-shape)),
                 "shape", hessian = FALSE)
    eval.d3 <- eval(d3)

    dl.dshape <-  attr(eval.d3, "gradient")

    c(w) * dl.dshape * dshape.deta
  }), list( .lshape = lshape, .eshape = eshape ))),
  weight = eval(substitute(expression({
    sd3 <- deriv3(~ log((2*ysim-1)^(-shape) - (2*ysim+1)^(-shape)),
                  "shape", hessian = FALSE)
    run.var <- 0
    for (ii in 1:( .nsimEIM )) {
      ysim <- rhzeta(n, shape = shape)
      eval.sd3 <- eval(sd3)
      dl.dshape <-  attr(eval.d3, "gradient")
      rm(ysim)
      temp3 <- dl.dshape
      run.var <- ((ii-1) * run.var + temp3^2) / ii
    }
    wz <- if (intercept.only)
        matrix(colMeans(cbind(run.var)),
               n, dimm(M), byrow = TRUE) else cbind(run.var)

    wz <- wz * dshape.deta^2
    c(w) * wz
  }), list( .nsimEIM = nsimEIM ))))
}  # hzeta




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

  if (!is.Numeric(shape, positive = TRUE))
    stop("'shape' must be numeric and have positive values")

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

  ox <- !is.finite(x)
  zero <- ox | round(x) != x | x < 1
  ans <- rep_len(0, nn)
  ans[!zero] <- (2*x[!zero]-1)^(-shape[!zero]) -
                (2*x[!zero]+1)^(-shape[!zero])
  if (log.arg) log(ans) else ans
}



 phzeta <- function(q, shape, log.p = FALSE) {


  nn <- max(length(q), length(shape))
  q     <- rep_len(q,     nn)
  shape <- rep_len(shape, nn)
  oq <- !is.finite(q)
  zero <- oq | q < 1
  q <- floor(q)
  ans <- 0 * q
  ans[!zero] <- 1 - (2*q[!zero]+1)^(-shape[!zero])

  ans[q == -Inf] <- 0  # 20141215 KaiH
  ans[q ==  Inf] <- 1  # 20141215 KaiH

  ans[shape <= 0] <- NaN
  if (log.p) log(ans) else ans
}  # phzeta



 qhzeta <- function(p, shape) {

  if (!is.Numeric(p, positive = TRUE) ||
      any(p >= 1))
    stop("argument 'p' must have values inside the interval (0,1)")

  nn <- max(length(p), length(shape))
  p     <- rep_len(p,     nn)
  shape <- rep_len(shape, nn)
  ans <- (((1 - p)^(-1/shape) - 1) / 2)  # p is in (0,1)
  ans[shape <= 0] <- NaN
  floor(ans + 1)
}  # qhzeta


 rhzeta <- function(n, shape) {


  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

  shape <- rep_len(shape, use.n)
  ans <- (runif(use.n)^(-1/shape) - 1) / 2
  ans[shape <= 0] <- NaN
  floor(ans + 1)
}  # rhzeta












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



  N <- max(length(x), length(shape1), length(shape2))
  if (length(x)      != N) x      <- rep_len(x,      N)
  if (length(shape1) != N) shape1 <- rep_len(shape1, N)
  if (length(shape2) != N) shape2 <- rep_len(shape2, N)

  logdensity <- rep_len(log(0), N)
  xok <- (0 <= x & x <= 1)
  logdensity[xok] <- log(shape1[xok]) + log(shape2[xok]) +
                     (shape1[xok] - 1) * log(x[xok]) +
                     (shape2[xok] - 1) * log1p(-x[xok]^shape1[xok])

  logdensity[shape1 <= 0] <- NaN
  logdensity[shape2 <= 0] <- NaN
  if (log.arg) logdensity else exp(logdensity)
}  # dkumar



rkumar <- function(n, shape1, shape2) {
  ans <- (1 - (runif(n))^(1/shape2))^(1/shape1)
  ans[(shape1 <= 0) | (shape2 <= 0)] <- NaN
  ans
}



qkumar <- function(p, shape1, shape2,
                   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 <- (-expm1((1/shape2) * log(-expm1(ln.p))))^(1/shape1)
      ans[ln.p > 0] <- NaN
    } else {
      ans <- (-expm1((1/shape2) * log1p(-p)))^(1/shape1)
      ans[p < 0] <- NaN
      ans[p == 0] <- 0
      ans[p == 1] <- 1
      ans[p > 1] <- NaN
    }
  } else {
    if (log.p) {
      ln.p <- p
      ans <- (-expm1(ln.p / shape2))^(1/shape1)
      ans[ln.p > 0] <- NaN
      ans
    } else {
      ans <- (-expm1((1/shape2) * log(p)))^(1/shape1)
      ans[p < 0] <- NaN
      ans[p == 0] <- 1
      ans[p == 1] <- 0
      ans[p > 1] <- NaN
    }
  }
  ans[(shape1 <= 0) | (shape2 <= 0)] = NaN
  ans
}  # qkumar



pkumar <- function(q, shape1, shape2,
                   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(shape2 * log1p(-q^shape1)))
      ans[q <= 0 ] <- -Inf
      ans[q >= 1] <- 0
    } else {
      ans <- -expm1(shape2 * log1p(-q^shape1))
      ans[q <= 0] <- 0
      ans[q >= 1] <- 1
    }
  } else {
    if (log.p) {
      ans <- shape2 * log1p(-q^shape1)
      ans[q <= 0] <- 0
      ans[q >= 1] <- -Inf
    } else {
      ans <- exp(shape2 * log1p(-q^shape1))
      ans[q <= 0] <- 1
      ans[q >= 1] <- 0
    }
  }

  ans[(shape1 <= 0) | (shape2 <= 0)] <- NaN
  ans
}  # pkumar











 kumar <-
  function(lshape1 = "loglink", lshape2 = "loglink",
           ishape1 = NULL,   ishape2 = NULL,
           gshape1 = exp(2*ppoints(5) - 1), tol12 = 1.0e-4,
           zero = NULL) {
  lshape1 <- as.list(substitute(lshape1))
  eshape1 <- link2list(lshape1)
  lshape1 <- attr(eshape1, "function.name")
  lshape2 <- as.list(substitute(lshape2))
  eshape2 <- link2list(lshape2)
  lshape2 <- attr(eshape2, "function.name")

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

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


  new("vglmff",
  blurb = c("Kumaraswamy distribution\n\n",
            "Links:    ",
            namesof("shape1", lshape1, eshape1, tag = FALSE), ", ",
            namesof("shape2", lshape2, eshape2, tag = FALSE), "\n",
            "Mean:     shape2 * beta(1 + 1 / shape1, shape2)"),
 constraints = eval(substitute(expression({
    constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
                     M = M, M1 = 2,
                     predictors.names = predictors.names)

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

  infos = eval(substitute(function(...) {
    list(M1 = 2,
         Q1 = 1,
         dpqrfun = "kumar",
         expected = TRUE,
         multipleResponses = TRUE,
         parameters.names = c("shape1", "shape2"),
         lshape1 = .lshape1 ,
         lshape2 = .lshape2 ,
         zero = .zero )
  }, list( .zero = zero,
           .lshape1 = lshape1, .lshape2 = lshape2 ))),

  rqresslot = eval(substitute(
    function(mu, y, w, eta, extra = NULL) {
    shape1 <- eta2theta(eta[, c(TRUE, FALSE)], .lshape1 , .eshape1 )
    shape2 <- eta2theta(eta[, c(FALSE, TRUE)], .lshape2 , .eshape2 )
    scrambleseed <- runif(1)  # To scramble the seed
    qnorm(pkumar(y, shape1, shape2))
  }, list( .lshape1 = lshape1, .lshape2 = lshape2,
           .eshape1 = eshape1, .eshape2 = eshape2 ))),

  initialize = eval(substitute(expression({
    checklist <- w.y.check(w = w, y = y, Is.positive.y = TRUE,
                           ncol.w.max = Inf, ncol.y.max = Inf,
                           out.wy = TRUE, colsyperw = 1,
                           maximize = TRUE)
    w <- checklist$w
    y <- checklist$y  # Now 'w' and 'y' have the same dimension.
    if (any((y <= 0) | (y >= 1)))
      stop("the response must be in (0, 1)")

    extra$ncoly <- ncoly <- ncol(y)
    extra$M1 <- M1 <- 2
    M <- M1 * ncoly
    mynames1 <- param.names("shape1", ncoly, skip1 = TRUE)
    mynames2 <- param.names("shape2", ncoly, skip1 = TRUE)
    predictors.names <-
        c(namesof(mynames1, .lshape1 , .eshape1 , tag = FALSE),
          namesof(mynames2, .lshape2 , .eshape2 , tag = FALSE))[
          interleave.VGAM(M, M1 = M1)]

    if (!length(etastart)) {
      kumar.Loglikfun <- function(shape1, y, x, w, extraargs) {
        mediany <- colSums(y * w) / colSums(w)
        shape2 <- log(0.5) / log1p(-(mediany^shape1))
        sum(c(w) * dkumar(y, shape1, shape2, log = TRUE))
      }

      shape1.grid <- as.vector( .gshape1 )
      shape1.init <- if (length( .ishape1 )) .ishape1 else
        grid.search(shape1.grid, objfun = kumar.Loglikfun,
                    y = y,  x = x, w = w)
      shape1.init <- matrix(shape1.init, n, ncoly, byrow = TRUE)

      mediany <- colSums(y * w) / colSums(w)
      shape2.init <- if (length( .ishape2 )) .ishape2 else
        log(0.5) / log1p(-(mediany^shape1.init))
      shape2.init <- matrix(shape2.init, n, ncoly, byrow = TRUE)

      etastart <-
        cbind(theta2eta(shape1.init, .lshape1 , earg = .eshape1 ),
              theta2eta(shape2.init, .lshape2 , earg = .eshape2 ))[,
                  interleave.VGAM(M, M1 = M1)]
    }
  }), list( .lshape1 = lshape1, .lshape2 = lshape2,
            .ishape1 = ishape1, .ishape2 = ishape2,
            .eshape1 = eshape1, .eshape2 = eshape2,
            .gshape1 = gshape1 ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    shape1 <- eta2theta(eta[, c(TRUE, FALSE)], .lshape1 , .eshape1 )
    shape2 <- eta2theta(eta[, c(FALSE, TRUE)], .lshape2 , .eshape2 )
    shape2 * (base::beta(1 + 1/shape1, shape2))
  }, list( .lshape1 = lshape1, .lshape2 = lshape2,
           .eshape1 = eshape1, .eshape2 = eshape2 ))),
  last = eval(substitute(expression({
    misc$link <- c(rep_len( .lshape1 , ncoly),
                   rep_len( .lshape2 , ncoly))[
        interleave.VGAM(M, M1 = M1)]
    temp.names <- c(mynames1, mynames2)[interleave.VGAM(M, M1 = M1)]
    names(misc$link) <- temp.names

    misc$earg <- vector("list", M)
    names(misc$earg) <- temp.names
    for (ii in 1:ncoly) {
      misc$earg[[M1*ii-1]] <- .eshape1
      misc$earg[[M1*ii  ]] <- .eshape2
    }
  }), list( .lshape1 = lshape1, .lshape2 = lshape2,
            .eshape1 = eshape1, .eshape2 = eshape2 ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL, summation = TRUE) {
    shape1 <- eta2theta(eta[, c(TRUE, FALSE)], .lshape1 , .eshape1 )
    shape2 <- eta2theta(eta[, c(FALSE, TRUE)], .lshape2 , .eshape2 )
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * dkumar(x = y, shape1, shape2, log = TRUE)
      if (summation) sum(ll.elts) else ll.elts
    }
  }, list( .lshape1 = lshape1, .lshape2 = lshape2,
           .eshape1 = eshape1, .eshape2 = eshape2 ))),
  vfamily = c("kumar"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    shape1 <- eta2theta(eta[, c(TRUE, FALSE)], .lshape1 , .eshape1 )
    shape2 <- eta2theta(eta[, c(FALSE, TRUE)], .lshape2 , .eshape2 )
    okay1 <- all(is.finite(shape1)) && all(0 < shape1) &&
             all(is.finite(shape2)) && all(0 < shape2)
    okay1
  }, list( .lshape1 = lshape1, .lshape2 = lshape2,
           .eshape1 = eshape1, .eshape2 = eshape2 ))),
  simslot = eval(substitute(
  function(object, nsim) {
    eta <- predict(object)
    shape1 <- eta2theta(eta[, c(TRUE, FALSE)], .lshape1 , .eshape1 )
    shape2 <- eta2theta(eta[, c(FALSE, TRUE)], .lshape2 , .eshape2 )
    rkumar(nsim * length(shape1), shape1 = shape1, shape2 = shape2)
  }, list( .lshape1 = lshape1, .lshape2 = lshape2,
           .eshape1 = eshape1, .eshape2 = eshape2 ))),
  deriv = eval(substitute(expression({
    shape1 <- eta2theta(eta[, c(TRUE, FALSE)], .lshape1 , .eshape1 )
    shape2 <- eta2theta(eta[, c(FALSE, TRUE)], .lshape2 , .eshape2 )
    dshape1.deta <- dtheta.deta(shape1, link = .lshape1 , .eshape1 )
    dshape2.deta <- dtheta.deta(shape2, link = .lshape2 , .eshape2 )
    dl.dshape1 <- 1 / shape1 + log(y) - (shape2 - 1) * log(y) *
                  (y^shape1) / (1 - y^shape1)
    dl.dshape2 <- 1 / shape2 + log1p(-y^shape1)
    dl.deta <- c(w) * cbind(dl.dshape1 * dshape1.deta,
                            dl.dshape2 * dshape2.deta)
    dl.deta[, interleave.VGAM(M, M1 = M1)]
  }), list( .lshape1 = lshape1, .lshape2 = lshape2,
            .eshape1 = eshape1, .eshape2 = eshape2 ))),
  weight = eval(substitute(expression({
    ned2l.dshape11 <- (1 + (shape2 / (shape2 - 2)) *
      ((digamma(shape2) -  digamma(2))^2 -
      (trigamma(shape2) - trigamma(2)))) / shape1^2
    ned2l.dshape22 <- 1 / shape2^2
    ned2l.dshape12 <-
       (digamma(2) - digamma(1 + shape2)) / ((shape2 - 1) * shape1)

    index1 <- (abs(shape2 - 1) < .tol12 )
    ned2l.dshape12[index1] <- -trigamma(2) / shape1[index1]
    index2 <- (abs(shape2 - 2) < .tol12 )
    ned2l.dshape11[index2] <-
      (1 - 2 * psigamma(2, deriv = 2)) / shape1[index2]^2

    wz <- array(c(c(w) * ned2l.dshape11 * dshape1.deta^2,
                  c(w) * ned2l.dshape22 * dshape2.deta^2,
                  c(w) * ned2l.dshape12 * dshape1.deta *
                                          dshape2.deta),
                dim = c(n, M / M1, 3))
    wz <- arwz2wz(wz, M = M, M1 = M1)
    wz
  }), list( .lshape1 = lshape1, .lshape2 = lshape2,
             .eshape1 = eshape1, .eshape2 = eshape2,
             .tol12 = tol12 ))))
}  # kumar






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




  N <- max(length(x), length(vee), length(sigma))
  if (length(x)      != N) x      <- rep_len(x,      N)
  if (length(vee)    != N) vee    <- rep_len(vee   , N)
  if (length(sigma ) != N) sigma  <- rep_len(sigma , N)

  logdensity <- rep_len(log(0), N)
  xok <- (x > 0)
  x.abs <- abs(x[xok] * vee[xok] / sigma[xok]^2)
  logdensity[xok] <- log(x[xok]) - 2 * log(sigma[xok]) +
                     (-(x[xok]^2+vee[xok]^2)/(2*sigma[xok]^2)) +
                log(besselI(x.abs, nu = 0, expon.scaled = TRUE)) +
                     x.abs
  logdensity[sigma <= 0] <- NaN
  logdensity[vee < 0] <- NaN

  logdensity[is.infinite(x)] <- -Inf  # 20141209 KaiH

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



rrice <- function(n, sigma, vee) {
  theta <- 1  # any number
  X <- rnorm(n, mean = vee * cos(theta), sd = sigma)
  Y <- rnorm(n, mean = vee * sin(theta), sd = sigma)
  sqrt(X^2 + Y^2)
}




marcumQ <- function(a, b, m = 1,
                    lower.tail = TRUE, log.p = FALSE, ... ) {
  pchisq(b^2, df = 2*m, ncp = a^2,
         lower.tail = lower.tail, log.p = log.p, ... )
}



price <- function(q, sigma, vee,
                  lower.tail = TRUE, log.p = FALSE, ...) {
  ans <- marcumQ(vee/sigma, q/sigma, m = 1,
                 lower.tail = lower.tail, log.p = log.p, ... )
  ans
}



qrice <- function(p, sigma, vee,
                  lower.tail = TRUE, log.p = FALSE, ... ) {
  sqrt(qchisq(p, df = 2, ncp = (vee/sigma)^2,
              lower.tail = lower.tail,
              log.p = log.p, ... )) * sigma
}









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




 riceff <-
  function(lsigma = "loglink", lvee = "loglink",
           isigma = NULL, ivee = NULL,
           nsimEIM = 100, zero = NULL, nowarning = FALSE) {



  lvee     <- as.list(substitute(lvee))
  evee     <- link2list(lvee)
  lvee     <- attr(evee, "function.name")


  lsigma <- as.list(substitute(lsigma))
  esigma <- link2list(lsigma)
  lsigma <- attr(esigma, "function.name")



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

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


  new("vglmff",
  blurb = c("Rice distribution\n\n",
            "Links:    ",
    namesof("sigma", lsigma, earg = esigma, tag = FALSE), ", ",
    namesof("vee",   lvee,   earg = evee,   tag = FALSE), "\n",
            "Mean:     ",
            "sigma*sqrt(pi/2)*exp(z/2)*((1-z)*",
            "besselI(-z/2, nu = 0) - z * besselI(-z/2, nu = 1)) ",
            "where z=-vee^2/(2*sigma^2)"),

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

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

  infos = eval(substitute(function(...) {
    list(M1 = 2,
         Q1 = 1,
         dpqrfun = "rice",
         expected = FALSE,
         multipleResponses = FALSE,
         parameters.names = c("sigma", "vee"),
         nsimEIM = .nsimEIM,
         lsigma = .lsigma ,
         lvee = .lvee ,
         zero = .zero )
  }, list( .zero = zero, .lsigma = lsigma, .lvee = lvee,
           .nsimEIM = nsimEIM ))),
  initialize = eval(substitute(expression({

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


    predictors.names <-
      c(namesof("sigma", .lsigma , earg = .esigma , tag = FALSE),
        namesof("vee",   .lvee   , earg = .evee   , tag = FALSE))




    if (!length(etastart)) {
      riceff.Loglikfun <- function(vee, y, x, w, extraargs) {
        sigma.init <- sd(rep(y, w))
        sum(c(w) * (log(y) - 2*log(sigma.init) +
                    log(besselI(y*vee/sigma.init^2, nu = 0)) -
                   (y^2 + vee^2) / (2*sigma.init^2)))
      }
    vee.grid <-
      seq(quantile(rep(y, w), probs = seq(0, 1, 0.2))["20%"],
          quantile(rep(y, w), probs = seq(0, 1, 0.2))["80%"],
          len = 11)
    vee.init <- if (length( .ivee )) .ivee else
                grid.search(vee.grid, objfun = riceff.Loglikfun,
                            y = y,  x = x, w = w)
      vee.init <- rep_len(vee.init, length(y))
      sigma.init <- if (length( .isigma )) .isigma else
          sqrt(max((weighted.mean(y^2, w) - vee.init^2)/2, 0.001))
      sigma.init <- rep_len(sigma.init, length(y))

      etastart <-
        cbind(theta2eta(sigma.init, .lsigma , earg = .esigma ),
              theta2eta(vee.init,   .lvee ,   earg = .evee ))
    }
  }), list( .lvee = lvee, .lsigma = lsigma,
            .ivee = ivee, .isigma = isigma,
            .evee = evee, .esigma = esigma ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    vee   <- eta2theta(eta[, 1], link = .lvee ,   earg = .evee )
    sigma <- eta2theta(eta[, 2], link = .lsigma , earg = .esigma )
    temp9 <- -vee^2 / (2*sigma^2)


      sigma * sqrt(pi/2) *
      ((1-temp9) * besselI(-temp9/2, nu = 0, expon = TRUE) -
          temp9  * besselI(-temp9/2, nu = 1, expon = TRUE))
  }, list( .lvee = lvee, .lsigma = lsigma,
           .evee = evee, .esigma = esigma ))),
  last = eval(substitute(expression({
    misc$link <-    c("sigma" = .lsigma , "vee" = .lvee )

    misc$earg <- list("sigma" = .esigma , "vee" = .evee )

    misc$expected <- TRUE
    misc$nsimEIM <- .nsimEIM
    misc$multipleResponses <- FALSE
  }),
  list( .lvee = lvee, .lsigma = lsigma,
        .evee = evee, .esigma = esigma,
        .nsimEIM = nsimEIM ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    sigma <- eta2theta(eta[, 1], .lsigma , earg = .esigma )
    vee   <- eta2theta(eta[, 2], .lvee   , earg = .evee )
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * drice(x = y, sigma = sigma,
                              vee = vee, log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .lvee = lvee, .lsigma = lsigma,
           .evee = evee, .esigma = esigma ))),
  vfamily = c("riceff"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    sigma <- eta2theta(eta[, 1], link = .lsigma , earg = .esigma )
    vee   <- eta2theta(eta[, 2], link = .lvee   , earg = .evee )
    okay1 <- all(is.finite(sigma)) && all(0 < sigma) &&
             all(is.finite(vee  )) && all(0 < vee  )
    okay1
  }, list( .lvee = lvee, .lsigma = lsigma,
           .evee = evee, .esigma = esigma ))),


  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)
    sigma <- eta2theta(eta[, 1], link = .lsigma , .esigma )
    vee   <- eta2theta(eta[, 2], link = .lvee   , .evee )
    rrice(nsim * length(vee),
          vee = vee, sigma = sigma)
  }, list( .lvee = lvee, .lsigma = lsigma,
           .evee = evee, .esigma = esigma ))),



  deriv = eval(substitute(expression({
    sigma <- eta2theta(eta[, 1], link = .lsigma , .esigma )
    vee   <- eta2theta(eta[, 2], link = .lvee   , .evee )

    dvee.deta <- dtheta.deta(vee, link = .lvee , earg = .evee )
    dsigma.deta <- dtheta.deta(sigma, link = .lsigma , .esigma )

    temp8 <- y * vee / sigma^2
    dl.dvee <- -vee/sigma^2 + (y/sigma^2) *
               besselI(temp8, nu = 1) / besselI(temp8, nu = 0)
    dl.dsigma <- -2/sigma + (y^2 + vee^2)/(sigma^3) -
                 (2 * temp8 / sigma) *
                 besselI(temp8, nu = 1) / besselI(temp8, nu = 0)

    c(w) * cbind(dl.dsigma * dsigma.deta,
                 dl.dvee   * dvee.deta)

  }), list( .lvee = lvee, .lsigma = lsigma,
            .evee = evee, .esigma = esigma,
            .nsimEIM = nsimEIM ))),
  weight = eval(substitute(expression({
    run.var <- run.cov <- 0
    for (ii in 1:( .nsimEIM )) {
      ysim <- rrice(n, vee = vee, sigma = sigma)
      temp8 <- ysim * vee / sigma^2
      dl.dvee <- -vee/sigma^2 + (ysim/sigma^2) *
                 besselI(temp8, nu = 1) / besselI(temp8, nu = 0)
      dl.dsigma <- -2/sigma + (ysim^2 + vee^2)/(sigma^3) -
                   (2 * temp8 / sigma) *
                   besselI(temp8, nu = 1) / besselI(temp8, nu = 0)

      rm(ysim)
      temp3 <- cbind(dl.dsigma, dl.dvee)
      run.var <- ((ii-1) * run.var + temp3^2) / ii
      run.cov <- ((ii-1) * run.cov + temp3[, 1] * temp3[, 2]) / ii
    }
    wz <- if (intercept.only)
        matrix(colMeans(cbind(run.var, run.cov)),
               n, dimm(M), byrow = TRUE) else
               cbind(run.var, run.cov)

    dtheta.detas <- cbind(dsigma.deta, dvee.deta)
    index0 <- iam(NA_real_, NA_real_, M = M,
                  both = TRUE, diag = TRUE)
    wz <- wz * dtheta.detas[, index0$row] *
               dtheta.detas[, index0$col]
    c(w) * wz
  }), list( .lvee = lvee, .lsigma = lsigma,
            .evee = evee, .esigma = esigma,
            .nsimEIM = nsimEIM ))))
}  # riceff





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


  L <- max(length(x), length(mu1), length(mu2))
  if (length(x)      != L) x      <- rep_len(x,      L)
  if (length(mu1)    != L) mu1    <- rep_len(mu1,    L)
  if (length(mu2)    != L) mu2    <- rep_len(mu2,    L)

  ok2 <- is.finite(mu1) & is.finite(mu2) & (mu1 >= 0) & (mu2 >= 0)
  ok3 <- (mu1 == 0) & (mu2 >  0)
  ok4 <- (mu1 >  0) & (mu2 == 0)
  ok5 <- (mu1 == 0) & (mu2 == 0)
    if (log.arg) {
      ans <- -mu1 - mu2 + 2 * sqrt(mu1*mu2) +
             0.5 * x * log(mu1) - 0.5 * x * log(mu2) +
             log(besselI(2 * sqrt(mu1*mu2),

                         nu = abs(x),

                         expon.scaled = TRUE))
      ans[ok3] <- dpois(-x[ok3], lambda = mu2[ok3], log = TRUE)
      ans[ok4] <- dpois(-x[ok4], lambda = mu1[ok4], log = TRUE)
      ans[ok5] <- dpois( x[ok5], lambda = 0.0,      log = TRUE)
      ans[x != round(x)] = log(0.0)
    } else {
      ans <- (mu1/mu2)^(x/2) * exp(-mu1-mu2 + 2 * sqrt(mu1*mu2)) *
             besselI(2 * sqrt(mu1*mu2),

                     nu = abs(x),

                     expon.scaled = TRUE)
      ans[ok3] <- dpois(x = -x[ok3], lambda = mu2[ok3])
      ans[ok4] <- dpois(x = -x[ok4], lambda = mu1[ok4])
      ans[ok5] <- dpois(x =  x[ok5], lambda = 0.0)
      ans[x != round(x)] <- 0.0
    }
    ans[!ok2] <- NaN
    ans
}  # dskellam 






rskellam <- function(n, mu1, mu2) {
  rpois(n, mu1) - rpois(n, mu2)
}



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



 skellam <-
  function(lmu1 = "loglink", lmu2 = "loglink",
           imu1 = NULL,   imu2 = NULL,
           nsimEIM = 100, parallel = FALSE, zero = NULL) {

  lmu1 <- as.list(substitute(lmu1))
  emu1 <- link2list(lmu1)
  lmu1 <- attr(emu1, "function.name")

  lmu2 <- as.list(substitute(lmu2))
  emu2 <- link2list(lmu2)
  lmu2 <- attr(emu2, "function.name")


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



  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("Skellam distribution\n\n",
         "Links:    ",
         namesof("mu1", lmu1, earg = emu1, tag = FALSE), ", ",
         namesof("mu2", lmu2, earg = emu2, tag = FALSE), "\n",
         "Mean:     mu1-mu2", "\n",
         "Variance: mu1+mu2"),
  constraints = eval(substitute(expression({
    constraints <- cm.VGAM(matrix(1, M, 1), x = x,
                           bool = .parallel ,
                           constraints = constraints,
                           apply.int = TRUE)
    constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
                     M = M, M1 = 2,
                     predictors.names = predictors.names)
  }), list( .parallel = parallel, .zero = zero ))),
  infos = eval(substitute(function(...) {
    list(M1 = 2,
         Q1 = 1,
         dpqrfun = "skellam",
         expected = FALSE,
         multipleResponses = FALSE,
         parameters.names = c("mu1", "mu2"),
         nsimEIM = .nsimEIM,
         lmu1 = .lmu1 ,
         lmu2 = .lmu2 ,
         zero = .zero )
  }, list( .zero = zero, .lmu1 = lmu1, .lmu2 = lmu2,
           .nsimEIM = nsimEIM ))),
  initialize = eval(substitute(expression({


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


    predictors.names <- c(
      namesof("mu1", .lmu1, earg = .emu1, tag = FALSE),
      namesof("mu2", .lmu2, earg = .emu2, tag = FALSE))


    if (!length(etastart)) {
      junk <- lm.wfit(x = x, y = c(y), w = c(w))
      var.y.est <- sum(c(w) * junk$resid^2) / junk$df.residual
      mean.init <- weighted.mean(y, w)

      mu1.init <- max((var.y.est + mean.init) / 2, 0.01)
      mu2.init <- max((var.y.est - mean.init) / 2, 0.01)
      mu1.init <- rep_len(if (length( .imu1 )) .imu1 else mu1.init, n)
      mu2.init <- rep_len(if (length( .imu2 )) .imu2 else mu2.init, n)

      etastart <- cbind(theta2eta(mu1.init, .lmu1, earg = .emu1 ),
                        theta2eta(mu2.init, .lmu2, earg = .emu2 ))
      }
  }), list( .lmu1 = lmu1, .lmu2 = lmu2,
            .imu1 = imu1, .imu2 = imu2,
            .emu1 = emu1, .emu2 = emu2 ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
      mu1 <- eta2theta(eta[, 1], link = .lmu1, earg = .emu1 )
      mu2 <- eta2theta(eta[, 2], link = .lmu2, earg = .emu2 )
      mu1 - mu2
  }, list( .lmu1 = lmu1, .lmu2 = lmu2,
           .emu1 = emu1, .emu2 = emu2 ))),
  last = eval(substitute(expression({
      misc$link <-    c("mu1" = .lmu1, "mu2" = .lmu2)

      misc$earg <- list("mu1" = .emu1, "mu2" = .emu2 )

      misc$expected <- TRUE
      misc$nsimEIM <- .nsimEIM
  }), list( .lmu1 = lmu1, .lmu2 = lmu2,
            .emu1 = emu1, .emu2 = emu2,
            .nsimEIM = nsimEIM ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    mu1 <- eta2theta(eta[, 1], link = .lmu1, earg = .emu1 )
    mu2 <- eta2theta(eta[, 2], link = .lmu2, earg = .emu2 )
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {

      ll.elts <-
        if ( is.logical( .parallel ) &&
             length( .parallel ) == 1 &&
             .parallel )
          c(w) * log(besselI(2*mu1, nu = y, expon = TRUE)) else
          c(w) * (-mu1 - mu2 +
                  0.5 * y * log(mu1) -
                  0.5 * y * log(mu2) +
                  2 * sqrt(mu1*mu2) +  # Use this when expon = TRUE
                  log(besselI(2 * sqrt(mu1*mu2),
                              nu = y, expon = TRUE)))
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .lmu1 = lmu1, .lmu2 = lmu2,
           .emu1 = emu1, .emu2 = emu2,
           .parallel = parallel ))),
  vfamily = c("skellam"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    mu1 <- eta2theta(eta[, 1], link = .lmu1, earg = .emu1 )
    mu2 <- eta2theta(eta[, 2], link = .lmu2, earg = .emu2 )
    okay1 <- all(is.finite(mu1)) && all(0 < mu1) &&
             all(is.finite(mu2)) && all(0 < mu2)
    okay1
  }, list( .lmu1 = lmu1, .lmu2 = lmu2,
           .emu1 = emu1, .emu2 = emu2 ))),




  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)
    mu1 <- eta2theta(eta[, 1], link = .lmu1, earg = .emu1 )
    mu2 <- eta2theta(eta[, 2], link = .lmu2, earg = .emu2 )
    rskellam(nsim * length(mu1), mu1, mu2)
  },
  list( .lmu1 = lmu1, .lmu2 = lmu2,
        .emu1 = emu1, .emu2 = emu2,
        .parallel = parallel ))),





  deriv = eval(substitute(expression({
    mu1 <- eta2theta(eta[, 1], link = .lmu1, earg = .emu1 )
    mu2 <- eta2theta(eta[, 2], link = .lmu2, earg = .emu2 )

    dmu1.deta <- dtheta.deta(mu1, link = .lmu1, earg = .emu1 )
    dmu2.deta <- dtheta.deta(mu2, link = .lmu2, earg = .emu2 )

    temp8 <- 2 * sqrt(mu1*mu2)
    temp9 <-  besselI(temp8, nu = y  , expon = TRUE)
    temp7 <- (besselI(temp8, nu = y-1, expon = TRUE) +
              besselI(temp8, nu = y+1, expon = TRUE)) / 2
    temp6 <- temp7 / temp9

    dl.dmu1 <- -1 + 0.5 * y / mu1 + sqrt(mu2/mu1) * temp6
    dl.dmu2 <- -1 - 0.5 * y / mu2 + sqrt(mu1/mu2) * temp6

    c(w) * cbind(dl.dmu1 * dmu1.deta,
                 dl.dmu2 * dmu2.deta)
  }), list( .lmu1 = lmu1, .lmu2 = lmu2,
            .emu1 = emu1, .emu2 = emu2,
            .nsimEIM = nsimEIM ))),
  weight = eval(substitute(expression({
    run.var <- run.cov <- 0
    for (ii in 1:( .nsimEIM )) {
      ysim <- rskellam(n, mu1=mu1, mu2=mu2)
      temp9 <-  besselI(temp8, nu = ysim,   expon = TRUE)
      temp7 <- (besselI(temp8, nu = ysim-1, expon = TRUE) +
                besselI(temp8, nu = ysim+1, expon = TRUE)) / 2
      temp6 <- temp7 / temp9
      dl.dmu1 <- -1 + 0.5 * ysim/mu1 + sqrt(mu2/mu1) * temp6
      dl.dmu2 <- -1 - 0.5 * ysim/mu2 + sqrt(mu1/mu2) * temp6
      rm(ysim)
      temp3 <- cbind(dl.dmu1, dl.dmu2)
      run.var <- ((ii-1) * run.var + temp3^2) / ii
      run.cov <- ((ii-1) * run.cov + temp3[, 1] * temp3[, 2]) / ii
    }
    wz <- if (intercept.only)
          matrix(colMeans(cbind(run.var, run.cov)),
                 n, dimm(M), byrow = TRUE) else
          cbind(run.var, run.cov)

    dtheta.detas <- cbind(dmu1.deta, dmu2.deta)
    index0 <- iam(NA_real_, NA_real_, M = M,
                  both = TRUE, diag = TRUE)
    wz <- wz * dtheta.detas[, index0$row] *
               dtheta.detas[, index0$col]
    c(w) * wz
  }),
  list( .lmu1 = lmu1, .lmu2 = lmu2,
        .emu1 = emu1, .emu2 = emu2,
        .nsimEIM = nsimEIM ))))
}  # skellam







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

  LLL <- max(length(x), length(shape))
  if (length(x)     != LLL) x     <- rep_len(x,     LLL)
  if (length(shape) != LLL) shape <- rep_len(shape, LLL)
  
  bad0 <- !is.finite(shape) | shape <= 0
  bad <- bad0 | !is.finite(x) | x < 1 | x != round(x)

  logpdf <- x + shape
  if (any(!bad)) {
    logpdf[!bad] <- log(shape[!bad]) +
                        lbeta(x[!bad], shape[!bad] + 1)
  }

  logpdf[!bad0 & is.infinite(x)] <- log(0)
  logpdf[!bad0 & x < 1         ] <- log(0)
  logpdf[!bad0 & x != round(x) ] <- log(0)
  logpdf[ bad0] <- NaN

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




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

  tq <- trunc(q)

  if (lower.tail) {
    ans <- 1 - tq * beta(abs(tq), shape+1)
    ans[q < 1] <- 0
    ans[is.infinite(q) & 0 < q] <- 1  # 20141215 KaiH
  } else {
    ans <-     tq * beta(abs(tq), shape+1)
    ans[q < 1] <- 1
    ans[is.infinite(q) & 0 < q] <- 0  # 20160713
  }

  ans[shape <= 0] <- NaN
  if (log.p) log(ans) else ans
  ans
}



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

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

  lo <- rep_len(1, LLL) - 0.5
  approx.ans <- lo  # True at lhs
  hi <- 2 * lo + 10.5
  dont.iterate <- bad
  done <- dont.iterate | p <= pyules(hi, shape)
  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] <= pyules(hi[!done],
                                       shape = shape[!done]))
    iter <- iter + 1
  }


  foo <- function(q, shape, p)
    pyules(q, shape) - p

  lhs <- dont.iterate | (p <= dyules(1, shape))

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

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

  ans[!bad0 & !is.na(p) & p == 0] <- 1
  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
}  # qyules



ryules <- function(n, shape) {

  rgeom(n, prob = exp(-rexp(n, rate = shape))) + 1
}





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



 yulesimon <- function(lshape = "loglink",
                       ishape = NULL, nsimEIM = 200,
                       zero = NULL) {

  if (length(ishape) &&
      !is.Numeric(ishape, positive = TRUE))
    stop("argument 'ishape' must be > 0")



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


  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("Yule-Simon distribution ",
            "f(y) = shape * beta(y, shape + 1), ",
            "shape > 0, y = 1, 2,..\n\n",
            "Link:    ",
            namesof("shape", lshape, earg = eshape), "\n\n",
            "Mean:     shape / (shape - 1), provided shape > 1\n",
            "Variance: shape^2 / ((shape - 1)^2 * (shape - 2)), ",
            "provided shape > 2"),
  constraints = eval(substitute(expression({
    constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
                     M = M, M1 = 1,
                     predictors.names = predictors.names)

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

  infos = eval(substitute(function(...) {
    list(M1 = 1,
         Q1 = 1,
         dpqrfun = "yules",
         expected = TRUE,
         multipleResponses = TRUE,
         nsimEIM = .nsimEIM ,
         parameters.names = c("shape"),
         zero = .zero )
  }, list( .zero = zero,
           .nsimEIM = nsimEIM ))),

  rqresslot = eval(substitute(
    function(mu, y, w, eta, extra = NULL) {
    Shape <- eta2theta(eta, .lshape , earg = .eshape )
    scrambleseed <- runif(1)  # To scramble the seed
    qnorm(runif(length(y),
                pyules(y - 1, Shape),
                pyules(y    , Shape)))
  }, list( .lshape = lshape, .eshape = eshape ))),

  initialize = eval(substitute(expression({


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



    ncoly <- ncol(y)

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


    mynames1  <- param.names("shape", ncoly, skip1 = TRUE)
    predictors.names <-
      namesof(mynames1, .lshape , earg = .eshape , tag = FALSE)

    if (!length(etastart)) {
      wmeany <- colSums(y * w) / colSums(w) + 1/8

      shape.init <- wmeany / (wmeany - 1)
      shape.init <- matrix(if (length( .ishape )) .ishape else
                           shape.init, n, M, byrow = TRUE)
      etastart <- theta2eta(shape.init, .lshape , .eshape )
    }
  }), list( .lshape = lshape, .eshape = eshape,
            .ishape = ishape ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    ans <- shape <- eta2theta(eta, .lshape , earg = .eshape )
    ans[shape >  1] <- shape / (shape - 1)
    ans[shape <= 1] <- NA
    ans
  }, list( .lshape = lshape, .eshape = eshape ))),
  last = eval(substitute(expression({
    M1 <- extra$M1
    misc$link <- c(rep_len( .lshape , ncoly))
    names(misc$link) <- mynames1

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

    misc$M1 <- M1
    misc$ishape <- .ishape
    misc$nsimEIM <- .nsimEIM
  }), list( .lshape = lshape, .eshape = eshape,
            .nsimEIM = nsimEIM,
            .ishape = ishape ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * dyules(x = y, shape = shape, log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .lshape = lshape, .eshape = eshape ))),
  vfamily = c("yulesimon"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    okay1 <- all(is.finite(shape)) && all(0 < shape)
    okay1
  }, list( .lshape = lshape, .eshape = eshape ))),





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

    pwts <- if (length(pwts <- object@prior.weights) > 0)
              pwts else weights(object, type = "prior")
    if (any(pwts != 1))
      warning("ignoring prior weights")
    eta <- predict(object)
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    ryules(nsim * length(shape), shape = shape)
  }, list( .lshape = lshape, .eshape = eshape ))),







  deriv = eval(substitute(expression({
    M1 <- 1
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    dl.dshape <- 1/shape + digamma(1+shape) - digamma(1+shape+y)
    dshape.deta <- dtheta.deta(shape, .lshape , earg = .eshape )
    c(w) * dl.dshape * dshape.deta
  }), list( .lshape = lshape, .eshape = eshape ))),
  weight = eval(substitute(expression({






    run.var <- 0
    for (ii in 1:( .nsimEIM )) {
      ysim <- ryules(n, shape <- shape)
      dl.dshape <- 1/shape + digamma(1+shape) - digamma(1+shape+ysim)
      rm(ysim)
      temp3 <- dl.dshape
      run.var <- ((ii-1) * run.var + temp3^2) / ii
    }
    wz <- if (intercept.only)
        matrix(colMeans(cbind(run.var)),
               n, M, byrow = TRUE) else cbind(run.var)

    wz <- wz * dshape.deta^2


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







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


  if ( log.arg ) {
    ans <- 2 * log(theta) + log1p(x) - theta * x - log1p(theta)
    ans[x < 0 | is.infinite(x)] <- log(0)  # 20141209 KaiH
  } else {
    ans <- theta^2 * (1 + x) * exp(-theta * x) / (1 + theta)
    ans[x < 0 | is.infinite(x)] <- 0  # 20141209 KaiH
  }
  ans[theta <= 0] <- NaN
  ans
}



plind <- function(q, theta, 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(-theta * q + log1p(q / (1 + 1/theta))))
      ans[q <= 0 ] <- -Inf
      ans[q == Inf] <- 0
    } else {
      ans <- -expm1(-theta * q + log1p(q / (1 + 1/theta)))
      ans[q <= 0] <- 0
      ans[q == Inf] <- 1
    }
  } else {
    if (log.p) {
      ans <- -theta * q + log1p(q / (1 + 1/theta))
      ans[q <= 0] <- 0
      ans[q == Inf] <- -Inf
    } else {
      ans <- exp(-theta * q + log1p(q / (1 + 1/theta)))
      ans[q <= 0] <- 1
      ans[q == Inf] <- 0
    }
  }
  ans[theta <= 0] <- NaN
  ans
}






rlind <- function(n, theta) {
  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



  ifelse(runif(use.n) < rep_len(1 / (1 + 1/theta), use.n),
         rexp(use.n, theta),
         rgamma(use.n, shape = 2, scale = 1 / theta))
}



 lindley <- function(link = "loglink",
                     itheta = NULL, zero = NULL) {


  if (length(itheta) &&
      !is.Numeric(itheta, positive = TRUE))
    stop("argument 'itheta' must be > 0")


  link <- as.list(substitute(link))
  earg <- link2list(link)
  link <- attr(earg, "function.name")




  new("vglmff",
  blurb = c("Lindley distribution f(y) = ",
        "theta^2 * (1 + y) * exp(-theta * y) / (1 + theta), ",
        "theta > 0, y > 0,\n\n",
        "Link:    ",
        namesof("theta", link, earg = earg), "\n\n",
        "Mean:     (theta + 2) / (theta * (theta + 1))\n",
  "Variance: (theta^2 + 4 * theta + 2) / (theta * (theta + 1))^2"),

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

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

  infos = eval(substitute(function(...) {
    list(M1 = 1,
         Q1 = 1,
         dpqrfun = "lind",
         expected = TRUE,
         hadof = TRUE,
         multipleResponses = TRUE,
         parameters.names = c("theta"),
         zero = .zero )
  }, list( .zero = zero ))),

  initialize = eval(substitute(expression({


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



    ncoly <- ncol(y)

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


    mynames1  <- param.names("theta", ncoly, skip1 = TRUE)
    predictors.names <- namesof(mynames1, .link , earg = .earg ,
                                tag = FALSE)

    if (!length(etastart)) {
      wmeany <- colSums(y * w) / colSums(w) + 1/8


      theta.init <- 1 / (wmeany + 1)
      theta.init <- matrix(if (length( .itheta )) .itheta else
                           theta.init, n, M, byrow = TRUE)
      etastart <- theta2eta(theta.init, .link , earg = .earg )
    }
  }), list( .link = link, .earg = earg, .itheta = itheta ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    theta <- eta2theta(eta, .link , earg = .earg )
    (theta + 2) / (theta * (theta + 1))
  }, list( .link = link, .earg = earg ))),
  last = eval(substitute(expression({
    M1 <- extra$M1
    misc$link <- c(rep_len( .link , ncoly))
    names(misc$link) <- mynames1

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

    misc$M1 <- M1
    misc$itheta <- .itheta
  }), list( .link = link, .earg = earg,
            .itheta = itheta ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    theta <- eta2theta(eta, .link , earg = .earg )
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * dlind(x = y, theta = theta, log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .link = link, .earg = earg ))),
  vfamily = c("lindley"),


  hadof = eval(substitute(
  function(eta, extra = list(), deriv = 1,
           linpred.index = 1, w = 1,
           dim.wz = c(NROW(eta), NCOL(eta) * (NCOL(eta)+1)/2),
           ...) {
    theta <- eta2theta(eta, .link , earg = .earg )
    numer <- theta^2 + 4 * theta + 2
    denom <- (theta * (1 + theta))^2
    ans <- c(w) *
    switch(as.character(deriv),
           "0" =  numer / denom,
           "1" =  (2 * theta + 4 - numer *
  2 * theta * (1 + theta) * (1 + 2 * theta) / denom) / denom,
           "2" = NA * theta,
           "3" = NA * theta,
           stop("argument 'deriv' must be 0, 1, 2 or 3"))
    if (deriv == 0) ans else
      retain.col(ans, linpred.index)  # Since M1 = 1
  }, list( .link = link, .earg = earg ))),



  validparams = eval(substitute(function(eta, y, extra = NULL) {
    theta <- eta2theta(eta, .link , earg = .earg )
    okay1 <- all(is.finite(theta)) && all(0 < theta)
    okay1
  }, list( .link = link, .earg = earg ))),




  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)
    theta <- eta2theta(eta, .link , earg = .earg )
    rlind(nsim * length(theta), theta = theta)
  }, list( .link = link, .earg = earg ))),



  deriv = eval(substitute(expression({
    M1 <- 1
    theta <- eta2theta(eta, .link , earg = .earg )

    dl.dtheta <- 2 / theta - 1 / (1 + theta) - y

    DTHETA.DETA <- dtheta.deta(theta, .link , earg = .earg )

    c(w) * dl.dtheta * DTHETA.DETA
  }), list( .link = link, .earg = earg ))),
  weight = eval(substitute(expression({
    ned2l.dtheta2 <- (theta^2 + 4 * theta + 2) / (theta *
                     (1 + theta))^2

    c(w) * ned2l.dtheta2 * DTHETA.DETA^2
  }), list( .zero = zero ))))
}  # lindley







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

  if ( log.arg ) {
    ans <- 2 * log(theta) + log(theta + 2 + x) -
           (x+3) * log1p(theta)
    ans[(x != round(x)) | (x < 0)] <- log(0)
  } else {
    ans <- theta^2 * (theta + 2 + x) / (theta + 1)^(x+3)
    ans[(x != round(x)) | (x < 0)] <- 0
  }
  ans[ # !is.finite(theta) |
     (theta <= 0)] <- NA
  ans
}







dslash <- function(x, mu = 0, sigma = 1, log = FALSE,
                   smallno = .Machine$double.eps * 1000) {
  if (!is.logical(log.arg <- log) || length(log) != 1)
    stop("bad input for argument 'log'")
  rm(log)

  if (!is.Numeric(sigma) || any(sigma <= 0))
    stop("argument 'sigma' must be positive")
  L <- max(length(x), length(mu), length(sigma))
  if (length(x)     != L) x     <- rep_len(x,     L)
  if (length(mu)    != L) mu    <- rep_len(mu,    L)
  if (length(sigma) != L) sigma <- rep_len(sigma, L)

  zedd <- (x-mu)/sigma
  if (log.arg) {
    ifelse(abs(zedd) < smallno,
           -log(2*sigma*sqrt(2*pi)),
           log1p(-exp(-zedd^2/2)) - log(sqrt(2*pi)*sigma*zedd^2))
  } else {
    ifelse(abs(zedd) < smallno,
           1/(2*sigma*sqrt(2*pi)),
           -expm1(-zedd^2/2)/(sqrt(2*pi)*sigma*zedd^2))
  }
}




pslash <-
  function(q, mu = 0, sigma = 1, very.negative = -10000,
           lower.tail = TRUE, log.p = FALSE) {
  if (anyNA(q))
    stop("argument 'q' must have non-missing values")
  if (!is.Numeric(mu))
    stop("argument 'mu' must have finite and non-missing values")
  if (!is.Numeric(sigma, positive = TRUE))
    stop("argument 'sigma' must have positive finite ",
         "non-missing values")
  if (!is.Numeric(very.negative, length.arg = 1) ||
     (very.negative >= 0))
    stop("argument 'very.negative' must be quite negative")

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

  L <- max(length(q), length(mu), length(sigma))
  if (length(q)     != L) q     <- rep_len(q,     L)
  if (length(mu)    != L) mu    <- rep_len(mu,    L)
  if (length(sigma) != L) sigma <- rep_len(sigma, L)

  zedd <- (q - mu)/sigma
  ans <- as.numeric(q * NA)
  extreme.q <- FALSE
  for (ii in 1:L) {
    use.trick <- (-abs(zedd[ii]) <= very.negative)
    if (use.trick) {
      ans[ii] <- ifelse(zedd[ii] < 0, 0.0, 1.0)
      extreme.q <- TRUE
    } else
    if ((zedd[ii] >= very.negative) &&
         zedd[ii] <= 0.0) {
      temp2 <- integrate(dslash, lower = q[ii], upper = mu[ii],
                         mu = mu[ii], sigma = sigma[ii])
      if (temp2$message != "OK")
        warning("integrate() failed on 'temp2'")
      ans[ii] <- 0.5 - temp2$value
    } else {
      temp1 <- integrate(dslash, lower = mu[ii], upper =  q[ii],
                         mu = mu[ii], sigma = sigma[ii])
      if (temp1$message != "OK")
        warning("integrate() failed")
      ans[ii] <- 0.5 + temp1$value
    }
  }
  if (extreme.q)
    warning("returning 0 or 1 values for extreme ",
            "values of argument 'q'")

  if (lower.tail) {
    if (log.p) log(ans) else ans
  } else {
    if (log.p) log1p(-ans) else -expm1(log(ans))
  }
}




rslash <- function (n, mu = 0, sigma = 1) {
  rnorm(n = n, mean = mu, sd = sigma) / runif(n = n)
}



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



 slash <- function(lmu = "identitylink", lsigma = "loglink",
                   imu = NULL, isigma = NULL,
                   gprobs.y = ppoints(8),
                   nsimEIM = 250, zero = NULL,
                   smallno = .Machine$double.eps * 1000) {

  lmu <- as.list(substitute(lmu))
  emu <- link2list(lmu)
  lmu <- attr(emu, "function.name")

  lsigma <- as.list(substitute(lsigma))
  esigma <- link2list(lsigma)
  lsigma <- attr(esigma, "function.name")


  if (length(isigma) &&
      !is.Numeric(isigma, positive = TRUE))
    stop("argument 'isigma' must be > 0")



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

  if (!is.Numeric(gprobs.y, positive = TRUE) ||
      max(gprobs.y) >= 1)
    stop("bad input for argument 'gprobs.y'")
  if (!is.Numeric(smallno, positive = TRUE) ||
      smallno > 0.1)
    stop("bad input for argument 'smallno'")


  new("vglmff",
  blurb = c("Slash distribution\n\n",
         "Links:    ",
     namesof("mu",    lmu,    earg = emu,    tag = FALSE), ", ",
     namesof("sigma", lsigma, earg = esigma, tag = FALSE), "\n",
         paste(
         "1-exp(-(((y-mu)/sigma)^2)/2))/(sqrt(2*pi)*",
         "sigma*((y-mu)/sigma)^2)",
         "\ty!=mu",
         "\n1/(2*sigma*sqrt(2*pi))",
         "\t\t\t\t\t\t\ty=mu\n")),

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


  infos = eval(substitute(function(...) {
    list(M1 = 2,
         Q1 = 1,
         dpqrfun = "slash",
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = c("mu", "sigma"),
         lmu    = .lmu ,
         lsigma = .lsigma ,
         zero = .zero )
  }, list( .zero = zero, .lmu = lmu, .lsigma = lsigma ))),


  rqresslot = eval(substitute(
    function(mu, y, w, eta, extra = NULL) {
    mu    <- eta2theta(eta[, 1], link = .lmu    , .emu    )
    sigma <- eta2theta(eta[, 2], link = .lsigma , .esigma )
    scrambleseed <- runif(1)  # To scramble the seed
    qnorm(pslash(y, mu = mu, sigma))
  }, list( .lmu = lmu, .lsigma = lsigma,
           .emu = emu, .esigma = esigma ))),


  initialize = eval(substitute(expression({

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


    predictors.names <- c(
        namesof("mu",    .lmu ,    earg = .emu,    tag = FALSE),
        namesof("sigma", .lsigma , earg = .esigma, tag = FALSE))


    if (!length(etastart)) {

      slash.Loglikfun <- function(mu, y, x, w, extraargs) {
          sigma <- if (is.Numeric(.isigma)) .isigma else
            max(0.01,
         ((quantile(rep(y, w), prob = 0.75)/2)-mu)/qnorm(0.75))
          zedd <- (y-mu)/sigma
          sum(c(w) * ifelse(abs(zedd)<.smallno,
                         -log(2*sigma*sqrt(2*pi)),
                         log1p(-exp(-zedd^2/2)) -
                         log(sqrt(2*pi) * sigma * zedd^2)))
      }
      gprobs.y <- .gprobs.y
      mu.grid <- quantile(rep(y, w), probs = gprobs.y)
      mu.grid <- seq(mu.grid[1], mu.grid[2], length=100)
      mu.init <- if (length( .imu )) .imu else
                 grid.search(mu.grid, objfun = slash.Loglikfun,
                             y = y,  x = x, w = w)
      sigma.init <- if (is.Numeric(.isigma)) .isigma else
        max(0.01,
           ((quantile(rep(y, w), prob = 0.75)/2) -
                      mu.init) / qnorm(0.75))
      mu.init <- rep_len(mu.init, length(y))
      etastart <- matrix(0, n, 2)
      etastart[, 1] <- theta2eta(mu.init, .lmu , earg = .emu )
      etastart[, 2] <- theta2eta(sigma.init, .lsigma , .esigma )
    }
  }), list( .lmu = lmu, .lsigma = lsigma,
            .imu = imu, .isigma = isigma,
            .emu = emu, .esigma = esigma,
            .gprobs.y = gprobs.y, .smallno = smallno))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
      NA * eta2theta(eta[, 1], link = .lmu , earg = .emu )
  }, list( .lmu = lmu, .emu = emu ))),
  last = eval(substitute(expression({
    misc$link <-    c("mu" = .lmu , "sigma" = .lsigma )

    misc$earg <- list("mu" = .emu , "sigma" = .esigma )

    misc$expected <- TRUE
    misc$nsimEIM <- .nsimEIM
  }), list( .lmu = lmu, .lsigma = lsigma,
            .emu = emu, .esigma = esigma, .nsimEIM = nsimEIM ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    mu    <- eta2theta(eta[, 1], link = .lmu    , earg = .emu )
    sigma <- eta2theta(eta[, 2], link = .lsigma , earg = .esigma )
    zedd <- (y - mu) / sigma
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * dslash(y, mu = mu, sigma, log = TRUE,
                               smallno = .smallno)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .lmu = lmu, .lsigma = lsigma,
           .emu = emu, .esigma = esigma, .smallno = smallno ))),
  vfamily = c("slash"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    mu    <- eta2theta(eta[, 1], link = .lmu    , .emu    )
    sigma <- eta2theta(eta[, 2], link = .lsigma , .esigma )

    okay1 <- all(is.finite(mu))    &&
             all(is.finite(sigma)) && all(0 < sigma)
    okay1
  }, list( .lmu = lmu, .lsigma = lsigma,
           .emu = emu, .esigma = esigma ))),






  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)
    mu    <- eta2theta(eta[, 1], link = .lmu    , earg = .emu )
    sigma <- eta2theta(eta[, 2], link = .lsigma , earg = .esigma )
    rslash(nsim * length(sigma), mu = mu, sigma = sigma)
  }, list( .lmu = lmu, .lsigma = lsigma,
           .emu = emu, .esigma = esigma, .smallno = smallno ))),




  deriv = eval(substitute(expression({
    mu    <- eta2theta(eta[, 1], link = .lmu    , .emu    )
    sigma <- eta2theta(eta[, 2], link = .lsigma , .esigma )

    dmu.deta    <- dtheta.deta(mu,    link = .lmu    , .emu    )
    dsigma.deta <- dtheta.deta(sigma, link = .lsigma , .esigma )

    zedd <- (y - mu) / sigma
    d3 <- deriv3(~ w * log(1 - exp(-(((y - mu) / sigma)^2) / 2)) -
                 log(sqrt(2 * pi) * sigma * ((y - mu) / sigma)^2),
                 c("mu", "sigma"))
    eval.d3 <- eval(d3)
    dl.dthetas <-  attr(eval.d3, "gradient")
    dl.dmu    <- dl.dthetas[, 1]
    dl.dsigma <- dl.dthetas[, 2]
    ind0 <- (abs(zedd) < .smallno)
    dl.dmu[ind0] <- 0
    dl.dsigma[ind0] <- -1 / sigma[ind0]
    c(w) * cbind(dl.dmu * dmu.deta, dl.dsigma * dsigma.deta)
  }), list( .lmu = lmu, .lsigma = lsigma,
            .emu = emu, .esigma = esigma, .smallno = smallno ))),
  weight = eval(substitute(expression({
    run.varcov <- 0
    ind1 <- iam(NA_real_, NA_real_, M = M,
                both = TRUE, diag = TRUE)
    sd3 <- deriv3(~ w * log(1 -
                            exp(-(((ysim - mu) / sigma)^2) / 2))-
             log(sqrt(2 * pi) * sigma * ((ysim - mu) / sigma)^2),
             c("mu", "sigma"))
    for (ii in 1:( .nsimEIM )) {
      ysim <- rslash(n, mu = mu, sigma = sigma)
      seval.d3 <- eval(sd3)

      dl.dthetas <-  attr(seval.d3, "gradient")
      dl.dmu    <- dl.dthetas[, 1]
      dl.dsigma <- dl.dthetas[, 2]

      temp3 <- cbind(dl.dmu, dl.dsigma)
      run.varcov <- run.varcov + temp3[, ind1$row] *
                                 temp3[, ind1$col]
    }
    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
    dthetas.detas <- cbind(dmu.deta, dsigma.deta)
    wz <- wz * dthetas.detas[, ind1$row] *
               dthetas.detas[, ind1$col]
    c(w) * wz
  }), list( .lmu = lmu, .lsigma = lsigma,
            .emu = emu, .esigma = esigma,
            .nsimEIM = nsimEIM, .smallno = smallno ))))
}  # slash





dnefghs <- function(x, tau, log = FALSE) {



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


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

  logdensity <- log(sin(pi*tau)) + (1-tau)*x -
                log(pi) - log1pexp(x)
  logdensity[tau < 0] <- NaN
  logdensity[tau > 1] <- NaN
  if (log.arg) logdensity else exp(logdensity)
}



 nefghs <- function(link = "logitlink",
                    itau = NULL, imethod = 1) {

  if (length(itau) &&
      !is.Numeric(itau, positive = TRUE) ||
      any(itau >= 1))
    stop("argument 'itau' must be in (0, 1)")

  link <- as.list(substitute(link))
  earg <- link2list(link)
  link <- attr(earg, "function.name")

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


  new("vglmff",
  blurb = c("Natural exponential family generalized hyperbolic ",
            "secant distribution\n",
            "f(y) = sin(pi*tau)*exp((1-tau)*y)/(pi*(1+exp(y))\n\n",
            "Link:    ",
            namesof("tau", link, earg = earg), "\n\n",
            "Mean:     pi / tan(pi * tau)\n"),
  infos = eval(substitute(function(...) {
    list(M1 = 1,
         Q1 = 1,
         dpqrfun = "nefgs",
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = c("tau"),
         ltau = .link )
  }, list( .link = link ))),

  initialize = eval(substitute(expression({
    temp5 <-
    w.y.check(w = w, y = y,
              ncol.w.max = 1,
              ncol.y.max = 1,
              out.wy = TRUE,
              maximize = TRUE)
    w <- temp5$w
    y <- temp5$y


    predictors.names <-
      namesof("tau", .link , earg = .earg , tag = FALSE)


    if (!length(etastart)) {
      wmeany <- if ( .imethod == 1) weighted.mean(y, w) else
                median(rep(y, w))
      if (abs(wmeany) < 0.01)
        wmeany <- 0.01
      tau.init <- atan(pi / wmeany) / pi + 0.5
      tau.init[tau.init < 0.03] <- 0.03
      tau.init[tau.init > 0.97] <- 0.97
      tau.init <- rep_len(if (length( .itau )) .itau else
                          tau.init, n)
      etastart <- theta2eta(tau.init, .link , earg = .earg )
    }
  }), list( .link = link, .earg = earg,
            .itau = itau,
            .imethod = imethod ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    tau <- eta2theta(eta, .link , earg = .earg )
    pi / tan(pi * tau)
  }, list( .link = link, .earg = earg ))),
  last = eval(substitute(expression({
    misc$link <-    c(tau = .link )

    misc$earg <- list(tau = .earg )

    misc$expected <- TRUE
    misc$imethod <- .imethod
  }), list( .link = link, .earg = earg,
            .imethod = imethod ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    tau <- eta2theta(eta, .link , earg = .earg )
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * dnefghs(x = y, tau = tau, log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .link = link, .earg = earg ))),
  vfamily = c("nefghs"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    tau <- eta2theta(eta, .link , earg = .earg )
    okay1 <- all(is.finite(tau)) && all(0 < tau)
    okay1
  }, list( .link = link, .earg = earg ))),


  deriv = eval(substitute(expression({
    tau <- eta2theta(eta, .link , earg = .earg )
    dl.dtau <- pi / tan(pi * tau) - y
    dtau.deta <- dtheta.deta(tau, .link , earg = .earg )
    w * dl.dtau * dtau.deta
  }), list( .link = link, .earg = earg ))),
  weight = eval(substitute(expression({
    ned2l.dtau2 <- (pi / sin(pi * tau))^2
    wz <- ned2l.dtau2 * dtau.deta^2
    c(w) * wz
  }), list( .link = link ))))
}  # nefghs





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




  logdensity <- shape1*x - lbeta(shape1, shape2) -
                (shape1 + shape2) * log1pexp(x)

  logdensity[is.infinite(x)] <- -Inf  # 20141209 KaiH

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




 logF <- function(lshape1 = "loglink", lshape2 = "loglink",
                  ishape1 = NULL, ishape2 = 1,
                  imethod = 1) {

  if (length(ishape1) &&
      !is.Numeric(ishape1, positive = TRUE))
    stop("argument 'ishape1' must be positive")
  if ( # length(ishape2) &&
      !is.Numeric(ishape2, positive = TRUE))
    stop("argument 'ishape2' must be positive")


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


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


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

  new("vglmff",
  blurb = c("log F distribution\n",
            "f(y) = exp(-shape2 * y) / (beta(shape1, shape2) * ",
            "(1 + exp(-y))^(shape1 + shape2))\n\n",
            "Link:    ",
            namesof("shape1", lshape1, earg = eshape1), ", ",
            namesof("shape2", lshape2, earg = eshape2), "\n\n",
            "Mean:     digamma(shape1) - digamma(shape2)"),
  infos = eval(substitute(function(...) {
    list(M1 = 2,
         Q1 = 1,
         dpqrfun = "logF",
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = c("shape1", "shape2"),
         lshape1 = .lshape1 ,
         lshape2 = .lshape2 ,
         imethod = .imethod )
  }, list( .lshape1 = lshape1, .imethod = imethod,
           .lshape2 = lshape2 ))),

  initialize = eval(substitute(expression({

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


    predictors.names <- c(
      namesof("shape1", .lshape1 , earg = .eshape1 , tag = FALSE),
      namesof("shape2", .lshape2 , earg = .eshape2 , tag = FALSE))


    if (!length(etastart)) {
      wmeany <- if ( .imethod == 1) weighted.mean(y, w) else
                median(rep(y, w))


      shape1.init <- shape2.init <- rep_len( .ishape2 , n)
      shape1.init <- if (length( .ishape1))
                       rep_len( .ishape1, n) else {
                index1 <- (y > wmeany)
                shape1.init[ index1] <- shape2.init[ index1] + 1/1
                shape1.init[!index1] <- shape2.init[!index1] - 1/1
                shape1.init <- pmax(shape1.init, 1/8)
                shape1.init
              }
      etastart <-
          cbind(theta2eta(shape1.init, .lshape1 , earg = .eshape1 ),
                theta2eta(shape2.init, .lshape2 , earg = .eshape2 ))
    }
  }), list( .lshape1 = lshape1, .lshape2 = lshape2,
            .eshape1 = eshape1, .eshape2 = eshape2,
            .ishape1 = ishape1, .ishape2 = ishape2,
            .imethod = imethod ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    shape1 <- eta2theta(eta[, 1], .lshape1 , earg = .eshape1 )
    shape2 <- eta2theta(eta[, 2], .lshape2 , earg = .eshape2 )
    digamma(shape1) - digamma(shape2)
  }, list( .lshape1 = lshape1, .lshape2 = lshape2,
           .eshape1 = eshape1, .eshape2 = eshape2 ))),
  last = eval(substitute(expression({
    misc$link <-    c(shape1 = .lshape1 , shape2 = .lshape2 )
    misc$earg <- list(shape1 = .eshape1 , shape2 = .eshape2 )



    extra$percentile <- numeric(ncol(y))
    locat <- cbind(digamma(shape1) - digamma(shape2))  # zz unsure
    for (ii in 1:ncol(y)) {
      y.use <- if (ncol(y) > 1) y[, ii] else y
      extra$percentile[ii] <-
          100 * weighted.mean(y.use <= locat[, ii],
                              w[, min(ii, ncol(w))])
    }




    misc$expected <- TRUE
    misc$imethod <- .imethod
  }), list( .lshape1 = lshape1, .lshape2 = lshape2,
            .eshape1 = eshape1, .eshape2 = eshape2,
            .imethod = imethod ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL, summation = TRUE) {
    shape1 <- eta2theta(eta[, 1], .lshape1 , earg = .eshape1 )
    shape2 <- eta2theta(eta[, 2], .lshape2 , earg = .eshape2 )
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * dlogF(x = y, shape1 = shape1,
                              shape2 = shape2, log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .lshape1 = lshape1, .lshape2 = lshape2,
           .eshape1 = eshape1, .eshape2 = eshape2 ))),
  vfamily = c("logF"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    shape1 <- eta2theta(eta[, 1], .lshape1 , earg = .eshape1 )
    shape2 <- eta2theta(eta[, 2], .lshape2 , earg = .eshape2 )
    okay1 <- all(is.finite(shape1)) && all(0 < shape1) &&
             all(is.finite(shape2)) && all(0 < shape2)
    okay1
  }, list( .lshape1 = lshape1, .lshape2 = lshape2,
           .eshape1 = eshape1, .eshape2 = eshape2 ))),













  deriv = eval(substitute(expression({
    shape1 <- eta2theta(eta[, 1], .lshape1 , earg = .eshape1 )
    shape2 <- eta2theta(eta[, 2], .lshape2 , earg = .eshape2 )

    tmp888 <- digamma(shape1 + shape2) - log1pexp(-y)
    dl.dshape1 <- tmp888 - digamma(shape1)
    dl.dshape2 <- tmp888 - digamma(shape2) - y

    dshape1.deta <- dtheta.deta(shape1, .lshape1 , .eshape1 )
    dshape2.deta <- dtheta.deta(shape2, .lshape2 , .eshape2 )

    c(w) * cbind(dl.dshape1 * dshape1.deta,
                 dl.dshape2 * dshape2.deta)
  }), list( .lshape1 = lshape1, .lshape2 = lshape2,
            .eshape1 = eshape1, .eshape2 = eshape2 ))),
  weight = eval(substitute(expression({
    tmp888 <- trigamma(shape1 + shape2)
    ned2l.dshape12 <- trigamma(shape1) - tmp888
    ned2l.dshape22 <- trigamma(shape2) - tmp888
    ned2l.dshape1shape2 <- -tmp888

    wz <- matrix(0, n, dimm(M))
    wz[, iam(1, 1, M = M)] <- ned2l.dshape12 * dshape1.deta^2
    wz[, iam(2, 2, M = M)] <- ned2l.dshape22 * dshape2.deta^2
    wz[, iam(1, 2, M = M)] <- ned2l.dshape1shape2 *
                              dshape1.deta * dshape2.deta

    c(w) * wz
  }), list( .lshape1 = lshape1, .lshape2 = lshape2,
            .eshape1 = eshape1, .eshape2 = eshape2 ))))
}  # logF







dbenf <- function(x, ndigits = 1, log = FALSE) {
  if (!is.Numeric(ndigits, length.arg = 1,
                  positive = TRUE, integer.valued = TRUE) ||
      ndigits > 2)
    stop("argument 'ndigits' must be 1 or 2")
  lowerlimit <- ifelse(ndigits == 1, 1, 10)
  upperlimit <- ifelse(ndigits == 1, 9, 99)

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


  ans <- x * NA
  indexTF <- is.finite(x) & (x >= lowerlimit)

  ans[indexTF] <- log10(1 + 1/x[indexTF])
  ans[!is.na(x) & !is.nan(x) &
     ((x < lowerlimit) |
      (x > upperlimit) |
      (x != round(x)))] <- 0.0
  if (log.arg) log(ans) else ans
}



rbenf <- function(n, ndigits = 1) {
  if (!is.Numeric(ndigits, length.arg = 1,
                  positive = TRUE, integer.valued = TRUE) ||
      ndigits > 2)
    stop("argument 'ndigits' must be 1 or 2")
  lowerlimit <- ifelse(ndigits == 1, 1, 10)
  upperlimit <- ifelse(ndigits == 1, 9, 99)
  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
  myrunif <- runif(use.n)

  ans <- rep_len(lowerlimit, use.n)
  for (ii in (lowerlimit+1):upperlimit) {
      indexTF <- (pbenf(ii-1, ndigits = ndigits) < myrunif) &
                 (myrunif <= pbenf(ii, ndigits = ndigits))
      ans[indexTF] <- ii
  }
  ans
}



 pbenf <-
  function(q, ndigits = 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.p'")

  if (!is.Numeric(ndigits, length.arg = 1,
                  positive = TRUE, integer.valued = TRUE) ||
      ndigits > 2)
    stop("argument 'ndigits' must be 1 or 2")
  lowerlimit <- ifelse(ndigits == 1, 1, 10)
  upperlimit <- ifelse(ndigits == 1, 9, 99)

  ans <- q * NA
  floorq <- floor(q)
  indexTF <- is.finite(q) & (floorq >= lowerlimit)

  if (ndigits == 1) {
    if (lower.tail) {
      if (log.p) {
        ans[indexTF] <- log(log10(1 + floorq[indexTF]))
        ans[q <  lowerlimit ] <- -Inf
        ans[q >= upperlimit] <- 0
      } else {
        ans[indexTF] <- log10(1 + floorq[indexTF])
        ans[q <  lowerlimit] <- 0
        ans[q >= upperlimit] <- 1
      }
    } else {
      if (log.p) {
        ans[indexTF] <- log1p(-log10(1 + floorq[indexTF]))
        ans[q <  lowerlimit] <- 0
        ans[q >= upperlimit] <- -Inf
      } else {
        ans[indexTF] <- log10(10 / (1 + floorq[indexTF]))
        ans[q <  lowerlimit] <- 1
        ans[q >= upperlimit] <- 0
      }
    }
  } else {
    if (lower.tail) {
      if (log.p) {
        ans[indexTF] <- log(log10((1 + floorq[indexTF])/10))
        ans[q <  lowerlimit ] <- -Inf
        ans[q >= upperlimit] <- 0
      } else {
        ans[indexTF] <- log10((1 + floorq[indexTF])/10)
        ans[q <  lowerlimit] <- 0
        ans[q >= upperlimit] <- 1
     }
    } else {
      if (log.p) {
        ans[indexTF] <- log(log10(100/(1 + floorq[indexTF])))
        ans[q <  lowerlimit] <- 0
        ans[q >= upperlimit] <- -Inf
      } else {
        ans[indexTF] <- log10(100/(1 + floorq[indexTF]))
        ans[q <  lowerlimit] <- 1
        ans[q >= upperlimit] <- 0
      }
    }
  }
  ans
}



if (FALSE)
qbenf <- function(p, ndigits = 1) {

  if (!is.Numeric(ndigits, length.arg = 1,
                  positive = TRUE, integer.valued = TRUE) ||
      ndigits > 2)
    stop("argument 'ndigits' must be 1 or 2")
  lowerlimit <- ifelse(ndigits == 1, 1, 10)
  upperlimit <- ifelse(ndigits == 1, 9, 99)
  bad <- !is.na(p) & !is.nan(p) & ((p < 0) | (p > 1))
  if (any(bad))
    stop("bad input for argument 'p'")

  ans <- rep_len(lowerlimit, length(p))
  for (ii in (lowerlimit+1):upperlimit) {
    indexTF <- is.finite(p) &
              (pbenf(ii-1, ndigits = ndigits) < p) &
              (p <= pbenf(ii, ndigits = ndigits))
    ans[indexTF] <- ii
  }

  ans[ is.na(p) |  is.nan(p)] <- NA
  ans[!is.na(p) & !is.nan(p) & (p == 0)] <- lowerlimit
  ans[!is.na(p) & !is.nan(p) & (p == 1)] <- upperlimit
  ans
}




qbenf <- function(p, ndigits = 1,
                  lower.tail = TRUE, log.p = FALSE) {
  if (!is.Numeric(ndigits, length.arg = 1,
                  positive = TRUE, integer.valued = TRUE) ||
      ndigits > 2)
    stop("argument 'ndigits' must be 1 or 2")

  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 (log.p) {
    bad <- ((p > 0) | is.na(p) | is.nan(p))
  } else {
    bad <- ((p < 0) | (p > 1) | is.na(p) | is.nan(p))
  }
  if (any(bad))
    stop("bad input for argument 'p'")

  lowerlimit <- ifelse(ndigits == 1, 1, 10)
  upperlimit <- ifelse(ndigits == 1, 9, 99)
  ans <- rep_len(lowerlimit, length(p))

  if (lower.tail) {
    for (ii in (lowerlimit+1):upperlimit) {
      indexTF <- is.finite(p) &
                 (pbenf(ii-1, ndigits = ndigits,
               lower.tail = lower.tail, log.p = log.p) < p) &
              (p <= pbenf(ii, ndigits = ndigits,
                    lower.tail = lower.tail, log.p = log.p))
      ans[indexTF] <- ii
    }
  } else {  ## when lower.tail = F, pbenf(ii-1) >= p & pben(ii) < p
    for (ii in (lowerlimit+1):upperlimit) {
      indexTF <- is.finite(p) &
                 (pbenf(ii-1, ndigits = ndigits,
                 lower.tail = lower.tail, log.p = log.p) >= p) &
                 (p > pbenf(ii, ndigits = ndigits,
                      lower.tail = lower.tail, log.p = log.p))
      ans[indexTF] <- ii
    }
  }

  if (lower.tail) {
    if (log.p) {
      ans[p > 0] <- NaN
      ans[p == -Inf] <- lowerlimit
    } else {
      ans[p < 0] <- NaN
      ans[p == 0] <- lowerlimit
      ans[p == 1] <- upperlimit
      ans[p > 1] <- NaN
    }
  } else {
    if (log.p) {
      ans[p > 0] <- NaN
      ans[p == -Inf] <- upperlimit
    } else {
      ans[p < 0] <- NaN
      ans[p == 0] <- upperlimit
      ans[p == 1] <- lowerlimit
      ans[p > 1] <- NaN
    }
  }
  ans
}


















 truncgeometric <-
  function(upper.limit = Inf,  # lower.limit = 1,  # Inclusive
           link = "logitlink", expected = TRUE,
           imethod = 1, iprob = NULL, zero = NULL) {

  if (is.finite(upper.limit) &&
      !is.Numeric(upper.limit, integer.valued = TRUE,
                  positive = TRUE))
    stop("bad input for argument 'upper.limit'")

  if (any(upper.limit < 0))
    stop("bad input for argument 'upper.limit'")



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


  link <- as.list(substitute(link))
  earg <- link2list(link)
  link <- attr(earg, "function.name")


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



  uu.ll <- min(upper.limit)


  new("vglmff",
  blurb = c("Truncated geometric distribution ",
            "(P[Y=y] =\n",
            "     ",
            "prob * (1 - prob)^y / [1-(1-prob)^",
             uu.ll+1, "], y = 0,1,...,",
             uu.ll, ")\n",
            "Link:     ",
            namesof("prob", link, earg = earg), "\n",
            "Mean:     mu = 1 / prob - 1 ",
            ifelse(is.finite(upper.limit),
                   paste("- (", upper.limit+1, ") * (1 - prob)^",
                         upper.limit+1, " / (1 - ",
                   "(1 - prob)^", upper.limit+1, ")", sep = ""),
                         "")),
  constraints = eval(substitute(expression({
    constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
                     M = M, M1 = 1,
                     predictors.names = predictors.names)
  }), list( .zero = zero ))),

  infos = eval(substitute(function(...) {
    list(M1 = 1,
         Q1 = 1,
         expected = .expected ,
         imethod = .imethod ,
         multipleResponses = TRUE,
         parameters.names = c("prob"),
         upper.limit = .upper.limit ,
         zero = .zero )
  }, list( .zero = zero,
           .expected = expected,
           .imethod = imethod,
           .upper.limit = upper.limit ))),

  initialize = eval(substitute(expression({

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


    ncoly <- ncol(y)
    M1 <- 1
    extra$ncoly <- ncoly
    extra$M1 <- M1
    M <- M1 * ncoly
    extra$upper.limit <- matrix( .upper.limit , n,
                                ncoly, byrow = TRUE)

    if (any(y > extra$upper.limit))
      stop("some response values greater than ",
           "argument 'upper.limit'")


    mynames1 <- param.names("prob", ncoly, skip1 = TRUE)
    predictors.names <-
      namesof(mynames1, .link , earg = .earg , tag = FALSE)


    if (!length(etastart)) {
      prob.init <- if ( .imethod == 2)
                      1 / (1 + y + 1/16) else
                  if ( .imethod == 3)
                      1 / (1 + apply(y, 2, median) + 1/16) else
                      1 / (1 + colSums(y * w) / colSums(w) + 1/16)

      if (!is.matrix(prob.init))
        prob.init <- matrix(prob.init, n, M, byrow = TRUE)


      if (length( .iprob ))
        prob.init <- matrix( .iprob , n, M, byrow = TRUE)


        etastart <- theta2eta(prob.init, .link , earg = .earg )
    }
  }), list( .link = link, .earg = earg,
            .upper.limit = upper.limit,
            .imethod = imethod, .iprob = iprob ))),

  linkinv = eval(substitute(function(eta, extra = NULL) {
    prob <- eta2theta(eta, .link , earg = .earg )
    QQQ <- 1 - prob
    upper.limit <- extra$upper.limit
    tmp1 <- QQQ^(upper.limit+1)
    answer <- 1 / prob - 1 - (upper.limit+1) * tmp1 / (1 - tmp1)
    answer[!is.finite(answer)] <- 1 / prob[!is.finite(answer)] - 1
    answer
  }, list( .link = link, .earg = earg ))),

  last = eval(substitute(expression({
    M1 <- extra$M1
    misc$link <- c(rep_len( .link , ncoly))
    names(misc$link) <- mynames1

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

    misc$M1 <- M1
    misc$multipleResponses <- TRUE
    misc$expected <- ( .expected )
    misc$imethod <- ( .imethod )
    misc$iprob <- ( .iprob )
  }), list( .link = link, .earg = earg,
            .iprob = iprob,
            .upper.limit = upper.limit,
            .expected = expected, .imethod = imethod ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    prob <- eta2theta(eta, .link , earg = .earg )
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      upper.limit <- extra$upper.limit
      ll.elts <- c(w) * (dgeom(x = y, prob = prob, log = TRUE) -
                         log1p(-(1.0 - prob)^(1 + upper.limit)))
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .link = link, .earg = earg ))),
  vfamily = c("truncgeometric"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    prob <- eta2theta(eta, .link , earg = .earg )
    okay1 <- all(is.finite(prob)) && all(0 < prob & prob < 1)
    okay1
  }, list( .link = link, .earg = earg ))),
  deriv = eval(substitute(expression({
    prob <- eta2theta(eta, .link , earg = .earg )
    sss <- upper.limit <- extra$upper.limit  # Is a matrix

    QQQ <- 1 - prob
    tmp1 <- QQQ^(upper.limit + 1)
    dl.dprob <- 1 / prob  + (0 - y) / (1 - prob) -
      (1 + upper.limit) * QQQ^(upper.limit - 0) / (1 - tmp1)
    dl.dprob[!is.finite(upper.limit)] <-
        1 / prob[!is.finite(upper.limit)] +
        (0 - y[!is.finite(upper.limit)]) / (1 -
             prob[!is.finite(upper.limit)])


    dprobdeta <- dtheta.deta(prob, .link , earg = .earg )
    c(w) * cbind(dl.dprob * dprobdeta)
  }), list( .link = link, .earg = earg,
            .upper.limit = upper.limit,
            .expected = expected ))),
  weight = eval(substitute(expression({

    eim.oim.fun <- function(mu.y, sss)
      ifelse(is.finite(sss),
             1/prob^2 + (0 + mu.y) / QQQ^2 - (1+sss) *
             ((sss-0) * QQQ^(sss-1) / (1 - tmp1) +
             (1+sss) * QQQ^(2*sss) / (1 - tmp1)^2),
             1 / (prob^2 * (1 - prob)))


    ned2l.dprob2 <- if ( .expected ) {
      eim.oim.fun(mu, sss)
    } else {
      eim.oim.fun(y, sss)
    }
    wz <- ned2l.dprob2 * dprobdeta^2
    if ( !( .expected ))
      wz <- wz - dl.dprob * d2theta.deta2(prob,
                                          .link , .earg )
    c(w) * wz
  }), list( .link = link, .earg = earg,
            .expected = expected ))))
}  # truncgeometric









 betaff <-
  function(A = 0, B = 1,
           lmu = "logitlink",
           lphi = "loglink",
           imu = NULL, iphi = NULL,  # imethod = 1,
           gprobs.y = ppoints(8),  # (1:9)/10,
           gphi  = exp(-3:5)/4,
           zero = NULL) {



  if (!is.Numeric(A, length.arg = 1) ||
      !is.Numeric(B, length.arg = 1) || A >= B)
    stop("A must be < B, and both must be of length one")

  stdbeta <- (A == 0 && B == 1)


  lmu <- as.list(substitute(lmu))
  emu <- link2list(lmu)
  lmu <- attr(emu, "function.name")



  lphi <- as.list(substitute(lphi))
  ephi <- link2list(lphi)
  lphi <- attr(ephi, "function.name")



  if (length(imu) && (!is.Numeric(imu, positive = TRUE) ||
     any(imu <= A) || any(imu >= B)))
    stop("bad input for argument 'imu'")
  if (length(iphi) && !is.Numeric(iphi, positive = TRUE))
    stop("bad input for argument 'iphi'")


  new("vglmff",
  blurb = c("Beta distribution parameterized by mu and a ",
            "precision parameter\n",
  if (stdbeta) paste("f(y) = y^(mu*phi-1) * (1-y)^((1-mu)*phi-1)",
            "/ beta(mu*phi,(1-mu)*phi),\n",
            "      0<y<1, 0<mu<1, phi>0\n\n") else
            paste("f(y) = (y-",A,")^(mu1*phi-1) * (",B,
  "-y)^(((1-mu1)*phi)-1) / \n(beta(mu1*phi,(1-mu1)*phi) * (",
            B, "-", A, ")^(phi-1)),\n",
            A," < y < ",B, ", ", A," < mu < ",B,
            ", mu = ", A, " + ", (B-A), " * mu1",
            ", phi > 0\n\n", sep = ""),
            "Links:    ",
            namesof("mu",  lmu,  earg = emu),  ", ",
            namesof("phi", lphi, earg = ephi)),
  constraints = eval(substitute(expression({
    constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
                     M = M, M1 = 2,
                     predictors.names = predictors.names)
  }), list( .zero = zero ))),
  infos = eval(substitute(function(...) {
    list(M1 = 2,
         Q1 = 1,
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = c("mu", "phi"),
         A = .A ,
         B = .B ,
         zero = .zero )
  }, list( .zero = zero,
           .A = A, .B = B ))),

  initialize = eval(substitute(expression({
    if (min(y) <= .A || max(y) >= .B)
      stop("data not within (A, B)")


    temp5 <-
    w.y.check(w = w, y = y,
              out.wy = TRUE,
              maximize = TRUE)
    w <- temp5$w
    y <- temp5$y
    extra$A <- .A  # Needed for @validparams
    extra$B <- .B


    predictors.names <-
        c(namesof("mu",  .lmu ,  .emu , short = TRUE),
          namesof("phi", .lphi , .ephi, short = TRUE))
    if (!length(etastart)) {
      NOS <- 1
      muu.init <-
      phi.init <- matrix(NA_real_, n, NOS)
      gprobs.y <- .gprobs.y
      gphi <- if (length( .iphi )) .iphi else .gphi

      betaff.Loglikfun <- function(muu, phi, y, x, w, extraargs) {
        zedd <- (y   - extraargs$A) / ( extraargs$B - extraargs$A)
        m1u  <- (muu - extraargs$A) / ( extraargs$B - extraargs$A)
        shape1 <- phi * m1u
        shape2 <- (1 - m1u) * phi
        sum(c(w) * (dbeta(x = zedd, shape1, shape2, log = TRUE) -
                    log(abs( extraargs$B - extraargs$A ))))
      }


      for (jay in 1:NOS) {  # For each response 'y_jay'... do:
        gmuu <- if (length( .imu )) .imu else
                quantile(y[, jay], probs = gprobs.y)


        try.this <-
          grid.search2(gmuu, gphi,
                       objfun = betaff.Loglikfun,
                       y = y[, jay],
                       w = w[, jay],
                       extraargs = list(A = .A , B = .B ),
                       ret.objfun = TRUE)  # Last value is \ell
        muu.init[, jay] <-  try.this["Value1"]
        phi.init[, jay] <-  try.this["Value2"]
      }  # for (jay ...)


if (FALSE) {
      mu.init <- if (is.Numeric( .imu )) .imu else {
                   if ( .imethod == 1) weighted.mean(y, w) else
                      (y + weighted.mean(y, w)) / 2
                 }
      mu1.init <- (mu.init - .A ) / ( .B - .A )  # In (0,1)
      phi.init <- if (is.Numeric( .iphi )) .iphi else
   max(0.01, -1 + ( .B - .A )^2 * mu1.init*(1-mu1.init)/var(y))
  }



      etastart <- matrix(0, n, 2)
      etastart[, 1] <- theta2eta(muu.init, .lmu  , earg = .emu  )
      etastart[, 2] <- theta2eta(phi.init, .lphi , earg = .ephi )
    }
  }), list( .lmu = lmu, .lphi = lphi, .imu = imu, .iphi = iphi,
            .A = A, .B = B, .emu = emu, .ephi = ephi,
            .gprobs.y = gprobs.y, .gphi = gphi  ))),

  linkinv = eval(substitute(function(eta, extra = NULL) {
     mu <- eta2theta(eta[, 1], .lmu , .emu )
     mu
  }, list( .lmu = lmu, .emu = emu, .A = A, .B = B))),
  last = eval(substitute(expression({
    misc$link <-    c(mu = .lmu , phi = .lphi )
    misc$earg <- list(mu = .emu , phi = .ephi )
    misc$limits <- c( .A , .B )
    misc$stdbeta <- .stdbeta
  }), list( .lmu = lmu, .lphi = lphi, .A = A, .B = B,
            .emu = emu, .ephi = ephi,
            .stdbeta = stdbeta ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    mu  <- eta2theta(eta[, 1], .lmu  , earg = .emu  )
    phi <- eta2theta(eta[, 2], .lphi , earg = .ephi )
    m1u <- if ( .stdbeta ) mu else (mu - .A ) / ( .B - .A )
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      shape1 <- phi * m1u
      shape2 <- (1 - m1u) * phi
      zedd <- (y - .A) / ( .B - .A)
      ll.elts <-
        c(w) * (dbeta(x = zedd, shape1 = shape1, shape2 = shape2,
                      log = TRUE) -
                log( abs( .B - .A )))
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .lmu = lmu, .lphi = lphi, .A = A, .B = B,
           .emu = emu, .ephi = ephi,
           .stdbeta = stdbeta ))),
  vfamily = "betaff",
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    mu  <- eta2theta(eta[, 1], .lmu  , .emu  )
    phi <- eta2theta(eta[, 2], .lphi , .ephi )
    okay1 <- all(is.finite(mu )) &&
             all(extra$A < mu & mu < extra$B) &&
             all(is.finite(phi)) && all(0 < phi)
    okay1
  }, list( .lmu = lmu, .lphi = lphi, .A = A, .B = B,
           .emu = emu, .ephi = ephi ))),


  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)
    mu  <- eta2theta(eta[, 1], .lmu  , earg = .emu  )
    phi <- eta2theta(eta[, 2], .lphi , earg = .ephi )
    m1u <- if ( .stdbeta ) mu else (mu - .A ) / ( .B - .A )
    shape1 <- phi * m1u
    shape2 <- (1 - m1u) * phi
    .A + ( .B - .A ) *
    rbeta(nsim * length(shape1), shape1 = shape1, shape2 = shape2)
  }, list( .lmu = lmu, .lphi = lphi, .A = A, .B = B,
           .emu = emu, .ephi = ephi,
           .stdbeta = stdbeta ))),





  deriv = eval(substitute(expression({
    mu  <- eta2theta(eta[, 1], .lmu  , .emu  )  # 20171222
    phi <- eta2theta(eta[, 2], .lphi , .ephi )
    m1u <- if ( .stdbeta ) mu else (mu - .A) / ( .B - .A)
    dmu.deta <- dtheta.deta(mu, .lmu , .emu )
    dmu1.dmu <- 1 / ( .B - .A )
    dphi.deta <- dtheta.deta(phi, .lphi , .ephi )
    temp1 <- m1u*phi
    temp2 <- (1-m1u)*phi
    if ( .stdbeta ) {
      dl.dmu1 <- phi * (digamma(temp2) -
                        digamma(temp1) + log(y) - log1p(-y))
      dl.dphi <- digamma(phi) - mu*digamma(temp1) -
                 (1-mu)*digamma(temp2) +
                 mu*log(y) + (1-mu)*log1p(-y)
    } else {
      dl.dmu1 <- phi*(digamma(temp2) - digamma(temp1) +
                     log(y-.A) - log( .B-y))
      dl.dphi <- digamma(phi) - m1u*digamma(temp1) -
                (1-m1u)*digamma(temp2) +
                m1u*log(y-.A) + (1-m1u)*log( .B-y) - log( .B -.A)
    }
      c(w) * cbind(dl.dmu1 * dmu1.dmu * dmu.deta,
                   dl.dphi * dphi.deta)
  }), list( .lmu = lmu, .lphi = lphi,
            .emu = emu, .ephi = ephi,
            .A = A, .B = B,
            .stdbeta = stdbeta ))),
  weight = eval(substitute(expression({
    ned2l.dmu12 <- (trigamma(temp1) + trigamma(temp2)) * phi^2
    ned2l.dphi2 <- -trigamma(phi) + trigamma(temp1) * m1u^2 +
                    trigamma(temp2) * (1-m1u)^2
    ned2l.dmu1phi <- temp1 * trigamma(temp1) -
                     temp2 * trigamma(temp2)
    wz <- matrix(NA_real_, n, dimm(M))
    wz[, iam(1, 1, M)] <- ned2l.dmu12 * dmu1.dmu^2 * dmu.deta^2
    wz[, iam(2, 2, M)] <- ned2l.dphi2 * dphi.deta^2
    wz[, iam(1, 2, M)] <- ned2l.dmu1phi * dmu1.dmu *
                                          dmu.deta * dphi.deta
    c(w) * wz
  }), list( .A = A, .B = B ))))
}  # betaff







 betaR <-
  function(lshape1 = "loglink", lshape2 = "loglink",
           i1 = NULL, i2 = NULL, trim = 0.05,
           A = 0, B = 1, parallel = FALSE, zero = NULL) {

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

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


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

  if (!is.Numeric(A, length.arg = 1) ||
     !is.Numeric(B, length.arg = 1) ||
     A >= B)
    stop("A must be < B, and both must be of length one")

  stdbeta <- (A == 0 && B == 1)  # stdbeta == T iff



  new("vglmff",
  blurb = c("Two-parameter Beta distribution ",
            "(shape parameters parameterization)\n",
            if (stdbeta)
  paste("y^(shape1-1) * (1-y)^(shape2-1) / B(shape1,shape2),",
            "0 <= y <= 1, shape1>0, shape2>0\n\n") else
            paste0("(y-", A, ")^(shape1-1) * (", B,
            "-y)^(shape2-1) / [B(shape1,shape2) * (",
            B, "-", A, ")^(shape1+shape2-1)], ",
             A, " <= y <= ",B ," shape1>0, shape2>0\n\n"),
            "Links:    ",
            namesof("shape1", lshape1, earg = eshape1),  ", ",
            namesof("shape2", lshape2, earg = eshape2)),
  constraints = eval(substitute(expression({
    constraints <- cm.VGAM(matrix(1, M, 1), x = x,
                           bool = .parallel ,
                           constraints, apply.int  = TRUE)
    constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
                     M = M, M1 = 2,
                     predictors.names = predictors.names)
  }), list( .parallel = parallel, .zero = zero ))),
  infos = eval(substitute(function(...) {
    list(M1 = 2,
         Q1 = 1,
         A  = .A,
         B  = .B,
         multipleResponses = FALSE,
         zero = .zero )
  }, list( .A = A, .B = B,
           .zero = zero ))),

  rqresslot = eval(substitute(
    function(mu, y, w, eta, extra = NULL) {
    shapes <- cbind(eta2theta(eta[, 1], .lshape1 , .eshape1 ),
                    eta2theta(eta[, 2], .lshape2 , .eshape2 ))
    zedd <- (y - .A ) / ( .B - .A )
    scrambleseed <- runif(1)  # To scramble the seed
    qnorm(pbeta(zedd, shape1, shape2))
  }, list( .lshape1 = lshape1, .lshape2 = lshape2,
           .A = A, .B = B,
           .eshape1 = eshape1, .eshape2 = eshape2 ))),

  initialize = eval(substitute(expression({
    if (min(y) <= .A || max(y) >= .B )
      stop("data not within (A, B)")

    if (NCOL(y) != 1)
      stop("response must be a vector or a one-column matrix")


    w.y.check(w = w, y = y)


    predictors.names <-
        c(namesof("shape1", .lshape1 , .eshape1 , short = TRUE),
          namesof("shape2", .lshape2 , .eshape2 , short = TRUE))

    if (!length(etastart)) {
      mu1d <- mean(y, trim = .trim )
      uu <- (mu1d - .A) / ( .B - .A)
      DD <- ( .B - .A)^2
      pinit <- max(0.01, uu^2 * (1 - uu) * DD / var(y) - uu)
      qinit <- max(0.01, pinit * (1 - uu) / uu)
      etastart <- matrix(0, n, 2)
      etastart[, 1] <- theta2eta( pinit, .lshape1 , .eshape1 )
      etastart[, 2] <- theta2eta( qinit, .lshape2 , .eshape2 )
    }
    if (is.Numeric( .i1 ))
      etastart[, 1] <- theta2eta( .i1 , .lshape1 , .eshape1 )
    if (is.Numeric( .i2 ))
      etastart[, 2] <- theta2eta( .i2 , .lshape2 , .eshape2 )
  }), list( .lshape1 = lshape1, .lshape2 = lshape2,
            .i1 = i1, .i2 = i2, .trim = trim, .A = A, .B = B,
            .eshape1 = eshape1, .eshape2 = eshape2 ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    shapes <- cbind(eta2theta(eta[, 1], .lshape1 , .eshape1 ),
                    eta2theta(eta[, 2], .lshape2 , .eshape2 ))
    .A + ( .B - .A ) * shapes[, 1] / (shapes[, 1] + shapes[, 2])
  }, list( .lshape1 = lshape1, .lshape2 = lshape2, .A = A, .B = B,
           .eshape1 = eshape1, .eshape2 = eshape2 ))),
  last = eval(substitute(expression({
    misc$link <-    c(shape1 = .lshape1 , shape2 = .lshape2 )
    misc$earg <- list(shape1 = .eshape1 , shape2 = .eshape2 )
    misc$limits <- c( .A , .B )
  }), list( .lshape1 = lshape1, .lshape2 = lshape2,
            .A = A, .B = B,
            .eshape1 = eshape1, .eshape2 = eshape2 ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta, extra = NULL,
             summation = TRUE) {
    shapes <- cbind(eta2theta(eta[, 1], .lshape1 , .eshape1 ),
                    eta2theta(eta[, 2], .lshape2 , .eshape2 ))
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      zedd <- (y - .A ) / ( .B - .A )
      ll.elts <-
        c(w) * (dbeta(zedd, shapes[, 1], shapes[, 2],
                      log = TRUE) - log( abs( .B - .A )))
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
    }, list( .lshape1 = lshape1, .lshape2 = lshape2,
             .A = A, .B = B,
             .eshape1 = eshape1, .eshape2 = eshape2 ))),
  vfamily = "betaR",
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    shapes <- cbind(eta2theta(eta[, 1], .lshape1 , .eshape1 ),
                    eta2theta(eta[, 2], .lshape2 , .eshape2 ))
    okay1 <- all(is.finite(shapes)) && all(0 < shapes)
    okay1
  }, list( .lshape1 = lshape1, .lshape2 = lshape2,
           .A = A, .B = B,
           .eshape1 = eshape1, .eshape2 = eshape2 ))),



  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)
    shapes <- cbind(eta2theta(eta[, 1], .lshape1 , .eshape1 ),
                    eta2theta(eta[, 2], .lshape2 , .eshape2 ))
    .A + ( .B - .A ) *
    rbeta(nsim * length(shapes[, 1]),
          shape1 = shapes[, 1], shape2 = shapes[, 2])
  }, list( .lshape1 = lshape1, .lshape2 = lshape2,
           .A = A, .B = B,
           .eshape1 = eshape1, .eshape2 = eshape2 ))),



  deriv = eval(substitute(expression({
    shapes <- cbind(eta2theta(eta[, 1], .lshape1 , .eshape1 ),
                    eta2theta(eta[, 2], .lshape2 , .eshape2 ))

    dshapes.deta <-
      cbind(dtheta.deta(shapes[, 1], .lshape1 , .eshape1 ),
            dtheta.deta(shapes[, 2], .lshape2 , .eshape2 ))

    dl.dshapes <- cbind(log(y - .A ), log( .B - y)) -
                  digamma(shapes) +
                  digamma(shapes[, 1] + shapes[, 2]) -
                  log( .B - .A )

    c(w) * dl.dshapes * dshapes.deta
  }), list( .lshape1 = lshape1, .lshape2 = lshape2,
            .A = A, .B = B,
            .eshape1 = eshape1, .eshape2 = eshape2 ))),
  weight = expression({
    trig.sum <- trigamma(shapes[, 1] + shapes[, 2])
    ned2l.dshape12 <- trigamma(shapes[, 1]) - trig.sum
    ned2l.dshape22 <- trigamma(shapes[, 2]) - trig.sum
    ned2l.dshape1shape2 <- -trig.sum
    wz <- matrix(NA_real_, n, dimm(M))  # dimm(M) == 3
    wz[, iam(1, 1, M)] <- ned2l.dshape12      * dshapes.deta[, 1]^2
    wz[, iam(2, 2, M)] <- ned2l.dshape22      * dshapes.deta[, 2]^2
    wz[, iam(1, 2, M)] <- ned2l.dshape1shape2 * dshapes.deta[, 1] *
                                                dshapes.deta[, 2]
    c(w) * wz
  }))
}  # betaR





 betaprime <-
   function(lshape = "loglink", ishape1 = 2, ishape2 = NULL,
            zero = NULL) {

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


  new("vglmff",
  blurb = c("Beta-prime distribution\n",
  "y^(shape1-1) * (1+y)^(-shape1-shape2) / Beta(shape1,shape2),",
            " y>0, shape1>0, shape2>0\n\n",
            "Links:    ",
            namesof("shape1", lshape, earg = eshape),  ", ",
            namesof("shape2", lshape, earg = eshape), "\n",
            "Mean:     shape1/(shape2-1) provided shape2>1"),
  constraints = eval(substitute(expression({
    constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
                     M = M, M1 = 2,
                     predictors.names = predictors.names)
  }), list( .zero = zero ))),
  infos = eval(substitute(function(...) {
    list(M1 = 2,
         Q1 = 1,
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = c("shape1", "shape2"),
         lshape1 = .lshape ,
         lshape2 = .lshape ,
         zero = .zero )
  }, list( .zero = zero, .lshape = lshape ))),

  initialize = eval(substitute(expression({

    w.y.check(w = w, y = y,
              Is.positive.y = TRUE,
              ncol.w.max = 1,
              ncol.y.max = 1)



    predictors.names <-
      c(namesof("shape1", .lshape , earg = .eshape , short = TRUE),
        namesof("shape2", .lshape , earg = .eshape , short = TRUE))
    if (is.numeric( .ishape1) && is.numeric( .ishape2 )) {
      vec <- c( .ishape1, .ishape2 )
      vec <- c(theta2eta(vec[1], .lshape , earg = .eshape ),
               theta2eta(vec[2], .lshape , earg = .eshape ))
      etastart <- matrix(vec, n, 2, byrow = TRUE)
    }
    if (!length(etastart)) {
      init1 <- if (length( .ishape1 ))
        rep_len( .ishape1 , n) else rep_len(1, n)
      init2 <- if (length( .ishape2 ))
        rep_len( .ishape2 , n) else 1 + init1 / (y + 0.1)
      etastart <-
        matrix(theta2eta(c(init1, init2), .lshape , .eshape ),
               n, 2, byrow = TRUE)
    }
  }), list( .lshape = lshape, .eshape = eshape,
            .ishape1 = ishape1, .ishape2 = ishape2 ))),

  linkinv = eval(substitute(function(eta, extra = NULL) {
    shapes <- eta2theta(eta, .lshape , earg = .eshape )
    ifelse(shapes[, 2] > 1, shapes[, 1] / (shapes[, 2] - 1), NA)
  }, list( .lshape = lshape, .eshape = eshape ))),
  last = eval(substitute(expression({
    misc$link <-    c(shape1 = .lshape , shape2 = .lshape )
    misc$earg <- list(shape1 = .eshape , shape2 = .eshape )
  }), list( .lshape = lshape, .eshape = eshape ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    shapes <- eta2theta(eta, .lshape , earg = .eshape )
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <-
        c(w) * ((shapes[, 1]-1) * log(y) -
                 lbeta(shapes[, 1], shapes[, 2]) -
                (shapes[, 2] + shapes[, 1]) * log1p(y))
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .lshape = lshape, .eshape = eshape ))),
  vfamily = "betaprime",
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    shapes <- eta2theta(eta, .lshape , earg = .eshape )
    okay1 <- all(is.finite(shapes)) && all(0 < shapes)
    okay1
  }, list( .lshape = lshape, .eshape = eshape ))),
  deriv = eval(substitute(expression({
    shapes <- eta2theta(eta, .lshape , earg = .eshape )
    dshapes.deta <- dtheta.deta(shapes, .lshape , earg = .eshape )
    dl.dshapes <- cbind(log(y) - log1p(y) - digamma(shapes[, 1]) +
                        digamma(shapes[, 1] + shapes[, 2]),
                        - log1p(y) - digamma(shapes[, 2]) +
                        digamma(shapes[, 1] + shapes[, 2]))
    c(w) * dl.dshapes * dshapes.deta
  }), list( .lshape = lshape, .eshape = eshape ))),
  weight = expression({
    temp2 <- trigamma(shapes[, 1] + shapes[, 2])
    ned2l.dshape12 <- trigamma(shapes[, 1]) - temp2
    ned2l.dshape22 <- trigamma(shapes[, 2]) - temp2
    ned2l.dshape1shape2 <- -temp2

    wz <- matrix(NA_real_, n, dimm(M))  #3=dimm(M)
    wz[, iam(1, 1, M)] <- ned2l.dshape12 * dshapes.deta[, 1]^2
    wz[, iam(2, 2, M)] <- ned2l.dshape22 * dshapes.deta[, 2]^2
    wz[, iam(1, 2, M)] <- ned2l.dshape1shape2 *
                          dshapes.deta[, 1] * dshapes.deta[, 2]

    c(w) * wz
  }))
}  # betaprime












 zoabetaR <-
  function(lshape1 = "loglink", lshape2 = "loglink",
           lpobs0 = "logitlink", lpobs1 = "logitlink",
           ishape1 = NULL, ishape2 = NULL, trim = 0.05,
           type.fitted = c("mean", "pobs0", "pobs1", "beta.mean"),
           parallel.shape = FALSE,
           parallel.pobs = FALSE,
           zero = NULL) {

  A <- 0
  B <- 1

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

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

  lprobb0 <- as.list(substitute(lpobs0))
  eprobb0 <- link2list(lprobb0)
  lprobb0 <- attr(eprobb0, "function.name")

  lprobb1 <- as.list(substitute(lpobs1))
  eprobb1 <- link2list(lprobb1)
  lprobb1 <- attr(eprobb1, "function.name")

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

  if (!is.Numeric(A, length.arg = 1) ||
      !is.Numeric(B, length.arg = 1) ||
     A >= B)
    stop("A must be < B, and both must be of length one")

  stdbeta <- (A == 0 && B == 1)  # stdbeta==TRUE iff



  type.fitted <- match.arg(type.fitted,
                   c("mean", "pobs0", "pobs1", "beta.mean"))[1]


  new("vglmff",
  blurb = c("Standard Beta distribution with 0- and \n",
            "1-inflation ",
            "(shape parameters parameterization)\n",
            if (stdbeta)
  paste("y^(shape1-1) * (1-y)^(shape2-1) / beta(shape1,shape2),",
            "0 <= y <= 1, shape1>0, shape2>0\n\n") else
            paste("(y-",A,")^(shape1-1) * (",B,
            "-y)^(shape2-1) / [beta(shape1,shape2) * (",
            B, "-", A, ")^(shape1+shape2-1)], ",
             A," <= y <= ",B," shape1>0, shape2>0, ",
            "0 < pobs0 < 1, 0 < pobs1 < 1 \n\n", sep = ""),
            "Links:    ",
            namesof("shape1", lshape1, earg = eshape1),  ", ",
            namesof("shape2", lshape2, earg = eshape2),  ", ",
            namesof("pobs0",  lprobb0, earg = eprobb0),  ", ",
            namesof("pobs1",  lprobb1, earg = eshape1)),



  constraints = eval(substitute(expression({


  if (is.null(constraints)) {
 

    constraints.orig <- constraints


    if (is.logical( .parallel.probb ) && .parallel.probb &&
        (cind0[1] + cind1[1] <= 1))
      warning("argument 'parallel.pobs' specified when there is",
              " only one of 'pobs0' and 'pobs1'")




    cmk.s <- kronecker(matrix(1, NOS, 1), rbind(1, 1, 0, 0))
    cmk.S <- kronecker(diag(NOS), rbind(diag(2), 0*diag(2)))
    con.s <- cm.VGAM(cmk.s, x = x,
                     bool = .parallel.shape ,  # Same as .parallel.b
                     constraints = constraints.orig,
                     apply.int = TRUE,
                     cm.default           = cmk.S,
                     cm.intercept.default = cmk.S)



    cmk.p <- kronecker(matrix(1, NOS, 1), rbind(0, 0, 1, 1))
    cmk.P <- kronecker(diag(NOS), rbind(0*diag(2), diag(2)))
    con.p <- cm.VGAM(cmk.p,
                     x = x,
                     bool = .parallel.probb ,  #
                     constraints = constraints.orig,
                     apply.int = TRUE,
                     cm.default           = cmk.P,
                     cm.intercept.default = cmk.P)

    con.use <- con.s
    for (klocal in seq_along(con.s)) {
      con.use[[klocal]] <-
        cbind(con.s[[klocal]], con.p[[klocal]])
 # Delete rows that are not needed:
      if (!cind0[1]) {
        vec.use <- rep_len(c(TRUE, TRUE, FALSE, TRUE),
                           nrow(con.use[[klocal]]))
        con.use[[klocal]] <-
       (con.use[[klocal]])[vec.use, ]
      }
      if (!cind1[1]) {
        vec.use <- rep_len(c(TRUE, TRUE, TRUE, FALSE),
                           nrow(con.use[[klocal]]))
        con.use[[klocal]] <-
       (con.use[[klocal]])[vec.use, ]
      }
      col.delete <- apply(con.use[[klocal]], 2,
                          function(HkCol) all(HkCol == 0))
      con.use[[klocal]] <- (con.use[[klocal]])[, !col.delete]
    }


    constraints <- con.use
  }  # if (is.null(constraints))




    constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
                     M = M, M1 = M1,
                     predictors.names = predictors.names)

  }), list( .parallel.shape = parallel.shape,
            .parallel.probb = parallel.pobs,
            .zero = zero ))),
  infos = eval(substitute(function(...) {
    list(M1 = NA,  # Either 3 or 4, data-dependent
         Q1 = 1,
         A  = .A ,
         B  = .B ,
         expected = TRUE,
         multipleResponses = TRUE,
         type.fitted  = .type.fitted ,
         zero = .zero )
  }, list( .A = A, .B = B,
           .type.fitted = type.fitted,
           .zero = zero ))),
  initialize = eval(substitute(expression({
    if (min(y) <  .A || max(y)  > .B)
      stop("data not within [A, B]")

    temp5 <-
    w.y.check(w = w, y = y,
              ncol.w.max = Inf,
              ncol.y.max = Inf,
              out.wy = TRUE,
              colsyperw = 1,
              maximize = TRUE)
    w <- temp5$w
    y <- temp5$y
    ncoly <- NOS <- ncol(y)
    if (ncoly > 1 && !( .stdbeta ))
      stop("can only input multiple responses with the ",
           "standard beta")

    cind0 <- colSums(ind0 <- y == 0) > 0
    cind1 <- colSums(ind1 <- y == 1) > 0
    if (!any(cind0 | cind1))
      stop("no 0s or 1s in the responses to perform 0- and/or ",
           "1-inflation! ",
           "Try using betaff() or betaR() instead.")

    if (ncoly > 1 && !all(cind0 == cind0[1]) &&  # FALSE &&
                     !all(cind0 == cind0[1]))
      stop("with multiple responses, cannot have 0-inflation in ",
           "some responses and 1-inflation in other responses")
    M1 <- 2 + cind0[1] + cind1[1]  # 4 when
    M <- M1 * NOS

    mynames1 <- param.names("shape1", ncoly, skip1 = TRUE)
    mynames2 <- param.names("shape2", ncoly, skip1 = TRUE)
    mynames3 <- param.names("pobs0",  ncoly, skip1 = TRUE)
    mynames4 <- param.names("pobs1",  ncoly, skip1 = TRUE)
    predictors.names <-
        c(namesof(mynames1, .lshape1 , earg = .eshape1 ,
                  short = TRUE),
          namesof(mynames2, .lshape2 , earg = .eshape2 ,
                  short = TRUE),
        if (cind0[1])
            namesof(mynames3, .lprobb0 , earg = .eprobb0 ,
                    short = TRUE) else
        NULL,
        if (cind1[1])
            namesof(mynames4, .lprobb1 , earg = .eprobb1 ,
                    short = TRUE) else
        NULL)[interleave.VGAM(M, M1 = M1)]

    extra$type.fitted <- .type.fitted
    extra$colnames.y  <- colnames(y)
    extra$M1          <- M1  # Determined from the data
    extra$cind0       <- cind0
    extra$cind1       <- cind1

    if (!length(etastart)) {

      p0init <- matrix(colMeans(ind0), n, ncoly, byrow = TRUE)
      p1init <- matrix(colMeans(ind1), n, ncoly, byrow = TRUE)


      mu1d <- matrix(NA_real_, n, NOS)
      for (jay in 1:ncoly) {
        yy <- y[, jay]
        yy <- yy[ .A < yy & yy < .B ]
        mu1d[, jay] <- weighted.mean(yy, trim = .trim )
      }
      uu <- (mu1d - .A ) / ( .B - .A )
      DD <- ( .B - .A )^2
      p.init <- if (is.Numeric( .ishape1 ))
        matrix( .ishape1 , n, ncoly, byrow = TRUE) else
        uu^2 * (1 - uu) * DD / var(yy) - uu
      p.init[p.init < 0.01] <- 0.01
      q.init <- if (is.Numeric( .ishape2 ))
        matrix( .ishape2 , n, ncoly, byrow = TRUE) else
        p.init * (1 - uu) / uu
      q.init[q.init < 0.01] <- 0.01
      etastart <- cbind(
        theta2eta(p.init, .lshape1 , earg = .eshape1 ),
        theta2eta(q.init, .lshape2 , earg = .eshape2 ),
        if (cind0[1])
        theta2eta(p0init, .lprobb0 , earg = .eprobb0 )
        else NULL,
        if (cind1[1])
        theta2eta(p1init, .lprobb1 , earg = .eprobb1 )
        else NULL)[,
          interleave.VGAM(M, M1 = M1)]
    }
  }), list( .lshape1 = lshape1, .lshape2 = lshape2,
            .eshape1 = eshape1, .eshape2 = eshape2,
            .lprobb0 = lprobb0, .lprobb1 = lprobb1,
            .eprobb0 = eprobb0, .eprobb1 = eprobb1,
            .ishape1 = ishape1, .ishape2 = ishape2,
            .trim = trim, .A = A, .B = B,
            .type.fitted = type.fitted,
            .stdbeta = stdbeta ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    M1 <- extra$M1
    cind0 <- extra$cind0
    cind1 <- extra$cind1
    NOS <- ncol(eta) / M1
      shape1 <- eta2theta(eta[, c(TRUE, rep(FALSE, M1 - 1)),
                              drop = FALSE],
                        .lshape1 , earg = .eshape1 )
    shape2 <- eta2theta(eta[, c(FALSE, TRUE,
                                rep(FALSE, M1 - 2)),
                            drop = FALSE],
                        .lshape2 , earg = .eshape2 )
    probb0 <- if (cind0[1])
      eta2theta(eta[, c(FALSE, FALSE, TRUE,
                        if (cind1[1]) FALSE else NULL),
                    drop = FALSE],
                .lprobb0 , earg = .eprobb0 ) else 0
    probb1 <- if (cind1[1])
      eta2theta(eta[, c(FALSE, FALSE,
                        if (cind0[1]) FALSE else NULL, TRUE),
                    drop = FALSE],
                .lprobb1 , earg = .eprobb1 ) else 0

    type.fitted <- match.arg(extra$type.fitted,
                     c("mean", "pobs0", "pobs1", "beta.mean"))[1]

    ans <-
      switch(type.fitted,
      "mean"      = (1 - probb0) * shape1 / (shape1 + shape2) +
                         probb1  * shape2 / (shape1 + shape2),
             "beta.mean" = shape1 / (shape1 + shape2),
             "pobs0"     = probb0,
             "pobs1"     = probb1)
    label.cols.y(ans, colnames.y = extra$colnames.y, NOS = NOS)
  },
  list( .lshape1 = lshape1, .lshape2 = lshape2, .A = A, .B = B,
        .eshape1 = eshape1, .eshape2 = eshape2,
        .lprobb0 = lprobb0, .lprobb1 = lprobb1,
        .eprobb0 = eprobb0, .eprobb1 = eprobb1 ))),
  last = eval(substitute(expression({
    misc$link <- rep_len( c( .lshape1 , .lshape2 ,
                             if (cind0[1]) .lprobb0 else NULL,
                             if (cind1[1]) .lprobb1 else NULL), M)
    names(misc$link) <- c(mynames1, mynames2,
                          if (cind0[1]) mynames3 else NULL,
                          if (cind1[1]) mynames4 else NULL)[
                        interleave.VGAM(M, M1 = M1)]

    misc$earg <- vector("list", M)
    names(misc$earg) <- names(misc$link)
    jay <- 1
    while (jay <= M) {
      misc$earg[[jay]] <- .eshape1
      jay <- jay + 1
      misc$earg[[jay]] <- .eshape2
      jay <- jay + 1
      if (cind0[1]) {
        misc$earg[[jay]] <- .eprobb0
        jay <- jay + 1
      }
      if (cind1[1]) {
        misc$earg[[jay]] <- .eprobb1
        jay <- jay + 1
      }
    }

    misc$supportlimits <- c( .A , .B )
  }), list( .lshape1 = lshape1, .lshape2 = lshape2,
            .eshape1 = eshape1, .eshape2 = eshape2,
            .lprobb0 = lprobb0, .lprobb1 = lprobb1,
            .eprobb0 = eprobb0, .eprobb1 = eprobb1,
            .A = A, .B = B ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta, extra = NULL,
             summation = TRUE) {
    M1 <- 4
    M1 <- extra$M1
    cind0 <- extra$cind0
    cind1 <- extra$cind1
    NOS <- ncol(eta) / M1
    shape1 <- eta2theta(eta[, c(TRUE, rep(FALSE, M1 - 1)),
                            drop = FALSE],
                        .lshape1 , earg = .eshape1 )
    shape2 <- eta2theta(eta[, c(FALSE, TRUE,
                                rep(FALSE, M1 - 2)),
                            drop = FALSE],
                        .lshape2 , earg = .eshape2 )
    probb0 <- if (cind0[1])
      eta2theta(eta[, c(FALSE, FALSE, TRUE,
                        if (cind1[1]) FALSE else NULL),
                    drop = FALSE],
                .lprobb0 , earg = .eprobb0 ) else 0
    probb1 <- if (cind1[1])
      eta2theta(eta[, c(FALSE, FALSE,
                        if (cind0[1]) FALSE else NULL, TRUE),
                    drop = FALSE],
                .lprobb1 , earg = .eprobb1 ) else 0

    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      zedd <- (y - .A ) / ( .B - .A )
      ll.elts <-
        c(w) * (dzoabeta(zedd, shape1 = shape1, shape2 = shape2,
                         pobs0 = probb0, pobs1 = probb1,
                      log = TRUE) - log( abs( .B - .A )))
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .lshape1 = lshape1, .lshape2 = lshape2, .A = A, .B = B,
           .eshape1 = eshape1, .eshape2 = eshape2,
           .lprobb0 = lprobb0, .lprobb1 = lprobb1,
           .eprobb0 = eprobb0, .eprobb1 = eprobb1 ))),
  vfamily = "zoabetaR",









  validparams = eval(substitute(function(eta, y, extra = NULL) {
    M1 <- 4
    M1 <- extra$M1
    cind0 <- extra$cind0
    cind1 <- extra$cind1
    NOS <- ncol(eta) / M1
    shape1 <- eta2theta(eta[, c(TRUE, rep(FALSE, M1 - 1)),
                            drop = FALSE],
                        .lshape1 , earg = .eshape1 )
    shape2 <- eta2theta(eta[, c(FALSE, TRUE,
                        rep(FALSE, M1 - 2)), drop = FALSE],
                        .lshape2 , earg = .eshape2 )
    probb0 <- if (cind0[1])
      eta2theta(eta[, c(FALSE, FALSE, TRUE,
                  if (cind1[1]) FALSE else NULL), drop = FALSE],
                .lprobb0 , earg = .eprobb0 ) else 0.5
    probb1 <- if (cind1[1])
      eta2theta(eta[, c(FALSE, FALSE,
            if (cind0[1]) FALSE else NULL, TRUE), drop = FALSE],
                .lprobb1 , earg = .eprobb1 ) else 0.5
    okay1 <- all(is.finite(shape1)) && all(0 < shape1) &&
       all(is.finite(shape2)) && all(0 < shape2) &&
       all(is.finite(probb0)) && all(0 < probb0 & probb0 < 1) &&
       all(is.finite(probb1)) && all(0 < probb1 & probb1 < 1)
    okay1
  }, list( .lshape1 = lshape1, .lshape2 = lshape2, .A = A, .B = B,
           .eshape1 = eshape1, .eshape2 = eshape2,
           .lprobb0 = lprobb0, .lprobb1 = lprobb1,
           .eprobb0 = eprobb0, .eprobb1 = eprobb1 ))),



  deriv = eval(substitute(expression({
    M1 <- 4
    M1 <- extra$M1
    cind0 <- extra$cind0
    cind1 <- extra$cind1
    NOS <- ncol(eta) / M1
    shape1 <- eta2theta(eta[, c(TRUE, rep(FALSE, M1 - 1)),
                            drop = FALSE],
                        .lshape1 , earg = .eshape1 )
    shape2 <- eta2theta(eta[, c(FALSE, TRUE,
                                rep(FALSE, M1 - 2)),
                            drop = FALSE],
                        .lshape2 , earg = .eshape2 )
    probb0 <- if (cind0[1])
      eta2theta(eta[, c(FALSE, FALSE, TRUE,
                        if (cind1[1]) FALSE else NULL),
                    drop = FALSE],
                .lprobb0 , earg = .eprobb0 ) else 0
    probb1 <- if (cind1[1])
      eta2theta(eta[, c(FALSE, FALSE,
                        if (cind0[1]) FALSE else NULL, TRUE),
                    drop = FALSE],
                .lprobb1 , earg = .eprobb1 ) else 0


    dshape1.deta <- dtheta.deta(shape1, .lshape1 , earg = .eshape1 )
    dshape2.deta <- dtheta.deta(shape2, .lshape2 , earg = .eshape2 )
    dprobb0.deta <- dtheta.deta(probb0, .lprobb0 , earg = .eprobb0 )
    dprobb1.deta <- dtheta.deta(probb1, .lprobb1 , earg = .eprobb1 )

    index0 <- y == 0
    index1 <- y == 1
    indexi <- !index0 & !index1  # In the interior, i.e., (0, 1)
    dig.sum <- digamma(shape1 + shape2)
    QQ <- 1 - probb0 - probb1

    if (cind0[1]) {
      dl.dprobb0 <- -1 / QQ
      dl.dprobb0[index0] <- 1 / probb0[index0]
      dl.dprobb0[index1] <- 0
    }

    if (cind1[1]) {
      dl.dprobb1 <- -1 / QQ
      dl.dprobb1[index0] <- 0
      dl.dprobb1[index1] <- 1 / probb1[index1]
    }

    dl.dshape1 <-    log(y) - digamma(shape1) + dig.sum
    dl.dshape2 <- log1p(-y) - digamma(shape2) + dig.sum
    dl.dshape1[!indexi] <- 0
    dl.dshape2[!indexi] <- 0

    myderiv <- c(w) *
      cbind(dl.dshape1 * dshape1.deta,
            dl.dshape2 * dshape2.deta,
            if (cind0[1]) dl.dprobb0 * dprobb0.deta else NULL,
            if (cind1[1]) dl.dprobb1 * dprobb1.deta else NULL)
    colnames(myderiv) <- NULL
    myderiv[, interleave.VGAM(M, M1 = M1)]
  }),
  list( .lshape1 = lshape1, .lshape2 = lshape2, .A = A, .B = B,
        .eshape1 = eshape1, .eshape2 = eshape2,
        .lprobb0 = lprobb0, .lprobb1 = lprobb1,
        .eprobb0 = eprobb0, .eprobb1 = eprobb1 ))),
  weight = expression({
    trig.sum <- trigamma(shape1 + shape2)


ned2l.dshape12 <- (trigamma(shape1) - trig.sum) * QQ
ned2l.dshape22 <- (trigamma(shape2) - trig.sum) * QQ
ned2l.dprobb02 <- (1 - probb1) / (probb0 * QQ)
ned2l.dprobb12 <- (1 - probb0) / (probb1 * QQ)
ned2l.dshape1shape2 <- -trig.sum * QQ  # (1 - probb0 - probb0) zz
ned2l.dshape2probb0 <- 0
ned2l.dprobb0probb1 <- 1 / QQ
ned2l.dshape1probb0 <- 0
ned2l.dshape2probb1 <- 0
ned2l.dshape1probb1 <- 0

ned2l.dshape1probb0 <- 0
    wz <- array(c(c(w) * ned2l.dshape12 * dshape1.deta^2,
                  c(w) * ned2l.dshape22 * dshape2.deta^2,
 if (cind0[1])  c(w) * ned2l.dprobb02 * dprobb0.deta^2 else NULL,
 if (cind1[1])  c(w) * ned2l.dprobb12 * dprobb1.deta^2 else NULL,
                c(w) * ned2l.dshape1shape2 * dshape1.deta *
                                             dshape2.deta,
 if (cind0[1])  c(w) * ned2l.dshape2probb0 * dshape2.deta *
                                             dprobb0.deta,
                c(w) * ned2l.dprobb0probb1 * dprobb0.deta *
                                             dprobb1.deta,
 if (cind0[1])  c(w) * ned2l.dshape1probb0 * dshape1.deta *
                                             dprobb0.deta,
 if (cind1[1])  c(w) * ned2l.dshape2probb1 * dshape2.deta *
                                             dprobb1.deta,
 if (cind1[1])  c(w) * ned2l.dshape1probb1 * dshape1.deta *
                                             dprobb1.deta),
                dim = c(n, M / M1, M1*(M1+1)/2))

    wz <- arwz2wz(wz, M = M, M1 = M1)  # tridiagonal
    wz
  }))
}  # zoabetaR






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

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

  logdensity <- rep_len(log(0), L)
  xok <- (0 <= x) & (x <= 1)
  logdensity[xok] <-
    log(2) + log(shape[xok]) + log1p(-x[xok]) +
    (shape[xok] - 1) * (log(x[xok]) + log(2) + log1p(-x[xok]/2))
  logdensity[shape >= 1] <- NaN
  if (log.arg) logdensity else exp(logdensity)
}



 ptopple <-
  function(q, 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 <- shape * (log(q) + log(2) + log1p(-q/2))
      ans[q <= 0 ] <- -Inf
      ans[q >= 1] <- 0
    } else {
      ans <- (q * (2 - q))^shape
      ans[q <= 0] <- 0
      ans[q >= 1] <- 1
    }
  } else {
    if (log.p) {
      ans <- log1p(-(q * (2 - q))^shape)
      ans[q <= 0] <- 0
      ans[q >= 1] <- -Inf
    } else {
      ans <- exp(log1p(-(q * (2 - q))^shape))
      ans[q <= 0] <- 1
      ans[q >= 1] <- 0
    }
  }
  ans[shape <= 0] <- NaN
  ans[shape >= 1] <- NaN
  ans
}



qtopple <- function(p, shape) {
  ans <- -expm1(0.5 * log1p(-p^(1/shape)))
  ans[shape <= 0] <- NaN
  ans[shape >= 1] <- NaN
  ans
}



rtopple <- function(n, shape) {
  qtopple(runif(n), shape)
}





 topple <-
  function(lshape = "logitlink", zero = NULL,
           gshape = ppoints(8),
           parallel = FALSE,
           type.fitted = c("mean", "percentiles", "Qlink"),
           percentiles = 50) {


  type.fitted <- match.arg(type.fitted,
                           c("mean", "percentiles", "Qlink"))[1]

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



  new("vglmff",
  blurb = c("Topp-Leone distribution ",
            "F(y; shape) = (y * (2 - y))^shape, ",
            "0 < y < 1, 0 < shape < 1\n",
            "Link:    ",
            namesof("shape", lshape, earg = eshape)),
  constraints = eval(substitute(expression({
    constraints <- cm.VGAM(matrix(1, M, 1), x = x,
                           bool = .parallel ,
                           constraints, apply.int = FALSE)
    constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
                     M = M, M1 = 1,
                     predictors.names = predictors.names)
  }), list( .parallel = parallel,
            .zero = zero ))),

  infos = eval(substitute(function(...) {
    list(M1 = 1,
         Q1 = 1,
         dpqrfun = "topple",
         expected = TRUE,
         hadof = TRUE,
         multipleResponses = TRUE,
         parallel = .parallel ,
         parameters.names = "shape",
         percentiles = .percentiles ,
         type.fitted = .type.fitted ,
         zero = .zero )
  }, list( .parallel = parallel,
           .percentiles = percentiles ,
           .type.fitted = type.fitted,
           .zero = zero ))),

  rqresslot = eval(substitute(
    function(mu, y, w, eta, extra = NULL) {
    Shape <- eta2theta(eta, .lshape , earg = .eshape )
    scrambleseed <- runif(1)  # To scramble the seed
    qnorm(runif(length(y),
                ptopple(y - 1, Shape),
                ptopple(y    , Shape)))
  }, list( .lshape = lshape, .eshape = eshape ))),


  initialize = eval(substitute(expression({
    temp5 <-
    w.y.check(w = w, y = y,
              Is.positive.y = TRUE,
              ncol.w.max = Inf,
              ncol.y.max = Inf,
              out.wy = TRUE,
              colsyperw = 1,
              maximize = TRUE)
    w <- temp5$w
    y <- temp5$y
    if (any(y >= 1))
      stop("response must be in (0, 1)")



    ncoly <- ncol(y)
    M1 <- 1
    M <- M1 * ncoly
    extra$ncoly <- ncoly
    extra$type.fitted <- .type.fitted
    extra$colnames.y  <- colnames(y)
    extra$percentiles <- .percentiles
    extra$M1 <- M1


    if ((NOS <- M / M1) > 1 && length( .percentiles ) > 1)
      stop("can only have one response when 'percentiles' is a ",
           "vector longer than unity")



    mynames1  <- param.names("shape", ncoly, skip1 = TRUE)
    predictors.names <-
      namesof(mynames1, .lshape , earg = .eshape , tag = FALSE)


    if (!length(etastart)) {
      shape.init <- matrix(0, nrow(x), ncoly)
      gshape <- .gshape
      topple.Loglikfun <- function(shape, y, x = NULL,
                                   w, extraargs = NULL) {
        sum(c(w) * dtopple(x = y, shape = shape, log = TRUE))
      }

      for (jay in 1:ncoly) {
        shape.init[, jay] <- grid.search(gshape,
                                         objfun = topple.Loglikfun,
                                         y = y[, jay], w = w[, jay])
      }
      etastart <- theta2eta(shape.init, .lshape , earg = .eshape )
    }
  }), list( .lshape = lshape, .gshape = gshape,
            .eshape = eshape,
            .percentiles = percentiles,
            .type.fitted = type.fitted
           ))),

  linkinv = eval(substitute(function(eta, extra = NULL) {
    type.fitted <-
      if (length(extra$type.fitted)) {
        extra$type.fitted
      } else {
        warning("cannot find 'type.fitted'. Returning the 'mean'.")
        "mean"
      }
    type.fitted <- match.arg(type.fitted,
                             c("mean", "percentiles", "Qlink"))[1]

    if (type.fitted == "Qlink") {
      eta2theta(eta, link = "logitlink")
    } else {
      shape <- eta2theta(eta, .lshape , earg = .eshape )
      pcent <- extra$percentiles
      perc.mat <- matrix(pcent, NROW(eta), length(pcent),
                         byrow = TRUE) / 100
      fv <-
        switch(type.fitted,
               "mean" = 1 - (gamma(1 + shape))^2 *
                            4^shape / gamma(2 * (1 + shape)),
               "percentiles" = qtopple(perc.mat,
                        shape = matrix(shape, nrow(perc.mat),
                                              ncol(perc.mat))))
      if (type.fitted == "percentiles")
        fv <- label.cols.y(fv, colnames.y = extra$colnames.y,
                           NOS = NCOL(eta), percentiles = pcent,
                           one.on.one = FALSE)
      fv
    }
  }, list( .lshape = lshape, .eshape = eshape ))),


  last = eval(substitute(expression({
    misc$earg <- vector("list", M)
    names(misc$earg) <- mynames1
    for (ilocal in 1:ncoly) {
      misc$earg[[ilocal]] <- .eshape
    }

    misc$link <- rep_len( .lshape , ncoly)
    names(misc$link) <- mynames1
  }), list( .lshape = lshape, .eshape = eshape ))),

  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * dtopple(x = y, shape = shape, log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .lshape = lshape, .eshape = eshape ))),
  vfamily = c("topple"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    okay1 <- all(is.finite(shape)) && all(0 < shape & shape < 1)
    okay1
  }, list( .lshape = lshape, .eshape = eshape ))),




  hadof = eval(substitute(
  function(eta, extra = list(), deriv = 1,
           linpred.index = 1, w = 1,
           dim.wz = c(NROW(eta), NCOL(eta) * (NCOL(eta)+1)/2),
           ...) {
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    ans <- c(w) *
    switch(as.character(deriv),
           "0" =   1 / shape^2,
           "1" =  -2 / shape^3,
           "2" =   6 / shape^4,
           "3" = -24 / shape^5,
           stop("argument 'deriv' must be 0, 1, 2 or 3"))
    if (deriv == 0)
      ans else retain.col(ans, linpred.index)  # M1 = 1
  }, list( .lshape = lshape, .eshape = eshape ))),




  simslot = eval(substitute(
  function(object, nsim) {
    pwts <- if (length(pwts <- object@prior.weights) > 0)
              pwts else weights(object, type = "prior")
    if (any(pwts != 1))
      warning("ignoring prior weights")
    eta <- predict(object)
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    rtopple(nsim * length(shape), shape = c(shape))
  }, list( .lshape = lshape,
           .eshape = eshape ))),


  deriv = eval(substitute(expression({
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    dl.dshape <- 1 / shape + log(y) + log(2) + log1p(-y/2)
    dshape.deta <- dtheta.deta(shape, .lshape , earg = .eshape )
    c(w) * dl.dshape * dshape.deta
  }), list( .lshape = lshape, .eshape = eshape ))),
  weight = eval(substitute(expression({
    ned2l.dshape2 <- 1 / shape^2
    wz <- c(w) * ned2l.dshape2 * dshape.deta^2
    wz
  }), list( .lshape = lshape, .eshape = eshape ))))
}  # topple





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



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

  ox <- !is.finite(x)
  zero <- ox | round(x) != x | x < 1
  ans <- rep_len(if (log.arg) log(0) else 0, LLL)
  if (any(!zero)) {
    if (log.arg) {
      ans[!zero] <- (-shape[!zero]-1) * log(x[!zero]) -
                    log(zeta(shape[!zero] + 1))
    } else {
      ans[!zero] <- x[!zero]^(-shape[!zero]-1) / zeta(
                                         shape[!zero]+1)
    }
  }
  if (any(ox))
    ans[ox] <- if (log.arg) log(0) else 0
  ans[shape <= 0] <- NaN  # Added 20160617
  ans
}



 pzeta <- function(q, shape, lower.tail = TRUE) {


  LLL <- max(lenq <- length(q), lens <- length(shape))
  if (length(q)     != LLL) q     <- rep_len(q,     LLL)
  if (length(shape) != LLL) shape <- rep_len(shape, LLL)
  ans <- rep_len(0, LLL)

  aa <- 12  # Same as Zeta.aux()
  qfloor <- floor(q)
  for (nn in 1:(aa-1))
    ans <- ans + as.numeric(nn <= qfloor) / nn^(shape+1)

  vecTF <- (aa-1 <= qfloor)
  if (lower.tail) {
    if (any(vecTF))
      ans[vecTF] <- zeta(shape[vecTF]+1) -
                    Zeta.aux(shape[vecTF]+1, qfloor[vecTF]  )
  } else {
    ans <- zeta(shape+1) - ans
    if (any(vecTF))
      ans[vecTF] <- Zeta.aux(shape[vecTF]+1, qfloor[vecTF]  )
  }
  ans / zeta(shape+1)
}  # pzeta



 qzeta <- function(p, shape) {

  LLL <- max(lenp <- length(p), lens <- length(shape))
  if (length(p)     != LLL) p     <- rep_len(p,     LLL)
  if (length(shape) != LLL) shape <- rep_len(shape, LLL)
  ans <- rep_len(0, LLL)

  # First, bracket the solution between 'lo' and 'hi'.
  lowsup <- 1
  lo <- rep_len(lowsup - 0.5, LLL)
  approx.ans <- lo  # True at lhs
  hi <- 2 * lo + 10.5
  dont.iterate <- p == 1 | shape <= 0
  done <- p <= pzeta(hi, shape) | dont.iterate
  while (!all(done)) {
    lo[!done] <- hi[!done]
    hi[!done] <- 2 * hi[!done] + 10.5  # Bug fixed
    done[!done] <- (p[!done] <= pzeta(hi[!done], shape[!done]))
  }

  foo <- function(q, shape, p)
    pzeta(q, shape) - p

  lhs <- (p <= dzeta(1, shape)) | dont.iterate

  approx.ans[!lhs] <-
    bisection.basic(foo, lo[!lhs], hi[!lhs], tol = 1/16,
                    shape = shape[!lhs], p = p[!lhs])
  faa <- floor(approx.ans)
   ans <- ifelse(pzeta(faa, shape) < p &
                 p <= pzeta(faa+1, shape), faa+1, faa)

  ans[p == 1] <- Inf
  ans[shape <= 0] <- NaN

  ans
}  # qzeta



rzeta <- function(n, shape) {
  qzeta(runif(n), shape)
}





 zetaff <-
  function(lshape = "loglink",
           ishape = NULL,
           gshape = 1 + exp(-seq(7)),
           zero = NULL) {


  if (length(ishape) && !is.Numeric(ishape, positive = TRUE))
    stop("argument 'ishape' must be > 0")

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


  new("vglmff",
  blurb = c("Zeta distribution ",
    "f(y; shape) = 1/(y^(shape+1) zeta(shape+1)), ",
    "shape>0, y = 1, 2,..\n\n",
    "Link:    ",
    namesof("shape", lshape, earg = eshape), "\n\n",
    "Mean:     zeta(shape) / zeta(shape+1), provided shape>1\n",
    "Variance: zeta(shape-1) / zeta(shape+1) - mean^2, if shape>2"),
  infos = eval(substitute(function(...) {
    list(M1 = 1,
         Q1 = 1,
         dpqrfun = "zeta",
         multipleResponses = TRUE,
         parameters.names = "shape",
         zero = .zero ,
         lshape = .lshape )
  }, list( .lshape = lshape,
           .zero = zero ))),

  rqresslot = eval(substitute(
    function(mu, y, w, eta, extra = NULL) {
    Shape <- eta2theta(eta, .lshape , earg = .eshape )
    scrambleseed <- runif(1)  # To scramble the seed
    qnorm(runif(length(y),
                pzeta(y - 1, Shape),
                pzeta(y    , Shape)))
  }, list( .lshape = lshape, .eshape = eshape ))),

  initialize = eval(substitute(expression({

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


    ncoly <- ncol(y)
    mynames1 <- param.names("shape", ncoly, skip1 = TRUE)
    predictors.names <-
      namesof(mynames1, .lshape , earg = .eshape , tag = FALSE)

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


    if (!length(etastart)) {
      zetaff.Loglikfun <- function(shape, y, x, w, extraargs) {
        sum(c(w) * dzeta(x = y, shape, log = TRUE))
      }


      gshape <- .gshape
      if (!length( .ishape )) {
        shape.init <- matrix(NA_real_, n, M, byrow = TRUE)
        for (jay in 1:ncoly) {
          shape.init[, jay] <- grid.search(gshape,
                               objfun = zetaff.Loglikfun,
                               y = y[, jay], x = x, w = w[, jay])
        }
      } else {
        shape.init <- matrix( .ishape , n, M, byrow = TRUE)
      }
      etastart <- theta2eta(shape.init, .lshape , earg = .eshape )
    }
  }), list( .lshape = lshape, .eshape = eshape,
            .ishape = ishape, .gshape = gshape ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    ans <- pp <- eta2theta(eta, .lshape , earg = .eshape )
    ans[pp > 1] <- zeta(pp[pp > 1]) / zeta(pp[pp > 1] + 1)
    ans[pp <= 1] <- NA
    ans
  }, list( .lshape = lshape, .eshape = eshape ))),
  last = eval(substitute(expression({
    misc$link <- rep_len( .lshape , ncoly)
    names(misc$link) <- mynames1

    misc$earg <- vector("list", M)
    names(misc$earg) <- mynames1
    for (jay in 1:ncoly) {
      misc$earg[[jay]] <- .eshape
    }

  }), list( .lshape = lshape, .eshape = eshape ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE,
             eta, extra = NULL, summation = TRUE) {
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * dzeta(x = y, shape, log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .lshape = lshape, .eshape = eshape ))),
  vfamily = c("zetaff"),


  hadof = eval(substitute(
  function(eta, extra = list(), deriv = 1,
           linpred.index = 1, w = 1,
           dim.wz = c(NROW(eta), NCOL(eta) * (NCOL(eta)+1)/2),
           ...) {
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    fred0 <- zeta(shape+1)
    fred1 <- zeta(shape+1, deriv = 1)
    fred2 <- zeta(shape+1, deriv = 2)

    ans <- c(w) *
    switch(as.character(deriv),
           "0" = fred2 / fred0 - (fred1/fred0)^2,
           "1" = (zeta(shape + 1, deriv = 3) -  # Curr. unavailable
            fred2 * fred1 / fred0) / fred0 -
            2 * (fred1 / fred0) * (
            fred2 /fred0 - (fred1/fred0)^2),
           "2" = NA * theta,
           "3" = NA * theta,
           stop("argument 'deriv' must be 0, 1, 2 or 3"))
    if (deriv == 0) ans else
      retain.col(ans, linpred.index)  # Since M1 = 1
  }, list( .lshape = lshape, .eshape = eshape ))),



  validparams = eval(substitute(function(eta, y, extra = NULL) {
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    okay1 <- all(is.finite(shape)) && all(0 < shape)
    okay1
  }, list( .lshape = lshape, .eshape = eshape ))),
  deriv = eval(substitute(expression({
    shape <- eta2theta(eta, .lshape , earg = .eshape )

    fred0 <- zeta(shape + 1)
    fred1 <- zeta(shape + 1, deriv = 1)
    dl.dshape <- -log(y) - fred1 / fred0
    dshape.deta <- dtheta.deta(shape, .lshape , earg = .eshape )
    c(w) * dl.dshape * dshape.deta
  }), list( .lshape = lshape, .eshape = eshape ))),
  weight = expression({
    NOS <- NCOL(y)
    ned2l.dshape2 <- zeta(shape + 1, deriv = 2) / fred0 -
                     (fred1 / fred0)^2
    wz <- ned2l.dshape2 * dshape.deta^2
    w.wz.merge(w = w, wz = wz, n = n, M = M, ndepy = NOS)
  }))
}  # zetaff



gharmonic2 <- function(n, shape = 1) {



  if (!is.Numeric(n, integer.valued = TRUE, positive = TRUE))
    stop("bad input for argument 'n'")

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

  aa <- 12
  ans <- rep_len(0, LLL)
  for (ii in 1:aa)
     ans <- ans + as.numeric(ii <= n) / ii^shape

  vecTF <- (aa < n)
  if (any(vecTF))
    ans[vecTF] <- zeta(shape[vecTF]) -
                  Zeta.aux(shape[vecTF], 1 + n[vecTF])
  ans
}  # gharmonic2



gharmonic <- function(n, shape = 1, deriv = 0) {





  if (!is.Numeric(n, integer.valued = TRUE, positive = TRUE))
      stop("bad input for argument 'n'")
  if (!is.Numeric(deriv, length.arg = 1, integer.valued = TRUE) ||
      deriv < 0)
    stop("bad input for argument 'deriv'")

  lognexponent <- deriv
  sign <- ifelse(deriv %% 2 == 0, 1, -1)

  ans <-
  if (length(n) == 1 && length(shape) == 1) {
      if (lognexponent != 0)
        sum(log(1:n)^lognexponent * (1:n)^(-shape)) else
      sum((1:n)^(-shape))
  } else {
    LEN <- max(length(n), length(shape))
    n <- rep_len(n, LEN)
    ans <- shape <- rep_len(shape, LEN)
    if (lognexponent != 0) {
      for (ii in 1:LEN)
        ans[ii] <- sum(log(1:n[ii])^lognexponent *
                      (1:n[ii])^(-shape[ii]))
    } else {
      for (ii in 1:LEN)
        ans[ii] <- sum((1:n[ii])^(-shape[ii]))
    }
    ans
  }
  sign * ans
}  # gharmonic






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


  if (!is.Numeric(x))
    stop("bad input for argument 'x'")
  if (!is.Numeric(N, integer.valued = TRUE, positive = TRUE))
    stop("bad input for argument 'N'")
  if (!is.Numeric(shape, positive = TRUE))
    stop("bad input for argument 'shape'")
  nn <- max(length(x), length(N), length(shape))
  if (length(x)     != nn) x     <- rep_len(x,     nn)
  if (length(N)     != nn) N     <- rep_len(N,     nn)
  if (length(shape) != nn) shape <- rep_len(shape, nn)

  ox <- !is.finite(x)
  zero <- ox | round(x) != x | x < 1 | x > N
  ans <- (if (log.arg) log(0) else 0) * x
  if (any(!zero))
    if (log.arg) {
      ans[!zero] <- (-shape[!zero]) * log(x[!zero]) -
                    log(gharmonic2(N[!zero], shape[!zero]))
    } else {
      ans[!zero] <- x[!zero]^(-shape[!zero]) / (
          gharmonic2(N[!zero], shape[!zero]))
    }
  ans
}  # dzipf



 pzipf <- function(q, N, shape, log.p = FALSE) {

  if (!is.Numeric(N, integer.valued = TRUE, positive = TRUE))
    stop("bad input for argument 'N'")

  nn <- max(length(q), length(N), length(shape))
  if (length(q)     != nn) q     <- rep_len(q,     nn)
  if (length(N)     != nn) N     <- rep_len(N,     nn)
  if (length(shape) != nn) shape <- rep_len(shape, nn)
    oq <- !is.finite(q)
  dont.iterate <- shape <= 0
    zeroOR1 <- oq | q < 1 | N <= q | dont.iterate
    floorq <- floor(q)
    ans <- 0 * floorq
    ans[oq | q >= N] <- 1
    if (any(!zeroOR1))
        ans[!zeroOR1] <-
          gharmonic2(floorq[!zeroOR1], shape[!zeroOR1]) /
          gharmonic2(     N[!zeroOR1], shape[!zeroOR1])

    ans[shape <= 0] <- NaN

    if (log.p) log(ans) else ans
}  # pzipf



qzipf <- function(p, N, shape) {
  if (!is.Numeric(p))
    stop("bad input for argument 'p'")
  if (!is.Numeric(N, integer.valued = TRUE, positive = TRUE))
    stop("bad input for argument 'N'")
  if (!is.Numeric(shape, positive = TRUE))
    stop("bad input for argument 'shape'")

  nn <- max(length(p), length(N), length(shape))
  if (length(p)     != nn) p     <- rep_len(p,     nn)
  if (length(N)     != nn) N     <- rep_len(N,     nn)
  if (length(shape) != nn) shape <- rep_len(shape, nn)

  a <- rep_len(1, nn)
  b <- rep_len(N, nn)
  approx.ans <- a  # True at lhs

  foo <- function(q, N, shape, p)
    pzipf(q, N, shape) - p

  dont.iterate <- p == 1 | shape <= 0
  lhs <- (p <= dzipf(1, N, shape)) | dont.iterate

  approx.ans[!lhs] <-
    bisection.basic(foo, a[!lhs], b[!lhs], shape = shape[!lhs],
                    tol = 1/16, p = p[!lhs], N = N[!lhs])
  faa <- floor(approx.ans)
  ans <- ifelse(pzipf(faa, N, shape) < p &
                p <= pzipf(faa+1, N, shape),
                faa+1, faa)

  ans[shape <= 0] <- NaN
  ans[p == 1] <- N

  ans
}  # qzipf



rzipf <- function(n, N, shape) {
  qzipf(runif(n), N, shape)
}



 zipf <-
  function(N = NULL, lshape = "loglink", ishape = NULL) {

  if (length(N) &&
    (!is.Numeric(N, positive = TRUE,
                 integer.valued = TRUE, length.arg = 1) ||
      N <= 1))
    stop("bad input for argument 'N'")
  enteredN <- length(N)

  if (length(ishape) && !is.Numeric(ishape, positive = TRUE))
      stop("argument 'ishape' must be > 0")

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


  new("vglmff",
  blurb = c("Zipf distribution f(y;s) = y^(-s) / sum((1:N)^(-s)),",
            " s > 0, y = 1, 2,...,N",
            ifelse(enteredN, paste(" = ", N, sep = ""), ""),
            "\n\n",
            "Link:    ",
            namesof("shape", lshape, earg = eshape),
            "\n\n",
            "Mean:    gharmonic(N, shape-1) / gharmonic(N, shape)"),
  infos = eval(substitute(function(...) {
    list(M1 = 1,
         Q1 = 1,
         dpqrfun = "zipf",
         multipleResponses = FALSE,
         parameters.names = "shape",
         N = enteredN,
         lshape = .lshape )
  }, list( .lshape = lshape,
           .enteredN = enteredN
         ))),

  rqresslot = eval(substitute(
    function(mu, y, w, eta, extra = NULL) {
    Shape <- eta2theta(eta, .lshape , earg = .eshape )
    scrambleseed <- runif(1)  # To scramble the seed
    qnorm(runif(length(y),
                pzipf(y - 1, Shape, N = extra$N),
                pzipf(y    , Shape, N = extra$N)))
  }, list( .lshape = lshape, .eshape = eshape ))),

  initialize = eval(substitute(expression({


    w.y.check(w = w, y = y,
              Is.integer.y = TRUE)


    predictors.names <- namesof("shape", .lshape ,
                                earg = .eshape , tag = FALSE)

    NN <- .N
    if (!is.Numeric(NN, length.arg = 1,
                    positive = TRUE, integer.valued = TRUE))
        NN <- max(y)
    if (max(y) > NN)
      stop("maximum of the response is greater than argument 'N'")
    if (any(y < 1))
      stop("all response values must be in 1, 2, ..., N( = ",
           NN, ")")
    extra$N <- NN

    if (!length(etastart)) {
      llfun <- function(shape, y, N, w) {
        sum(c(w) * dzipf(x = y, N = extra$N,
                         shape = shape, log = TRUE))
      }
      shape.init <- if (length( .ishape )) .ishape else
        getInitVals(gvals = seq(0.1, 3, length.out = 19),
                    llfun = llfun,
                    y = y, N = extra$N, w = w)
      shape.init <- rep_len(shape.init, length(y))
      if ( .lshape == "logloglink")
        shape.init[shape.init <= 1] <- 1.2
      etastart <- theta2eta(shape.init, .lshape , earg = .eshape )
    }
  }), list( .lshape = lshape, .eshape = eshape,
            .ishape = ishape, .N = N ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    gharmonic2(extra$N,
               shape = shape - 1) / gharmonic2(extra$N,
                                               shape = shape)
  }, list( .lshape = lshape, .eshape = eshape ))),
  last = eval(substitute(expression({
    misc$expected <- FALSE
    misc$link <-    c(shape = .lshape)
    misc$earg <- list(shape = .eshape )
    misc$N <- extra$N
  }), list( .lshape = lshape, .eshape = eshape ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * dzipf(x = y, N = extra$N, shape = shape,
                              log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .lshape = lshape, .eshape = eshape ))),
  vfamily = c("zipf"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    okay1 <- all(is.finite(shape)) && all(0 < shape)
    okay1
  }, list( .lshape = lshape, .eshape = eshape ))),




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

    pwts <- if (length(pwts <- object@prior.weights) > 0)
              pwts else weights(object, type = "prior")
    if (any(pwts != 1))
      warning("ignoring prior weights")
    eta <- predict(object)
    extra <- object@extra
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    rzipf(nsim * length(shape), N = extra$N, shape = shape)
  }, list( .lshape = lshape, .eshape = eshape ))),



  deriv = eval(substitute(expression({
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    fred1 <- gharmonic(extra$N, shape, deriv = 1)

    fred0 <- gharmonic2(extra$N, shape)

    dl.dshape <- -log(y) - fred1 / fred0
    dshape.deta <- dtheta.deta(shape, .lshape , earg = .eshape )
    d2shape.deta2 <- d2theta.deta2(shape, .lshape , .eshape )
    c(w) * dl.dshape * dshape.deta
  }), list( .lshape = lshape, .eshape = eshape ))),
  weight = expression({
    d2l.dshape <- gharmonic(extra$N, shape, deriv = 2) / fred0 -
                  (fred1/fred0)^2
    wz <- c(w) * (dshape.deta^2 * d2l.dshape -
                  d2shape.deta2 * dl.dshape)
    wz
  }))
}  # zipf






 ddiffzeta <- 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(shape), length(x), length(start))
  if (length(x)     != LLL) x     <- rep_len(x,     LLL)
  if (length(shape) != LLL) shape <- rep_len(shape, LLL)
  if (length(start) != LLL) start <- rep_len(start, LLL)

  ox <- !is.finite(x)
  zero <- ox | round(x) != x | x < start
  ans <- rep_len(if (log.arg) log(0) else 0, LLL)
  if (any(!zero)) {
    ans[!zero] <- (start[!zero] /      x[!zero]) ^(shape[!zero]) -
                  (start[!zero] / (1 + x[!zero]))^(shape[!zero])
    if (log.arg)
      ans[!zero] <- log(ans[!zero])
  }
  if (any(ox))
    ans[ox] <- if (log.arg) log(0) else 0
  ans[shape <= 0] <- NaN
  ans[start != round(start) | start < 1] <- NaN
  ans
}  # ddiffzeta



 pdiffzeta <-
  function(q, shape, start = 1, lower.tail = TRUE) {

  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)

  if (lower.tail) {
    ans <- 1 - (start / floor(1 + q))^shape
  } else {
    ans <-     (start / floor(1 + q))^shape
  }
  ans[q < start] <- if (lower.tail) 0 else 1
  ans[shape <= 0] <- NaN
  ans[start != round(start) | start < 1] <- NaN
  ans
}  # pdiffzeta




 qdiffzeta <- 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)

  lo <- rep_len(start, LLL)
  approx.ans <- lo  # True at lhs
  hi <- 2 * lo + 10.5
  dont.iterate <- p == 1 | shape <= 0 |
                  start != round(start) | start < 1
  done <- p <= pdiffzeta(hi, shape, start = start) | dont.iterate
  max.iter <- 100
  iter <- 0
  while (!all(done) && iter < max.iter) {
    lo[!done] <- hi[!done]
    hi[!done] <- 2 * hi[!done] + 10.5  # Bug fixed
    done[!done] <- is.infinite(hi[!done]) |
                   (p[!done] <= pdiffzeta(hi[!done], shape[!done],
                                          start[!done]))
    iter <- iter + 1
  }

  foo <- function(q, shape, start, p)
    pdiffzeta(q, shape, start) - p

   lhs <- (p <= ddiffzeta(start, shape, start = start)) |
           dont.iterate

  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)
  ans <- ifelse(pdiffzeta(faa  , shape, start = start) <  p &
                p <= pdiffzeta(faa+1, shape, start = start),
                faa+1, faa)

  ans[p == 1] <- Inf
  ans[shape <= 0] <- NaN
  ans[start != round(start) | start < 1] <- NaN

  ans
}  # qdiffzeta




rdiffzeta <- function(n, shape, start = 1) {
  rr <- runif(n)
  qdiffzeta(rr, shape, start = start)
}



 diffzeta <-
  function(start = 1, lshape = "loglink", ishape = NULL) {


  if (!is.Numeric(start, positive = TRUE,
                  integer.valued = TRUE, length.arg = 1))
    stop("bad input for argument 'start'")
  enteredstart <- length(start)

  if (length(ishape) && !is.Numeric(ishape, positive = TRUE))
    stop("argument 'ishape' must be > 0")

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


  new("vglmff",
  blurb = c("Difference in 2 Zeta distributions ",
        "f(y; shape) = y^(-shape) / sum((1:start)^(-shape)), ",
        "shape > 0; y = start, start+1,...",
        ifelse(enteredstart, paste0("start = ", start), ""),
        "\n\n",
        "Link:    ",
        namesof("shape", lshape, earg = eshape),
        "\n\n",
        "Mean:    gharmonic(start, shape-1) / ",
                 "gharmonic(start, shape)"),
  infos = eval(substitute(function(...) {
    list(M1 = 1,
         Q1 = 1,
         expected = TRUE,
         multipleResponses = TRUE,
         start = .start ,
         parameters.names = "shape")
  }, list( .start = start ))),

  rqresslot = eval(substitute(
    function(mu, y, w, eta, extra = NULL) {
    Shape <- eta2theta(eta, .lshape , earg = .eshape )
    scrambleseed <- runif(1)  # To scramble the seed
    qnorm(runif(length(y),
                pdiffzeta(y - 1, Shape, start = extra$start),
                pdiffzeta(y    , Shape, start = extra$start)))
  }, list( .lshape = lshape, .eshape = eshape ))),


  initialize = eval(substitute(expression({
  start <- .start
  temp5 <-
    w.y.check(w = w, y = y,
              ncol.w.max = Inf,
              ncol.y.max = Inf,
              Is.integer.y = TRUE,
              Is.positive.y = TRUE,
              out.wy = TRUE,
              colsyperw = 1,
              maximize = TRUE)
    w <- temp5$w
    y <- temp5$y
    if (any(y < start))
      stop("some response values less than 'start'")


    predictors.names <- namesof("shape", .lshape , earg = .eshape ,
                                tag = FALSE)

    extra$start <- start
    if (!length(etastart)) {
      llfun <- function(shape, y, start, w) {
        sum(c(w) * ddiffzeta(x = y, start = extra$start,
                             shape = shape, log = TRUE))
      }
      shape.init <- if (length( .ishape )) .ishape else
        getInitVals(gvals = seq(0.1, 3.0, length.out = 19),
                    llfun = llfun,
                    y = y, start = extra$start, w = w)
      shape.init <- rep_len(shape.init, length(y))
      if ( .lshape == "logloglink")
        shape.init[shape.init <= 1] <- 1.2
      etastart <- theta2eta(shape.init, .lshape , .eshape )
    }
  }),
  list( .lshape = lshape,
        .eshape = eshape, .ishape = ishape, .start = start ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    aa <- extra$start
    if (length(aa) != 1 || aa < 1 || round(aa) != aa)
      stop("the 'start' variable must be of unit length")
    if (aa == 1)
      return(zeta(shape))
    mymat <- matrix(1:aa, NROW(eta), aa, byrow = TRUE)
    temp1 <- rowSums(1 / mymat^shape)
    (aa^shape) * (zeta(shape) - temp1 + 1 / aa^(shape-1))
  }, list( .lshape = lshape, .eshape = eshape ))),
  last = eval(substitute(expression({
    misc$expected <- FALSE
    misc$link <-    c(shape = .lshape )
    misc$earg <- list(shape = .eshape )
    misc$start <- extra$start
  }), list( .lshape = lshape, .eshape = eshape ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * ddiffzeta(x = y, start = extra$start,
                                  shape = shape, log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .lshape = lshape, .eshape = eshape ))),
  vfamily = c("diffzeta"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    okay1 <- all(is.finite(shape)) && all(0 < shape)
    okay1
  }, list( .lshape = lshape, .eshape = eshape ))),






  deriv = eval(substitute(expression({
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    temp1 <- extra$start /  y
    temp2 <- extra$start / (y+1)
    AA <- temp1^shape - temp2^shape
    Aprime <- log(temp1) * temp1^shape -
              log(temp2) * temp2^shape

    dl.dshape <- Aprime / AA
    dshape.deta <- dtheta.deta(shape, .lshape , earg = .eshape )
    c(w) * dl.dshape * dshape.deta
  }), list( .lshape = lshape, .eshape = eshape ))),
  weight = expression({
    ned2l.dshape <- (Aprime / AA)^2  # Not quite FS. Half FS.
    wz <- c(w) * ned2l.dshape * dshape.deta^2
    wz
  }))
}  # diffzeta






ddiffzeta <- 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(shape), length(x), length(start))
  if (length(x)     != LLL) x     <- rep_len(x,     LLL)
  if (length(shape) != LLL) shape <- rep_len(shape, LLL)
  if (length(start) != LLL) start <- rep_len(start, LLL)

  ox <- !is.finite(x)
  zero <- ox | round(x) != x | x < start
  ans <- rep_len(if (log.arg) log(0) else 0, LLL)
  if (any(!zero)) {
    ans[!zero] <- (start[!zero] /      x[!zero]) ^(shape[!zero]) -
                  (start[!zero] / (1 + x[!zero]))^(shape[!zero])
    if (log.arg)
      ans[!zero] <- log(ans[!zero])
  }
  if (any(ox))
    ans[ox] <- if (log.arg) log(0) else 0
  ans[shape <= 0] <- NaN
  ans[start != round(start) | start < 1] <- NaN
  ans
}



 pdiffzeta <- function(q, shape, start = 1, lower.tail = TRUE) {

  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)

  if (lower.tail) {
    ans <- 1 - (start / floor(1 + q))^shape
  } else {
    ans <-     (start / floor(1 + q))^shape
  }
  ans[q < start] <- if (lower.tail) 0 else 1
  ans[shape <= 0] <- NaN
  ans[start != round(start) | start < 1] <- NaN
  ans
}  # pdiffzeta




 qdiffzeta <- 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)

  lo <- rep_len(start, LLL)
  approx.ans <- lo  # True at lhs
  hi <- 2 * lo + 10.5
  dont.iterate <- p == 1 | shape <= 0 |
                  start != round(start) | start < 1
  done <- p <= pdiffzeta(hi, shape, start = start) | dont.iterate
  max.iter <- 100
  iter <- 0
  while (!all(done) && iter < max.iter) {
    lo[!done] <- hi[!done]
    hi[!done] <- 2 * hi[!done] + 10.5  # Bug fixed
    done[!done] <- is.infinite(hi[!done]) |
                   (p[!done] <= pdiffzeta(hi[!done], shape[!done],
                                          start[!done]))
    iter <- iter + 1
  }

  foo <- function(q, shape, start, p)
    pdiffzeta(q, shape, start) - p

   lhs <- (p <= ddiffzeta(start, shape, start = start)) |
           dont.iterate

  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)
  ans <- ifelse(pdiffzeta(faa  , shape, start = start) <  p &
                p <= pdiffzeta(faa+1, shape, start = start),
                faa+1, faa)

  ans[p == 1] <- Inf
  ans[shape <= 0] <- NaN
  ans[start != round(start) | start < 1] <- NaN

  ans
}  # qdiffzeta




rdiffzeta <- function(n, shape, start = 1) {
  rr <- runif(n)
  qdiffzeta(rr, shape, start = start)
}



 diffzeta <-
  function(start = 1, lshape = "loglink", ishape = NULL) {


  if (!is.Numeric(start, positive = TRUE,
                  integer.valued = TRUE, length.arg = 1))
    stop("bad input for argument 'start'")
  enteredstart <- length(start)

  if (length(ishape) && !is.Numeric(ishape, positive = TRUE))
    stop("argument 'ishape' must be > 0")

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


  new("vglmff",
  blurb = c("Difference in 2 Zeta distributions ",
        "f(y; shape) = y^(-shape) / sum((1:start)^(-shape)), ",
        "shape > 0, start, start+1,...",
        ifelse(enteredstart, paste0("start = ", start), ""),
        "\n\n",
        "Link:    ",
        namesof("shape", lshape, earg = eshape),
        "\n\n",
        "Mean:    gharmonic(start, shape-1) / ",
                 "gharmonic(start, shape)"),
  infos = eval(substitute(function(...) {
    list(M1 = 1,
         Q1 = 1,
         dpqrfun = "diffzipf",
         expected = TRUE,
         multipleResponses = TRUE,
         start = .start ,
         parameters.names = "shape")
  }, list( .start = start ))),

  rqresslot = eval(substitute(
    function(mu, y, w, eta, extra = NULL) {
    Shape <- eta2theta(eta, .lshape , earg = .eshape )
    scrambleseed <- runif(1)  # To scramble the seed
    qnorm(runif(length(y),
                pdiffzeta(y - 1, Shape, start = extra$start),
                pdiffzeta(y    , Shape, start = extra$start)))
  }, list( .lshape = lshape, .eshape = eshape ))),

  initialize = eval(substitute(expression({
  start <- .start
  temp5 <-
    w.y.check(w = w, y = y,
              ncol.w.max = Inf,
              ncol.y.max = Inf,
              Is.integer.y = TRUE,
              Is.positive.y = TRUE,
              out.wy = TRUE,
              colsyperw = 1,
              maximize = TRUE)
    w <- temp5$w
    y <- temp5$y
    if (any(y < start))
      stop("some response values less than 'start'")


    predictors.names <- namesof("shape", .lshape , earg = .eshape ,
                                tag = FALSE)

    extra$start <- start
    if (!length(etastart)) {
      llfun <- function(shape, y, start, w) {
        sum(c(w) * ddiffzeta(x = y, start = extra$start,
                             shape = shape, log = TRUE))
      }
      shape.init <- if (length( .ishape )) .ishape else
        getInitVals(gvals = seq(0.1, 3.0, length.out = 19),
                    llfun = llfun,
                    y = y, start = extra$start, w = w)
      shape.init <- rep_len(shape.init, length(y))
      if ( .lshape == "logloglink")
        shape.init[shape.init <= 1] <- 1.2
      etastart <- theta2eta(shape.init, .lshape , earg = .eshape )
    }
  }), list( .lshape = lshape,
            .eshape = eshape, .ishape = ishape, .start = start ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    aa <- extra$start
    if (length(aa) != 1 || aa < 1 || round(aa) != aa)
      stop("the 'start' variable must be of unit length")
    if (aa == 1)
      return(zeta(shape))
    mymat <- matrix(1:aa, NROW(eta), aa, byrow = TRUE)
    temp1 <- rowSums(1 / mymat^shape)
    (aa^shape) * (zeta(shape) - temp1 + 1 / aa^(shape-1))
  }, list( .lshape = lshape, .eshape = eshape ))),
  last = eval(substitute(expression({
    misc$expected <- FALSE
    misc$link <-    c(shape = .lshape )
    misc$earg <- list(shape = .eshape )
    misc$start <- extra$start
  }), list( .lshape = lshape, .eshape = eshape ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * ddiffzeta(x = y, start = extra$start,
                                  shape = shape, log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .lshape = lshape, .eshape = eshape ))),
  vfamily = c("diffzeta"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    okay1 <- all(is.finite(shape)) && all(0 < shape)
    okay1
  }, list( .lshape = lshape, .eshape = eshape ))),






  deriv = eval(substitute(expression({
    shape <- eta2theta(eta, .lshape , earg = .eshape )
    temp1 <- extra$start /  y
    temp2 <- extra$start / (y+1)
    AA <- temp1^shape - temp2^shape
    Aprime <- log(temp1) * temp1^shape -
              log(temp2) * temp2^shape

    dl.dshape <- Aprime / AA
    dshape.deta <- dtheta.deta(shape, .lshape , earg = .eshape )
    c(w) * dl.dshape * dshape.deta
  }), list( .lshape = lshape, .eshape = eshape ))),
  weight = expression({
    ned2l.dshape <- (Aprime / AA)^2  # Not quite FS. Half FS.
    wz <- c(w) * ned2l.dshape * dshape.deta^2
    wz
  }))
}  # diffzeta

Try the VGAM package in your browser

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

VGAM documentation built on Sept. 19, 2023, 9:06 a.m.