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 gammaff.mm dgamma.mm 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 bigumbelIexp bilogistic binormal binormalcop binormalcop biplackettcop bistudentt dbiamhcop dbiclaytoncop dbiclaytoncop dbifgmcop dbifrankcop dbilogis dbinorm dbinorm dbinormcop dbinormcop dbiplackcop dbistudentt dgamma.mm dtrinorm freund61 gammaff.mm 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-2024 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)
}  # dtrinorm



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) <- paste0("X", 1:3)
  ans
}  # rtrinorm





trinormal.control <-
  function(summary.HDEtest = FALSE,
           ...) {
  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 dim nx3
    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)
}  # dbiclaytoncop



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)
}  # rbiclaytoncop



 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 +ve
          } 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)
}  # dbistudentt




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, NOS = NOS,
                 colnames.y = extra$colnames.y)
  }  ,
  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 ))))
}  # bistudentt








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)
}  # dbinormcop






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(qnorm(q1), qnorm(q2), cov12 = rho)
}  # pbinormcop



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

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






 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 ))))
}  # binormalcop















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 ))))
}  # bilogistic






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)
}  # dbilogis



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
}  # pbilogis



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
}  # rbilogis






 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 ))))
}  # freund61









 dgamma.mm <-
     function(x, shape, scale = 1, log = FALSE,
              sh.byrow = TRUE) {
  if (!is.matrix(shape))
    shape <- as.matrix(shape)
  if (!is.matrix(x))
    x <- as.matrix(x)
  if (ncol(x) < 2)
    stop("argument 'x' must have at least two columns")
  if (!is.logical(log.arg <- log) || length(log) != 1)
    stop("bad input for argument 'log'")
  rm(log)

  n <- max(nrow(x), length(scale), nrow(shape))
  Q <- max(ncol(x), ncol(shape))
  x <- matrix(c(x), n, Q)   # , byrow = xs.byrow
  shape <- matrix(c(shape), n, Q, byrow = sh.byrow)
  scale <- rep_len(scale, n)
  x.Q <- x[, ncol(x)]  # Last coln of x
  x <- x - cbind(0, x[, -ncol(x), drop = FALSE])  # n x Q
  logdensity <-
    rowSums((shape - 1) * log(x)) - x.Q / scale -
    rowSums(lgamma(shape)) -
    rowSums(shape) * log(scale)

  if (log.arg) logdensity else exp(logdensity)
}  # dgamma.mm










 gammaff.mm <-
    function(lscale = "loglink",
             lshape = "loglink",
             iscale = NULL,    # 1,   # NULL,
             ishape = NULL,    # 2.1,   #NULL,
             imethod = 1,
             eq.shapes = FALSE,
             sh.byrow = TRUE,
             zero = "shape") {
  lscale <- as.list(substitute(lscale))
  escale <- link2list(lscale)
  lscale <- attr(escale, "function.name")

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

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

  if (!is.null(iscale))
    if (!is.Numeric(iscale, positive = TRUE))
      stop("argument 'iscale' must be positive or NULL")
  if (!is.null(ishape))
    if (!is.Numeric(ishape, positive = TRUE))
      stop("argument 'ishape' must be positive or NULL")

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



  new("vglmff",
  blurb = c("Multivariate gamma distribution: ",
            "Mathai and Moschopoulos (1992)\n",
            "Links:    ",
        namesof("scale",    lscale, earg = escale ), ", ",
        namesof("shape1",   lshape, earg = eshape), ", ..., ",
        namesof("shapeM-1", lshape, earg = eshape), "\n\n"),
  constraints = eval(substitute(expression({
    constraints.orig <- constraints
    Msub1 <- M - 1  # >= 2, aka Q
    NOS <- 1

    cmk.s <- cbind(c(1, rep_len(0, Msub1)),
                   c(0, rep_len(1, Msub1))) 
    con.s <- cm.VGAM(cmk.s,
                     x = x,
                     bool = .eq.shapes ,
                     constraints = constraints.orig,
                     apply.int = TRUE,
                     cm.default           = diag(M),
                     cm.intercept.default = diag(M))
    constraints <- con.s
    constraints <- cm.zero.VGAM(constraints, x = x,
                                .zero , M = M,
                      predictors.names = predictors.names,
                                M1 = NA)
  }), list( .zero = zero, .eq.shapes = eq.shapes ))),

  infos = eval(substitute(function(...) {
    list(M1 = NA,
         Q1 = NA,   # >= 2 is required
         eq.shapes = .eq.shapes ,
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = c("scale", "shape"),
         lscale = .lscale ,
         lshape = .lshape ,
         zero = .zero )
  },
  list( .zero = zero, .eq.shapes = eq.shapes ,
        .lscale = lscale , .lshape = lshape ))),


  initialize = eval(substitute(expression({
    if (NCOL(y) < 2 || !is.matrix(y))
      stop("the response must be a 2 column matrix")

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

    extra$colnames.y  <- colnames(y)

    if (any(y[, -1, drop = FALSE] -
            y[, -ncol(y), drop = FALSE] <= 0))
      stop("each row must be strictly increasing ",
           "from the first column to the last")

    mynames1 <- param.names("shape", Msub1, skip1 = TRUE)
    predictors.names <-
        c(namesof("scale",  .lscale , .escale , short = TRUE),
          namesof(mynames1, .lshape , .eshape , short = TRUE))
    if (!length(etastart)) {
      if ( .imethod == 1 && (
           length( .iscale ) == 0 ||
           length( .ishape ) == 0)) {
        all.mean <- colMeans(y)
        all.vars <- apply(y, 2, var)
        sc.mme <- tail(all.vars / all.mean, 1)  # Last one
        sh.mme <-  all.mean *
          (all.mean - c(0, head(all.mean, -1))) / all.vars
        sc.mme <- as.vector(sc.mme)
        sh.mme <- as.vector(sh.mme)
      }

      use.iscale <-
        if (is.numeric( .iscale )) .iscale else sc.mme
      use.ishape <-
        if (is.numeric( .ishape )) .ishape else sh.mme


      etastart <-
        cbind(theta2eta(use.iscale , .lscale , .escale ),
              theta2eta(matrix(use.ishape , n, Msub1,
                               byrow = .sh.byrow ),
                        .lshape , earg = .eshape ))
    }
  }),
  list( .lscale = lscale, .lshape = lshape,
        .escale = escale, .eshape = eshape,
        .iscale = iscale, .ishape = ishape,
        .imethod = imethod, .sh.byrow = sh.byrow ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    FF <- FALSE
    NOS <- 1  # NCOL(eta) / c(M1 = 3)
    SC <- eta2theta(eta[,  1],            .lscale , .escale )
    SH <- eta2theta(eta[, -1, drop = FF], .lshape , .eshape )
    fv.mat <-  matrix(SC * SH[, 1], NROW(eta), ncol(eta) - 1)
    colnames(fv.mat) <- paste0("y", 1:ncol(fv.mat))
    for (jay in 2:ncol(fv.mat)) {
      SH[, jay] <- SH[, jay] + SH[, jay - 1]
      fv.mat[, jay] <- SC * SH[, jay]
    }
      
    label.cols.y(fv.mat, colnames.y = extra$colnames.y,
                 NOS = NOS)
  },
  list( .lscale = lscale, .lshape = lshape, 
        .escale = escale, .eshape = eshape ))),
  last = eval(substitute(expression({
    misc$link <- c( .lscale , rep_len( .lshape , Msub1))
    names(misc$link) <- c("scale", mynames1)
    misc$earg <- vector("list", M)
    names(misc$earg) <- names(misc$link)
    misc$earg[[1]] <- ( .escale )
    for (ii in 1:Msub1)
      misc$earg[[ii + 1]] <- ( .eshape )

    misc$iscale <- ( .iscale )
    misc$ishape <- ( .ishape )
    misc$expected <- TRUE
    misc$multipleResponses <- FALSE
  }),
  list( .lscale = lscale, .lshape = lshape,
        .escale = escale, .eshape = eshape,
        .iscale = iscale, .ishape = ishape,
        .imethod = imethod ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    FF <- FALSE
    SC <- eta2theta(eta[,  1],            .lscale , .escale )
    SH <- eta2theta(eta[, -1, drop = FF], .lshape , .eshape )

    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * 
        dgamma.mm(y, shape = SH, scale = SC, log = TRUE,
                  sh.byrow = FALSE )  # Internally
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
    },
  list( .lscale = lscale, .lshape = lshape,
        .escale = escale, .eshape = eshape ))),
  vfamily = c("gammaff.mm"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    FF <- FALSE
    SC <- eta2theta(eta[,  1],            .lscale , .escale )
    SH <- eta2theta(eta[, -1, drop = FF], .lshape , .eshape )
    okay1 <- all(is.finite(SC)) && all(0 < SC) &&
             all(is.finite(SH)) && all(0 < SH)
    okay1
  },
  list( .lscale = lscale, .lshape = lshape,
        .escale = escale, .eshape = eshape ))),
  deriv = eval(substitute(expression({
    FF <- FALSE
    SC <- eta2theta(eta[,  1],            .lscale , .escale )
    SH <- eta2theta(eta[, -1, drop = FF], .lshape , .eshape )
    colnames(SH) <- NULL  # More tidy
    sumshapes <- rowSums(SH)
    dSC.deta <- dtheta.deta(SC, .lscale , earg = .escale )
    dSH.deta <- dtheta.deta(SH, .lshape , earg = .eshape )

    dl.dscale  <- (y[, ncol(y)] / SC  - sumshapes) / SC
    dl.dshapes <- -digamma(SH) - log(SC)
    dl.dshapes[,  1] <- dl.dshapes[,  1] + log(y[, 1])
    dl.dshapes[, -1] <- dl.dshapes[, -1] +
                        log(y[, -1] - y[, -ncol(y)])

    c(w) * cbind(dl.dscale  * dSC.deta,
                 dl.dshapes * dSH.deta)
  }),
  list( .lscale = lscale, .lshape = lshape,
        .escale = escale, .eshape = eshape ))),
  weight = eval(substitute(expression({
    wz <- matrix(0, n, dimm(M))  # Most elts are 0 identically
    wz[, iam(1, 1, M)] <- (sumshapes / SC^2) * dSC.deta^2
    for (jay in 1:(M - 1))  # Diagonals
      wz[, iam(1 + jay, 1 + jay, M)] <-
        trigamma(SH[, jay]) * (dSH.deta[, jay])^2
    for (jay in 1:(M - 1))  # Sides
      wz[, iam(1, 1 + jay, M)] <-
        dSH.deta[, jay] * dSC.deta / SC


    c(w) * wz
  }),
  list( .lscale = lscale, .lshape = lshape ))))
}  # gammaff.mm





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 pdf 4 apar == 1
    ans[index, 2] <- runif(sum(index))
  }
  ans
}  # rbifrankcop



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
}  # pbifrankcop





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)
}  # dbifrank




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
  }
}  # dbifrankcop




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 , .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 , .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 ))))
}  # bifrankcop







 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 = 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("theta", .ltheta , .etheta , short = TRUE))

    if (!length(etastart)) {
      theta.init <- if (length( .itheta)) {
        rep_len( .itheta , n)
      } else {
        1 / (y[, 2] - 1 + 0.01)
      }
      etastart <-
        cbind(theta2eta(theta.init, .ltheta , .etheta ))
    }
  }),
  list(
      .ltheta = ltheta, .etheta = etheta, .itheta = itheta
  ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    NOS <- NCOL(eta) / c(M1 = 1)
    theta <- eta2theta(eta, .ltheta , .etheta )
    fv.mat <- cbind(theta * exp(theta), 1 + 1 / theta)
    label.cols.y(fv.mat, colnames.y = extra$colnames.y,
                 NOS = NOS)
  }, list( .ltheta = ltheta, .etheta = etheta ))),
  last = eval(substitute(expression({
    misc$link <-    c("theta" = .ltheta )

    misc$earg <- list("theta" = .etheta )

    misc$expected <- .expected
    misc$multipleResponses <- FALSE
  }), list( .ltheta = ltheta,
            .etheta = etheta, .expected = expected ))),

  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    theta <- eta2theta(eta, .ltheta , .etheta )
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <- c(w) * (-exp(-theta) * y[, 1] / theta -
                         theta * y[, 2])
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .ltheta = ltheta, .etheta = etheta ))),
  vfamily = c("gammahyperbola"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    theta <- eta2theta(eta, .ltheta , .etheta )
    okay1 <- all(is.finite(theta)) && all(0 < theta)
    okay1
  }, list( .ltheta = ltheta, .etheta = etheta ))),
  deriv = eval(substitute(expression({
    theta <- eta2theta(eta, .ltheta , .etheta )
    Dl.dtheta <- exp(-theta) * y[, 1] *
        (1 + theta) / theta^2 - y[, 2]
    DTHETA.deta <- dtheta.deta(theta, .ltheta , .etheta )
    c(w) * Dl.dtheta * DTHETA.deta
  }), list( .ltheta = ltheta, .etheta = etheta ))),
  weight = eval(substitute(expression({
    temp300 <- 2 + theta * (2 + theta)
    if ( .expected ) {
      D2l.dtheta2 <- temp300 / theta^2
      wz <- c(w) * DTHETA.deta^2 * D2l.dtheta2
    } else {
      D2l.dtheta2 <- temp300 * y[, 1] * exp(-theta) / theta^3
      D2theta.deta2 <- d2theta.deta2(theta, .ltheta )
      wz <- c(w) * (DTHETA.deta^2 * D2l.dtheta2 -
                    Dl.dtheta * D2theta.deta2)
    }
    wz
  }),
  list( .ltheta = ltheta,
        .etheta = etheta, .expected = expected ))))
}  # gammahyperbola




 bifgmexp <-
  function(lapar = "rhobitlink",
           iapar = NULL, tola0 = 0.01,
           imethod = 1) {
  lapar <- as.list(substitute(lapar))
  earg  <- link2list(lapar)
  lapar <- attr(earg, "function.name")

  if (length(iapar) &&
     (!is.Numeric(iapar, length.arg = 1) ||
      abs(iapar) >= 1))
    stop("argument 'iapar' must be a single number in (-1, 1)")

  if (!is.Numeric(tola0, length.arg = 1, positive = TRUE))
      stop("argument 'tola0' must be a single positive number")

  if (length(iapar) && abs(iapar) <= tola0)
      stop("argument 'iapar' must not be between ",
           "-tola0 and tola0")
  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 Farlie-Gumbel-Morgenstern ",
            "exponential distribution\n",  # Morgenstern's
            "Links:    ",
            namesof("apar", lapar, earg = earg )),
  initialize = eval(substitute(expression({
    temp5 <-
    w.y.check(w = w, y = y,
              Is.nonnegative.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 = .earg , short = TRUE))

    extra$colnames.y  <- colnames(y)

    if (!length(etastart)) {
      ainit  <- if (length(.iapar))
                  rep_len( .iapar , n) else {
                  mean1 <- if ( .imethod == 1)
                           median(y[, 1]) else
                           mean(y[, 1])
        mean2 <- if ( .imethod == 1)
          median(y[, 2]) else mean(y[, 2])
        Finit <- 0.01 + mean(y[, 1] <= mean1 & y[, 2] <= mean2)
        ((Finit + expm1(-mean1) +
           exp(-mean2)) / exp(-mean1 - mean2) - 1) / (
           expm1(-mean1) * expm1(-mean2))
          }
        etastart <-
          theta2eta(rep_len(ainit, n), .lapar , earg = .earg )
      }
  }), list( .iapar = iapar, .lapar = lapar, .earg = earg,
            .imethod = imethod ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    NOS <- NCOL(eta) / c(M1 = 1)
    Q1 <- 2
    fv.mat <- matrix(1, NROW(eta), NOS * Q1)
    label.cols.y(fv.mat, colnames.y = extra$colnames.y,
                 NOS = NOS)
  }, list( .lapar = lapar, .earg = earg ))),
  last = eval(substitute(expression({
    misc$link <-    c("apar" = .lapar )

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

    misc$expected <- FALSE
    misc$pooled.weight <- pooled.weight
    misc$multipleResponses <- FALSE
  }), list( .lapar = lapar, .earg = earg ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
      alpha  <- eta2theta(eta, .lapar , earg = .earg )
      alpha[abs(alpha) < .tola0 ] <- .tola0
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      denom <- (1 + alpha - 2 * alpha * (exp(-y[, 1]) +
                                         exp(-y[, 2])) +
               4 * alpha * exp(-y[, 1] - y[, 2]))
      ll.elts <- c(w) * (-y[, 1] - y[, 2] + log(denom))
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .lapar = lapar, .earg = earg, .tola0 = tola0 ))),
  vfamily = c("bifgmexp"),  # morgenstern
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    alpha <- eta2theta(eta, .lapar , earg = .earg )
    okay1 <- all(is.finite(alpha)) && all(abs(alpha) < 1)
    okay1
  }, list( .lapar = lapar, .earg = earg, .tola0 = tola0 ))),
  deriv = eval(substitute(expression({
    alpha <- eta2theta(eta, .lapar , earg = .earg )
    alpha[abs(alpha) < .tola0 ] <- .tola0
    numerator <- 1 - 2 * (exp(-y[, 1]) + exp(-y[, 2])) +
                 4 * exp(-y[, 1] - y[, 2])
    denom <- (1 + alpha - 2 * alpha * (exp(-y[, 1]) +
                                       exp(-y[, 2])) +
             4 * alpha * exp(-y[, 1] - y[, 2]))
    dl.dalpha <- numerator / denom

    dalpha.deta <- dtheta.deta(alpha,  .lapar , earg = .earg )

    c(w) * cbind(dl.dalpha * dalpha.deta)
  }), list( .lapar = lapar, .earg = earg, .tola0 = tola0 ))),
  weight = eval(substitute(expression({
    d2l.dalpha2 <- dl.dalpha^2
    d2alpha.deta2 <- d2theta.deta2(alpha, .lapar ,
                                   earg = .earg )
    wz <- c(w) * (dalpha.deta^2 * d2l.dalpha2 -
                  d2alpha.deta2 * dl.dalpha)
    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, .earg = earg ))))
}  # bifgmexp





rbifgmcop <- 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))
    stop("bad input for argument 'apar'")
  if (any(abs(apar) > 1))
    stop("argument 'apar' has values out of range")

  y1 <- V1 <- runif(use.n)
  V2 <- runif(use.n)
  temp <- 2*y1 - 1
  A <- apar * temp - 1
  B <- sqrt(1 - 2 * apar * temp + (apar*temp)^2 +
            4 * apar * V2 * temp)
  y2 <- 2 * V2 / (B - A)
  matrix(c(y1, y2), nrow = use.n, ncol = 2)
}  # rbifgmcop



dbifgmcop <- 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(apar))
    stop("bad input for 'apar'")
  if (any(abs(apar) > 1))
    stop("'apar' values out of range")
  if ( !is.logical( log.arg ) ||
       length( log.arg ) != 1 )
    stop("bad input for argument 'log'")

  L <- max(length(x1), length(x2), length(apar))
  if (length(x1)    != L)  x1   <- rep_len(x1,   L)
  if (length(x2)    != L)  x2   <- rep_len(x2,   L)
  if (length(apar)  != L)  apar <- rep_len(apar, L)
  ans <- 0 * x1
  xnok <- (x1 <= 0) | (x1 >= 1) | (x2 <= 0) | (x2 >= 1)
  if ( log.arg ) {
    ans[!xnok] <- log1p(apar[!xnok] * (1 - 2 * x1[!xnok]) *
                       (1 - 2 * x2[!xnok]))
    ans[xnok] <- log(0)
  } else {
    ans[!xnok] <- 1 + apar[!xnok] * (1 - 2 * x1[!xnok]) *
                  (1 - 2 * x2[!xnok])
    ans[xnok] <- 0
    if (any(ans < 0))
      stop("negative values in the density ",
           "(apar out of range)")
  }
  ans
}  # dbifgmcop



pbifgmcop <- 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))  stop("bad input for 'apar'")
  if (any(abs(apar) > 1)) stop("'apar' values out of range")

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

  x <- q1
  y <- q2
  index <- (x >= 1 & y <  1) |
           (y >= 1 & x <  1) |
           (x <= 0 | y <= 0) |
           (x >= 1 & y >= 1)
  ans <- as.numeric(index)
  if (any(!index)) {
    ans[!index] <-
          q1[!index]  *      q2[!index] * (1 + apar[!index] *
     (1 - q1[!index]) * (1 - q2[!index]))
  }
  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
}  # pbifgmcop






 bifgmcop <-
   function(lapar = "rhobitlink", iapar = NULL,
            imethod = 1) {

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


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

  if (length(iapar) &&
     (abs(iapar) >= 1))
    stop("'iapar' should be less than 1 in absolute value")


  new("vglmff",
  blurb = c("Farlie-Gumbel-Morgenstern copula \n",
            "Links:    ",
            namesof("apar", lapar, earg = earg )),
  initialize = eval(substitute(expression({
    if (any(y < 0) || any(y > 1))
      stop("the response must have values in the unit square")

    temp5 <-
    w.y.check(w = w, y = y,
              Is.nonnegative.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 <-
      namesof("apar", .lapar , earg = .earg , short = TRUE)

    extra$colnames.y  <- colnames(y)

    if (!length(etastart)) {
      ainit  <- if (length( .iapar ))  .iapar else {


      if ( .imethod == 1) {
        3 * cor(y[, 1], y[, 2], method = "spearman")
      } else if ( .imethod == 2) {
        9 * kendall.tau(y[, 1], y[, 2]) / 2
      } else {
        mean1 <- if ( .imethod == 1)
                 weighted.mean(y[, 1], w) else
                 median(y[, 1])
        mean2 <- if ( .imethod == 1)
                 weighted.mean(y[, 2], w) else
                 median(y[, 2])
        Finit <- weighted.mean(y[, 1] <= mean1 &
                               y[, 2] <= mean2, w)
        (Finit / (mean1 * mean2) - 1) / (
        (1 - mean1) * (1 - mean2))
      }
    }

    ainit <- min(0.95, max(ainit, -0.95))
    etastart <- theta2eta(rep_len(ainit, n), .lapar ,
                          earg = .earg )
    }
  }),
  list( .iapar = iapar, .lapar = lapar, .earg = earg,
        .imethod = imethod ))),
  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, .earg = earg ))),
  last = eval(substitute(expression({
    misc$link <-    c("apar" = .lapar )
    misc$earg <- list("apar" = .earg  )

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


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

    pwts <- if (length(pwts <- object@prior.weights) > 0)
              pwts else weights(object, type = "prior")
    if (any(pwts != 1))
      warning("ignoring prior weights")
    eta <- predict(object)
    alpha <- eta2theta(eta, .lapar , earg = .earg )
    rbifgmcop(nsim * length(alpha), apar = c(alpha))
  }, list( .lapar = lapar, .earg = earg ))),



  deriv = eval(substitute(expression({
    alpha <- eta2theta(eta, .lapar , earg = .earg )

    dalpha.deta <- dtheta.deta(alpha, .lapar , earg = .earg )

    numerator <- (1 - 2 * y[, 1])  * (1 - 2 * y[, 2])
    denom <- 1 + alpha * numerator

    mytolerance <- .Machine$double.eps
    bad <- (denom <= mytolerance)   # Range violation
    if (any(bad)) {
      cat("There are some range violations in @deriv\n")
      flush.console()
      denom[bad] <- 2 * mytolerance
    }
    dl.dalpha <- numerator / denom
    c(w) * cbind(dl.dalpha * dalpha.deta)
  }), list( .lapar = lapar, .earg = earg))),

  weight = eval(substitute(expression({
  wz <- lerch(alpha^2, 2, 1.5) / 4  # Checked and correct
  wz <- wz * dalpha.deta^2
    c(w) * wz
  }), list( .lapar = lapar, .earg = earg))))
}  # bifgmcop







 bigumbelIexp <-
  function(lapar = "identitylink", iapar = NULL, imethod = 1) {

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


  if (length(iapar) &&
      !is.Numeric(iapar, length.arg = 1))
    stop("'iapar' must be a single number")
  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("Gumbel's Type I bivariate exponential ",
            "distribution\n",
            "Links:    ",
            namesof("apar", lapar, earg = earg )),
  initialize = eval(substitute(expression({

    temp5 <-
    w.y.check(w = w, y = y,
              Is.nonnegative.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

    extra$colnames.y  <- colnames(y)



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

    if (!length(etastart)) {
      ainit  <- if (length( .iapar ))
                  rep_len( .iapar, n) else {
        mean1 <- if ( .imethod == 1)
                 median(y[, 1]) else mean(y[, 1])
        mean2 <- if ( .imethod == 1)
                 median(y[, 2]) else mean(y[, 2])
        Finit <- 0.01 + mean(y[, 1] <= mean1 &
                             y[, 2] <= mean2)
        (log(Finit + expm1(-mean1) + exp(-mean2)) +
         mean1 + mean2) / (mean1 * mean2)
      }
      etastart <-
        theta2eta(rep_len(ainit,  n), .lapar , earg = .earg )
      }
  }), list( .iapar = iapar, .lapar = lapar, .earg = earg,
            .imethod = imethod ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    NOS <- NCOL(eta) / c(M1 = 1)
    Q1 <- 2
    alpha <- eta2theta(eta, .lapar , earg = .earg )
    fv.mat <- matrix(1, NROW(eta), NOS * Q1)
    label.cols.y(fv.mat, colnames.y = extra$colnames.y,
                 NOS = NOS)
  }, list( .lapar = lapar, .earg = earg ))),
  last = eval(substitute(expression({
    misc$link <-    c("apar" = .lapar )

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

    misc$expected <- FALSE
    misc$pooled.weight <- pooled.weight
    misc$multipleResponses <- FALSE
  }), list( .lapar = lapar, .earg = earg ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    alpha  <- eta2theta(eta, .lapar , earg = .earg )
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      denom <- (alpha*y[, 1] - 1) * (alpha*y[, 2] - 1) + alpha
      mytolerance <- .Machine$double.xmin
      bad <- (denom <= mytolerance)  # Range violation
      if (any(bad)) {
        cat("There are some range violations in @deriv\n")
        flush.console()
      }




      if (summation) {
      sum(bad) * (-1.0e10) +
      sum(w[!bad] * (-y[!bad, 1] - y[!bad, 2] +
                     alpha[!bad] * y[!bad, 1] * y[!bad, 2] +
                     log(denom[!bad])))
      } else {
        stop("argument 'summation = FALSE' does not work yet")
      }
    }
  }, list( .lapar = lapar, .earg = earg ))),
  vfamily = c("bigumbelIexp"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    alpha <- eta2theta(eta, .lapar , earg = .earg )
    okay1 <- all(is.finite(alpha))
    okay1
  }, list( .lapar = lapar, .earg = earg ))),
  deriv = eval(substitute(expression({
    alpha <- eta2theta(eta, .lapar , earg = .earg )
    numerator <- (alpha * y[, 1] - 1) * y[, 2] +
                 (alpha * y[, 2] - 1) * y[, 1] + 1
    denom <- (alpha * y[, 1] - 1) * (alpha * y[, 2] - 1) + alpha
    denom <- abs(denom)

    dl.dalpha <- numerator / denom + y[, 1] * y[, 2]

    dalpha.deta <- dtheta.deta(alpha,  .lapar , earg = .earg )

    c(w) * cbind(dl.dalpha * dalpha.deta)
  }), list( .lapar = lapar, .earg = earg ))),
  weight = eval(substitute(expression({
    d2l.dalpha2 <- (numerator/denom)^2 - 2*y[, 1]*y[, 2] / denom
    d2alpha.deta2 <- d2theta.deta2(alpha, .lapar , earg = .earg )
    wz <- c(w) * (dalpha.deta^2 * d2l.dalpha2 -
                  d2alpha.deta2 * dl.dalpha)
    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, .earg = earg ))))
}  # bigumbelIexp







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

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

  x <- q1; y <- q2
  index <- (x >= 1 & y <  1) | (y >= 1 & x <  1) |
           (x <= 0 | y <= 0) | (x >= 1 & y >= 1) |
           (abs(oratio - 1) < 1.0e-6)  #  .Machine$double.eps
  ans <- as.numeric(index)
  if (any(!index)) {
    temp1 <- 1 + (oratio[!index]  - 1) *
             (q1[!index] + q2[!index])
    temp2 <- temp1 - sqrt(temp1^2 - 4 * oratio[!index] *
             (oratio[!index] - 1) * q1[!index] * q2[!index])
    ans[!index] <- 0.5 * temp2 / (oratio[!index] - 1)
  }

  ind2 <- (abs(oratio - 1) < 1.0e-6)  # .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
}



rbiplackcop <- function(n, oratio) {
  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


  y1 <- U <- runif(use.n)
  V <- runif(use.n)
  Z <- V * (1-V)
  y2 <- (2*Z*(y1*oratio^2 + 1 - y1) + oratio * (1 - 2 * Z) -
  (1 - 2 * V) *
  sqrt(oratio * (oratio + 4*Z*y1*(1-y1)*
                 (1-oratio)^2))) / (oratio +
  Z*(1-oratio)^2)
  matrix(c(y1, 0.5 * y2), nrow = use.n, ncol = 2)
}



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


  ans <- log(oratio) + log1p((oratio - 1) *
         (x1+x2 - 2*x1*x2)) - 1.5 *
         log((1 + (x1+x2)*(oratio - 1))^2 -
             4 * oratio * (oratio - 1)*x1*x2)
  ans[ # !is.na(x1) & !is.na(x2) & !is.na(oratio) &
     ((x1 < 0) | (x1 > 1) | (x2 < 0) | (x2 > 1))] <- log(0)


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



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



 biplackettcop <-
    function(link = "loglink", ioratio = NULL,
                      imethod = 1, nsimEIM = 200) {

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


  if (length(ioratio) &&
     (!is.Numeric(ioratio, positive = TRUE)))
    stop("'ioratio' must be positive")

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


  new("vglmff",
  blurb = c("Plackett distribution (bivariate copula)\n",
            "Links:    ",
            namesof("oratio", link, earg = earg )),
  initialize = eval(substitute(expression({
    if (any(y < 0) || any(y > 1))
      stop("the response must have values in the unit square")

    temp5 <-
    w.y.check(w = w, y = y,
              Is.nonnegative.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 <-
      namesof("oratio", .link , earg = .earg, short = TRUE)

    extra$colnames.y  <- colnames(y)

    if (!length(etastart)) {
      orinit <- if (length( .ioratio ))  .ioratio else {
          if ( .imethod == 2) {
            scorp <- cor(y)[1, 2]
            if (abs(scorp) <= 0.1) 1 else
            if (abs(scorp) <= 0.3) 3^sign(scorp) else
            if (abs(scorp) <= 0.6) 5^sign(scorp) else
            if (abs(scorp) <= 0.8) 20^sign(scorp) else
                                   40^sign(scorp)
          } else {
            y10 <- weighted.mean(y[, 1], w)
            y20 <- weighted.mean(y[, 2], w)
         (0.5 + sum(w[(y[, 1] <  y10) & (y[, 2] <  y20)])) *
         (0.5 + sum(w[(y[, 1] >= y10) & (y[, 2] >= y20)])) / (
        ((0.5 + sum(w[(y[, 1] <  y10) & (y[, 2] >= y20)])) *
         (0.5 + sum(w[(y[, 1] >= y10) & (y[, 2] <  y20)]))))
          }
        }
        etastart <- theta2eta(rep_len(orinit, n), .link , .earg )
    }
  }), list( .ioratio = ioratio, .link = link, .earg = earg,
            .imethod = imethod ))),
  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( .link = link, .earg = earg ))),
  last = eval(substitute(expression({
    misc$link <-    c(oratio = .link)

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

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


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

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



  deriv = eval(substitute(expression({
    oratio <- eta2theta(eta, .link , earg = .earg )
    doratio.deta <- dtheta.deta(oratio, .link , .earg )
    y1 <- y[, 1]
    y2 <- y[, 2]
    de3 <- deriv3(~ (log(oratio) + log(1+(oratio - 1) *
                 (y1+y2-2*y1*y2)) - 1.5 *
                 log((1 + (y1+y2)*(oratio - 1))^2 -
                 4 * oratio * (oratio - 1)*y1*y2)),
                 name = "oratio", hessian = FALSE)
    eval.de3 <- eval(de3)

    dl.doratio <-  attr(eval.de3, "gradient")

    c(w) * dl.doratio * doratio.deta
  }), list( .link = link, .earg = earg ))),
  weight = eval(substitute(expression({
    sd3 <- deriv3(~ (log(oratio) + log(1+(oratio - 1) *
          (y1sim+y2sim-2*y1sim*y2sim)) - 1.5 *
          log((1 + (y1sim+y2sim)*(oratio - 1))^2 -
          4 * oratio * (oratio - 1)*y1sim*y2sim)),
                    name = "oratio", hessian = FALSE)
    run.var <- 0
    for (ii in 1:( .nsimEIM )) {
      ysim <- rbiplackcop(n, oratio = oratio)
      y1sim <- ysim[, 1]
      y2sim <- ysim[, 1]
        eval.sd3 <- eval(sd3)
        dl.doratio <-  attr(eval.sd3, "gradient")
        rm(ysim, y1sim, y2sim)
        temp3 <- dl.doratio
        run.var <- ((ii - 1) * run.var + temp3^2) / ii
    }
    wz <- if (intercept.only)
        matrix(colMeans(cbind(run.var)),
               n, dimm(M), byrow = TRUE) else cbind(run.var)

    wz <- wz * doratio.deta^2
    c(w) * wz
  }),
  list( .link = link, .earg = earg, .nsimEIM = nsimEIM ))))
}  # biplackettcop







dbiamhcop <- function(x1, x2, apar, 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(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)
  temp <- 1 - apar*(1-x1)*(1-x2)

  if (log.arg) {
    ans <- log1p(-apar+2*apar*x1*x2/temp) - 2*log(temp)
    ans[(x1 <= 0) | (x1 >= 1) | (x2 <= 0) | (x2 >= 1)] <-
        log(0)
  } else {
    ans <- (1-apar+2*apar*x1*x2/temp) / (temp^2)
    ans[(x1 <= 0) | (x1 >= 1) | (x2 <= 0) | (x2 >= 1)] <- 0
  }
  ans[abs(apar) > 1] <- NA
  ans
}  # dbiamhcop



pbiamhcop <- 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)) stop("bad input for 'apar'")

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

  x <- q1
  y <- q2
  index <- (x >= 1 & y < 1) | (y >= 1 & x <  1) |
           (x <= 0 | y<= 0) | (x >= 1 & y >= 1)
  ans <- as.numeric(index)
  if (any(!index)) {
    ans[!index] <- (q1[!index] * q2[!index]) / (1 -
       apar[!index] * (1-q1[!index]) * (1-q2[!index]))
  }
  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[abs(apar) > 1] <- NA
  ans
}  # pbiamhcop


rbiamhcop <- 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 (any(abs(apar) > 1))
    stop("'apar' values out of range")

  U1 <- V1 <- runif(use.n)
  V2 <- runif(use.n)
  b <- 1-V1
  A <- -apar*(2*b*V2+1)+2*apar^2*b^2*V2+1
  B <- apar^2*(4*b^2*V2-4*b*V2+1)+apar*(4*V2-4*b*V2-2)+1
  U2 <- (2*V2*(apar*b - 1)^2)/(A+sqrt(B))
  matrix(c(U1, U2), nrow = use.n, ncol = 2)
}  # rbiamhcop


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


 biamhcop <-
    function(lapar = "rhobitlink", iapar = NULL,
             imethod = 1, nsimEIM = 250) {
  lapar <- as.list(substitute(lapar))
  eapar <- link2list(lapar)
  lapar <- attr(eapar, "function.name")



  if (length(iapar) && (abs(iapar) > 1))
    stop("'iapar' should be <= 1 in absolute value")
  if (!is.Numeric(imethod, length.arg = 1,
                  integer.valued = TRUE, positive = TRUE) ||
    imethod > 2)
    stop("imethod must be 1 or 2")

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


  new("vglmff",
  blurb = c("Ali-Mikhail-Haq distribution\n",
            "Links:    ",
            namesof("apar", lapar, earg = eapar )),
  initialize = eval(substitute(expression({
    if (any(y < 0) || any(y > 1))
        stop("the response must have values in ",
             "the unit square")

    temp5 <-
    w.y.check(w = w, y = y,
              Is.nonnegative.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)) {
      ainit  <- if (length( .iapar ))  .iapar else {
          mean1 <- if ( .imethod == 1)
                   weighted.mean(y[, 1], w) else
                   median(y[, 1])
          mean2 <- if ( .imethod == 1)
                   weighted.mean(y[, 2], w) else
                   median(y[, 2])
          Finit <- weighted.mean(y[, 1] <= mean1 &
                                 y[, 2] <= mean2, w)
          (1 - (mean1 * mean2 / Finit)) / (
          (1-mean1) * (1-mean2))
      }
      ainit <- min(0.95, max(ainit, -0.95))
      etastart <- theta2eta(rep_len(ainit, n), .lapar , .eapar )
    }
  }), list( .lapar = lapar, .eapar = eapar, .iapar = iapar,
            .imethod = imethod))),
  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$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) * dbiamhcop(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("biamhcop"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    apar <- eta2theta(eta, .lapar, earg = .eapar )
    okay1 <- all(is.finite(apar)) && all(abs(apar) < 1)
    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 )
    rbiamhcop(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 )

    y1 <- y[, 1]
    y2 <- y[, 2]
    de3 <- deriv3(~ (log(1 - apar+
                   (2 * apar*y1*y2/(1-apar*(1-y1)*(1-y2)))) -
                    2 * log(1 - apar*(1-y1)*(1-y2))) ,
                    name = "apar", hessian = FALSE)
    eval.de3 <- eval(de3)

    dl.dapar <-  attr(eval.de3, "gradient")

    c(w) * dl.dapar * dapar.deta
  }), list( .lapar = lapar, .eapar = eapar ))),
  weight = eval(substitute(expression({
    sd3 <- deriv3(~ (log(1 - apar +
                    (2 * apar * y1sim * y2sim / (1 - apar *
                    (1 - y1sim) * (1-y2sim)))) -
                     2 * log(1-apar*(1-y1sim)*(1-y2sim))),
                     name = "apar", hessian = FALSE)
    run.var <- 0
    for (ii in 1:( .nsimEIM )) {
      ysim <- rbiamhcop(n, apar = apar)
      y1sim <- ysim[, 1]
      y2sim <- ysim[, 1]
      eval.sd3 <- eval(sd3)
      dl.apar <-  attr(eval.sd3, "gradient")
      rm(ysim, y1sim, y2sim)
      temp3 <- dl.dapar
      run.var <- ((ii - 1) * run.var + temp3^2) / ii
    }

    wz <- if (intercept.only)
        matrix(colMeans(cbind(run.var)),
               n, dimm(M), byrow = TRUE) else cbind(run.var)

    wz <- wz * dapar.deta^2

    c(w) * wz
  }), list( .lapar = lapar,
            .eapar = eapar, .nsimEIM = nsimEIM ))))
}  # biamhcop








dbinorm <-
  function(x1, x2, mean1 = 0, mean2 = 0,
           var1 = 1, var2 = 1, cov12 = 0,
           log = FALSE) {
  if (!is.logical(log.arg <- log) || length(log) != 1)
    stop("bad input for argument 'log'")
  rm(log)


  sd1 <- sqrt(var1)
  sd2 <- sqrt(var2)
  rho <- cov12 / (sd1 * sd2)


  temp5 <- 1 - rho^2
  zedd1 <- (x1 - mean1) / sd1
  zedd2 <- (x2 - mean2) / sd2
  logpdf <- -log(2 * pi) - log(sd1) - log(sd2) -
     0.5 * log1p(-rho^2) +
   -(0.5 / temp5)  * (zedd1^2 +
                     (-2 * rho * zedd1 + zedd2) * zedd2)

  logpdf[is.infinite(x1) | is.infinite(x2)] <- log(0)
  if (log.arg) logpdf else exp(logpdf)
}  # dbinorm



rbinorm <- function(n, mean1 = 0, mean2 = 0,
                    var1 = 1, var2 = 1, cov12 = 0) {

  Y1 <- rnorm(n)
  Y2 <- rnorm(n)
  X1 <- sqrt(var1) * Y1 + mean1
  delta <- sqrt(var2 - (cov12^2) / var1)
  X2 <- cov12 * Y1 / sqrt(var1) + delta * Y2 + mean2

  ans <- cbind(X1, X2)
  ans[is.na(delta), ] <- NA

  ans
}  # rbinorm




 binormal <-
  function(lmean1 = "identitylink",
           lmean2 = "identitylink",
           lsd1   = "loglink",
           lsd2   = "loglink",
           lrho   = "rhobitlink",
           imean1 = NULL,       imean2 = NULL,
           isd1   = NULL,       isd2   = NULL,
           irho   = NULL,       imethod = 1,
           eq.mean = FALSE,     eq.sd = FALSE,
           zero = c("sd", "rho"),
           rho.arg = NA  # 20210923; possibly a known value
           ) {


  if (length(rho.arg) != 1)
    stop("argument 'rho.arg' must be scalar")
  est.rho <- is.na(rho.arg)  # Estimate rho?
  if (!est.rho && (!is.Numeric(rho.arg) ||
      rho.arg <= -1 || 1 <= rho.arg))
    stop("bad input for argument 'rho.arg'")

  if (!est.rho && is.character(zero) && any(zero == "rho")) {
    zero <- zero[zero != "rho"]
    if (length(zero) == 0)
      zero <- NULL  # Make sure
  }

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

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

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




  trivial1 <- is.logical(eq.mean) &&
              length(eq.mean) == 1 && !eq.mean
  trivial2 <- is.logical(eq.sd  ) &&
              length(eq.sd  ) == 1 && !eq.sd

  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 normal distribution\n",
            "Links:    ",
            namesof("mean1", lmean1, earg = emean1 ), ", ",
            namesof("mean2", lmean2, earg = emean2 ), ", ",
            namesof("sd1",   lsd1,   earg = esd1   ), ", ",
            namesof("sd2",   lsd2,   earg = esd2   ),
            if (est.rho) ", ",
            if (est.rho) namesof("rho", lrho, earg = erho )),
  constraints = eval(substitute(expression({

    constraints.orig <- constraints

    if (is.null(constraints.orig)) {

    M1.use <- M1 <- ifelse( .est.rho , 5, 4)
    NOS <- M / M1

            

    cm1.m <-
    cmk.m <- kronecker(diag(NOS),
                       rbind(diag(2),
                       matrix(0, ifelse( .est.rho , 3, 2), 2)))
    con.m <- cm.VGAM(kronecker(diag(NOS),
                       rbind(1, 1, 0, 0,
                             if ( .est.rho ) 0 else NULL)),
                     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, 2, 2),
                             diag(2),
              if ( .est.rho ) matrix(0, 1, 2) else NULL))
    con.s <- cm.VGAM(kronecker(diag(NOS),
                               rbind(0, 0, 1, 1,
                               if ( .est.rho ) 0 else NULL)),
                     x = x,
                     bool = .eq.sd ,  #
                     constraints = constraints.orig,
                     apply.int = TRUE,
                     cm.default           = cmk.s,
                     cm.intercept.default = cm1.s)


    con.use <- con.m
    for (klocal in seq_along(con.m)) {
      con.use[[klocal]] <-
        cbind(con.m[[klocal]],
              con.s[[klocal]],
              if ( .est.rho )
          kronecker(matrix(1, NOS, 1), diag(5)[, 5]) else NULL)

    }
  constraints <-
    cm.zero.VGAM(con.use,   # constraints,  # Prior to 20210923
                 x = x, .zero , M = M,
                 predictors.names = predictors.names,
                 M1 = M1.use)
  }  #  if (is.null(constraints.orig))



  }),
  list( .zero = zero, .est.rho = est.rho, .rho.arg = rho.arg,
        .eq.sd   = eq.sd,
        .eq.mean = eq.mean ))),

  infos = eval(substitute(function(...) {
    list(M1 = ifelse( .est.rho , 5, 4),
         Q1 = 2,
         expected = TRUE,
         multipleResponses = FALSE,
         parameters.names = c("mean1", "mean2", "sd1", "sd2",
                              if ( .est.rho ) "rho" else NULL),
         eq.mean = .eq.mean ,
         eq.sd   = .eq.sd   ,
         zero    = .zero )
    }, list( .zero    = zero, .est.rho = est.rho,
             .eq.mean = eq.mean,
             .eq.sd   = eq.sd    ))),

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

    temp5 <-
    w.y.check(w = w, y = y,
              ncol.w.max = 1,
              ncol.y.max = Q1,
              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("sd1",   .lsd1 ,   earg = .esd1 ,   short = TRUE),
      namesof("sd2",   .lsd2 ,   earg = .esd2 ,   short = TRUE),
      if ( .est.rho )
      namesof("rho",   .lrho ,   .erho ,   short = TRUE) else
      NULL)

    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)
      isd1 <- rep_len(if (length( .isd1 )) .isd1 else
                      sd(y[, 1]), n)
      isd2 <- rep_len(if (length( .isd2 )) .isd2 else
                      sd(y[, 2]), n)
      irho <- rep_len(if (length( .irho )) .irho else
                      cor(y[, 1], y[, 2]), n)

      if ( .imethod == 2) {
        imean1 <- abs(imean1) + 0.01
        imean2 <- abs(imean2) + 0.01
      }
      etastart <-
        cbind(theta2eta(imean1, .lmean1 , earg = .emean1 ),
              theta2eta(imean2, .lmean2 , earg = .emean2 ),
              theta2eta(isd1,   .lsd1 ,   earg = .esd1 ),
              theta2eta(isd2,   .lsd2 ,   earg = .esd2 ),
              if ( .est.rho )
              theta2eta(irho,   .lrho ,   earg = .erho ) else
              NULL)
    }
  }),
  list( .lmean1 = lmean1, .lmean2 = lmean2,
        .emean1 = emean1, .emean2 = emean2,
        .lsd1   = lsd1  , .lsd2   = lsd2  , .lrho = lrho,
        .esd1   = esd1  , .esd2   = esd2  , .erho = erho,
        .imethod = imethod,
        .est.rho = est.rho, .rho.arg = rho.arg,
        .imean1 = imean1, .imean2 = imean2,
        .isd1   = isd1,   .isd2   = isd2,
        .irho   = irho ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    NOS <- ncol(eta) / c(M1 = 5)
    mean1 <- eta2theta(eta[, 1], .lmean1 , earg = .emean1 )
    mean2 <- eta2theta(eta[, 2], .lmean2 , earg = .emean2 )
    fv.mat <- cbind(mean1, mean2)
    label.cols.y(fv.mat, colnames.y = extra$colnames.y,
                 NOS = NOS)
  }  ,
  list( .lmean1 = lmean1, .lmean2 = lmean2,
        .est.rho = est.rho, .rho.arg = rho.arg,
        .emean1 = emean1, .emean2 = emean2,
        .lsd1   = lsd1  , .lsd2   = lsd2  , .lrho = lrho,
        .esd1   = esd1  , .esd2   = esd2  , .erho = erho ))),

  last = eval(substitute(expression({
    misc$link <-
      c("mean1" = .lmean1 ,
        "mean2" = .lmean2 ,
        "sd1"   = .lsd1 ,
        "sd2"   = .lsd2 ,
        if ( .est.rho ) c("rho" = .lrho ) else NULL)

    if ( .est.rho ) {
    misc$earg <- list("mean1" = .emean1 ,
                      "mean2" = .emean2 ,
                      "sd1"   = .esd1 ,
                      "sd2"   = .esd2 ,
                      "rho"   = .erho )
    } else {
    misc$earg <- list("mean1" = .emean1 ,
                      "mean2" = .emean2 ,
                      "sd1"   = .esd1 ,
                      "sd2"   = .esd2 )
    }

    misc$expected <- TRUE
    misc$multipleResponses <- FALSE
  }) ,
  list( .lmean1 = lmean1, .lmean2 = lmean2,
        .emean1 = emean1, .emean2 = emean2,
        .est.rho = est.rho, .rho.arg = rho.arg,
        .lsd1   = lsd1  , .lsd2   = lsd2  , .lrho = lrho,
        .esd1   = esd1  , .esd2   = esd2  , .erho = erho ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    mean1 <- eta2theta(eta[, 1], .lmean1 , .emean1 )
    mean2 <- eta2theta(eta[, 2], .lmean2 , .emean2 )
    sd1   <- eta2theta(eta[, 3], .lsd1   , .esd1   )
    sd2   <- eta2theta(eta[, 4], .lsd2   , .esd2   )
    Rho   <- if ( .est.rho )
             eta2theta(eta[, 5], .lrho   , .erho ) else
             rep( .rho.arg , length = nrow(eta)) 

    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      ll.elts <-
        c(w) * dbinorm(x1 = y[, 1], x2 = y[, 2],
                       mean1 = mean1, mean2 = mean2,
                       var1 = sd1^2, var2 = sd2^2,
                       cov12 = Rho * sd1 * sd2,
                       log = TRUE)
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  } ,
  list( .lmean1 = lmean1, .lmean2 = lmean2,
        .emean1 = emean1, .emean2 = emean2,
        .est.rho = est.rho, .rho.arg = rho.arg,
        .lsd1   = lsd1  , .lsd2   = lsd2  , .lrho = lrho,
        .esd1   = esd1  , .esd2   = esd2  , .erho = erho,
        .imethod = imethod ))),
  vfamily = c("binormal"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    mean1 <- eta2theta(eta[, 1], .lmean1, earg = .emean1)
    mean2 <- eta2theta(eta[, 2], .lmean2, earg = .emean2)
    sd1   <- eta2theta(eta[, 3], .lsd1  , earg = .esd1  )
    sd2   <- eta2theta(eta[, 4], .lsd2  , earg = .esd2  )
    Rho   <- if ( .est.rho )
             eta2theta(eta[, 5], .lrho  , earg = .erho  ) else
             rep( .rho.arg , length = nrow(eta)) 
    okay1 <- all(is.finite(mean1)) &&
             all(is.finite(mean2)) &&
             all(is.finite(sd1  )) && all(0 < sd1) &&
             all(is.finite(sd2  )) && all(0 < sd2) &&
             all(is.finite(Rho  )) && all(abs(Rho) < 1)
    okay1
  } ,
  list( .lmean1 = lmean1, .lmean2 = lmean2,
        .emean1 = emean1, .emean2 = emean2,
        .est.rho = est.rho, .rho.arg = rho.arg,
        .lsd1   = lsd1  , .lsd2   = lsd2  , .lrho = lrho,
        .esd1   = esd1  , .esd2   = esd2  , .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)
    mean1 <- eta2theta(eta[, 1], .lmean1 , earg = .emean1 )
    mean2 <- eta2theta(eta[, 2], .lmean2 , earg = .emean2 )
    sd1   <- eta2theta(eta[, 3], .lsd1   , earg = .esd1   )
    sd2   <- eta2theta(eta[, 4], .lsd2   , earg = .esd2   )
    Rho   <- if ( .est.rho )
             eta2theta(eta[, 5], .lrho   , earg = .erho   ) else
             rep( .rho.arg , length = nrow(eta)) 
    rbinorm(nsim * length(sd1),
            mean1 = mean1, mean2 = mean2,
            var1 = sd1^2, var2 = sd2^2,
            cov12 = Rho * sd1 * sd2)
  } ,
  list( .lmean1 = lmean1, .lmean2 = lmean2,
        .emean1 = emean1, .emean2 = emean2,
        .est.rho = est.rho, .rho.arg = rho.arg,
        .lsd1   = lsd1  , .lsd2   = lsd2  , .lrho = lrho,
        .esd1   = esd1  , .esd2   = esd2  , .erho = erho ))),




  deriv = eval(substitute(expression({
    mean1 <- eta2theta(eta[, 1], .lmean1, earg = .emean1)
    mean2 <- eta2theta(eta[, 2], .lmean2, earg = .emean2)
    sd1   <- eta2theta(eta[, 3], .lsd1  , earg = .esd1  )
    sd2   <- eta2theta(eta[, 4], .lsd2  , earg = .esd2  )
    Rho   <- if ( .est.rho )
             eta2theta(eta[, 5], .lrho  , earg = .erho ) else
             rep( .rho.arg , length = nrow(eta)) 

    zedd1 <- (y[, 1] - mean1) / sd1
    zedd2 <- (y[, 2] - mean2) / sd2
    temp5 <- 1 - Rho^2

    SigmaInv <- matrix(0, n, dimm(2))
    SigmaInv[, iam(1, 1, M = 2)] <- 1 / ((sd1^2) * temp5)
    SigmaInv[, iam(2, 2, M = 2)] <- 1 / ((sd2^2) * temp5)
    SigmaInv[, iam(1, 2, M = 2)] <- -Rho / (sd1 * sd2 * temp5)
    dl.dmeans <- mux22(t(SigmaInv), y - cbind(mean1, mean2),
                       M = 2, as.matrix = TRUE)
    dl.dsd1   <- -1 / sd1 +
                 zedd1 * (zedd1 - Rho * zedd2) / (sd1 * temp5)
    dl.dsd2   <- -1 / sd2 +
                 zedd2 * (zedd2 - Rho * zedd1) / (sd2 * temp5)
    dl.drho   <- -Rho * (zedd1^2 - 2 * Rho * zedd1 * zedd2 +
                         zedd2^2) / temp5^2 +
                zedd1 * zedd2 / temp5 +
                Rho / temp5

    dmean1.deta <- dtheta.deta(mean1, .lmean1)
    dmean2.deta <- dtheta.deta(mean2, .lmean2)
    dsd1.deta   <- dtheta.deta(sd1  , .lsd1  )
    dsd2.deta   <- dtheta.deta(sd2  , .lsd2  )
    drho.deta   <- dtheta.deta(Rho  , .lrho  )
    dthetas.detas  <- cbind(dmean1.deta,
                            dmean2.deta,
                            dsd1.deta,
                            dsd2.deta,
                            if ( .est.rho ) drho.deta else NULL)

    c(w) * cbind(dl.dmeans[, 1],
                 dl.dmeans[, 2],
                 dl.dsd1,
                 dl.dsd2,
                 if ( .est.rho ) dl.drho else NULL) *
    dthetas.detas
  }),
  list( .lmean1 = lmean1, .lmean2 = lmean2,
        .emean1 = emean1, .emean2 = emean2,
        .est.rho = est.rho, .rho.arg = rho.arg,
        .lsd1   = lsd1  , .lsd2   = lsd2  , .lrho = lrho,
        .esd1   = esd1  , .esd2   = esd2  , .erho = erho,
        .imethod = imethod ))),
  weight = eval(substitute(expression({
    wz <- matrix(0.0, n, dimm(M))
    wz[, iam(1, 1, M)] <- SigmaInv[, iam(1, 1, M = 2)]
    wz[, iam(2, 2, M)] <- SigmaInv[, iam(2, 2, M = 2)]
    wz[, iam(1, 2, M)] <- SigmaInv[, iam(1, 2, M = 2)]
    wz[, iam(3, 3, M)] <- (1 + 1 / temp5) / sd1^2
    wz[, iam(4, 4, M)] <- (1 + 1 / temp5) / sd2^2
    wz[, iam(3, 4, M)] <- -(Rho^2) / (temp5 * sd1 * sd2)
    if ( .est.rho ) {
      wz[, iam(5, 5, M)] <- (1 + Rho^2) / temp5^2
      wz[, iam(3, 5, M)] <- -Rho / (sd1 * temp5)
      wz[, iam(4, 5, M)] <- -Rho / (sd2 * temp5)
    }
    for (ilocal in 1:M)
      for (jlocal in ilocal:M)
        wz[, iam(ilocal, jlocal, M)] <-
        wz[, iam(ilocal, jlocal, M)] * dthetas.detas[, ilocal] *
                                       dthetas.detas[, jlocal]
      c(w) * wz
  }),
  list( .lmean1 = lmean1, .lmean2 = lmean2,
        .emean1 = emean1, .emean2 = emean2,
        .est.rho = est.rho, .rho.arg = rho.arg,
        .lsd1   = lsd1  , .lsd2   = lsd2  , .lrho = lrho,
        .esd1   = esd1  , .esd2   = esd2  , .erho = erho,
        .imethod = imethod ))))
}  # binormal









gumbelI <-
    function(la = "identitylink", earg = list(),
             ia = NULL, imethod = 1) {

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



  if (length(ia) && !is.Numeric(ia, length.arg = 1))
      stop("'ia' must be a single number")

  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("Gumbel's Type I Bivariate Distribution\n",
          "Links:    ",
          namesof("a", la, earg =  earg )),
  initialize = eval(substitute(expression({
    if (!is.matrix(y) || ncol(y) != 2)
        stop("the response must be a 2 column matrix")

    if (any(y < 0))
        stop("the response must have non-negative values only")


    extra$colnames.y  <- colnames(y)

    predictors.names <-
      c(namesof("a", .la, earg =  .earg , short = TRUE))
    if (!length(etastart)) {
      ainit  <- if (length( .ia ))  rep_len( .ia , n) else {
        mean1 <- if ( .imethod == 1)
                 median(y[,1]) else mean(y[,1])
        mean2 <- if ( .imethod == 1)
                 median(y[,2]) else mean(y[,2])
        Finit <- 0.01 + mean(y[,1] <= mean1 & y[,2] <= mean2)
        (log(Finit+expm1(-mean1)+exp(-mean2))+
         mean1+mean2)/(mean1*mean2)
      }
      etastart <- theta2eta(rep_len(ainit, n), .la , .earg )
      }
  }), list( .ia = ia,
            .la = la, .earg = earg, .imethod = imethod ))),
  linkinv = eval(substitute(function(eta, extra = NULL) {
    NOS <- NCOL(eta) / c(M1 = 1)
    Q1 <- 2
    fv.mat <- matrix(1, NROW(eta), NOS * Q1)
    label.cols.y(fv.mat, colnames.y = extra$colnames.y,
                 NOS = NOS)
  }, list( .la = la ))),
  last = eval(substitute(expression({
    misc$link <-    c("a" = .la )
    misc$earg <- list("a" = .earg )

    misc$expected <- FALSE
    misc$pooled.weight <- pooled.weight
  }), list( .la = la, .earg = earg ))),
  loglikelihood = eval(substitute(
    function(mu, y, w, residuals = FALSE, eta,
             extra = NULL,
             summation = TRUE) {
    alpha  <- eta2theta(eta, .la, earg =  .earg )
    if (residuals) {
      stop("loglikelihood residuals not implemented yet")
    } else {
      denom <- (alpha*y[,1] - 1) * (alpha*y[,2] - 1) + alpha
      mytolerance <- .Machine$double.xmin
      bad <- (denom <= mytolerance)   # Range violation
      if (any(bad)) {
        cat("There are some range violations in @deriv\n")
        flush.console()
          denom[bad] <- 2 * mytolerance
      }
      ll.elts <- c(w) * (-y[,1] - y[,2] +
                         alpha*y[,1]*y[,2] + log(denom))
      if (summation) {
        sum(ll.elts)
      } else {
        ll.elts
      }
    }
  }, list( .la = la, .earg = earg ))),
  vfamily = c("gumbelI"),
  validparams = eval(substitute(function(eta, y, extra = NULL) {
    alpha  <- eta2theta(eta, .la , earg = .earg )
    okay1 <- all(is.finite(alpha))
    okay1
  } , list( .la = la, .earg = earg ))),
  deriv = eval(substitute(expression({
      alpha  <- eta2theta(eta, .la, earg =  .earg )
      numerator <- (alpha*y[,1] - 1) * y[,2] +
                   (alpha*y[,2] - 1)*y[,1] + 1
      denom <- (alpha*y[,1] - 1) * (alpha*y[,2] - 1) + alpha
      denom <- abs(denom)
      dl.dalpha <- numerator / denom + y[,1]*y[,2]
      dalpha.deta <- dtheta.deta(alpha,  .la, earg =  .earg )
      c(w) * cbind(dl.dalpha * dalpha.deta)
  }), list( .la = la, .earg = earg ))),
  weight = eval(substitute(expression({
    d2l.dalpha2 <- (numerator/denom)^2 - 2*y[,1]*y[,2] / denom
    d2alpha.deta2 <- d2theta.deta2(alpha, .la, earg =  .earg )
    wz <- w * (dalpha.deta^2 * d2l.dalpha2 -
               d2alpha.deta2 * dl.dalpha)
    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( .la = la, .earg = earg ))))
}  # gumbelI






kendall.tau <-
    function(x, y, exact = FALSE, max.n = 3000) {

  if ((N <- length(x)) != length(y))
    stop("arguments 'x' and 'y' do not have equal lengths")

  NN <- if (!exact && N > max.n) {
    cindex <- sample.int(n = N, size = max.n, replace = FALSE)
    x <- x[cindex]
    y <- y[cindex]
    max.n
  } else {
    N
  }


  ans3 <-
    c( .C("VGAM_C_kend_tau",
         as.double(x), as.double(y),
         as.integer(NN), ans = double(3),
         NAOK = TRUE)$ans)

  con <- ans3[1] + ans3[2] / 2  # Ties put half and half
  dis <- ans3[3] + ans3[2] / 2
  (con - dis) / (con + dis)
}  # kendall.tau




if (FALSE)
kendall.tau <- function(x, y, exact = TRUE, max.n = 1000) {

  if ((N <- length(x)) != length(y))
    stop("arguments 'x' and 'y' do not have equal lengths")
  index <- iam(NA, NA, M = N, both = TRUE)

  index$row.index <- index$row.index[-(1:N)]
  index$col.index <- index$col.index[-(1:N)]

  NN <- if (!exact && N > max.n) {
    cindex <- sample.int(n = N, size = max.n, replace = FALSE)
    index$row.index <- index$row.index[cindex]
    index$col.index <- index$col.index[cindex]
    max.n
  } else{
    choose(N, 2)
  }

  con <- sum((x[index$row.index] - x[index$col.index]) *
             (y[index$row.index] - y[index$col.index]) > 0)
  dis <- NN - con
  (con - dis) / (con + dis)
}  # kendall.tau





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

  u1 <- qt(x1, df = df)
  u2 <- qt(x2, df = df)

  logdensity <-
    -(df/2 + 1) * log1p(
    (u1^2 + u2^2 - 2 * rho * u1 * u2) / (df * (1 - rho^2))) -
    log(2*pi) - 0.5 * log1p(-rho^2) -
  dt(u1, df = df, log = TRUE) -
  dt(u2, df = df, log = TRUE)

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

Try the VGAM package in your browser

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

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