R/family.bivariate.R

Defines functions dbistudenttcop kendall.tau kendall.tau gumbelI binormal rbinorm dbinorm biamhcop biamhcop.control rbiamhcop pbiamhcop dbiamhcop biplackettcop biplackettcop.control dbiplackcop rbiplackcop pbiplackcop bigumbelIexp bifgmcop pbifgmcop dbifgmcop rbifgmcop bifgmexp gammahyperbola bifrankcop bifrankcop.control dbifrankcop dbifrank pbifrankcop rbifrankcop bigamma.mckay freund61 rbilogis pbilogis dbilogis bilogistic bilogistic.control binormalcop rbinormcop pbinormcop dbinormcop bistudentt bistudent.deriv.dof dbistudentt biclaytoncop rbiclaytoncop dbiclaytoncop trinormal trinormal.control rtrinorm dtrinorm

Documented in biamhcop biclaytoncop biclaytoncop bifgmcop bifgmexp bifrankcop bigamma.mckay bigumbelIexp bilogistic binormal binormalcop binormalcop biplackettcop bistudentt dbiamhcop dbiclaytoncop dbiclaytoncop dbifgmcop dbifrankcop dbilogis dbinorm dbinorm dbinormcop dbinormcop dbiplackcop dbistudentt dtrinorm freund61 gammahyperbola kendall.tau kendall.tau pbiamhcop pbifgmcop pbifrankcop pbilogis pbinormcop pbinormcop pbiplackcop rbiamhcop rbiclaytoncop rbiclaytoncop rbifgmcop rbifrankcop rbilogis rbinorm rbinormcop rbinormcop rbiplackcop rtrinorm trinormal

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












 dtrinorm <-
  function(x1, x2, x3, mean1 = 0, mean2 = 0, mean3 = 0,
           var1 = 1, var2 = 1, var3 = 1,
           cov12 = 0, cov23 = 0, cov13 = 0,
           log = FALSE) {

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

  M <- 3
  n <- max(length(x1), length(x2), length(x3),
           length(mean1), length(mean2), length(mean3),
           length(var1 ), length(var2 ), length(var3 ),
           length(cov12), length(cov23), length(cov13))

  sd1 <- sqrt(var1)
  sd2 <- sqrt(var2)
  sd3 <- sqrt(var3)
  rho12 <- cov12 / (sd1 * sd2)
  rho13 <- cov13 / (sd1 * sd3)
  rho23 <- cov23 / (sd2 * sd3)

  bbb <- 1 - rho12^2 - rho13^2 - rho23^2 + 2 * rho12 * rho13 * rho23
  logdet <- 2 * (log(sd1) + log(sd2) + log(sd3)) + log(bbb)

  Sigmainv <- matrix(0, n, dimm(3))  # sum(3:1)
  Sigmainv[, iam(1, 1, M = M)] <- (1 - rho23^2) / (bbb * sd1^2)      
  Sigmainv[, iam(2, 2, M = M)] <- (1 - rho13^2) / (bbb * sd2^2)      
  Sigmainv[, iam(3, 3, M = M)] <- (1 - rho12^2) / (bbb * sd3^2)      
  Sigmainv[, iam(1, 2, M = M)] <- (rho13 * rho23 - rho12) / (
                                   sd1 * sd2 * bbb)
  Sigmainv[, iam(2, 3, M = M)] <- (rho12 * rho13 - rho23) / (
                                   sd2 * sd3 * bbb)
  Sigmainv[, iam(1, 3, M = M)] <- (rho12 * rho23 - rho13) / (
                                   sd1 * sd3 * bbb)

  ymatt <- rbind(x1 - mean1, x2 - mean2, x3 - mean3)
  dim(ymatt) <- c(nrow(ymatt), 1, ncol(ymatt))  # For mux5()
  qform <- mux5(x = ymatt, cc = Sigmainv, M = 3, matrix.arg = TRUE)

  logpdf <- -1.5 * log(2 * pi) - 0.5 * logdet - 0.5 * c(qform)
  logpdf[is.infinite(x1) | is.infinite(x2) | is.infinite(x3)] <- log(0)

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



rtrinorm <- function(n, mean1 = 0, mean2 = 0, mean3 = 0,
                     var1 = 1, var2 = 1, var3 = 1,
                     cov12 = 0, cov23 = 0, cov13 = 0) {

  Y1 <- rnorm(n, mean1, sqrt(var1))
  ans2 <- rbinorm(n,
                  mean1 = mean2 + cov12 * (Y1 - mean1) / var1,
                  mean2 = mean3 + cov13 * (Y1 - mean1) / var1,
                  var1  = var2  - cov12 * cov12 / var1,
                  var2  = var3  - cov13 * cov13 / var1,
                  cov12 = cov23 - cov12 * cov13 / var1)

  ans <- cbind(Y1, ans2)
  colnames(ans) <- paste("X", 1:3, sep = "")
  ans
}





trinormal.control <-
  function(summary.HDEtest = FALSE,  # Overwrites the summary() default.
           ...) {
  list(summary.HDEtest = summary.HDEtest)
}



 trinormal <-
   function(
           zero = c("sd", "rho"),
           eq.mean = FALSE,
           eq.sd = FALSE,
           eq.cor = FALSE,
           lmean1 = "identitylink",
           lmean2 = "identitylink",
           lmean3 = "identitylink",
           lsd1   = "loglink",
           lsd2   = "loglink",
           lsd3   = "loglink",
           lrho12 = "rhobitlink",
           lrho23 = "rhobitlink",
           lrho13 = "rhobitlink",
           imean1 = NULL,       imean2 = NULL,       imean3 = NULL,
           isd1   = NULL,       isd2   = NULL,       isd3   = NULL,
           irho12 = NULL,       irho23 = NULL,       irho13 = NULL,
           imethod = 1) {


  lmean1 <- as.list(substitute(lmean1))
  emean1 <- link2list(lmean1)
  lmean1 <- attr(emean1, "function.name")

  lmean2 <- as.list(substitute(lmean2))
  emean2 <- link2list(lmean2)
  lmean2 <- attr(emean2, "function.name")

  lmean3 <- as.list(substitute(lmean3))
  emean3 <- link2list(lmean3)
  lmean3 <- attr(emean3, "function.name")

  lsd1 <- as.list(substitute(lsd1))
  esd1 <- link2list(lsd1)
  lsd1 <- attr(esd1, "function.name")

  lsd2 <- as.list(substitute(lsd2))
  esd2 <- link2list(lsd2)
  lsd2 <- attr(esd2, "function.name")

  lsd3 <- as.list(substitute(lsd3))
  esd3 <- link2list(lsd3)
  lsd3 <- attr(esd3, "function.name")

  lrho12 <- as.list(substitute(lrho12))
  erho12 <- link2list(lrho12)
  lrho12 <- attr(erho12, "function.name")

  lrho23 <- as.list(substitute(lrho23))
  erho23 <- link2list(lrho23)
  lrho23 <- attr(erho23, "function.name")

  lrho13 <- as.list(substitute(lrho13))
  erho13 <- link2list(lrho13)
  lrho13 <- attr(erho13, "function.name")


  if (!is.logical(eq.mean) || length(eq.mean) != 1)
    stop("argument 'eq.mean' must be a single logical")
  if (!is.logical(eq.sd) || length(eq.sd) != 1)
    stop("argument 'eq.sd' must be a single logical")
  if (!is.logical(eq.cor) || length(eq.cor) != 1)
    stop("argument 'eq.cor' must be a single logical")


  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("Trivariate normal distribution\n",
            "Links:    ",
            namesof("mean1", lmean1, earg = emean1 ), ", ",
            namesof("mean2", lmean2, earg = emean2 ), ", ",
            namesof("mean3", lmean3, earg = emean3 ), ", ",
            namesof("sd1",   lsd1,   earg = esd1   ), ", ",
            namesof("sd2",   lsd2,   earg = esd2   ), ", ",
            namesof("sd3",   lsd3,   earg = esd3   ), ",\n",
            "          ",
            namesof("rho12", lrho12, earg = erho12 ), ", ",
            namesof("rho23", lrho23, earg = erho23 ), ", ",
            namesof("rho13", lrho13, earg = erho13 )),
  constraints = eval(substitute(expression({

    constraints.orig <- constraints
    M1 <- 9
    NOS <- M / M1

    cm1.m <-
    cmk.m <- kronecker(diag(NOS), rbind(diag(3), matrix(0, 6, 3)))
    con.m <- cm.VGAM(kronecker(diag(NOS), eijfun(1:3, 9)),
                     x = x,
                     bool = .eq.mean ,  #
                     constraints = constraints.orig,
                     apply.int = TRUE,
                     cm.default           = cmk.m,
                     cm.intercept.default = cm1.m)


    cm1.s <-
    cmk.s <- kronecker(diag(NOS),
                       rbind(matrix(0, 3, 3), diag(3), matrix(0, 3, 3)))
    con.s <- cm.VGAM(kronecker(diag(NOS), eijfun(4:6, 9)),
                     x = x,
                     bool = .eq.sd ,  #
                     constraints = constraints.orig,
                     apply.int = TRUE,
                     cm.default           = cmk.s,
                     cm.intercept.default = cm1.s)



    cm1.r <-
    cmk.r <- kronecker(diag(NOS),
                       rbind(matrix(0, 3, 3), matrix(0, 3, 3), diag(3)))
    con.r <- cm.VGAM(kronecker(diag(NOS), eijfun(7:9, 9)),
                     x = x,
                     bool = .eq.cor ,  #
                     constraints = constraints.orig,
                     apply.int = TRUE,
                     cm.default           = cmk.r,
                     cm.intercept.default = cm1.r)



    con.use <- con.m
    for (klocal in seq_along(con.m)) {
      con.use[[klocal]] <-
        cbind(con.m[[klocal]],
              con.s[[klocal]],
              con.r[[klocal]])
    }




    constraints <- con.use
    constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
                                predictors.names = predictors.names,
                                M1 = M1)
  }), list( .zero = zero,
            .eq.sd   = eq.sd,
            .eq.mean = eq.mean,
            .eq.cor  = eq.cor  ))),

  infos = eval(substitute(function(...) {
    list(M1 = 9,
         Q1 = 3,
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = c("mean1", "mean2", "mean3",
                              "sd1",   "sd2",   "sd3",
                              "rho12", "rho13", "rho23"),
         eq.cor  = .eq.cor ,
         eq.mean = .eq.mean ,
         eq.sd   = .eq.sd   ,
         zero    = .zero )
    }, list( .zero    = zero,
             .eq.cor  = eq.cor,
             .eq.mean = eq.mean,
             .eq.sd   = eq.sd    ))),

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



    predictors.names <- c(
      namesof("mean1", .lmean1 , earg = .emean1 , short = TRUE),
      namesof("mean2", .lmean2 , earg = .emean2 , short = TRUE),
      namesof("mean3", .lmean3 , earg = .emean3 , short = TRUE),
      namesof("sd1",   .lsd1 ,   earg = .esd1 ,   short = TRUE),
      namesof("sd2",   .lsd2 ,   earg = .esd2 ,   short = TRUE),
      namesof("sd3",   .lsd3 ,   earg = .esd3 ,   short = TRUE),
      namesof("rho12", .lrho12 , earg = .erho12 , short = TRUE),
      namesof("rho23", .lrho23 , earg = .erho23 , short = TRUE),
      namesof("rho13", .lrho13 , earg = .erho13 , short = TRUE))

    extra$colnames.y  <- colnames(y)

    if (!length(etastart)) {
      imean1 <- rep_len(if (length( .imean1 )) .imean1 else
                   weighted.mean(y[, 1], w = w), n)
      imean2 <- rep_len(if (length( .imean2 )) .imean2 else
                   weighted.mean(y[, 2], w = w), n)
      imean3 <- rep_len(if (length( .imean3 )) .imean3 else
                   weighted.mean(y[, 3], w = w), n)
      isd1 <- rep_len(if (length( .isd1 )) .isd1 else  sd(y[, 1]), n)
      isd2 <- rep_len(if (length( .isd2 )) .isd2 else  sd(y[, 2]), n)
      isd3 <- rep_len(if (length( .isd3 )) .isd3 else  sd(y[, 3]), n)
      irho12 <- rep_len(if (length( .irho12 )) .irho12 else
                        cor(y[, 1], y[, 2]), n)
      irho23 <- rep_len(if (length( .irho23 )) .irho23 else
                        cor(y[, 2], y[, 3]), n)
      irho13 <- rep_len(if (length( .irho13 )) .irho13 else
                        cor(y[, 1], y[, 3]), n)

      if ( .imethod == 2) {
        imean1 <- abs(imean1) + 0.01
        imean2 <- abs(imean2) + 0.01
        imean3 <- abs(imean3) + 0.01
      }
      etastart <-
        cbind(theta2eta(imean1, .lmean1 , earg = .emean1 ),
              theta2eta(imean2, .lmean2 , earg = .emean2 ),
              theta2eta(imean3, .lmean3 , earg = .emean3 ),
              theta2eta(isd1,   .lsd1 ,   earg = .esd1 ),
              theta2eta(isd2,   .lsd2 ,   earg = .esd2 ),
              theta2eta(isd3,   .lsd3 ,   earg = .esd3 ),
              theta2eta(irho12, .lrho12 , earg = .erho12 ),
              theta2eta(irho23, .lrho23 , earg = .erho23 ),
              theta2eta(irho13, .lrho13 , earg = .erho13 ))
    }
  }), list( .lmean1 = lmean1, .lmean2 = lmean2, .lmean3 = lmean3,
            .emean1 = emean1, .emean2 = emean2, .emean3 = emean3,
            .imean1 = imean1, .imean2 = imean2, .imean3 = imean3,
            .lsd1   = lsd1  , .lsd2   = lsd2  , .lsd3   = lsd3  ,
            .esd1   = esd1  , .esd2   = esd2  , .esd3   = esd3  ,
            .isd1   = isd1,   .isd2   = isd2,   .isd3   = isd3,
            .lrho12 = lrho12, .lrho13 = lrho13, .lrho23 = lrho23,
            .erho12 = erho12, .erho13 = erho13, .erho23 = erho23,
            .irho12 = irho12, .irho13 = irho13, .irho23 = irho23,
            .imethod = imethod ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    NOS <- ncol(eta) / c(M1 = 9)
    mean1 <- eta2theta(eta[, 1], .lmean1 , earg = .emean1 )
    mean2 <- eta2theta(eta[, 2], .lmean2 , earg = .emean2 )
    mean3 <- eta2theta(eta[, 3], .lmean3 , earg = .emean3 )
    fv.mat <- cbind(mean1, mean2, mean3)
    label.cols.y(fv.mat, colnames.y = extra$colnames.y, NOS = NOS)
  } , list( .lmean1 = lmean1, .lmean2 = lmean2, .lmean3 = lmean3,
            .emean1 = emean1, .emean2 = emean2, .emean3 = emean3,
            .lsd1   = lsd1  , .lsd2   = lsd2  , .lsd3   = lsd3  ,
            .esd1   = esd1  , .esd2   = esd2  , .esd3   = esd3  ,
            .erho12 = erho12, .erho13 = erho13, .erho23 = erho23 ))),

  last = eval(substitute(expression({
    misc$link <-    c("mean1" = .lmean1 ,
                      "mean2" = .lmean2 ,
                      "mean3" = .lmean3 ,
                      "sd1"   = .lsd1 ,
                      "sd2"   = .lsd2 ,
                      "sd3"   = .lsd3 ,
                      "rho12" = .lrho12 ,
                      "rho23" = .lrho23 ,
                      "rho13" = .lrho13 )

    misc$earg <- list("mean1" = .emean1 ,
                      "mean2" = .emean2 ,
                      "mean3" = .emean3 ,
                      "sd1"   = .esd1 ,
                      "sd2"   = .esd2 ,
                      "sd3"   = .esd3 ,
                      "rho12" = .erho12 ,
                      "rho23" = .erho23 ,
                      "rho13" = .erho13 )
  }), list( .lmean1 = lmean1, .lmean2 = lmean2, .lmean3 = lmean3,
            .emean1 = emean1, .emean2 = emean2, .emean3 = emean3,
            .lsd1   = lsd1  , .lsd2   = lsd2  , .lsd3   = lsd3  ,
            .esd1   = esd1  , .esd2   = esd2  , .esd3   = esd3  ,
            .lrho12 = lrho12, .lrho13 = lrho13, .lrho23 = lrho23,
            .erho12 = erho12, .erho13 = erho13, .erho23 = erho23 ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    mean1 <- eta2theta(eta[, 1], .lmean1 , earg = .emean1 )
    mean2 <- eta2theta(eta[, 2], .lmean2 , earg = .emean2 )
    mean3 <- eta2theta(eta[, 3], .lmean3 , earg = .emean3 )
    sd1   <- eta2theta(eta[, 4], .lsd1   , earg = .esd1   )
    sd2   <- eta2theta(eta[, 5], .lsd2   , earg = .esd2   )
    sd3   <- eta2theta(eta[, 6], .lsd3   , earg = .esd3   )
    Rho12 <- eta2theta(eta[, 7], .lrho12 , earg = .erho12 )
    Rho23 <- eta2theta(eta[, 8], .lrho23 , earg = .erho23 )
    Rho13 <- eta2theta(eta[, 9], .lrho13 , earg = .erho13 )

    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <-
        c(w) * dtrinorm(x1 = y[, 1], x2 = y[, 2], x3 = y[, 3],
                       mean1 = mean1, mean2 = mean2, mean3 = mean3,
                       var1 = sd1^2, var2 = sd2^2, var3 = sd3^2, 
                       cov12 = Rho12 * sd1 * sd2,
                       cov23 = Rho23 * sd2 * sd3,
                       cov13 = Rho13 * sd1 * sd3,
                       log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  } , list( .lmean1 = lmean1, .lmean2 = lmean2, .lmean3 = lmean3,
            .emean1 = emean1, .emean2 = emean2, .emean3 = emean3,
            .lsd1   = lsd1  , .lsd2   = lsd2  , .lsd3   = lsd3  ,
            .esd1   = esd1  , .esd2   = esd2  , .esd3   = esd3  ,
            .lrho12 = lrho12, .lrho13 = lrho13, .lrho23 = lrho23,
            .erho12 = erho12, .erho13 = erho13, .erho23 = erho23 ))),
  vfamily = c("trinormal"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    mean1 <- eta2theta(eta[, 1], .lmean1 , earg = .emean1 )
    mean2 <- eta2theta(eta[, 2], .lmean2 , earg = .emean2 )
    mean3 <- eta2theta(eta[, 3], .lmean3 , earg = .emean3 )
    sd1   <- eta2theta(eta[, 4], .lsd1   , earg = .esd1   )
    sd2   <- eta2theta(eta[, 5], .lsd2   , earg = .esd2   )
    sd3   <- eta2theta(eta[, 6], .lsd3   , earg = .esd3   )
    Rho12 <- eta2theta(eta[, 7], .lrho12 , earg = .erho12 )
    Rho23 <- eta2theta(eta[, 8], .lrho23 , earg = .erho23 )
    Rho13 <- eta2theta(eta[, 9], .lrho13 , earg = .erho13 )
    okay1 <- all(is.finite(mean1)) &&
             all(is.finite(mean2)) &&
             all(is.finite(mean3)) &&
             all(is.finite(sd1  )) && all(0 < sd1)        &&
             all(is.finite(sd2  )) && all(0 < sd2)        &&
             all(is.finite(sd3  )) && all(0 < sd3)        &&
             all(is.finite(Rho12)) && all(abs(Rho12) < 1) &&
             all(is.finite(Rho23)) && all(abs(Rho23) < 1) &&
             all(is.finite(Rho13)) && all(abs(Rho13) < 1) &&
             all(0 < 1 - Rho12^2 - Rho13^2 - Rho23^2 +
                     2 * Rho12 * Rho13 * Rho23)
    okay1
  } , list( .lmean1 = lmean1, .lmean2 = lmean2, .lmean3 = lmean3,
            .emean1 = emean1, .emean2 = emean2, .emean3 = emean3,
            .lsd1   = lsd1  , .lsd2   = lsd2  , .lsd3   = lsd3  ,
            .esd1   = esd1  , .esd2   = esd2  , .esd3   = esd3  ,
            .lrho12 = lrho12, .lrho13 = lrho13, .lrho23 = lrho23,
            .erho12 = erho12, .erho13 = erho13, .erho23 = erho23 ))),



  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)
    mean1 <- eta2theta(eta[, 1], .lmean1 , earg = .emean1 )
    mean2 <- eta2theta(eta[, 2], .lmean2 , earg = .emean2 )
    mean3 <- eta2theta(eta[, 3], .lmean3 , earg = .emean3 )
    sd1   <- eta2theta(eta[, 4], .lsd1   , earg = .esd1   )
    sd2   <- eta2theta(eta[, 5], .lsd2   , earg = .esd2   )
    sd3   <- eta2theta(eta[, 6], .lsd3   , earg = .esd3   )
    Rho12 <- eta2theta(eta[, 7], .lrho12 , earg = .erho12 )
    Rho23 <- eta2theta(eta[, 8], .lrho23 , earg = .erho23 )
    Rho13 <- eta2theta(eta[, 9], .lrho13 , earg = .erho13 )
    rtrinorm(nsim * length(sd1),
             mean1 = mean1, mean2 = mean2, mean3 = mean3,
             var1 = sd1^2, var2 = sd2^2, var3 = sd3^2,
             cov12 = Rho12 * sd1 * sd2,
             cov23 = Rho23 * sd2 * sd3,
             cov13 = Rho13 * sd1 * sd3)
  } , list( .lmean1 = lmean1, .lmean2 = lmean2, .lmean3 = lmean3,
            .emean1 = emean1, .emean2 = emean2, .emean3 = emean3,
            .lsd1   = lsd1  , .lsd2   = lsd2  , .lsd3   = lsd3  ,
            .esd1   = esd1  , .esd2   = esd2  , .esd3   = esd3  ,
            .lrho12 = lrho12, .lrho13 = lrho13, .lrho23 = lrho23,
            .erho12 = erho12, .erho13 = erho13, .erho23 = erho23 ))),




  deriv = eval(substitute(expression({
    mean1 <- eta2theta(eta[, 1], .lmean1 , earg = .emean1 )
    mean2 <- eta2theta(eta[, 2], .lmean2 , earg = .emean2 )
    mean3 <- eta2theta(eta[, 3], .lmean3 , earg = .emean3 )
    sd1   <- eta2theta(eta[, 4], .lsd1   , earg = .esd1   )
    sd2   <- eta2theta(eta[, 5], .lsd2   , earg = .esd2   )
    sd3   <- eta2theta(eta[, 6], .lsd3   , earg = .esd3   )
    rho12 <- eta2theta(eta[, 7], .lrho12 , earg = .erho12 )
    rho23 <- eta2theta(eta[, 8], .lrho23 , earg = .erho23 )
    rho13 <- eta2theta(eta[, 9], .lrho13 , earg = .erho13 )
    bbb <- 1 - rho12^2 - rho13^2 - rho23^2 + 2 * rho12 * rho13 * rho23



  Sigmainv <- matrix(0, n, dimm(3))  # sum(3:1)
  Sigmainv[, iam(1, 1, M = 3)] <- (1 - rho23^2) / (bbb * sd1^2)      
  Sigmainv[, iam(2, 2, M = 3)] <- (1 - rho13^2) / (bbb * sd2^2)      
  Sigmainv[, iam(3, 3, M = 3)] <- (1 - rho12^2) / (bbb * sd3^2)      
  Sigmainv[, iam(1, 2, M = 3)] <- (rho13 * rho23 - rho12) / (
                                   sd1 * sd2 * bbb)
  Sigmainv[, iam(2, 3, M = 3)] <- (rho12 * rho13 - rho23) / (
                                   sd2 * sd3 * bbb)
  Sigmainv[, iam(1, 3, M = 3)] <- (rho12 * rho23 - rho13) / (
                                   sd1 * sd3 * bbb)

 
    dem <- bbb * (sd1 * sd2 * sd3)^2
    ymat.cen <- y - cbind(mean1, mean2, mean3)  # Usual dimensions n x 3
    ymatt.cen <- t(ymat.cen)
    dim(ymatt.cen) <- c(nrow(ymatt.cen), 1, ncol(ymatt.cen))  # 4 mux5()
    dl.dmeans <- mux22(t(Sigmainv), ymat.cen, M = 3, as.matrix = TRUE)



  SI.sd1 <- Sigmainv * 0
  SI.sd1[, iam(1, 1, M = 3)] <- -2 * Sigmainv[, iam(1, 1, M = 3)] / sd1
  SI.sd1[, iam(2, 2, M = 3)] <- 0
  SI.sd1[, iam(3, 3, M = 3)] <- 0
  SI.sd1[, iam(1, 2, M = 3)] <- -1 * Sigmainv[, iam(1, 2, M = 3)] / sd1
  SI.sd1[, iam(2, 3, M = 3)] <- 0
  SI.sd1[, iam(1, 3, M = 3)] <- -1 * Sigmainv[, iam(1, 3, M = 3)] / sd1

  SI.sd2 <- Sigmainv * 0
  SI.sd2[, iam(2, 2, M = 3)] <- -2 * Sigmainv[, iam(2, 2, M = 3)] / sd2
  SI.sd2[, iam(1, 1, M = 3)] <- 0
  SI.sd2[, iam(3, 3, M = 3)] <- 0
  SI.sd2[, iam(1, 2, M = 3)] <- -1 * Sigmainv[, iam(1, 2, M = 3)] / sd2
  SI.sd2[, iam(1, 3, M = 3)] <- 0
  SI.sd2[, iam(2, 3, M = 3)] <- -1 * Sigmainv[, iam(2, 3, M = 3)] / sd2

  SI.sd3 <- Sigmainv * 0
  SI.sd3[, iam(3, 3, M = 3)] <- -2 * Sigmainv[, iam(3, 3, M = 3)] / sd3
  SI.sd3[, iam(2, 2, M = 3)] <- 0
  SI.sd3[, iam(1, 1, M = 3)] <- 0
  SI.sd3[, iam(1, 3, M = 3)] <- -1 * Sigmainv[, iam(1, 3, M = 3)] / sd3
  SI.sd3[, iam(1, 2, M = 3)] <- 0
  SI.sd3[, iam(2, 3, M = 3)] <- -1 * Sigmainv[, iam(2, 3, M = 3)] / sd3


    dl.dsd1   <- -1 / sd1 - 0.5 *
      c(mux5(x = ymatt.cen, cc = SI.sd1, M = 3, matrix.arg = TRUE))
    dl.dsd2   <- -1 / sd2 - 0.5 *
      c(mux5(x = ymatt.cen, cc = SI.sd2, M = 3, matrix.arg = TRUE))
    dl.dsd3   <- -1 / sd3 - 0.5 *
      c(mux5(x = ymatt.cen, cc = SI.sd3, M = 3, matrix.arg = TRUE))


  dbbb.drho12 <- 2 * (rho13 * rho23 - rho12)
  dbbb.drho23 <- 2 * (rho12 * rho13 - rho23)
  dbbb.drho13 <- 2 * (rho12 * rho23 - rho13)
  SI.rho12 <- Sigmainv * 0
  SI.rho12[, iam(1, 1, M = 3)] <-
    -1 * Sigmainv[, iam(1, 1, M = 3)] * dbbb.drho12 / bbb
  SI.rho12[, iam(2, 2, M = 3)] <-
    -1 * Sigmainv[, iam(2, 2, M = 3)] * dbbb.drho12 / bbb
  SI.rho12[, iam(3, 3, M = 3)] <-
    (-2 * rho12 - (1 - rho12^2) * dbbb.drho12 / bbb) / (bbb * sd3^2)
  SI.rho12[, iam(1, 2, M = 3)] <-
    (-1 - (rho13 * rho23 - rho12) * dbbb.drho12 / bbb) / (
     bbb * sd1 * sd2)
  SI.rho12[, iam(2, 3, M = 3)] <-
    (rho13 - (rho12 * rho13 - rho23) * dbbb.drho12 / bbb) / (
     bbb * sd2 * sd3)
  SI.rho12[, iam(1, 3, M = 3)] <-
    (rho23 - (rho12 * rho23 - rho13) * dbbb.drho12 / bbb) / (
     bbb * sd1 * sd3)


  SI.rho23 <- Sigmainv * 0
  SI.rho23[, iam(1, 1, M = 3)] <-
    (-2 * rho23 - (1 - rho23^2) * dbbb.drho23 / bbb) / (bbb * sd1^2)
  SI.rho23[, iam(2, 2, M = 3)] <-
    -1 * Sigmainv[, iam(2, 2, M = 3)] * dbbb.drho23 / bbb
  SI.rho23[, iam(3, 3, M = 3)] <-
    -1 * Sigmainv[, iam(3, 3, M = 3)] * dbbb.drho23 / bbb
  SI.rho23[, iam(1, 2, M = 3)] <-
    (rho13 - (rho13 * rho23 - rho12) * dbbb.drho23 / bbb) / (
     bbb * sd1 * sd2)
  SI.rho23[, iam(2, 3, M = 3)] <-
    (-1 - (rho12 * rho13 - rho23) * dbbb.drho23 / bbb) / (
     bbb * sd2 * sd3)
  SI.rho23[, iam(1, 3, M = 3)] <-
    (rho12 - (rho12 * rho23 - rho13) * dbbb.drho23 / bbb) / (
     bbb * sd1 * sd3)

  SI.rho13 <- Sigmainv * 0
  SI.rho13[, iam(1, 1, M = 3)] <-
    -1 * Sigmainv[, iam(1, 1, M = 3)] * dbbb.drho13 / bbb
  SI.rho13[, iam(2, 2, M = 3)] <-
    (-2 * rho13 - (1 - rho13^2) * dbbb.drho13 / bbb) / (bbb * sd2^2)
  SI.rho13[, iam(3, 3, M = 3)] <-
    -1 * Sigmainv[, iam(3, 3, M = 3)] * dbbb.drho13 / bbb
  SI.rho13[, iam(1, 2, M = 3)] <-
    (rho23 - (rho13 * rho23 - rho12) * dbbb.drho13 / bbb) / (
     bbb * sd1 * sd2)
  SI.rho13[, iam(2, 3, M = 3)] <-
    (rho12 - (rho12 * rho13 - rho23) * dbbb.drho13 / bbb) / (
     bbb * sd2 * sd3)
  SI.rho13[, iam(1, 3, M = 3)] <-
    (-1 - (rho12 * rho23 - rho13) * dbbb.drho13 / bbb) / (
     bbb * sd1 * sd3)

    dl.drho12 <- -0.5 * dbbb.drho12 / bbb - 0.5 *
      c(mux5(x = ymatt.cen, cc = SI.rho12, M = 3, matrix.arg = TRUE))
    dl.drho23 <- -0.5 * dbbb.drho23 / bbb - 0.5 *
      c(mux5(x = ymatt.cen, cc = SI.rho23, M = 3, matrix.arg = TRUE))
    dl.drho13 <- -0.5 * dbbb.drho13 / bbb - 0.5 *
      c(mux5(x = ymatt.cen, cc = SI.rho13, M = 3, matrix.arg = TRUE))


    dmean1.deta <- dtheta.deta(mean1, .lmean1 )
    dmean2.deta <- dtheta.deta(mean2, .lmean2 )
    dmean3.deta <- dtheta.deta(mean3, .lmean3 )
    dsd1.deta   <- dtheta.deta(sd1  , .lsd1   )
    dsd2.deta   <- dtheta.deta(sd2  , .lsd2   )
    dsd3.deta   <- dtheta.deta(sd3  , .lsd3   )
    drho12.deta <- dtheta.deta(rho12, .lrho12 )
    drho23.deta <- dtheta.deta(rho23, .lrho23 )
    drho13.deta <- dtheta.deta(rho13, .lrho13 )
    dThetas.detas  <- cbind(dmean1.deta,
                            dmean2.deta,
                            dmean3.deta,
                            dsd1.deta,
                            dsd2.deta,
                            dsd3.deta,
                            drho12.deta,
                            drho23.deta,
                            drho13.deta)
    c(w) * cbind(dl.dmeans,  # dl.dmeans[, 1:3],
                 dl.dsd1,
                 dl.dsd2,
                 dl.dsd3,
                 dl.drho12,
                 dl.drho23,
                 dl.drho13) * dThetas.detas
  }), list( .lmean1 = lmean1, .lmean2 = lmean2, .lmean3 = lmean3,
            .emean1 = emean1, .emean2 = emean2, .emean3 = emean3,
            .lsd1   = lsd1  , .lsd2   = lsd2  , .lsd3   = lsd3  ,
            .esd1   = esd1  , .esd2   = esd2  , .esd3   = esd3  ,
            .lrho12 = lrho12, .lrho13 = lrho13, .lrho23 = lrho23,
            .erho12 = erho12, .erho13 = erho13, .erho23 = erho23 ))),

  weight = eval(substitute(expression({
    wz <- matrix(0, n, dimm(M))
    wz[, iam(1, 1, M)] <- Sigmainv[, iam(1, 1, M = 3)]
    wz[, iam(2, 2, M)] <- Sigmainv[, iam(2, 2, M = 3)]
    wz[, iam(3, 3, M)] <- Sigmainv[, iam(3, 3, M = 3)]
    wz[, iam(1, 2, M)] <- Sigmainv[, iam(1, 2, M = 3)]
    wz[, iam(2, 3, M)] <- Sigmainv[, iam(2, 3, M = 3)]
    wz[, iam(1, 3, M)] <- Sigmainv[, iam(1, 3, M = 3)]


if (FALSE) {
    wz[, iam(4, 4, M)] <- -1 / sd1^2 + 
      (1 - rho23^2 + 2 * bbb) / (bbb * sd1^2)
    wz[, iam(5, 5, M)] <- -1 / sd2^2 + 
      (1 - rho13^2 + 2 * bbb) / (bbb * sd2^2)
    wz[, iam(6, 6, M)] <- -1 / sd3^2 +
      (1 - rho12^2 + 2 * bbb) / (bbb * sd3^2)
    wz[, iam(4, 5, M)] <- 0 -
      rho12 * (rho13 * rho23 - rho12) / (sd1 * sd2 * bbb)
    wz[, iam(5, 6, M)] <- 0 -
      rho23 * (rho12 * rho13 - rho23) / (sd2 * sd3 * bbb)
    wz[, iam(4, 6, M)] <- 0 -
      rho13 * (rho12 * rho23 - rho13) / (sd1 * sd3 * bbb)
}

if (FALSE) {
    d2bbb.drho12.12 <- -2
    d2bbb.drho23.23 <- -2
    d2bbb.drho13.13 <- -2
    d2bbb.drho12.13 <-  2 * rho23
    d2bbb.drho12.23 <-  2 * rho13
    d2bbb.drho13.23 <-  2 * rho12
    wz[, iam(7, 7, M)] <-
      0.5 * (d2bbb.drho12.12 - dbbb.drho12 * dbbb.drho12 / bbb) / bbb
    wz[, iam(8, 8, M)] <-
      0.5 * (d2bbb.drho23.23 - dbbb.drho23 * dbbb.drho23 / bbb) / bbb
    wz[, iam(9, 9, M)] <-
      0.5 * (d2bbb.drho13.13 - dbbb.drho13 * dbbb.drho13 / bbb) / bbb
    wz[, iam(7, 8, M)] <-
      0.5 * (d2bbb.drho12.23 - dbbb.drho12 * dbbb.drho23 / bbb) / bbb
    wz[, iam(7, 9, M)] <-
      0.5 * (d2bbb.drho12.13 - dbbb.drho12 * dbbb.drho13 / bbb) / bbb
    wz[, iam(8, 9, M)] <-
      0.5 * (d2bbb.drho13.23 - dbbb.drho13 * dbbb.drho23 / bbb) / bbb
}



  mux43mat <- function(A, B, C, D, aa, bb) {

    s <- rep(0, length(A[, 1]))
    for (i1 in 1:3)
      for (i2 in 1:3)
        for (i3 in 1:3)
          s <- s + A[, iam(aa, i1, M = 3)] *
                   B[, iam(i1, i2, M = 3)] *
                   C[, iam(i2, i3, M = 3)] *
                   D[, iam(i3, bb, M = 3)]
    s
  }  # mux43mat



  Sigma <- matrix(0, n, dimm(3))  # sum(3:1)
  Sigma[, iam(1, 1, M = 3)] <- sd1^2      
  Sigma[, iam(2, 2, M = 3)] <- sd2^2      
  Sigma[, iam(3, 3, M = 3)] <- sd3^2      
  Sigma[, iam(1, 2, M = 3)] <- rho12 * sd1 * sd2
  Sigma[, iam(2, 3, M = 3)] <- rho23 * sd2 * sd3
  Sigma[, iam(1, 3, M = 3)] <- rho13 * sd1 * sd3



  for (ii in 1:3)
    wz[, iam(4, 4, M)] <-
    wz[, iam(4, 4, M)] +
    0.5 * mux43mat(Sigma, SI.sd1, Sigma, SI.sd1, ii, ii)

  for (ii in 1:3)
    wz[, iam(5, 5, M)] <-
    wz[, iam(5, 5, M)] +
    0.5 * mux43mat(Sigma, SI.sd2, Sigma, SI.sd2, ii, ii)

  for (ii in 1:3)
    wz[, iam(6, 6, M)] <-
    wz[, iam(6, 6, M)] +
    0.5 * mux43mat(Sigma, SI.sd3, Sigma, SI.sd3, ii, ii)

  for (ii in 1:3)
    wz[, iam(4, 5, M)] <-
    wz[, iam(4, 5, M)] +
    0.5 * mux43mat(Sigma, SI.sd1, Sigma, SI.sd2, ii, ii)

  for (ii in 1:3)
    wz[, iam(5, 6, M)] <-
    wz[, iam(5, 6, M)] +
    0.5 * mux43mat(Sigma, SI.sd2, Sigma, SI.sd3, ii, ii)

  for (ii in 1:3)
    wz[, iam(4, 6, M)] <-
    wz[, iam(4, 6, M)] +
    0.5 * mux43mat(Sigma, SI.sd1, Sigma, SI.sd3, ii, ii)






  for (ii in 1:3)
    wz[, iam(4, 7, M)] <-
    wz[, iam(4, 7, M)] +
    0.5 * mux43mat(Sigma, SI.sd1, Sigma, SI.rho12, ii, ii)

  for (ii in 1:3)
    wz[, iam(4, 8, M)] <-
    wz[, iam(4, 8, M)] +
    0.5 * mux43mat(Sigma, SI.sd1, Sigma, SI.rho23, ii, ii)

  for (ii in 1:3)
    wz[, iam(4, 9, M)] <-
    wz[, iam(4, 9, M)] +
    0.5 * mux43mat(Sigma, SI.sd1, Sigma, SI.rho13, ii, ii)

  for (ii in 1:3)
    wz[, iam(5, 7, M)] <-
    wz[, iam(5, 7, M)] +
    0.5 * mux43mat(Sigma, SI.sd2, Sigma, SI.rho12, ii, ii)

  for (ii in 1:3)
    wz[, iam(5, 8, M)] <-
    wz[, iam(5, 8, M)] +
    0.5 * mux43mat(Sigma, SI.sd2, Sigma, SI.rho23, ii, ii)

  for (ii in 1:3)
    wz[, iam(5, 9, M)] <-
    wz[, iam(5, 9, M)] +
    0.5 * mux43mat(Sigma, SI.sd2, Sigma, SI.rho13, ii, ii)

    for (ii in 1:3)
    wz[, iam(6, 7, M)] <-
    wz[, iam(6, 7, M)] +
    0.5 * mux43mat(Sigma, SI.sd3, Sigma, SI.rho12, ii, ii)

  for (ii in 1:3)
    wz[, iam(6, 8, M)] <-
    wz[, iam(6, 8, M)] +
    0.5 * mux43mat(Sigma, SI.sd3, Sigma, SI.rho23, ii, ii)

  for (ii in 1:3)
    wz[, iam(6, 9, M)] <-
    wz[, iam(6, 9, M)] +
    0.5 * mux43mat(Sigma, SI.sd3, Sigma, SI.rho13, ii, ii)




  for (ii in 1:3)
    wz[, iam(7, 7, M)] <-
    wz[, iam(7, 7, M)] +
    0.5 * mux43mat(Sigma, SI.rho12, Sigma, SI.rho12, ii, ii)

  for (ii in 1:3)
    wz[, iam(8, 8, M)] <-
    wz[, iam(8, 8, M)] +
    0.5 * mux43mat(Sigma, SI.rho23, Sigma, SI.rho23, ii, ii)

  for (ii in 1:3)
    wz[, iam(9, 9, M)] <-
    wz[, iam(9, 9, M)] +
    0.5 * mux43mat(Sigma, SI.rho13, Sigma, SI.rho13, ii, ii)

  for (ii in 1:3)
    wz[, iam(7, 8, M)] <-
    wz[, iam(7, 8, M)] +
    0.5 * mux43mat(Sigma, SI.rho12, Sigma, SI.rho23, ii, ii)

  for (ii in 1:3)
    wz[, iam(8, 9, M)] <-
    wz[, iam(8, 9, M)] +
    0.5 * mux43mat(Sigma, SI.rho23, Sigma, SI.rho13, ii, ii)

  for (ii in 1:3)
    wz[, iam(7, 9, M)] <-
    wz[, iam(7, 9, M)] +
    0.5 * mux43mat(Sigma, SI.rho12, Sigma, SI.rho13, ii, ii)



  ind5 <- iam(NA, NA, M = M, both = TRUE, diag = TRUE)
  wz <- wz * dThetas.detas[, ind5$row.index] *
             dThetas.detas[, ind5$col.index]
  c(w) * wz
  }), list( .lmean1 = lmean1, .lmean2 = lmean2, .lmean3 = lmean3,
            .emean1 = emean1, .emean2 = emean2, .emean3 = emean3,
            .lsd1   = lsd1  , .lsd2   = lsd2  , .lsd3   = lsd3  ,
            .esd1   = esd1  , .esd2   = esd2  , .esd3   = esd3  ,
            .lrho12 = lrho12, .lrho13 = lrho13, .lrho23 = lrho23,
            .erho12 = erho12, .erho13 = erho13, .erho23 = erho23 ))))
}  # trinormal















dbiclaytoncop <- function(x1, x2, apar = 0, log = FALSE) {
  if (!is.logical(log.arg <- log) || length(log) != 1)
    stop("bad input for argument 'log'")
  rm(log)

  A <- x1^(-apar) + x2^(-apar) - 1
  logdensity <- log1p(apar) -
                (1 + apar) * (log(x1) + log(x2)) -
                (2 + 1 / apar) * log(abs(A))  # Avoid warning

  out.square <- (x1 < 0) | (x1 > 1) | (x2 < 0) | (x2 > 1)
  logdensity[out.square] <- log(0.0)


  index0 <- (rep_len(apar, length(A)) < sqrt(.Machine$double.eps))
  if (any(index0))
    logdensity[index0] <- log(1.0)


  index1 <- (rep_len(apar, length(A)) < 0.0) | (A < 0.0)
  if (any(index1))
    logdensity[index1] <- NaN






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



rbiclaytoncop <- function(n, apar = 0) {
  if (any(apar < 0))
    stop("argument 'apar' must be greater or equal to 0")

  u1 <- runif(n = n)
  v2 <- runif(n = n)

  u2 <- (u1^(-apar) *
        (v2^(-apar / (1 + apar)) - 1) + 1)^(-1 / apar)


  index0 <- (rep_len(apar, length(u1)) < sqrt(.Machine$double.eps))
  if (any(index0))
    u2[index0] <- runif(sum(index0))

  cbind(u1, u2)
}  # dbiclaytoncop



 biclaytoncop <- function(lapar    = "loglink",
                          iapar    = NULL,
                          imethod   = 1,
                          parallel  = FALSE,
                          zero = NULL) {

  apply.parint <- TRUE


  lapar <- as.list(substitute(lapar))
  eapar <- link2list(lapar)
  lapar <- attr(eapar, "function.name")


  if (length(iapar) && any(iapar <= 0))
    stop("argument 'iapar' must have values in (0, Inf)")



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



  new("vglmff",
  blurb = c("Bivariate Clayton copula distribution)\n",
            "Links:    ", namesof("apar", lapar, earg = eapar)),

  constraints = eval(substitute(expression({
    constraints <- cm.VGAM(matrix(1, M, 1), x = x,
                           bool = .parallel ,
                           constraints = constraints,
                           apply.int = .apply.parint )

    constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
                                predictors.names = predictors.names,
                                M1 = 1)
  }), list( .zero = zero,
            .apply.parint = apply.parint,
            .parallel = parallel ))),

  infos = eval(substitute(function(...) {
    list(M1 = 1,
         Q1 = 2,
         apply.parint = .apply.parint ,
         parameters.names = c("apar"),
         lapar = .lapar ,
         parallel = .parallel ,
         zero = .zero )
    }, list( .zero = zero,
             .apply.parint = apply.parint,
             .lapar = lapar,
             .parallel = parallel ))),

  initialize = eval(substitute(expression({
    M1 <- 1
    Q1 <- 2

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

    w <- temp5$w
    y <- temp5$y


    ncoly <- ncol(y)
    extra$ncoly <- ncoly
    extra$M1 <- M1
    extra$Q1 <- Q1
    M <- M1 * (ncoly / Q1)
    mynames1 <- param.names("apar", M / M1, skip1 = TRUE)
    predictors.names <-
      namesof(mynames1, .lapar , earg = .eapar , short = TRUE)

    extra$colnames.y  <- colnames(y)

    if (!length(etastart)) {

      apar.init <- matrix(if (length( .iapar )) .iapar else NA_real_,
                          n, M / M1, byrow = TRUE)

      if (!length( .iapar ))
        for (spp. in 1:(M / M1)) {
          ymatj <- y[, (Q1 * spp. - 1):(Q1 * spp.)]


          apar.init0 <- if ( .imethod == 1) {
            k.tau <- kendall.tau(ymatj[, 1], ymatj[, 2], exact = FALSE,
                                 max.n = 500)

            max(0.1, 2 * k.tau / (1 - k.tau))  # Must be positive
          } else if ( .imethod == 2) {
            spearman.rho <-  max(0.05, cor(ymatj[, 1],
                                           ymatj[, 2], meth = "spearman"))
            rhobitlink(spearman.rho)
          } else {
            pearson.rho <- max(0.05, cor(ymatj[, 1], ymatj[, 2]))
            rhobitlink(pearson.rho)
          }

          if (anyNA(apar.init[, spp.]))
            apar.init[, spp.] <- apar.init0
        }

      etastart <- theta2eta(apar.init, .lapar , earg = .eapar )
    }
  }), list( .imethod = imethod,
            .lapar = lapar,
            .eapar = eapar,
            .iapar = iapar ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    NOS <- NCOL(eta) / c(M1 = 1)
    Q1 <- 2
    fv.mat <- matrix(0.5, NROW(eta), NOS * Q1)
    label.cols.y(fv.mat, colnames.y = extra$colnames.y, NOS = NOS)
  }  , list( .lapar = lapar,
             .eapar = eapar ))),

  last = eval(substitute(expression({
    M1 <- extra$M1
    Q1 <- extra$Q1
    misc$link <- rep_len( .lapar , M)
    temp.names <- mynames1
    names(misc$link) <- temp.names

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

    misc$M1 <- M1
    misc$Q1 <- Q1
    misc$imethod <- .imethod
    misc$expected <- TRUE
    misc$parallel  <- .parallel
    misc$apply.parint <- .apply.parint
    misc$multipleResponses <- TRUE
  }) , list( .imethod = imethod,
             .parallel = parallel, .apply.parint = apply.parint,
             .lapar = lapar,
             .eapar = eapar ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    Alpha <- eta2theta(eta, .lapar , earg = .eapar )

    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {

      ll.elts <-
        c(w) * dbiclaytoncop(x1  = c(y[, c(TRUE, FALSE)]),
                             x2  = c(y[, c(FALSE, TRUE)]),
                             apar = c(Alpha), log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  } , list( .lapar = lapar, .eapar = eapar, .imethod = imethod ))),
  vfamily = c("biclaytoncop"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    Alpha <- eta2theta(eta, .lapar , earg = .eapar )
    okay1 <- all(is.finite(Alpha)) && all(0 < Alpha)
    okay1
  } , list( .lapar = lapar, .eapar = eapar, .imethod = imethod ))),

  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)
    Alpha <- eta2theta(eta, .lapar , earg = .eapar )
    rbiclaytoncop(nsim * length(Alpha),
                  apar = c(Alpha))
  }  , list( .lapar = lapar,
             .eapar = eapar ))),


  deriv = eval(substitute(expression({
    Alpha <- eta2theta(eta, .lapar , earg = .eapar )
    Yindex1 <- extra$Q1 * (1:(extra$ncoly/extra$Q1)) - 1
    Yindex2 <- extra$Q1 * (1:(extra$ncoly/extra$Q1))





    AA <- y[, Yindex1]^(-Alpha) + y[, Yindex2]^(-Alpha) - 1
    dAA.dapar <- -y[, Yindex1]^(-Alpha) * log(y[, Yindex1]) -
                  y[, Yindex2]^(-Alpha) * log(y[, Yindex2])
    dl.dapar <- 1 / (1 + Alpha) - log(y[, Yindex1] * y[, Yindex2]) -
                dAA.dapar / AA * (2 + 1 / Alpha ) + log(AA) / Alpha^2



    dapar.deta <- dtheta.deta(Alpha, .lapar , earg = .eapar )

    dl.deta <- c(w) * cbind(dl.dapar) * dapar.deta
    dl.deta
  }), list( .lapar = lapar,
            .eapar = eapar,
            .imethod = imethod ))),

  weight = eval(substitute(expression({
    par <- Alpha + 1
    denom1 <- (3 * par -2) * (2 * par - 1)
    denom2 <- 2 * (par - 1)
    v1 <- trigamma(1 / denom2)
    v2 <- trigamma(par / denom2)
    v3 <- trigamma((2 * par - 1) / denom2)
    Rho. <- (1 + par  * (v1 - v2) / denom2 +
                        (v2 - v3) / denom2) / denom1

    ned2l.dapar  <- 1 / par^2 + 2 / (par * (par - 1) * (2 * par - 1)) +
                    4 * par / (3 * par - 2) -
                    2 * (2 * par - 1) * Rho. / (par - 1)

    wz <- ned2l.dapar * dapar.deta^2
    c(w) * wz
  }), list( .lapar = lapar,
            .eapar = eapar,
            .imethod = imethod ))))
}  # biclaytoncop

















dbistudentt <- function(x1, x2, df, rho = 0, log = FALSE) {




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

  logdensity <-
    -(df/2 + 1) * log1p(
    (x1^2 + x2^2 - 2 * rho * x1 * x2) / (df * (1 - rho^2))) -
    log(2 * pi) - 0.5 * log1p(-rho^2)  # -

  logdensity[df <= 0] <- NaN  # Not picked up by dt().

  logdensity[is.infinite(x1) | is.infinite(x2)] <- log(0)

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




if (FALSE)
bistudent.deriv.dof <-  function(u, v, nu, rho) {


  t1 <- qt(u, nu, 1, 0)
  t2 <- qt(v, nu, 1, 0)
  t3 <- -(nu + 2.0) / 2.0
  t10 <- nu * (1.0 - rho * rho)
  t4 <- -2.0 * t1 * t2 / t10
  t11 <- (t1 * t1 + t2 * t2 - 2.0 * rho * t1 * t2)
  t5 <- 2.0 * t11 * rho / t10 / (1.0 - rho * rho)
  t6 <- 1.0 + (t11 / t10)
  t7 <- rho / (1.0 - rho * rho)
  out <- (t3 * (t4 + t5) / t6  +  t7)
}







 bistudentt <-
   function(ldf     = "logloglink",
            lrho    = "rhobitlink",
            idf     = NULL,
            irho    = NULL,
            imethod = 1,
            parallel = FALSE,
            zero = "rho") {




  apply.parint <- TRUE

  ldof <- as.list(substitute(ldf))
  edof <- link2list(ldof)
  ldof <- attr(edof, "function.name")

  lrho <- as.list(substitute(lrho))
  erho <- link2list(lrho)
  lrho <- attr(erho, "function.name")


  idof <- idf
  if (length(idof) &&
      any(idof <= 1))
    stop("argument 'idf' must have values in (1,Inf)")


  if (length(irho) &&
      any(abs(irho) >= 1))
    stop("argument 'irho' must have values in (-1,1)")



  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("Bivariate student-t distribution\n",
            "Links:    ",
            namesof("df",  ldof, earg = edof), ", ",
            namesof("rho", lrho, earg = erho)),

  constraints = eval(substitute(expression({
    constraints <- cm.VGAM(matrix(1, M, 1), x = x,
                           bool = .parallel ,
                           constraints = constraints,
                           apply.int = .apply.parint )

    constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
                                predictors.names = predictors.names,
                                M1 = 2)
  }), list( .zero = zero,
            .apply.parint = apply.parint,
            .parallel = parallel ))),

  infos = eval(substitute(function(...) {
    list(M1 = 2,
         Q1 = 2,
         parameters.names = c("df", "rho"),
         apply.parint = .apply.parint ,
         parallel = .parallel ,
         zero = .zero )
  }, list( .zero = zero,
           .apply.parint = apply.parint,
           .parallel = parallel ))),

  initialize = eval(substitute(expression({
    M1 <- 2
    Q1 <- 2

    temp5 <-
    w.y.check(w = w, y = y,
              ncol.w.max = Inf,
              ncol.y.max = Inf,
              ncol.y.min = Q1,
              out.wy = TRUE,
              colsyperw = Q1,
              maximize = TRUE)

    w <- temp5$w
    y <- temp5$y


    ncoly <- ncol(y)
    extra$ncoly <- ncoly
    extra$M1 <- M1
    extra$Q1 <- Q1
    M <- M1 * (ncoly / Q1)
    mynames1 <- param.names("df",  M / M1, skip1 = TRUE)
    mynames2 <- param.names("rho", M / M1, skip1 = TRUE)
    predictors.names <- c(
      namesof(mynames1, .ldof , earg = .edof , short = TRUE),
      namesof(mynames2, .lrho , earg = .erho , short = TRUE))[
              interleave.VGAM(M, M1 = M1)]


    extra$colnames.y  <- colnames(y)

    if (!length(etastart)) {

      dof.init <- matrix(if (length( .idof )) .idof else 0 + NA,
                         n, M / M1, byrow = TRUE)
      rho.init <- matrix(if (length( .irho )) .irho else 0 + NA,
                         n, M / M1, byrow = TRUE)

      if (!length( .idof ) || !length( .irho ))
      for (spp. in 1:(M / M1)) {
        ymatj <- y[, (M1 * spp. - 1):(M1 * spp.)]


        dof.init0 <- if ( .imethod == 1) {


          2 + rexp(n = 1, rate = 0.1)
        } else {
          10
        }

        if (anyNA(dof.init[, spp.]))
          dof.init[, spp.] <- dof.init0


        rho.init0 <- if ( .imethod == 2) {
          runif(n, min = -1 + 0.1, max = 1 - 0.1)
        } else {
          cor(ymatj[, 1], ymatj[, 2])
        }

        if (anyNA(rho.init[, spp.]))
          rho.init[, spp.] <- rho.init0

      }

      etastart <-
        cbind(theta2eta(dof.init, .ldof , earg = .edof ),
              theta2eta(rho.init, .lrho , earg = .erho ))

      etastart <- etastart[, interleave.VGAM(M, M1 = M1)]

    }
  }), list( .imethod = imethod,
            .lrho = lrho, .ldof = ldof,
            .erho = erho, .edof = edof,
            .idof = idof, .irho = irho ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    NOS <- ncol(eta) / c(M1 = 2)
    Q1 <- 2
    fv.mat <- matrix(0, nrow(eta), Q1 * NOS)
    label.cols.y(fv.mat, colnames.y = extra$colnames.y, NOS = NOS)
  }  , list( .lrho = lrho, .ldof = ldof,
             .erho = erho, .edof = edof ))),

  last = eval(substitute(expression({

    M1 <- extra$M1
    Q1 <- extra$Q1
    misc$link <-
      c(rep_len( .ldof , M / M1),
        rep_len( .lrho , M / M1))[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:(M / M1)) {
      misc$earg[[M1*ii-1]] <- .edof
      misc$earg[[M1*ii  ]] <- .erho
    }

    misc$M1 <- M1
    misc$Q1 <- Q1
    misc$imethod <- .imethod
    misc$expected <- TRUE
    misc$parallel  <- .parallel
    misc$apply.parint <- .apply.parint
    misc$multipleResponses <- TRUE

  }) , list( .imethod = imethod,
             .parallel = parallel,
             .apply.parint = apply.parint,
             .lrho = lrho, .ldof = ldof,
             .erho = erho, .edof = edof ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    Dof <- eta2theta(eta[, c(TRUE, FALSE), drop = FALSE],
                     .ldof , earg = .edof )
    Rho <- eta2theta(eta[, c(FALSE, TRUE), drop = FALSE],
                     .lrho , earg = .erho )

    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      Yindex1 <- extra$Q1 * (1:(extra$ncoly/extra$Q1)) - 1
      Yindex2 <- extra$Q1 * (1:(extra$ncoly/extra$Q1))
      ll.elts <-
        c(w) * dbistudentt(x1  = y[, Yindex1, drop = FALSE],
                           x2  = y[, Yindex2, drop = FALSE],
                           df  = Dof,
                           rho = Rho, log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .lrho = lrho, .ldof = ldof,
           .erho = erho, .edof = edof,
           .imethod = imethod ))),
  vfamily = c("bistudentt"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    Dof <- eta2theta(eta[, c(TRUE, FALSE), drop = FALSE],
                     .ldof , earg = .edof )
    Rho <- eta2theta(eta[, c(FALSE, TRUE), drop = FALSE],
                     .lrho , earg = .erho )
    okay1 <- all(is.finite(Dof)) && all(0 < Dof) &&
             all(is.finite(Rho)) && all(abs(Rho) < 1)
    okay1
  }, list( .lrho = lrho, .ldof = ldof,
           .erho = erho, .edof = edof,
           .imethod = imethod ))),
  deriv = eval(substitute(expression({
    M1 <- Q1 <- 2
    Dof <- eta2theta(eta[, c(TRUE, FALSE), drop = FALSE],
                     .ldof , earg = .edof )
    Rho <- eta2theta(eta[, c(FALSE, TRUE), drop = FALSE],
                     .lrho , earg = .erho )
    Yindex1 <- extra$Q1 * (1:(extra$ncoly/extra$Q1)) - 1
    Yindex2 <- extra$Q1 * (1:(extra$ncoly/extra$Q1))


    x1 <- c(y[, Yindex1])  # Convert into a vector
    x2 <- c(y[, Yindex2])

    dee3 <- deriv3( ~
        -(Dof/2 + 1) * log(1 +
        (x1^2 + x2^2 - 2 * Rho * x1 * x2) / (Dof * (1 - Rho^2))) -
        log(2 * pi) - 0.5 * log(1 - Rho^2),
        namevec = c("Dof", "Rho"), hessian = FALSE)
    eval.d3 <- eval(dee3)

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

    dl.ddof <- matrix(dl.dthetas[, "Dof"], n, length(Yindex1))
    dl.drho <- matrix(dl.dthetas[, "Rho"], n, length(Yindex2))


  if (FALSE) {
    dd <- cbind(y, Rho, Dof)
    pp <- apply(dd, 1, function(x)
                BiCopPDF(x[1], x[2], family = 2, x[3], x[4]))
    alt.dl.ddof <- apply(dd, 1, function(x)
                     BiCopDeriv(x[1], x[2], family = 2,
                                x[3], x[4], "par2")) / pp
    alt.dl.drho <- apply(dd, 1, function(x)
                     BiCopDeriv(x[1], x[2], family = 2,
                                x[3], x[4], "par")) / pp



  }





    ddof.deta <- dtheta.deta(Dof, .ldof , earg = .edof )
    drho.deta <- dtheta.deta(Rho, .lrho , earg = .erho )

    ans <- c(w) * cbind(dl.ddof * ddof.deta,
                        dl.drho * drho.deta)
    ans <- ans[, interleave.VGAM(M, M1 = M1)]
    ans
  }), list( .lrho = lrho, .ldof = ldof,
            .erho = erho, .edof = edof,
            .imethod = imethod ))),

  weight = eval(substitute(expression({
    wz11 <- beta(2, Dof / 2) / Dof -
            beta(3, Dof / 2) * (Dof + 2) / (4 * Dof)
    wz12 <- -Rho / (2 * (1 - Rho^2)) * (beta(2, Dof / 2) -
            beta(3, Dof / 2) * (Dof + 2) / 2)
    wz22 <- (1 + Rho^2) / (1 - Rho^2)^2 +
            (Dof^2 + 2 * Dof) * Rho^2 *
             beta(3, Dof / 2) / (4 * (1 - Rho^2)^2)
    wz22 <- wz22 + (Dof^2 + 2 * Dof) * (2 - 3 * Rho^2 + Rho^6) *
            beta(3, Dof / 2) / (16 * (1 - Rho^2)^4)
    wz22 <- wz22 + (Dof^2 + 2 * Dof) * (1 + Rho^2) *  # Replace - by + ???
            beta(2, Dof / 2) / (4 * (1 - Rho^2)^2)  # denom == 4 or 2 ???
    ned2l.ddof2   <- wz11
    ned2l.ddofrho <- wz12
    ned2l.drho2   <- wz22

    wz <- array(c(c(w) * ned2l.ddof2 * ddof.deta^2,
                  c(w) * ned2l.drho2 * drho.deta^2,
                  c(w) * ned2l.ddofrho * ddof.deta * drho.deta),
                dim = c(n, M / M1, 3))
    wz <- arwz2wz(wz, M = M, M1 = M1)
    wz
  }), list( .lrho = lrho, .ldof = ldof,
            .erho = erho, .edof = edof,
            .imethod = imethod ))))
}







dbinormcop <- function(x1, x2, rho = 0, log = FALSE) {
  if (!is.logical(log.arg <- log) || length(log) != 1)
    stop("bad input for argument 'log'")
  rm(log)

  x1 <- qnorm(x1)
  x2 <- qnorm(x2)

  logdensity <- (2 * rho * x1 * x2 -
                 rho^2 * (x1^2 + x2^2)) / (2 * (1 - rho^2)) -
                0.5 * log1p(-rho^2)

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




pbinormcop <- function(q1, q2, rho = 0) {

  if (!is.Numeric(q1, positive = TRUE) ||
      any(q1 >= 1))
    stop("bad input for argument 'q1'")
  if (!is.Numeric(q2, positive = TRUE) ||
      any(q2 >= 1))
    stop("bad input for argument 'q2'")
  if (!is.Numeric(rho) ||
      any(abs(rho) >= 1))
    stop("bad input for argument 'rho'")

  pbinorm(q1 = qnorm(q1), q2 = qnorm(q2), cov12 = rho)
}


rbinormcop <- function(n, rho = 0  #, inverse = FALSE
                      ) {

  inverse <- FALSE
  ymat <- rbinorm(n = n, cov12 = rho)
  if (inverse) {
    ymat
  } else {
    cbind(y1 = pnorm(ymat[, 1]),
          y2 = pnorm(ymat[, 2]))
  }
}





 binormalcop <- function(lrho    = "rhobitlink",
                         irho    = NULL,
                         imethod = 1,
                         parallel = FALSE,
                         zero = NULL) {



  apply.parint <- TRUE


  lrho <- as.list(substitute(lrho))
  erho <- link2list(lrho)
  lrho <- attr(erho, "function.name")


  if (length(irho) &&
      any(abs(irho) >= 1))
    stop("argument 'irho' must have values in (-1,1)")



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

  new("vglmff",
  blurb = c("Gaussian copula ",
            "(based on the bivariate normal distribution)\n",
            "Links:    ",
            namesof("rho", lrho, earg = erho)),

  constraints = eval(substitute(expression({
    constraints <- cm.VGAM(matrix(1, M, 1), x = x,
                           bool = .parallel ,
                           constraints = constraints,
                           apply.int = .apply.parint )

    constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
                                predictors.names = predictors.names,
                                M1 = 1)
  }), list( .zero = zero,
            .apply.parint = apply.parint,
            .parallel = parallel ))),

  infos = eval(substitute(function(...) {
    list(M1 = 1,
         Q1 = 2,
         parameters.names = c("rho"),
         apply.parint = .apply.parint ,
         parallel = .parallel ,
         zero = .zero )
  }, list( .zero = zero,
           .apply.parint = apply.parint,
           .parallel = parallel ))),

  initialize = eval(substitute(expression({
    M1 <- 1
    Q1 <- 2

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

    w <- temp5$w
    y <- temp5$y


    ncoly <- ncol(y)
    extra$ncoly <- ncoly
    extra$M1 <- M1
    extra$Q1 <- Q1
    M <- M1 * (ncoly / Q1)
    mynames1 <- param.names("rho", M / M1, skip1 = TRUE)
    predictors.names <- c(
      namesof(mynames1, .lrho , earg = .erho , short = TRUE))


    extra$colnames.y  <- colnames(y)

    if (!length(etastart)) {

      rho.init <- matrix(if (length( .irho )) .irho else 0 + NA,
                         n, M / M1, byrow = TRUE)

      if (!length( .irho ))
      for (spp. in 1:(M / M1)) {
        ymatj <- y[, (Q1 * spp. - 1):(Q1 * spp.)]


        rho.init0 <- if ( .imethod == 1) {
          sin(kendall.tau(ymatj[, 1], ymatj[, 2],
                          exact = FALSE,
                          max.n = 200) * pi / 2)
        } else if ( .imethod == 2) {
          sin(cor(ymatj[, 1], ymatj[, 2],
                  method = "spearman") * pi / 6) * 2
        } else {
          cor(ymatj[, 1], ymatj[, 2])
        }





        if (anyNA(rho.init[, spp.]))
          rho.init[, spp.] <- rho.init0
      }

      etastart <- theta2eta(rho.init, .lrho , earg = .erho )
    }
  }), list( .imethod = imethod,
            .lrho = lrho,
            .erho = erho,
            .irho = irho ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    NOS <- NCOL(eta) / c(M1 = 1)
    Q1 <- 2
    fv.mat <- matrix(0.5, NROW(eta), NOS * Q1)
    label.cols.y(fv.mat, colnames.y = extra$colnames.y, NOS = NOS)
  }  , list( .lrho = lrho,
             .erho = erho ))),

  last = eval(substitute(expression({
    M1 <- extra$M1
    Q1 <- extra$Q1
    misc$link <- rep_len( .lrho , M)
    temp.names <- mynames1
    names(misc$link) <- temp.names

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

    misc$M1 <- M1
    misc$Q1 <- Q1
    misc$imethod <- .imethod
    misc$expected <- TRUE
    misc$parallel  <- .parallel
    misc$apply.parint <- .apply.parint
    misc$multipleResponses <- TRUE

  }) , list( .imethod = imethod,
             .parallel = parallel,
             .apply.parint = apply.parint,
             .lrho = lrho,
             .erho = erho ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    Rho <- eta2theta(eta, .lrho , earg = .erho )

    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      Yindex1 <- extra$Q1 * (1:(extra$ncoly/extra$Q1)) - 1
      Yindex2 <- extra$Q1 * (1:(extra$ncoly/extra$Q1))
      ll.elts <-
        c(w) * dbinormcop(x1  = y[, Yindex1, drop = FALSE],
                          x2  = y[, Yindex2, drop = FALSE],
                          rho = Rho, log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  } , list( .lrho = lrho,
            .erho = erho,
            .imethod = imethod ))),
  vfamily = c("binormalcop"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    Rho <- eta2theta(eta, .lrho , earg = .erho )
    okay1 <- all(is.finite(Rho)) && all(abs(Rho) < 1)
    okay1
  }, list( .lrho = lrho, .erho = erho, .imethod = imethod ))),



  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)
    Rho <- eta2theta(eta, .lrho , earg = .erho )
    rbinormcop(nsim * length(Rho),
               rho = c(Rho))
  }  , list( .lrho = lrho,
             .erho = erho ))),



  deriv = eval(substitute(expression({
    Rho <- eta2theta(eta, .lrho , earg = .erho )
    Yindex1 <- extra$Q1 * (1:(extra$ncoly/extra$Q1)) - 1
    Yindex2 <- extra$Q1 * (1:(extra$ncoly/extra$Q1))

    temp7 <- 1 - Rho^2
    q.y <- qnorm(y)

    dl.drho <- ((1 + Rho^2) * q.y[, Yindex1] * q.y[, Yindex2] -
                Rho * (q.y[, Yindex1]^2 + q.y[, Yindex2]^2)) / temp7^2 +
                Rho / temp7

    drho.deta <- dtheta.deta(Rho, .lrho , earg = .erho )

    c(w) * cbind(dl.drho) * drho.deta
  }), list( .lrho = lrho,
            .erho = erho,
            .imethod = imethod ))),

  weight = eval(substitute(expression({
    ned2l.drho  <- (1 + Rho^2) / temp7^2
    wz <- ned2l.drho * drho.deta^2
    c(w) * wz
  }), list( .lrho = lrho,
            .erho = erho,
            .imethod = imethod ))))
}











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


 bilogistic  <- function(llocation = "identitylink",
                         lscale = "loglink",
                         iloc1 = NULL, iscale1 = NULL,
                         iloc2 = NULL, iscale2 = NULL,
                         imethod = 1,
                         nsimEIM = 250,
                         zero = NULL) {

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

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



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



  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("Bivariate logistic distribution\n\n",
            "Link:    ",
            namesof("location1", llocat, elocat), ", ",
            namesof("scale1",    lscale, escale), ", ",
            namesof("location2", llocat, elocat), ", ",
            namesof("scale2",    lscale, escale), "\n", "\n",
            "Means:     location1, location2"),
  constraints = eval(substitute(expression({
    constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
                                predictors.names = predictors.names,
                                M1 = 4)
  }), list( .zero = zero))),


  infos = eval(substitute(function(...) {
    list(M1 = 4,
         Q1 = 2,
         expected = FALSE,
         parameters.names =
             c("location1", "scale1", "location2", "scale2"),
         multipleResponses = FALSE,
         zero = .zero )
  }, list( .zero = zero
         ))),


  initialize = eval(substitute(expression({

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

    extra$colnames.y  <- colnames(y)


    predictors.names <-
      c(namesof("location1", .llocat, .elocat , tag = FALSE),
        namesof("scale1",    .lscale, .escale , tag = FALSE),
        namesof("location2", .llocat, .elocat , tag = FALSE),
        namesof("scale2",    .lscale, .escale , tag = FALSE))

    if (!length(etastart)) {
      if ( .imethod == 1) {
        locat.init1 <- y[, 1]
        scale.init1 <- sqrt(3) * sd(y[, 1]) / pi
        locat.init2 <- y[, 2]
        scale.init2 <- sqrt(3) * sd(y[, 2]) / pi
      } else {
        locat.init1 <- median(rep(y[, 1], w))
        locat.init2 <- median(rep(y[, 2], w))
        const4 <- sqrt(3) / (sum(w) * pi)
        scale.init1 <- const4 * sum(c(w) *(y[, 1] - locat.init1)^2)
        scale.init2 <- const4 * sum(c(w) *(y[, 2] - locat.init2)^2)
      }
      loc1.init <- if (length( .iloc1 )) rep_len( .iloc1 , n) else
                                         rep_len(locat.init1, n)
      loc2.init <- if (length( .iloc2 )) rep_len( .iloc2 , n) else
                                         rep_len(locat.init2, n)
      scale1.init <- if (length( .iscale1 )) rep_len( .iscale1 , n) else
                                             rep_len(1, n)
      scale2.init <- if (length( .iscale2 )) rep_len( .iscale2 , n) else
                                             rep_len(1, n)

      if ( .llocat == "loglink")
        locat.init1 <- abs(locat.init1) + 0.001
      if ( .llocat == "loglink")
        locat.init2 <- abs(locat.init2) + 0.001

      etastart <-
        cbind(theta2eta(locat.init1, .llocat , .elocat ),
              theta2eta(scale1.init, .lscale , .escale ),
              theta2eta(locat.init2, .llocat , .elocat ),
              theta2eta(scale2.init, .lscale , .escale ))
    }
  }), list(.imethod = imethod,
           .iloc1 = iloc1, .iloc2 = iloc2,
           .llocat = llocat, .lscale = lscale,
           .elocat = elocat, .escale = escale,
           .iscale1 = iscale1, .iscale2 = iscale2))),
  linkinv = function(eta, extra = NULL) {
    NOS <- NCOL(eta) / c(M1 = 4)
    fv.mat <- eta[, 1:2]
    label.cols.y(fv.mat, colnames.y = extra$colnames.y, NOS = NOS)
  },
  last = eval(substitute(expression({
    misc$link <-    c(location1 = .llocat , scale1 = .lscale ,
                      location2 = .llocat , scale2 = .lscale )

    misc$earg <- list(location1 = .elocat , scale1 = .escale ,
                      location2 = .elocat , scale2 = .escale )

  }), list( .llocat = llocat, .lscale = lscale,
            .elocat = elocat, .escale = escale ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    locat1 <- eta2theta(eta[, 1], .llocat , .elocat )
    Scale1 <- eta2theta(eta[, 2], .lscale , .escale )
    locat2 <- eta2theta(eta[, 3], .llocat , .elocat )
    Scale2 <- eta2theta(eta[, 4], .lscale , .escale )

    zedd1 <- (y[, 1]-locat1) / Scale1
    zedd2 <- (y[, 2]-locat2) / Scale2

    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <-
        c(w) * (-zedd1 - zedd2 -
                3 * log1p(exp(-zedd1) + exp(-zedd2)) -
                log(Scale1) - log(Scale2))
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .llocat = llocat, .lscale = lscale,
           .elocat = elocat, .escale = escale ))),
  vfamily = c("bilogistic"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    locat1 <- eta2theta(eta[, 1], .llocat , .elocat )
    Scale1 <- eta2theta(eta[, 2], .lscale , .escale )
    locat2 <- eta2theta(eta[, 3], .llocat , .elocat )
    Scale2 <- eta2theta(eta[, 4], .lscale , .escale )
    okay1 <- all(is.finite(locat1)) &&
             all(is.finite(Scale1)) && all(0 < Scale1) &&
             all(is.finite(locat2)) &&
             all(is.finite(Scale2)) && all(0 < Scale2)
    okay1
  }, list( .llocat = llocat, .lscale = lscale,
           .elocat = elocat, .escale = escale ))),


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

    pwts <- if (length(pwts <- object@prior.weights) > 0)
              pwts else weights(object, type = "prior")
    if (any(pwts != 1))
      warning("ignoring prior weights")
    eta <- predict(object)
    locat1 <- eta2theta(eta[, 1], .llocat , .elocat )
    Scale1 <- eta2theta(eta[, 2], .lscale , .escale )
    locat2 <- eta2theta(eta[, 3], .llocat , .elocat )
    Scale2 <- eta2theta(eta[, 4], .lscale , .escale )
    rbilogis(nsim * length(locat1),
             loc1 = locat1, scale1 = Scale1,
             loc2 = locat2, scale2 = Scale2)
  }, list( .llocat = llocat, .lscale = lscale,
           .elocat = elocat, .escale = escale ))),




  deriv = eval(substitute(expression({
    locat1 <- eta2theta(eta[, 1], .llocat , .elocat )
    Scale1 <- eta2theta(eta[, 2], .lscale , .escale )
    locat2 <- eta2theta(eta[, 3], .llocat , .elocat )
    Scale2 <- eta2theta(eta[, 4], .lscale , .escale )

    zedd1 <- (y[, 1] - locat1) / Scale1
    zedd2 <- (y[, 2] - locat2) / Scale2
    ezedd1 <- exp(-zedd1)
    ezedd2 <- exp(-zedd2)
    denom <- 1 + ezedd1 + ezedd2

    dl.dlocat1 <- (1 - 3 * ezedd1 / denom) / Scale1
    dl.dlocat2 <- (1 - 3 * ezedd2 / denom) / Scale2
    dl.dscale1 <- (zedd1 - 1 - 3 * ezedd1 * zedd1 / denom) / Scale1
    dl.dscale2 <- (zedd2 - 1 - 3 * ezedd2 * zedd2 / denom) / Scale2

    dlocat1.deta <- dtheta.deta(locat1, .llocat , .elocat )
    dlocat2.deta <- dtheta.deta(locat2, .llocat , .elocat )
    dscale1.deta <- dtheta.deta(Scale1, .lscale , .escale )
    dscale2.deta <- dtheta.deta(Scale2, .lscale , .escale )

    derivnew <- c(w) * cbind(dl.dlocat1 * dlocat1.deta,
                             dl.dscale1 * dscale1.deta,
                             dl.dlocat2 * dlocat2.deta,
                             dl.dscale2 * dscale2.deta)
    derivnew
  }), list( .llocat = llocat, .lscale = lscale,
            .elocat = elocat, .escale = escale ))),
  weight = eval(substitute(expression({
    run.varcov <- 0
    ind1 <- iam(NA_real_, NA_real_, M = M, both = TRUE, diag = TRUE)
    for (ii in 1:( .nsimEIM )) {
      ysim <- rbilogis( .nsimEIM * length(locat1),
                       loc1 = locat1, scale1 = Scale1,
                       loc2 = locat2, scale2 = Scale2)

    zedd1 <- (ysim[, 1] - locat1) / Scale1
    zedd2 <- (ysim[, 2] - locat2) / Scale2
    ezedd1 <- exp(-zedd1)
    ezedd2 <- exp(-zedd2)
    denom <- 1 + ezedd1 + ezedd2

    dl.dlocat1 <- (1 - 3 * ezedd1 / denom) / Scale1
    dl.dlocat2 <- (1 - 3 * ezedd2 / denom) / Scale2
    dl.dscale1 <- (zedd1 - 1 - 3 * ezedd1 * zedd1 / denom) / Scale1
    dl.dscale2 <- (zedd2 - 1 - 3 * ezedd2 * zedd2 / denom) / Scale2


      rm(ysim)
      temp3 <- cbind(dl.dlocat1,
                     dl.dscale1,
                     dl.dlocat2,
                     dl.dscale2)
      run.varcov <- run.varcov + temp3[, ind1$row] * temp3[, ind1$col]
    }  # ii
    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(dlocat1.deta,
                           dscale1.deta,
                           dlocat2.deta,
                           dscale2.deta)
    wz <- wz * dthetas.detas[, ind1$row] * dthetas.detas[, ind1$col]
    c(w) * wz
  }), list( .lscale = lscale,
            .escale = escale,
            .llocat = llocat,
            .nsimEIM = nsimEIM ))))
}






dbilogis <- function(x1, x2, loc1 = 0, scale1 = 1,
                     loc2 = 0, scale2 = 1, log = FALSE) {
  if (!is.logical(log.arg <- log) || length(log) != 1)
    stop("bad input for argument 'log'")
  rm(log)




  L <- max(length(x1), length(x2),
           length(loc1), length(loc2),
           length(scale1), length(scale2))
  if (length(x1    ) != L) x1     <- rep_len(x1,     L)
  if (length(x2    ) != L) x2     <- rep_len(x2,     L)
  if (length(loc1  ) != L) loc1   <- rep_len(loc1,   L)
  if (length(loc2  ) != L) loc2   <- rep_len(loc2,   L)
  if (length(scale1) != L) scale1 <- rep_len(scale1, L)
  if (length(scale2) != L) scale2 <- rep_len(scale2, L)
  zedd1 <- (x1 - loc1) / scale1
  zedd2 <- (x2 - loc2) / scale2




  logdensity <- log(2) - zedd1 - zedd2 - log(scale1) -
                log(scale1) - 3 * log1p(exp(-zedd1) + exp(-zedd2))


  logdensity[x1 == -Inf | x2 == -Inf] <- log(0)  # 20141216 KaiH


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



pbilogis <-
  function(q1, q2, loc1 = 0, scale1 = 1, loc2 = 0, scale2 = 1) {

  ans <- 1 / (1 + exp(-(q1-loc1)/scale1) + exp(-(q2-loc2)/scale2))
  ans[scale1 <= 0] <- NA
  ans[scale2 <= 0] <- NA
  ans
}



rbilogis <- function(n, loc1 = 0, scale1 = 1, loc2 = 0, scale2 = 1) {


  y1 <- rlogis(n = n, location = loc1, scale = scale1)
  ezedd1 <- exp(-(y1-loc1)/scale1)
  y2 <- loc2 - scale2 *
        log(1/sqrt(runif(n) / (1 + ezedd1)^2) - 1 - ezedd1)
  ans <- cbind(y1, y2)
  ans[scale2 <= 0, ] <- NA
  ans
}




 freund61 <-
  function(la  = "loglink",
           lap = "loglink",
           lb  = "loglink",
           lbp = "loglink",
           ia = NULL, iap = NULL, ib = NULL, ibp = NULL,
           independent = FALSE,
           zero = NULL) {
  la <- as.list(substitute(la))
  ea <- link2list(la)
  la <- attr(ea, "function.name")

  lap <- as.list(substitute(lap))
  eap <- link2list(lap)
  lap <- attr(eap, "function.name")

  lb <- as.list(substitute(lb))
  eb <- link2list(lb)
  lb <- attr(eb, "function.name")


  lbp <- as.list(substitute(lbp))
  ebp <- link2list(lbp)
  lbp <- attr(ebp, "function.name")



  new("vglmff",
  blurb = c("Freund (1961) bivariate exponential distribution\n",
            "Links:    ",
            namesof("a",  la,  earg = ea ), ", ",
            namesof("ap", lap, earg = eap), ", ",
            namesof("b",  lb,  earg = eb ), ", ",
            namesof("bp", lbp, earg = ebp)),
  constraints = eval(substitute(expression({
    M1 <- 4
    Q1 <- 2
    constraints <- cm.VGAM(matrix(c(1, 1,0,0, 0,0, 1, 1), M, 2), x = x,
                           bool = .independent ,
                           constraints = constraints,
                           apply.int = TRUE)
    constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
                                predictors.names = predictors.names,
                                M1 = 4)
  }), list( .independent = independent, .zero = zero))),



  infos = eval(substitute(function(...) {
    list(M1 = 4,
         Q1 = 2,
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = c("a", "ap", "b", "bp"),
         la    = .la  ,
         lap   = .lap ,
         lb    = .lb  ,
         lbp   = .lbp ,
         independent = .independent ,
         zero = .zero )
    }, list( .zero = zero,
             .la    = la  ,
             .lap   = lap ,
             .lb    = lb  ,
             .lbp   = lbp ,
             .independent = independent ))),


  initialize = eval(substitute(expression({

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


    predictors.names <-
      c(namesof("a",  .la  , earg = .ea  , short = TRUE),
        namesof("ap", .lap , earg = .eap , short = TRUE),
        namesof("b",  .lb  , earg = .eb  , short = TRUE),
        namesof("bp", .lbp , earg = .ebp , short = TRUE))
    extra$y1.lt.y2 = y[, 1] < y[, 2]

    if (!(arr <- sum(extra$y1.lt.y2)) || arr == n)
      stop("identifiability problem: either all y1<y2 or y2<y1")

    if (!length(etastart)) {
      sumx  <- sum(y[ extra$y1.lt.y2, 1]);
      sumxp <- sum(y[!extra$y1.lt.y2, 1])
      sumy  <- sum(y[ extra$y1.lt.y2, 2]);
      sumyp <- sum(y[!extra$y1.lt.y2, 2])

      if (FALSE) { # Noise:
        arr <- min(arr + n/10, n*0.95)
        sumx <- sumx * 1.1; sumxp <- sumxp * 1.2;
        sumy <- sumy * 1.2; sumyp <- sumyp * 1.3;
      }
      ainit  <- if (length( .ia  )) rep_len( .ia  , n) else
            arr  / (sumx  + sumyp)
      apinit <- if (length( .iap )) rep_len( .iap , n) else
         (n-arr) / (sumxp - sumyp)
      binit  <- if (length( .ib  )) rep_len( .ib  , n) else
         (n-arr) / (sumx  + sumyp)
      bpinit <- if (length( .ibp )) rep_len( .ibp , n) else
            arr  / (sumy - sumx)

      etastart <-
        cbind(theta2eta(rep_len(ainit,  n), .la  , earg = .ea  ),
              theta2eta(rep_len(apinit, n), .lap , earg = .eap ),
              theta2eta(rep_len(binit,  n), .lb  , earg = .eb  ),
              theta2eta(rep_len(bpinit, n), .lbp , earg = .ebp ))
    }
  }), list( .la = la, .lap = lap, .lb = lb, .lbp = lbp,
            .ea = ea, .eap = eap, .eb = eb, .ebp = ebp,
            .ia = ia, .iap = iap, .ib = ib, .ibp = ibp))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    NOS <- NCOL(eta) / c(M1 = 4)
    alpha  <- eta2theta(eta[, 1], .la,  earg = .ea  )
    alphap <- eta2theta(eta[, 2], .lap, earg = .eap )
    beta   <- eta2theta(eta[, 3], .lb,  earg = .eb  )
    betap  <- eta2theta(eta[, 4], .lbp, earg = .ebp )
    fv.mat <- cbind((alphap + beta) / (alphap * (alpha + beta)),
                    (alpha + betap) / (betap * (alpha + beta)))
    label.cols.y(fv.mat, colnames.y = extra$colnames.y, NOS = NOS)
  }, list( .la = la, .lap = lap, .lb = lb, .lbp = lbp,
           .ea = ea, .eap = eap, .eb = eb, .ebp = ebp ))),
  last = eval(substitute(expression({
    misc$link <-    c("a" = .la , "ap" = .lap , "b" = .lb , "bp" = .lbp )
    misc$earg <- list("a" = .ea , "ap" = .eap , "b" = .eb , "bp" = .ebp )

    misc$multipleResponses <- FALSE
  }), list( .la = la, .lap = lap, .lb = lb, .lbp = lbp,
            .ea = ea, .eap = eap, .eb = eb, .ebp = ebp ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    alpha  <- eta2theta(eta[, 1], .la,  earg = .ea  )
    alphap <- eta2theta(eta[, 2], .lap, earg = .eap )
    beta   <- eta2theta(eta[, 3], .lb,  earg = .eb  )
    betap  <- eta2theta(eta[, 4], .lbp, earg = .ebp )
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      tmp88 <- extra$y1.lt.y2
      ell1 <- log(alpha[tmp88]) + log(betap[tmp88]) -
             betap[tmp88] * y[tmp88, 2] -
             (alpha+beta-betap)[tmp88] * y[tmp88, 1]
      ell2 <- log(beta[!tmp88]) + log(alphap[!tmp88]) -
             alphap[!tmp88] * y[!tmp88, 1] -
             (alpha+beta-alphap)[!tmp88] * y[!tmp88, 2]
      all.vec <- alpha
      all.vec[ tmp88] <- ell1
      all.vec[!tmp88] <- ell2
      ll.elts <- c(w) * all.vec
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .la = la, .lap = lap, .lb = lb, .lbp = lbp,
           .ea = ea, .eap = eap, .eb = eb, .ebp = ebp ))),
  vfamily = c("freund61"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    alpha  <- eta2theta(eta[, 1], .la,  earg = .ea  )
    alphap <- eta2theta(eta[, 2], .lap, earg = .eap )
    beta   <- eta2theta(eta[, 3], .lb,  earg = .eb  )
    betap  <- eta2theta(eta[, 4], .lbp, earg = .ebp )
    okay1 <- all(is.finite(alpha )) && all(0 < alpha ) &&
             all(is.finite(alphap)) && all(0 < alphap) &&
             all(is.finite(beta  )) && all(0 < beta  ) &&
             all(is.finite(betap )) && all(0 < betap )
    okay1
  }, list( .la = la, .lap = lap, .lb = lb, .lbp = lbp,
           .ea = ea, .eap = eap, .eb = eb, .ebp = ebp ))),
  deriv = eval(substitute(expression({
    tmp88  <- extra$y1.lt.y2
    alpha  <- eta2theta(eta[, 1], .la,  earg = .ea  )
    alphap <- eta2theta(eta[, 2], .lap, earg = .eap )
    beta   <- eta2theta(eta[, 3], .lb,  earg = .eb  )
    betap  <- eta2theta(eta[, 4], .lbp, earg = .ebp )

    dalpha.deta  <- dtheta.deta(alpha,  .la,  earg = .ea  )
    dalphap.deta <- dtheta.deta(alphap, .lap, earg = .eap )
    dbeta.deta   <- dtheta.deta(beta,   .lb,  earg = .eb  )
    dbetap.deta  <- dtheta.deta(betap,  .lbp, earg = .ebp )

    d1 <- 1/alpha - y[, 1]
    d1[!tmp88] <- -y[!tmp88, 2]
    d2 <- 0 * alphap
    d2[!tmp88] <- 1/alphap[!tmp88] - y[!tmp88, 1] + y[!tmp88, 2]
    d3 <- -y[, 1]
    d3[!tmp88] <- 1/beta[!tmp88] - y[!tmp88, 2]
    d4 <- 1/betap - y[, 2] + y[, 1]
    d4[!tmp88] <- 0

    c(w) * cbind(d1 * dalpha.deta,
                 d2 * dalphap.deta,
                 d3 * dbeta.deta,
                 d4 * dbetap.deta)
  }), list( .la = la, .lap = lap, .lb = lb, .lbp = lbp,
            .ea = ea, .eap = eap, .eb = eb, .ebp = ebp ))),
  weight = eval(substitute(expression({
    py1.lt.y2 <- alpha / (alpha+beta)
    d11 <- py1.lt.y2 / alpha^2
    d22 <- (1-py1.lt.y2) / alphap^2
    d33 <- (1-py1.lt.y2) / beta^2
    d44 <- py1.lt.y2 / betap^2

    wz <- matrix(0, n, M)  # diagonal
    wz[, iam(1, 1, M)] <- dalpha.deta^2  * d11
    wz[, iam(2, 2, M)] <- dalphap.deta^2 * d22
    wz[, iam(3, 3, M)] <- dbeta.deta^2   * d33
    wz[, iam(4, 4, M)] <- dbetap.deta^2  * d44

    c(w) * wz
  }), list( .la = la, .lap = lap, .lb = lb, .lbp = lbp,
            .ea = ea, .eap = eap, .eb = eb, .ebp = ebp ))))
}








 bigamma.mckay <- function(lscale = "loglink",
                           lshape1 = "loglink",
                           lshape2 = "loglink",
                           iscale = NULL,
                           ishape1 = NULL,
                           ishape2 = NULL,
                           imethod = 1,
                           zero = "shape") {
  lscale <- as.list(substitute(lscale))
  escale <- link2list(lscale)
  lscale <- attr(escale, "function.name")

  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.null(iscale))
    if (!is.Numeric(iscale, positive = TRUE))
      stop("argument 'iscale' must be positive or NULL")
  if (!is.null(ishape1))
    if (!is.Numeric(ishape1, positive = TRUE))
      stop("argument 'ishape1' must be positive or NULL")
  if (!is.null(ishape2))
    if (!is.Numeric(ishape2, positive = TRUE))
      stop("argument 'ishape2' must be positive or NULL")

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



  new("vglmff",
  blurb = c("Bivariate gamma: McKay's distribution\n",
            "Links:    ",
            namesof("scale",  lscale,  earg = escale ), ", ",
            namesof("shape1", lshape1, earg = eshape1), ", ",
            namesof("shape2", lshape2, earg = eshape2)),
  constraints = eval(substitute(expression({
    constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
                                predictors.names = predictors.names,
                                M1 = 3)
  }), list( .zero = zero ))),



  infos = eval(substitute(function(...) {
    list(M1 = 3,
         Q1 = 2,
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = c("scale", "shape1", "shape2"),
         lscale  = .lscale  ,
         lshape1 = .lshape1 ,
         lshape2 = .lshape2 ,
         zero = .zero )
    }, list( .zero = zero,
             .lscale  = lscale ,
             .lshape1 = lshape1,
             .lshape2 = lshape2 ))),


  initialize = eval(substitute(expression({

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


    extra$colnames.y  <- colnames(y)

    if (any(y[, 1] >= y[, 2]))
      stop("the second column minus the first column must be a vector ",
           "of positive values")


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

    if (!length(etastart)) {
      momentsY <- if ( .imethod == 1) {
        cbind(median(y[, 1]),  # This may not be monotonic
              median(y[, 2])) + 0.01
      } else {
        cbind(weighted.mean(y[, 1], w),
              weighted.mean(y[, 2], w))
      }

      mcg2.loglik <- function(thetaval, y, x, w, extraargs) {
        ainit <- a <- thetaval
        momentsY <- extraargs$momentsY
          p <- (1/a) * abs(momentsY[1]) + 0.01
          q <- (1/a) * abs(momentsY[2] - momentsY[1]) + 0.01
            sum(c(w) * (-(p+q)*log(a) - lgamma(p) - lgamma(q) +
                 (p - 1)*log(y[, 1]) +
                 (q - 1)*log(y[, 2]-y[, 1]) - y[, 2] / a ))
      }

      a.grid <- if (length( .iscale )) c( .iscale ) else
         c(0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50, 100)
      extraargs <- list(momentsY = momentsY)
      ainit <- grid.search(a.grid, objfun = mcg2.loglik,
                           y = y, x = x, w = w, maximize = TRUE,
                           extraargs = extraargs)
      ainit <- rep_len(if (is.Numeric( .iscale )) .iscale else ainit, n)
      pinit <- (1/ainit) * abs(momentsY[1]) + 0.01
      qinit <- (1/ainit) * abs(momentsY[2] - momentsY[1]) + 0.01

      pinit <- rep_len(if (is.Numeric( .ishape1 )) .ishape1 else pinit, n)
      qinit <- rep_len(if (is.Numeric( .ishape2 )) .ishape2 else qinit, n)

      etastart <-
        cbind(theta2eta(ainit, .lscale),
              theta2eta(pinit, .lshape1),
              theta2eta(qinit, .lshape2))
    }
  }), list( .lscale = lscale, .lshape1 = lshape1, .lshape2 = lshape2,
            .escale = escale, .eshape1 = eshape1, .eshape2 = eshape2,
            .iscale = iscale, .ishape1 = ishape1, .ishape2 = ishape2,
            .imethod = imethod ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    NOS <- NCOL(eta) / c(M1 = 3)
    a <- eta2theta(eta[, 1], .lscale  ,  .escale )
    p <- eta2theta(eta[, 2], .lshape1 , .eshape1 )
    q <- eta2theta(eta[, 3], .lshape2 , .eshape2 )
    fv.mat <-  cbind("y1" = p*a,
                     "y2" = (p+q)*a)  # Overwrite the colnames:
    label.cols.y(fv.mat, colnames.y = extra$colnames.y, NOS = NOS)
  }, list( .lscale = lscale, .lshape1 = lshape1, .lshape2 = lshape2,
           .escale = escale, .eshape1 = eshape1, .eshape2 = eshape2 ))),
  last = eval(substitute(expression({
    misc$link <-    c("scale"  = .lscale ,
                      "shape1" = .lshape1 ,
                      "shape2" = .lshape2 )

    misc$earg <- list("scale"  = .escale ,
                      "shape1" = .eshape1 ,
                      "shape2" = .eshape2 )

    misc$ishape1 <- .ishape1
    misc$ishape2 <- .ishape2
    misc$iscale <- .iscale
    misc$expected <- TRUE
    misc$multipleResponses <- FALSE
  }), list( .lscale = lscale, .lshape1 = lshape1, .lshape2 = lshape2,
            .escale = escale, .eshape1 = eshape1, .eshape2 = eshape2,
            .iscale = iscale, .ishape1 = ishape1, .ishape2 = ishape2,
            .imethod = imethod ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    a <- eta2theta(eta[, 1], .lscale  ,  .escale )
    p <- eta2theta(eta[, 2], .lshape1 , .eshape1 )
    q <- eta2theta(eta[, 3], .lshape2 , .eshape2 )

    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <-
        c(w) * (-(p+q)*log(a) - lgamma(p) - lgamma(q) +
               (p - 1)*log(y[, 1]) + (q - 1)*log(y[, 2]-y[, 1]) -
               y[, 2] / a)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .lscale = lscale, .lshape1 = lshape1, .lshape2 = lshape2,
           .escale = escale, .eshape1 = eshape1, .eshape2 = eshape2 ))),
  vfamily = c("bigamma.mckay"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    aparam <- eta2theta(eta[, 1], .lscale  ,  .escale )
    shape1 <- eta2theta(eta[, 2], .lshape1 , .eshape1 )
    shape2 <- eta2theta(eta[, 3], .lshape2 , .eshape2 )
    okay1 <- all(is.finite(aparam)) && all(0 < aparam) &&
             all(is.finite(shape1)) && all(0 < shape1) &&
             all(is.finite(shape2)) && all(0 < shape2)
    okay1
  }, list( .lscale = lscale, .lshape1 = lshape1, .lshape2 = lshape2,
           .escale = escale, .eshape1 = eshape1, .eshape2 = eshape2 ))),
  deriv = eval(substitute(expression({
    aparam <- eta2theta(eta[, 1], .lscale  ,  .escale )
    shape1 <- eta2theta(eta[, 2], .lshape1 , .eshape1 )
    shape2 <- eta2theta(eta[, 3], .lshape2 , .eshape2 )

    dl.da <- (-(shape1+shape2) + y[, 2] / aparam) / aparam
    dl.dshape1 <- -log(aparam) - digamma(shape1) + log(y[, 1])
    dl.dshape2 <- -log(aparam) - digamma(shape2) + log(y[, 2]-y[, 1])

    c(w) * cbind(dl.da      * dtheta.deta(aparam, .lscale),
                 dl.dshape1 * dtheta.deta(shape1, .lshape1),
                 dl.dshape2 * dtheta.deta(shape2, .lshape2))
  }), list( .lscale = lscale, .lshape1 = lshape1, .lshape2 = lshape2,
            .escale = escale, .eshape1 = eshape1, .eshape2 = eshape2 ))),
  weight = eval(substitute(expression({
    d11 <- (shape1+shape2) / aparam^2
    d22 <- trigamma(shape1)
    d33 <- trigamma(shape2)
    d12 <- 1 / aparam
    d13 <- 1 / aparam
    d23 <- 0

    wz <- matrix(0, n, dimm(M))
    wz[, iam(1, 1, M)] <- dtheta.deta(aparam, .lscale  )^2 * d11
    wz[, iam(2, 2, M)] <- dtheta.deta(shape1, .lshape1 )^2 * d22
    wz[, iam(3, 3, M)] <- dtheta.deta(shape2, .lshape2 )^2 * d33
    wz[, iam(1, 2, M)] <- dtheta.deta(aparam, .lscale  ) *
                          dtheta.deta(shape1, .lshape1 ) * d12
    wz[, iam(1, 3, M)] <- dtheta.deta(aparam, .lscale  ) *
                          dtheta.deta(shape2, .lshape2 ) * d13
    wz[, iam(2, 3, M)] <- dtheta.deta(shape1, .lshape1 ) *
                          dtheta.deta(shape2, .lshape2 ) * d23

    c(w) * wz
  }), list( .lscale = lscale, .lshape1 = lshape1,
                              .lshape2 = lshape2 ))))
}











rbifrankcop <- function(n, apar) {
  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
  if (!is.Numeric(apar, positive = TRUE))
    stop("bad input for argument 'apar'")
  if (length(apar) != use.n)
    apar <- rep_len(apar, use.n)
  U <- runif(use.n)
  V <- runif(use.n)

  T <- apar^U + (apar - apar^U) * V
  X <- U
  index <- (abs(apar - 1) < .Machine$double.eps)
  Y <- U
  if (any(!index))
    Y[!index] <- logb(T[!index] / (T[!index] +
                      (1 - apar[!index]) * V[!index]),
                      base = apar[!index])
  ans <- matrix(c(X, Y), nrow = use.n, ncol = 2)
  if (any(index)) {
    ans[index, 1] <- runif(sum(index))  # Uniform density for apar == 1
    ans[index, 2] <- runif(sum(index))
  }
  ans
}


pbifrankcop <- function(q1, q2, apar) {
  if (!is.Numeric(q1))                     stop("bad input for 'q1'")
  if (!is.Numeric(q2))                     stop("bad input for 'q2'")
  if (!is.Numeric(apar, positive = TRUE)) stop("bad input for 'apar'")

  L <- max(length(q1), length(q2), length(apar))
  if (length(apar ) != L) apar  <- rep_len(apar, L)
  if (length(q1   ) != L) q1    <- rep_len(q1,   L)
  if (length(q2   ) != L) q2    <- rep_len(q2,   L)

  x <- q1; y <- q2
  index <- (x >= 1 & y <  1) | (y >= 1 & x <  1) |
           (x <= 0 | y <= 0) | (x >= 1 & y >= 1) |
           (abs(apar - 1) < .Machine$double.eps)
  ans <- as.numeric(index)
  if (any(!index))
  ans[!index] <- logb(1 + ((apar[!index])^(x[!index]) - 1)*
                 ((apar[!index])^(y[!index]) - 1)/(apar[!index] - 1),
                 base = apar[!index])
  ind2 <- (abs(apar - 1) < .Machine$double.eps)
  ans[ind2] <- x[ind2] * y[ind2]
  ans[x >= 1 & y <  1] <- y[x >= 1 & y < 1]  # P(Y2 < q2) = q2
  ans[y >= 1 & x <  1] <- x[y >= 1 & x < 1]  # P(Y1 < q1) = q1
  ans[x <= 0 | y <= 0] <- 0
  ans[x >= 1 & y >= 1] <- 1
  ans
}





if (FALSE)
dbifrank <- function(x1, x2, apar, log = FALSE) {
  if (!is.logical(log.arg <- log) || length(log) != 1)
    stop("bad input for argument 'log'")
  rm(log)
  logdens <- (x1+x2)*log(apar) + log(apar-1) + log(log(apar)) -
             2 * log(apar - 1 + (apar^x1 - 1) * (apar^x2 - 1))

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




dbifrankcop <- function(x1, x2, apar, log = FALSE) {
  if (!is.logical(log.arg <- log) || length(log) != 1)
    stop("bad input for argument 'log'")
  rm(log)


  if (!is.Numeric(x1))                     stop("bad input for 'x1'")
  if (!is.Numeric(x2))                     stop("bad input for 'x2'")
  if (!is.Numeric(apar, positive = TRUE)) stop("bad input for 'apar'")

  L <- max(length(x1), length(x2), length(apar))
  if (length(apar ) != L) apar  <- rep_len(apar, L)
  if (length(x1   ) != L) x1    <- rep_len(x1,   L)
  if (length(x2   ) != L) x2    <- rep_len(x2,   L)

  if (log.arg) {
    denom <- apar-1 + (apar^x1  - 1) * (apar^x2  - 1)
    denom <- abs(denom)
    log((apar - 1) * log(apar)) + (x1+x2)*log(apar) - 2 * log(denom)
  } else {
    temp <- (apar - 1) + (apar^x1 - 1) * (apar^x2 - 1)
    index <- (abs(apar - 1) < .Machine$double.eps)
    ans <- x1
    if (any(!index))
      ans[!index] <- (apar[!index] - 1) * log(apar[!index]) *
                     (apar[!index])^(x1[!index] +
                                     x2[!index]) / (temp[!index])^2
    ans[x1 <= 0 | x2 <= 0 | x1 >= 1 | x2 >= 1] <- 0
    ans[index] <- 1
    ans
  }
}




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






 bifrankcop <- function(lapar = "loglink", iapar = 2, nsimEIM = 250) {

  lapar <- as.list(substitute(lapar))
  eapar <- link2list(lapar)
  lapar <- attr(eapar, "function.name")


  if (!is.Numeric(iapar, positive = TRUE))
    stop("argument 'iapar' must be positive")


  if (length(nsimEIM) &&
     (!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("Frank's bivariate copula\n",
            "Links:    ",
            namesof("apar", lapar, earg = eapar )),
  initialize = eval(substitute(expression({

    if (any(y <= 0) || any(y >= 1))
      stop("the response must have values between 0 and 1")

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


    predictors.names <-
      c(namesof("apar", .lapar , earg = .eapar, short = TRUE))

    extra$colnames.y  <- colnames(y)

    if (!length(etastart)) {
      apar.init <- rep_len(.iapar , n)
      etastart <- cbind(theta2eta(apar.init, .lapar , earg = .eapar ))
    }
  }), list( .lapar = lapar, .eapar = eapar, .iapar = iapar))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    NOS <- NCOL(eta) / c(M1 = 1)
    Q1 <- 2
    fv.mat <- matrix(0.5, NROW(eta), NOS * Q1)
    label.cols.y(fv.mat, colnames.y = extra$colnames.y, NOS = NOS)
  }, list( .lapar = lapar, .eapar = eapar ))),
  last = eval(substitute(expression({
    misc$link <-    c("apar" = .lapar )

    misc$earg <- list("apar" = .eapar )

    misc$expected <- TRUE
    misc$nsimEIM <- .nsimEIM
    misc$pooled.weight <- pooled.weight
    misc$multipleResponses <- FALSE
  }), list( .lapar = lapar, .eapar = eapar, .nsimEIM = nsimEIM ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    apar <- eta2theta(eta, .lapar , earg = .eapar )
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * dbifrankcop(x1 = y[, 1], x2 = y[, 2],
                                    apar = apar, log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .lapar = lapar, .eapar = eapar ))),
  vfamily = c("bifrankcop"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    apar <- eta2theta(eta, .lapar , earg = .eapar )
    okay1 <- all(is.finite(apar)) && all(0 < apar)
    okay1
  }, list( .lapar = lapar, .eapar = eapar ))),



  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)
    apar <- eta2theta(eta, .lapar , earg = .eapar )
    rbifrankcop(nsim * length(apar), apar = c(apar))
  }, list( .lapar = lapar, .eapar = eapar ))),




  deriv = eval(substitute(expression({
    apar <- eta2theta(eta, .lapar , earg = .eapar )
    dapar.deta <- dtheta.deta(apar, .lapar , earg = .eapar )

    de3 <- deriv3(~ (log((apar - 1) * log(apar)) + (y1+y2)*log(apar) -
                      2 * log(apar-1 + (apar^y1  - 1) * (apar^y2  - 1))),
                    name = "apar", hessian = TRUE)

    denom <- apar-1 + (apar^y[, 1]  - 1) * (apar^y[, 2]  - 1)
    tmp700 <- 2*apar^(y[, 1]+y[, 2]) - apar^y[, 1] - apar^y[, 2]
    numerator <- 1 + y[, 1] * apar^(y[, 1] - 1) * (apar^y[, 2]  - 1) +
                     y[, 2] * apar^(y[, 2] - 1) * (apar^y[, 1]  - 1)
    Dl.dapar <- 1/(apar - 1) + 1/(apar*log(apar)) +
                (y[, 1]+y[, 2])/apar - 2 * numerator / denom
    c(w) * Dl.dapar * dapar.deta
  }), list( .lapar = lapar,
            .eapar = eapar, .nsimEIM = nsimEIM ))),
  weight = eval(substitute(expression({
  if ( is.Numeric( .nsimEIM)) {

    pooled.weight <- FALSE  # For @last


    run.mean <- 0
    for (ii in 1:( .nsimEIM )) {
      ysim <- rbifrankcop(n, apar = apar)
        y1 <- ysim[, 1]; y2 <- ysim[, 2];
        eval.de3 <- eval(de3)
        d2l.dthetas2 <-  attr(eval.de3, "hessian")
        rm(ysim)
        temp3 <- -d2l.dthetas2[, 1, 1]   # M = 1
        run.mean <- ((ii - 1) * run.mean + temp3) / ii
    }
    wz <- if (intercept.only)
        matrix(mean(run.mean), n, dimm(M)) else run.mean

    wz <- wz * dapar.deta^2
    c(w) * wz
  } else {
      nump <- apar^(y[, 1]+y[, 2]-2) * (2 * y[, 1] * y[, 2] +
                    y[, 1]*(y[, 1] - 1) + y[, 2]*(y[, 2] - 1)) -
                    y[, 1]*(y[, 1] - 1) * apar^(y[, 1]-2) -
                    y[, 2]*(y[, 2] - 1) * apar^(y[, 2]-2)
      D2l.dapar2 <- 1/(apar - 1)^2 + (1+log(apar))/(apar*log(apar))^2 +
                    (y[, 1]+y[, 2])/apar^2 + 2 *
                    (nump / denom - (numerator/denom)^2)
      d2apar.deta2 <- d2theta.deta2(apar, .lapar , earg = .eapar )
      wz <- c(w) * (dapar.deta^2 * D2l.dapar2 - Dl.dapar * d2apar.deta2)
      if (TRUE && intercept.only) {
        wz <- cbind(wz)
        sumw <- sum(w)
        for (iii in 1:ncol(wz))
          wz[,iii] <- sum(wz[, iii]) / sumw
        pooled.weight <- TRUE
        wz <- c(w) * wz   # Put back the weights
      } else {
        pooled.weight <- FALSE
      }
    wz
  }
  }), list( .lapar = lapar,
            .eapar = eapar, .nsimEIM = nsimEIM ))))
}





 gammahyperbola <-
  function(ltheta = "loglink", itheta = NULL, expected = FALSE) {

  ltheta <- as.list(substitute(ltheta))
  etheta <- link2list(ltheta)
  ltheta <- attr(etheta, "function.name")

  if (!is.logical(expected) || length(expected) != 1)
      stop("argument 'expected' must be a single logical")


  new("vglmff",
  blurb = c("Gamma hyperbola bivariate distribution\n",
            "Links:    ",
            namesof("theta", ltheta, etheta)),
  initialize = eval(substitute(expression({
    if (any(y[, 1] <= 0) || any(y[, 2] <= 1))
      stop("the response has values that are out of range")

    temp5 <-
    w.y.check(w = w, y = y,
              Is.positive.y = TRUE,
              ncol.w.max