Running examples in 'bzinb-Ex.R' failed
The error most likely occurred in:

> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: bnb
> ### Title: The bivariate negative binomial distribution
> ### Aliases: bnb lik.bnb rbnb
> 
> ### ** Examples
> 
> # generating a pair of random vectors
> set.seed(1)
> data1 <- rbnb(n = 100, a0 = 2, a1 = 1, a2 = 1, 
+                 b1 = 1, b2 = 1)
> 
> lik.bnb(xvec = data1[, 1], yvec = data1[ ,2], 
+           a0 = 1, a1 = 1, a2 = 1, b1 = 1, b2 = 1) 
[1] -447.4564
> 
> bnb(xvec = data1[,1], yvec = data1[,2], showFlag = FALSE)
 ----------- FAILURE REPORT -------------- 
 --- failure: the condition has length > 1 ---
 --- srcref --- 
: 
 --- package (from environment) --- 
bzinb
 --- call from context --- 
bzinb.base(xvec, yvec, initial = initial, tol = tol, maxiter = maxiter, 
    showFlag = showFlag, vcov = vcov, bnb = 1)
 --- call from argument --- 
if (class(cov.mat) == "try-error") {
    std.param = setNames(rep(NA, dim.std), c(par.names, "rho", 
        "logit.rho"))
    cov.mat <- NA
} else {
    if (!bnb) 
        var.p4 <- sum(cov.mat[6:8, 6:8])
    d.g <- rho * c(1/param[1] - 1/{
        2 * (param[1] + param[2])
    } - 1/{
        2 * (param[1] + param[3])
    }, -1/{
        2 * (param[1] + param[2])
    }, -1/{
        2 * (param[1] + param[3])
    }, 1/{
        2 * param[4] * (param[4] + 1)
    }, 1/{
        2 * param[5] * (param[5] + 1)
    })
    var.rho <- t(d.g) %*% cov.mat[1:5, 1:5] %*% d.g
    var.logit.rho <- var.rho/rho^2/(1 - rho)^2
    if (!bnb) {
        std.param = sqrt(c(diag(cov.mat), p4 = var.p4, rho = var.rho, 
            logit.rho = var.logit.rho))
    }
    else {
        std.param = sqrt(c(diag(cov.mat), rho = var.rho, logit.rho = var.logit.rho))
    }
}
 --- R stacktrace ---
where 1: bzinb.base(xvec, yvec, initial = initial, tol = tol, maxiter = maxiter, 
    showFlag = showFlag, vcov = vcov, bnb = 1)
where 2: doTryCatch(return(expr), name, parentenv, handler)
where 3: tryCatchOne(expr, names, parentenv, handlers[[1L]])
where 4: tryCatchList(expr, classes, parentenv, handlers)
where 5: tryCatch(expr, error = function(e) {
    call <- conditionCall(e)
    if (!is.null(call)) {
        if (identical(call[[1L]], quote(doTryCatch))) 
            call <- sys.call(-4L)
        dcall <- deparse(call)[1L]
        prefix <- paste("Error in", dcall, ": ")
        LONG <- 75L
        sm <- strsplit(conditionMessage(e), "\n")[[1L]]
        w <- 14L + nchar(dcall, type = "w") + nchar(sm[1L], type = "w")
        if (is.na(w)) 
            w <- 14L + nchar(dcall, type = "b") + nchar(sm[1L], 
                type = "b")
        if (w > LONG) 
            prefix <- paste0(prefix, "\n  ")
    }
    else prefix <- "Error : "
    msg <- paste0(prefix, conditionMessage(e), "\n")
    .Internal(seterrmessage(msg[1L]))
    if (!silent && isTRUE(getOption("show.error.messages"))) {
        cat(msg, file = outFile)
        .Internal(printDeferredWarnings())
    }
    invisible(structure(msg, class = "try-error", condition = e))
})
where 6: try(bzinb.base(xvec, yvec, initial = initial, tol = tol, maxiter = maxiter, 
    showFlag = showFlag, vcov = vcov, bnb = 1))
where 7: bnb(xvec = data1[, 1], yvec = data1[, 2], showFlag = FALSE)

 --- value of length: 2 type: logical ---
[1] FALSE FALSE
 --- function from context --- 
function (xvec, yvec, initial = NULL, tol = 1e-08, maxiter = 50000, 
    showFlag = FALSE, vcov = FALSE, bnb = 0) 
{
    se = TRUE
    xy.reduced <- as.data.frame(table(xvec, yvec))
    names(xy.reduced) <- c("x", "y", "freq")
    xy.reduced <- xy.reduced[xy.reduced$freq != 0, ]
    xy.reduced$x <- as.numeric(as.character(xy.reduced$x))
    xy.reduced$y <- as.numeric(as.character(xy.reduced$y))
    xy.reduced$freq <- as.numeric(as.character(xy.reduced$freq))
    n <- sum(xy.reduced$freq)
    n.reduced <- as.integer(length(xy.reduced$freq))
    if (max(xvec) == 0 & max(yvec) == 0) {
        return(c(rep(1e-10, 5), 1, 0, 0, 0, 0, 1, 0, if (se) {
            rep(NA, 11)
        }))
    }
    info <- if (se) {
        matrix(0, ncol = 8, nrow = 8, dimnames = list(abp.names[-9], 
            abp.names[-9]))
    }
    else {
        0
    }
    if (is.null(initial)) {
        xbar <- mean(xvec)
        ybar <- mean(yvec)
        xybar <- mean(c(xbar, ybar))
        s2.x <- var(xvec)
        s2.y <- var(yvec)
        if (is.na(s2.x) | is.na(s2.y)) {
            s2.x <- s2.y <- 1
        }
        cor.xy <- cor(xvec, yvec)
        if (is.na(cor.xy)) {
            cor.xy <- 0
        }
        zero <- sum(xvec == 0 & yvec == 0)/n
        initial <- rep(NA, 9)
        names(initial) <- c("a0", "a1", "a2", "b1", "b2", "p1", 
            "p2", "p3", "p4")
        initial[4] <- s2.x/ifelse(xbar == 0, 1e-04, xbar)
        initial[5] <- s2.y/ifelse(ybar == 0, 1e-04, ybar)
        initial[2:3] <- c(xbar, ybar)/pmax(initial[4:5], c(0.1, 
            0.1))
        initial[1] <- min(initial[2:3]) * abs(cor.xy)
        initial[2:3] <- initial[2:3] - initial[1]
        initial[6:9] <- bin.profile(xvec, yvec)
        initial[6:9] <- initial[6:9]/sum(initial[6:9])
        initial <- pmax(initial, 1e-05)
        if (is.na(sum(initial))) {
            initial[is.na(initial)] <- 1
        }
    }
    else {
        names(initial) <- abp.names
    }
    param = initial
    lik = -Inf
    lik.vec = rep(0, maxiter + 1)
    nonconv = 0L
    em.out <- em(param2 = param, xvec = xy.reduced$x, yvec = xy.reduced$y, 
        freq = xy.reduced$freq, n = n.reduced, se = as.integer(se), 
        maxiter = as.integer(maxiter), tol = as.double(tol), 
        showFlag = as.integer(showFlag), bnb = bnb)
    names(em.out) <- c("param2", "expt", "info", "iter", "nonconv", 
        "trajectory")
    param = setNames(em.out$param2, abp.names)
    expt = setNames(em.out$expt, expt.names)
    info = if (se) {
        matrix(em.out$info, 8, 8, dimnames = list(abp.names[-9], 
            abp.names[-9]))
    }
    else 0
    iter = em.out$iter
    nonconv = em.out$nonconv
    trajectory = em.out$trajectory
    if (nonconv == 1) 
        warning("The iteration exited before reaching convergence.")
    rho <- param[1]/sqrt((param[1] + param[2]) * (param[1] + 
        param[3])) * sqrt(param[4] * param[5]/(param[4] + 1)/(param[5] + 
        1))
    logit.rho <- qlogis(rho)
    if (bnb) {
        info = info[1:5, 1:5]
    }
    if (se) {
        par.names <- if (bnb) {
            abp.names[1:5]
        }
        else {
            abp.names
        }
        dim.std <- if (bnb) {
            7
        }
        else {
            11
        }
        fullrank <- if (bnb) {
            5
        }
        else {
            8
        }
        qr.info <- try(qr(info))
        if (class(qr.info) == "try-error") {
            warning("The information matrix has NA/NaN/Inf and thus the standard error is not properly estimatd.")
            std.param = setNames(rep(NA, dim.std), c(par.names, 
                "rho", "logit.rho"))
            cov.mat <- NA
        }
        else if (qr(info)$rank < fullrank) {
            warning("The information matrix is (essentially) not full rank, and thus the standard error is not reliable.")
            std.param = setNames(rep(NA, dim.std), c(par.names, 
                "rho", "logit.rho"))
            cov.mat <- NA
        }
        else {
            cov.mat <- try(solve(info))
            if (class(cov.mat) == "try-error") {
                std.param = setNames(rep(NA, dim.std), c(par.names, 
                  "rho", "logit.rho"))
                cov.mat <- NA
            }
            else {
                if (!bnb) 
                  var.p4 <- sum(cov.mat[6:8, 6:8])
                d.g <- rho * c(1/param[1] - 1/{
                  2 * (param[1] + param[2])
                } - 1/{
                  2 * (param[1] + param[3])
                }, -1/{
                  2 * (param[1] + param[2])
                }, -1/{
                  2 * (param[1] + param[3])
                }, 1/{
                  2 * param[4] * (param[4] + 1)
                }, 1/{
                  2 * param[5] * (param[5] + 1)
                })
                var.rho <- t(d.g) %*% cov.mat[1:5, 1:5] %*% d.g
                var.logit.rho <- var.rho/rho^2/(1 - rho)^2
                if (!bnb) {
                  std.param = sqrt(c(diag(cov.mat), p4 = var.p4, 
                    rho = var.rho, logit.rho = var.logit.rho))
                }
                else {
                  std.param = sqrt(c(diag(cov.mat), rho = var.rho, 
                    logit.rho = var.logit.rho))
                }
            }
        }
    }
    result <- list(rho = matrix(c(rho, logit.rho, if (se) std.param[c("rho", 
        "logit.rho")] else rep(NA, 2)), ncol = 2, dimnames = list(c("rho", 
        "logit.rho"), c("Estimate", "Std.err"))), coefficients = matrix(c(param, 
        if (se) std.param[1:9] else rep(NA, 9)), ncol = 2, dimnames = list(abp.names, 
        c("Estimate", "Std.err"))), lik = expt[1], iter = iter)
    if (se & vcov) {
        result$info = info
        result$vcov = cov.mat
    }
    return(result)
}
<bytecode: 0x25ce9b0>
<environment: namespace:bzinb>
 --- function search by body ---
Function bzinb.base in namespace bzinb has this body.
 ----------- END OF FAILURE REPORT -------------- 
Fatal error: the condition has length > 1