CRAN Package Check Results for Package mvProbit

Last updated on 2020-02-01 08:49:10 CET.

Flavor Version Tinstall Tcheck Ttotal Status Flags
r-devel-linux-x86_64-debian-clang 0.1-8 4.40 190.33 194.73 ERROR
r-devel-linux-x86_64-debian-gcc 0.1-8 3.35 235.52 238.87 OK
r-devel-linux-x86_64-fedora-clang 0.1-8 234.94 ERROR
r-devel-linux-x86_64-fedora-gcc 0.1-8 232.12 ERROR
r-devel-windows-ix86+x86_64 0.1-8 8.00 415.00 423.00 OK
r-devel-windows-ix86+x86_64-gcc8 0.1-8 13.00 504.00 517.00 OK
r-patched-linux-x86_64 0.1-8 3.90 301.12 305.02 OK
r-patched-solaris-x86 0.1-8 429.50 OK
r-release-linux-x86_64 0.1-8 3.53 296.16 299.69 OK
r-release-windows-ix86+x86_64 0.1-8 10.00 493.00 503.00 OK
r-release-osx-x86_64 0.1-8 OK
r-oldrel-windows-ix86+x86_64 0.1-8 5.00 389.00 394.00 OK
r-oldrel-osx-x86_64 0.1-8 OK

Additional issues

ATLAS

Check Details

Version: 0.1-8
Check: examples
Result: ERROR
    Running examples in 'mvProbit-Ex.R' failed
    The error most likely occurred in:
    
    > base::assign(".ptime", proc.time(), pos = "CheckExEnv")
    > ### Name: mvProbit
    > ### Title: Estimation of Multivariate Probit Models
    > ### Aliases: mvProbit print.mvProbit
    > ### Keywords: models regression
    >
    > ### ** Examples
    >
    > ## generate a simulated data set
    > set.seed( 123 )
    > # number of observations
    > nObs <- 50
    >
    > # generate explanatory variables
    > xMat <- cbind(
    + const = rep( 1, nObs ),
    + x1 = as.numeric( rnorm( nObs ) > 0 ),
    + x2 = rnorm( nObs ) )
    >
    > # model coefficients
    > beta <- cbind( c( 0.8, 1.2, -0.8 ),
    + c( -0.6, 1.0, -1.6 ),
    + c( 0.5, -0.6, 1.2 ) )
    >
    > # covariance matrix of error terms
    > library( miscTools )
    > sigma <- symMatrix( c( 1, 0.2, 0.4, 1, -0.1, 1 ) )
    >
    > # generate dependent variables
    > yMatLin <- xMat %*% beta
    > yMat <- ( yMatLin + rmvnorm( nObs, sigma = sigma ) ) > 0
    > colnames( yMat ) <- paste( "y", 1:3, sep = "" )
    >
    > # estimation (BHHH optimizer and GHK algorithm)
    > estResult <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
    + data = as.data.frame( cbind( xMat, yMat ) ), iterlim = 1, nGHK = 50 )
    Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
     ----------- FAILURE REPORT --------------
     --- failure: the condition has length > 1 ---
     --- srcref ---
    :
     --- package (from environment) ---
    mvProbit
     --- call from context ---
    pmvnormWrap(upper = xBetaTmp, sigma = sigmaTmp, algorithm = algorithm,
     nGHK = nGHK, random.seed = randomSeed, ...)
     --- call from argument ---
    if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
    }
     --- R stacktrace ---
    where 1: pmvnormWrap(upper = xBetaTmp, sigma = sigmaTmp, algorithm = algorithm,
     nGHK = nGHK, random.seed = randomSeed, ...)
    where 2: mvProbitLogLikInternal(yMat = yMat, xMat = xMat, coef = param,
     sigma = NULL, algorithm = llAlgorithm, nGHK = llNGHK, returnGrad = llIntGrad,
     oneSidedGrad = llOneSidedGrad, eps = llEps, randomSeed = llRandom.seed,
     ...)
    where 3: fnOrig(theta, ...)
    where 4: fn(start1, fixed = fixed, sumObs = TRUE, returnHessian = returnHessian,
     ...)
    where 5: maxNRCompute(fn = function (theta, fnOrig, gradOrig = NULL, hessOrig = NULL,
     fixed, sumObs = FALSE, returnHessian = TRUE, ...)
    {
     nParam <- length(theta)
     f <- fnOrig(theta, ...)
     if (any(is.na(f))) {
     attr(f, "gradient") <- NA
     attr(f, "hessian") <- NA
     return(f)
     }
     gr <- attr(f, "gradient")
     if (is.null(gr)) {
     if (!is.null(gradOrig)) {
     gr <- gradOrig(theta, ...)
     }
     else {
     gr <- numericGradient(f = fnOrig, t0 = theta, fixed = fixed,
     ...)
     }
     }
     if (is.matrix(gr)) {
     if (ncol(gr) != length(theta)) {
     stop(paste0("if gradient is a matrix, it must have length(parameter) colums (currently ",
     length(theta), "), not ", ncol(gr)))
     }
     activeGr <- gr[, !fixed]
     }
     else {
     activeGr <- gr[!fixed]
     }
     if (any(is.na(activeGr))) {
     attr(f, "gradient") <- gr
     attr(f, "hessian") <- NA
     return(f)
     }
     if (observationGradient(gr, length(theta))) {
     gr <- as.matrix(gr)
     }
     if (is.null(dim(gr))) {
     gr[fixed] <- NA
     }
     else {
     gr[, fixed] <- NA
     }
     if (isTRUE(returnHessian)) {
     h <- attr(f, "hessian")
     if (is.null(h)) {
     if (!is.null(hessOrig)) {
     h <- as.matrix(hessOrig(theta, ...))
     }
     else {
     llFunc <- function(theta, ...) {
     return(sum(fnOrig(theta, ...)))
     }
     if (!is.null(attr(f, "gradient"))) {
     gradFunc <- function(theta, ...) {
     return(sumGradients(attr(fnOrig(theta, ...),
     "gradient"), nParam))
     }
     }
     else if (!is.null(gradOrig)) {
     gradFunc <- function(theta, ...) {
     return(sumGradients(gradOrig(theta, ...),
     nParam))
     }
     }
     else {
     gradFunc <- NULL
     }
     h <- numericHessian(f = llFunc, grad = gradFunc,
     t0 = theta, fixed = fixed, ...)
     }
     }
     if ((dim(h)[1] != nParam) | (dim(h)[2] != nParam)) {
     stop("Wrong hessian dimension. Needed ", nParam,
     "x", nParam, " but supplied ", dim(h)[1], "x",
     dim(h)[2])
     }
     else {
     h[fixed, ] <- NA
     h[, fixed] <- NA
     }
     }
     else if (tolower(returnHessian) == "bhhh") {
     h <- NULL
     if (is.null(dim(gr)) & any(is.na(gr[!fixed]))) {
     h <- NA
     }
     else if (is.matrix(gr)) {
     if (any(is.na(gr[, !fixed]))) {
     h <- NA
     }
     }
     if (is.null(h)) {
     checkBhhhGrad(g = gr, theta = theta, analytic = (!is.null(attr(f,
     "gradient")) || !is.null(gradOrig)), fixed = fixed)
     h <- -crossprod(gr)
     }
     attr(h, "type") = "BHHH"
     }
     else {
     h <- NULL
     }
     if (sumObs) {
     f <- sumKeepAttr(f)
     }
     if (sumObs) {
     gr <- sumGradients(gr, nParam)
     }
     if (!is.null(gradOrig) && !is.null(attr(f, "gradient"))) {
     attr(f, "gradBoth") <- TRUE
     }
     if (!is.null(hessOrig) && !is.null(attr(f, "hessian"))) {
     attr(f, "hessBoth") <- TRUE
     }
     attr(f, "gradient") <- gr
     attr(f, "hessian") <- h
     return(f)
    }, fnOrig = function (param, yMat, xMat, llAlgorithm, llNGHK,
     llIntGrad, llOneSidedGrad, llEps, llRandom.seed, ...)
    {
     logLikVal <- mvProbitLogLikInternal(yMat = yMat, xMat = xMat,
     coef = param, sigma = NULL, algorithm = llAlgorithm,
     nGHK = llNGHK, returnGrad = llIntGrad, oneSidedGrad = llOneSidedGrad,
     eps = llEps, randomSeed = llRandom.seed, ...)
     return(logLikVal)
    }, start = c(b_1_0 = 0.598378200054935, b_1_1 = 6.336646267785,
    b_1_2 = -1.01231147815159, b_2_0 = -1.03564793824529, b_2_1 = 1.26302210717712,
    b_2_2 = -1.70638777681512, b_3_0 = 0.410893287564513, b_3_1 = -0.373151392833251,
    b_3_2 = 2.92074017047352, R_1_2 = 0.172617843070694, R_1_3 = -0.0160393443285098,
    R_2_3 = -0.0704265674853818), finalHessian = "BHHH", bhhhHessian = TRUE,
     fixed = c(b_1_0 = FALSE, b_1_1 = FALSE, b_1_2 = FALSE, b_2_0 = FALSE,
     b_2_1 = FALSE, b_2_2 = FALSE, b_3_0 = FALSE, b_3_1 = FALSE,
     b_3_2 = FALSE, R_1_2 = FALSE, R_1_3 = FALSE, R_2_3 = FALSE
     ), control = new("MaxControl", tol = 1e-08, reltol = 1.49011611938477e-08,
     gradtol = 1e-06, steptol = 1e-10, lambdatol = 1e-06,
     qrtol = 1e-10, qac = "stephalving", marquardt_lambda0 = 0.01,
     marquardt_lambdaStep = 2, marquardt_maxLambda = 1e+12,
     nm_alpha = 1, nm_beta = 0.5, nm_gamma = 2, sann_cand = NULL,
     sann_temp = 10, sann_tmax = 10L, sann_randomSeed = 123L,
     iterlim = 1L, printLevel = 0L), yMat = c(0, 1, 1, 1,
     1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1,
     1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1,
     1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1,
     0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1,
     1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1,
     1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0,
     1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1,
     1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0), xMat = c(1, 1, 1,
     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0,
     1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0,
     1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0,
     1, 0, 0.253318513994755, -0.028546755348703, -0.0428704572913161,
     1.36860228401446, -0.225770985659268, 1.51647060442954, -1.54875280423022,
     0.584613749636069, 0.123854243844614, 0.215941568743973,
     0.379639482759882, -0.502323453109302, -0.33320738366942,
     -1.01857538310709, -1.07179122647558, 0.303528641404258,
     0.448209778629426, 0.0530042267305041, 0.922267467879738,
     2.05008468562714, -0.491031166056535, -2.30916887564081,
     1.00573852446226, -0.709200762582393, -0.688008616467358,
     1.0255713696967, -0.284773007051009, -1.22071771225454, 0.18130347974915,
     -0.138891362439045, 0.00576418589988693, 0.38528040112633,
     -0.370660031792409, 0.644376548518833, -0.220486561818751,
     0.331781963915697, 1.09683901314935, 0.435181490833803, -0.325931585531227,
     1.14880761845109, 0.993503855962119, 0.54839695950807, 0.238731735111441,
     -0.627906076039371, 1.36065244853001, -0.600259587147127,
     2.18733299301658, 1.53261062618519, -0.235700359100477, -1.02642090030678
     ), llAlgorithm = "GHK", llNGHK = 50, llIntGrad = TRUE, llOneSidedGrad = FALSE,
     llEps = 1e-06, llRandom.seed = 123)
    where 6: eval(as.call(cl))
    where 7: eval(as.call(cl))
    where 8: maxNR(fn = fn, grad = grad, hess = hess, start = start, finalHessian = finalHessian,
     bhhhHessian = TRUE, ...)
    where 9: maxRoutine(fn = logLik, grad = grad, hess = hess, start = start,
     constraints = constraints, ...)
    where 10: maxLik(logLik = logLik, start = start, method = method, finalHessian = finalHessian,
     yMat = yMat, xMat = xMat, llAlgorithm = algorithm, llNGHK = nGHK,
     llIntGrad = intGrad, llOneSidedGrad = oneSidedGrad, llEps = eps,
     llRandom.seed = random.seed, ...)
    where 11: mvProbit(cbind(y1, y2, y3) ~ x1 + x2, data = as.data.frame(cbind(xMat,
     yMat)), iterlim = 1, nGHK = 50)
    
     --- value of length: 2 type: logical ---
    [1] FALSE FALSE
     --- function from context ---
    function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
     nGHK = NULL, ...)
    {
     if (!is.matrix(sigma)) {
     stop("argument 'sigma' must be a matrix")
     }
     else if (nrow(sigma) != ncol(sigma)) {
     stop("argument 'sigma' must be a square matrix")
     }
     else if (!all(is.numeric(sigma))) {
     stop("argument 'sigma must be a numeric matrix")
     }
     if (!all(is.numeric(lower))) {
     stop("argument 'lower' must be numeric")
     }
     else if (length(lower) == 1) {
     lower <- rep(lower, nrow(sigma))
     }
     else if (length(lower) != nrow(sigma)) {
     stop("argument 'lower' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     if (!all(is.numeric(upper))) {
     stop("argument 'upper' must be numeric")
     }
     else if (length(upper) == 1) {
     upper <- rep(upper, nrow(sigma))
     }
     else if (length(upper) != nrow(sigma)) {
     stop("argument 'upper' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     algOkay <- TRUE
     ghk <- FALSE
     if (!is.list(algorithm) && length(algorithm) != 1) {
     stop("argument 'algorithm' must be a single function",
     " or a single character string")
     }
     else if (is.function(algorithm)) {
     algResult <- do.call(algorithm, list())
     if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
     "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (is.character(algorithm)) {
     if (tolower(algorithm) == "ghk") {
     ghk <- TRUE
     }
     else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     if (!algOkay) {
     stop("argument 'algorithm' must be either one of the functions",
     " 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
     " \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
     }
     if (length(random.seed) != 1) {
     stop("argument 'random.seed' must be a single numerical value")
     }
     else if (!is.numeric(random.seed)) {
     stop("argument 'random.seed' must be numerical")
     }
     if (exists(".Random.seed")) {
     savedSeed <- .Random.seed
     }
     set.seed(random.seed)
     if (exists("savedSeed")) {
     on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
     }
     else {
     on.exit(rm(.Random.seed, envir = sys.frame()))
     }
     if (ghk) {
     if (is.null(nGHK)) {
     stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
     }
     else if (length(nGHK) != 1) {
     stop("argument 'nGHK' must be a single integer value")
     }
     else if (!is.numeric(nGHK)) {
     stop("argument 'nGHK' must be numeric")
     }
     else if (nGHK <= 0) {
     stop("argument 'nGHK' must be positive")
     }
     L <- try(t(chol(sigma)), silent = TRUE)
     if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
     }
     trunpt <- rep(NA, length(lower))
     above <- rep(NA, length(lower))
     for (i in 1:length(lower)) {
     if (lower[i] == -Inf) {
     trunpt[i] <- upper[i]
     above[i] <- 1
     }
     else if (upper[i] == Inf) {
     trunpt[i] <- lower[i]
     above[i] <- 0
     }
     else {
     stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
     " or the upper truncation point must be 'Inf'")
     }
     }
     sink(tempfile())
     on.exit(sink(), add = TRUE)
     result <- ghkvec(L = L, trunpt = trunpt, above = above,
     r = nGHK)
     result <- drop(result)
     }
     else {
     result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
     algorithm = algorithm, ...)
     }
     return(result)
    }
    <bytecode: 0x2590cb0>
    <environment: namespace:mvProbit>
     --- function search by body ---
    Function pmvnormWrap in namespace mvProbit has this body.
     ----------- END OF FAILURE REPORT --------------
    Error in if (class(L) == "try-error") { : the condition has length > 1
    Calls: mvProbit ... fn -> fnOrig -> mvProbitLogLikInternal -> pmvnormWrap
    Execution halted
Flavor: r-devel-linux-x86_64-debian-clang

Version: 0.1-8
Check: tests
Result: ERROR
     Running 'mvProbitEst.R' [29s/33s]
     Running 'mvProbitMargEff2.R' [2s/2s]
     Running 'mvProbitTest.R' [3s/3s]
     Running 'mvProbitTest5.R' [119s/130s]
     Comparing 'mvProbitTest5.Rout' to 'mvProbitTest5.Rout.save' ...1604c1604
    < 5 0.00 -0.01 0 0.00 0.00 0.00 0.00
    ---
    > 5 0.00 0.00 0 0.00 0.00 0.00 0.00
    1618,1620c1618,1620
    < 1 0 0.00 0 0 0 0 0
    < 5 0 -0.01 0 0 0 0 0
    < 10 0 0.00 0 0 0 0 0
    ---
    > 1 0 0 0 0 0 0 0
    > 5 0 0 0 0 0 0 0
    > 10 0 0 0 0 0 0 0
    1728c1728
    < 5 -0.01
    ---
    > 5 0.00
     Running 'pmvnormWrapTest.R' [2s/3s]
    Running the tests in 'tests/mvProbitEst.R' failed.
    Complete output:
     > library( "mvProbit" )
     Loading required package: mvtnorm
     Loading required package: maxLik
     Loading required package: miscTools
    
     Please cite the 'maxLik' package as:
     Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
    
     If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
     https://r-forge.r-project.org/projects/maxlik/
     Loading required package: abind
     > options( digits = 4 )
     >
     > ## generate a simulated data set
     > set.seed( 123 )
     > # number of observations
     > nObs <- 50
     >
     > # generate explanatory variables
     > xMat <- cbind(
     + const = rep( 1, nObs ),
     + x1 = as.numeric( rnorm( nObs ) > 0 ),
     + x2 = rnorm( nObs ) )
     >
     > # model coefficients
     > beta <- cbind( c( 0.8, 1.2, -0.8 ),
     + c( -0.6, 1.0, -1.6 ),
     + c( 0.5, -0.6, 1.2 ) )
     >
     > # covariance matrix of error terms
     > sigma <- miscTools::symMatrix( c( 1, 0.2, 0.4, 1, -0.1, 1 ) )
     >
     > # generate dependent variables
     > yMatLin <- xMat %*% beta
     > yMat <- ( yMatLin + rmvnorm( nObs, sigma = sigma, pre0.9_9994 = TRUE ) ) > 0
     > colnames( yMat ) <- paste( "y", 1:3, sep = "" )
     >
     > # create data frame
     > dat <- as.data.frame( cbind( xMat, yMat ) )
     >
     > # estimation with the BHHH algorithm, two-sided gradients
     > estResultBHHH <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + start = c( beta ), startSigma = sigma,
     + data = dat, tol = 0.5,
     + algorithm = GenzBretz() )
     > print( estResultBHHH )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, algorithm = GenzBretz(), tol = 0.5)
    
     Coefficients:
     b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
     1.02508 1.30536 -1.02809 -0.78248 1.14318 -1.66938 0.58229 -0.69667
     b_3_2 R_1_2 R_1_3 R_2_3
     1.23283 0.86322 0.36022 -0.08406
    
     > summary( estResultBHHH )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, algorithm = GenzBretz(), tol = 0.5)
    
     Coefficients:
     Estimate Std. error t value Pr(> t)
     b_1_0 1.02508 0.44548 2.301 0.0214 *
     b_1_1 1.30536 0.73511 1.776 0.0758 .
     b_1_2 -1.02809 0.68851 -1.493 0.1354
     b_2_0 -0.78248 0.38798 -2.017 0.0437 *
     b_2_1 1.14318 0.53026 2.156 0.0311 *
     b_2_2 -1.66938 0.81308 -2.053 0.0401 *
     b_3_0 0.58229 0.64107 0.908 0.3637
     b_3_1 -0.69667 0.89818 -0.776 0.4380
     b_3_2 1.23283 0.62452 1.974 0.0484 *
     R_1_2 0.86322 18.60071 0.046 0.9630
     R_1_3 0.36022 0.64712 0.557 0.5778
     R_2_3 -0.08406 0.34909 -0.241 0.8097
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     BHHH maximisation, 1 iterations
     Return code 2: successive function values within tolerance limit
     Log-likelihood: -56.91 on 12 Df
    
     > logLik( estResultBHHH )
     'log Lik.' -56.91 (df=12)
     > estResultBHHHA <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + start = c( beta, sigma[ lower.tri( sigma ) ] ),
     + data = dat, tol = 0.5,
     + algorithm = GenzBretz() )
     > all.equal( estResultBHHH, estResultBHHHA )
     [1] "Component \"call\": target, current do not match when deparsed"
     >
     > # estimation with the BHHH algorithm, one-sided gradients
     > estResultBHHH1 <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + start = c( beta ), startSigma = sigma,
     + data = dat, tol = 0.5,
     + algorithm = GenzBretz(), oneSidedGrad = TRUE )
     > print( estResultBHHH1 )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, algorithm = GenzBretz(), oneSidedGrad = TRUE,
     tol = 0.5)
    
     Coefficients:
     b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
     1.02508 1.30536 -1.02809 -0.78248 1.14318 -1.66938 0.58229 -0.69667
     b_3_2 R_1_2 R_1_3 R_2_3
     1.23283 0.86322 0.36022 -0.08406
    
     > summary( estResultBHHH1 )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, algorithm = GenzBretz(), oneSidedGrad = TRUE,
     tol = 0.5)
    
     Coefficients:
     Estimate Std. error t value Pr(> t)
     b_1_0 1.02508 0.44548 2.301 0.0214 *
     b_1_1 1.30536 0.73511 1.776 0.0758 .
     b_1_2 -1.02809 0.68851 -1.493 0.1354
     b_2_0 -0.78248 0.38797 -2.017 0.0437 *
     b_2_1 1.14318 0.53026 2.156 0.0311 *
     b_2_2 -1.66938 0.81307 -2.053 0.0401 *
     b_3_0 0.58229 0.64106 0.908 0.3637
     b_3_1 -0.69667 0.89817 -0.776 0.4380
     b_3_2 1.23283 0.62452 1.974 0.0484 *
     R_1_2 0.86322 18.60041 0.046 0.9630
     R_1_3 0.36022 0.64712 0.557 0.5778
     R_2_3 -0.08406 0.34909 -0.241 0.8097
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     BHHH maximisation, 1 iterations
     Return code 2: successive function values within tolerance limit
     Log-likelihood: -56.91 on 12 Df
    
     > logLik( estResultBHHH1 )
     'log Lik.' -56.91 (df=12)
     > all.equal( estResultBHHH, estResultBHHH1, tol = 1e-5 )
     [1] "Component \"call\": target, current do not match when deparsed"
     >
     > # estimation with the BFGS algorithm, two-sided gradients
     > estResultBFGS <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + start = c( beta ), startSigma = sigma, method = "BFGS",
     + data = dat,
     + reltol = 0.5, algorithm = GenzBretz() )
     > print( estResultBFGS )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, method = "BFGS", algorithm = GenzBretz(),
     reltol = 0.5)
    
     Coefficients:
     b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
     0.53607 1.16101 -1.14212 -0.44952 1.05496 -1.57624 0.40751 -0.85400
     b_3_2 R_1_2 R_1_3 R_2_3
     1.25958 0.29483 0.32822 -0.02382
    
     > summary( estResultBFGS )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, method = "BFGS", algorithm = GenzBretz(),
     reltol = 0.5)
    
     Coefficients:
     Estimate Std. error t value Pr(> t)
     b_1_0 0.53607 0.82709 0.648 0.5169
     b_1_1 1.16101 0.91515 1.269 0.2046
     b_1_2 -1.14212 0.93847 -1.217 0.2236
     b_2_0 -0.44952 0.79215 -0.567 0.5704
     b_2_1 1.05496 0.88398 1.193 0.2327
     b_2_2 -1.57624 0.73187 -2.154 0.0313 *
     b_3_0 0.40751 0.48897 0.833 0.4046
     b_3_1 -0.85400 0.84871 -1.006 0.3143
     b_3_2 1.25958 0.68040 1.851 0.0641 .
     R_1_2 0.29483 1.85817 0.159 0.8739
     R_1_3 0.32822 0.70923 0.463 0.6435
     R_2_3 -0.02382 0.34466 -0.069 0.9449
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     BFGS maximization, 4 iterations
     Return code 0: successful convergence
     Log-likelihood: -52.94 on 12 Df
    
     > logLik( estResultBFGS )
     'log Lik.' -52.94 (df=12)
     >
     > # estimation with the BFGS algorithm, one-sided gradients
     > estResultBFGS1 <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + start = c( beta ), startSigma = sigma, method = "BFGS",
     + data = dat,
     + reltol = 0.5, algorithm = GenzBretz(), oneSidedGrad = TRUE )
     > print( estResultBFGS1 )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, method = "BFGS", algorithm = GenzBretz(),
     oneSidedGrad = TRUE, reltol = 0.5)
    
     Coefficients:
     b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
     0.53607 1.16101 -1.14213 -0.44952 1.05496 -1.57624 0.40751 -0.85400
     b_3_2 R_1_2 R_1_3 R_2_3
     1.25958 0.29483 0.32822 -0.02382
    
     > summary( estResultBFGS1 )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, method = "BFGS", algorithm = GenzBretz(),
     oneSidedGrad = TRUE, reltol = 0.5)
    
     Coefficients:
     Estimate Std. error t value Pr(> t)
     b_1_0 0.53607 0.82709 0.648 0.5169
     b_1_1 1.16101 0.91515 1.269 0.2046
     b_1_2 -1.14213 0.93847 -1.217 0.2236
     b_2_0 -0.44952 0.79215 -0.567 0.5704
     b_2_1 1.05496 0.88398 1.193 0.2327
     b_2_2 -1.57624 0.73187 -2.154 0.0313 *
     b_3_0 0.40751 0.48897 0.833 0.4046
     b_3_1 -0.85400 0.84871 -1.006 0.3143
     b_3_2 1.25958 0.68040 1.851 0.0641 .
     R_1_2 0.29483 1.85817 0.159 0.8739
     R_1_3 0.32822 0.70923 0.463 0.6435
     R_2_3 -0.02382 0.34466 -0.069 0.9449
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     BFGS maximization, 4 iterations
     Return code 0: successful convergence
     Log-likelihood: -52.94 on 12 Df
    
     > logLik( estResultBFGS1 )
     'log Lik.' -52.94 (df=12)
     > all.equal( estResultBFGS, estResultBFGS1, tol = 1e-5 )
     [1] "Component \"call\": target, current do not match when deparsed"
     >
     > # estimation with the BFGS algorithm, one-sided gradients, no starting values
     > estResultBFGS1a <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + data = dat, method = "BFGS",
     + reltol = 0.5, algorithm = GenzBretz(), oneSidedGrad = TRUE )
     > print( estResultBFGS1a )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, method = "BFGS",
     algorithm = GenzBretz(), oneSidedGrad = TRUE, reltol = 0.5)
    
     Coefficients:
     b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
     0.79353 1.88367 -1.74922 -0.40036 0.91577 -1.52106 0.97790 -1.69433
     b_3_2 R_1_2 R_1_3 R_2_3
     1.57625 0.75270 0.23721 -0.03453
    
     > summary( estResultBFGS1a )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, method = "BFGS",
     algorithm = GenzBretz(), oneSidedGrad = TRUE, reltol = 0.5)
    
     Coefficients:
     Estimate Std. error t value Pr(> t)
     b_1_0 0.79353 0.90412 0.878 0.3801
     b_1_1 1.88367 0.90668 2.078 0.0378 *
     b_1_2 -1.74922 0.97400 -1.796 0.0725 .
     b_2_0 -0.40036 0.69056 -0.580 0.5621
     b_2_1 0.91577 0.97998 0.934 0.3501
     b_2_2 -1.52106 0.72006 -2.112 0.0347 *
     b_3_0 0.97790 0.64497 1.516 0.1295
     b_3_1 -1.69433 0.94788 -1.787 0.0739 .
     b_3_2 1.57625 0.73399 2.148 0.0318 *
     R_1_2 0.75270 4.32729 0.174 0.8619
     R_1_3 0.23721 0.90535 0.262 0.7933
     R_2_3 -0.03453 0.35889 -0.096 0.9234
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     BFGS maximization, 3 iterations
     Return code 0: successful convergence
     Log-likelihood: -49.36 on 12 Df
    
     > logLik( estResultBFGS1a )
     'log Lik.' -49.36 (df=12)
     >
     > # estimation with the BFGS algorithm, one-sided gradients, no starting values for beta
     > estResultBFGS1b <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + startSigma = sigma, data = dat,
     + method = "BFGS", reltol = 0.5, algorithm = GenzBretz(), oneSidedGrad = TRUE )
     > print( estResultBFGS1b )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, startSigma = sigma,
     method = "BFGS", algorithm = GenzBretz(), oneSidedGrad = TRUE,
     reltol = 0.5)
    
     Coefficients:
     b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1 b_3_2 R_1_2
     0.8110 1.9060 -1.7627 -0.3770 0.9441 -1.5332 0.9372 -1.7498 1.6024 0.6817
     R_1_3 R_2_3
     0.4044 0.1663
    
     > summary( estResultBFGS1b )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, startSigma = sigma,
     method = "BFGS", algorithm = GenzBretz(), oneSidedGrad = TRUE,
     reltol = 0.5)
    
     Coefficients:
     Estimate Std. error t value Pr(> t)
     b_1_0 0.8110 0.9297 0.872 0.3830
     b_1_1 1.9060 0.9160 2.081 0.0375 *
     b_1_2 -1.7627 1.0131 -1.740 0.0819 .
     b_2_0 -0.3770 0.7049 -0.535 0.5928
     b_2_1 0.9441 1.0345 0.913 0.3615
     b_2_2 -1.5332 0.7517 -2.040 0.0414 *
     b_3_0 0.9372 0.6291 1.490 0.1363
     b_3_1 -1.7498 0.9996 -1.751 0.0800 .
     b_3_2 1.6024 0.7848 2.042 0.0412 *
     R_1_2 0.6817 3.4094 0.200 0.8415
     R_1_3 0.4044 0.8024 0.504 0.6143
     R_2_3 0.1663 0.3389 0.491 0.6235
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     BFGS maximization, 3 iterations
     Return code 0: successful convergence
     Log-likelihood: -49.57 on 12 Df
    
     > logLik( estResultBFGS1b )
     'log Lik.' -49.57 (df=12)
     >
     > # estimation with the BFGS algorithm, one-sided gradients, no starting values for sigma
     > estResultBFGS1s <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + start = c( beta ), data = dat,
     + method = "BFGS", reltol = 0.5, algorithm = GenzBretz(), oneSidedGrad = TRUE )
     > print( estResultBFGS1s )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     method = "BFGS", algorithm = GenzBretz(), oneSidedGrad = TRUE,
     reltol = 0.5)
    
     Coefficients:
     b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
     0.55814 1.14958 -1.12747 -0.44426 1.07325 -1.57226 0.37255 -0.85470
     b_3_2 R_1_2 R_1_3 R_2_3
     1.23011 0.34981 -0.01177 0.09435
    
     > summary( estResultBFGS1s )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     method = "BFGS", algorithm = GenzBretz(), oneSidedGrad = TRUE,
     reltol = 0.5)
    
     Coefficients:
     Estimate Std. error t value Pr(> t)
     b_1_0 0.55814 0.86615 0.644 0.5193
     b_1_1 1.14958 0.90607 1.269 0.2045
     b_1_2 -1.12747 0.94191 -1.197 0.2313
     b_2_0 -0.44426 0.85526 -0.519 0.6035
     b_2_1 1.07325 0.99314 1.081 0.2798
     b_2_2 -1.57226 0.78557 -2.001 0.0453 *
     b_3_0 0.37255 0.56172 0.663 0.5072
     b_3_1 -0.85470 0.95214 -0.898 0.3694
     b_3_2 1.23011 0.71506 1.720 0.0854 .
     R_1_2 0.34981 2.02289 0.173 0.8627
     R_1_3 -0.01177 0.70776 -0.017 0.9867
     R_2_3 0.09435 0.33920 0.278 0.7809
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     BFGS maximization, 4 iterations
     Return code 0: successful convergence
     Log-likelihood: -52.78 on 12 Df
    
     > logLik( estResultBFGS1s )
     'log Lik.' -52.78 (df=12)
     >
     > # estimation with the BFGS algorithm, Miwa algorithm for obtaining integrals
     > estResultBFGSm <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + start = c( beta ), startSigma = sigma,
     + data = dat, method = "BFGS",
     + reltol = 0.5, algorithm = Miwa( steps = 64 ) )
     There were 13 warnings (use warnings() to see them)
     > print( estResultBFGSm )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, method = "BFGS", algorithm = Miwa(steps = 64),
     reltol = 0.5)
    
     Coefficients:
     b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
     0.53609 1.16102 -1.14213 -0.44952 1.05496 -1.57624 0.40753 -0.85400
     b_3_2 R_1_2 R_1_3 R_2_3
     1.25960 0.29491 0.32819 -0.02386
    
     > summary( estResultBFGSm )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, method = "BFGS", algorithm = Miwa(steps = 64),
     reltol = 0.5)
    
     Coefficients:
     Estimate Std. error t value Pr(> t)
     b_1_0 0.53609 0.82702 0.648 0.5168
     b_1_1 1.16102 0.91511 1.269 0.2045
     b_1_2 -1.14213 0.93837 -1.217 0.2235
     b_2_0 -0.44952 0.79215 -0.567 0.5704
     b_2_1 1.05496 0.88386 1.194 0.2326
     b_2_2 -1.57624 0.73173 -2.154 0.0312 *
     b_3_0 0.40753 0.48894 0.834 0.4046
     b_3_1 -0.85400 0.84860 -1.006 0.3142
     b_3_2 1.25960 0.68024 1.852 0.0641 .
     R_1_2 0.29491 1.85789 0.159 0.8739
     R_1_3 0.32819 0.70903 0.463 0.6435
     R_2_3 -0.02386 0.34468 -0.069 0.9448
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     BFGS maximization, 4 iterations
     Return code 0: successful convergence
     Log-likelihood: -52.94 on 12 Df
    
     > logLik( estResultBFGSm )
     'log Lik.' -52.94 (df=12)
     > all.equal( estResultBFGS, estResultBFGSm, tol = 1e-3 )
     [1] "Component \"call\": target, current do not match when deparsed"
     >
     > # estimation with the BFGS algorithm, GHK algorithm for obtaining integrals
     > estResultBFGSg <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + start = c( beta ), startSigma = sigma,
     + data = dat, method = "BFGS",
     + reltol = 0.5 )
     ----------- FAILURE REPORT --------------
     --- failure: the condition has length > 1 ---
     --- srcref ---
     :
     --- package (from environment) ---
     mvProbit
     --- call from context ---
     pmvnormWrap(upper = xBetaTmp, sigma = sigmaTmp, algorithm = algorithm,
     nGHK = nGHK, random.seed = randomSeed, ...)
     --- call from argument ---
     if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
     }
     --- R stacktrace ---
     where 1: pmvnormWrap(upper = xBetaTmp, sigma = sigmaTmp, algorithm = algorithm,
     nGHK = nGHK, random.seed = randomSeed, ...)
     where 2: mvProbitLogLikInternal(yMat = yMat, xMat = xMat, coef = param,
     sigma = NULL, algorithm = llAlgorithm, nGHK = llNGHK, returnGrad = llIntGrad,
     oneSidedGrad = llOneSidedGrad, eps = llEps, randomSeed = llRandom.seed,
     ...)
     where 3: fnOrig(theta, ...)
     where 4: logLikFunc(theta, fnOrig = function (param, yMat, xMat, llAlgorithm,
     llNGHK, llIntGrad, llOneSidedGrad, llEps, llRandom.seed,
     ...)
     {
     logLikVal <- mvProbitLogLikInternal(yMat = yMat, xMat = xMat,
     coef = param, sigma = NULL, algorithm = llAlgorithm,
     nGHK = llNGHK, returnGrad = llIntGrad, oneSidedGrad = llOneSidedGrad,
     eps = llEps, randomSeed = llRandom.seed, ...)
     return(logLikVal)
     }, gradOrig = ..2, hessOrig = ..3, yMat = c(1, 1, 1, 1, 1, 1,
     1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1,
     1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0,
     1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0,
     0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0,
     0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1,
     0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1,
     0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0), xMat = c(1,
     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1,
     0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1,
     1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0.253318513994755,
     -0.028546755348703, -0.0428704572913161, 1.36860228401446, -0.225770985659268,
     1.51647060442954, -1.54875280423022, 0.584613749636069, 0.123854243844614,
     0.215941568743973, 0.379639482759882, -0.502323453109302, -0.33320738366942,
     -1.01857538310709, -1.07179122647558, 0.303528641404258, 0.448209778629426,
     0.0530042267305041, 0.922267467879738, 2.05008468562714, -0.491031166056535,
     -2.30916887564081, 1.00573852446226, -0.709200762582393, -0.688008616467358,
     1.0255713696967, -0.284773007051009, -1.22071771225454, 0.18130347974915,
     -0.138891362439045, 0.00576418589988693, 0.38528040112633, -0.370660031792409,
     0.644376548518833, -0.220486561818751, 0.331781963915697, 1.09683901314935,
     0.435181490833803, -0.325931585531227, 1.14880761845109, 0.993503855962119,
     0.54839695950807, 0.238731735111441, -0.627906076039371, 1.36065244853001,
     -0.600259587147127, 2.18733299301658, 1.53261062618519, -0.235700359100477,
     -1.02642090030678), llAlgorithm = "GHK", llNGHK = 1000, llIntGrad = TRUE,
     llOneSidedGrad = FALSE, llEps = 1e-06, llRandom.seed = 123)
     where 5: eval(f, sys.frame(sys.parent()))
     where 6: eval(f, sys.frame(sys.parent()))
     where 7: callWithoutArgs(theta, fName = fName, args = names(formals(sumt)),
     ...)
     where 8: (function (theta, fName, ...)
     {
     return(callWithoutArgs(theta, fName = fName, args = names(formals(sumt)),
     ...))
     })(theta = c(b_1_0 = 0.8, b_1_1 = 1.2, b_1_2 = -0.8, b_2_0 = -0.6,
     b_2_1 = 1, b_2_2 = -1.6, b_3_0 = 0.5, b_3_1 = -0.6, b_3_2 = 1.2,
     R_1_2 = 0.2, R_1_3 = 0.4, R_2_3 = -0.1), fName = "logLikFunc",
     fnOrig = function (param, yMat, xMat, llAlgorithm, llNGHK,
     llIntGrad, llOneSidedGrad, llEps, llRandom.seed, ...)
     {
     logLikVal <- mvProbitLogLikInternal(yMat = yMat, xMat = xMat,
     coef = param, sigma = NULL, algorithm = llAlgorithm,
     nGHK = llNGHK, returnGrad = llIntGrad, oneSidedGrad = llOneSidedGrad,
     eps = llEps, randomSeed = llRandom.seed, ...)
     return(logLikVal)
     }, gradOrig = NULL, hessOrig = NULL, yMat = c(1, 1, 1, 1,
     1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0,
     1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0,
     0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1,
     1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1,
     0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1,
     0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0,
     0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1,
     1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0), xMat = c(1, 1, 1,
     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0,
     1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0,
     1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0,
     1, 0, 0.253318513994755, -0.028546755348703, -0.0428704572913161,
     1.36860228401446, -0.225770985659268, 1.51647060442954, -1.54875280423022,
     0.584613749636069, 0.123854243844614, 0.215941568743973,
     0.379639482759882, -0.502323453109302, -0.33320738366942,
     -1.01857538310709, -1.07179122647558, 0.303528641404258,
     0.448209778629426, 0.0530042267305041, 0.922267467879738,
     2.05008468562714, -0.491031166056535, -2.30916887564081,
     1.00573852446226, -0.709200762582393, -0.688008616467358,
     1.0255713696967, -0.284773007051009, -1.22071771225454, 0.18130347974915,
     -0.138891362439045, 0.00576418589988693, 0.38528040112633,
     -0.370660031792409, 0.644376548518833, -0.220486561818751,
     0.331781963915697, 1.09683901314935, 0.435181490833803, -0.325931585531227,
     1.14880761845109, 0.993503855962119, 0.54839695950807, 0.238731735111441,
     -0.627906076039371, 1.36065244853001, -0.600259587147127,
     2.18733299301658, 1.53261062618519, -0.235700359100477, -1.02642090030678
     ), llAlgorithm = "GHK", llNGHK = 1000, llIntGrad = TRUE,
     llOneSidedGrad = FALSE, llEps = 1e-06, llRandom.seed = 123)
     where 9: do.call(callWithoutSumt, argList)
     where 10: maxOptim(fn = fn, grad = grad, hess = hess, start = start, method = "BFGS",
     fixed = fixed, constraints = constraints, finalHessian = finalHessian,
     parscale = parscale, control = mControl, ...)
     where 11: maxRoutine(fn = logLik, grad = grad, hess = hess, start = start,
     constraints = constraints, ...)
     where 12: maxLik(logLik = logLik, start = start, method = method, finalHessian = finalHessian,
     yMat = yMat, xMat = xMat, llAlgorithm = algorithm, llNGHK = nGHK,
     llIntGrad = intGrad, llOneSidedGrad = oneSidedGrad, llEps = eps,
     llRandom.seed = random.seed, ...)
     where 13: mvProbit(cbind(y1, y2, y3) ~ x1 + x2, start = c(beta), startSigma = sigma,
     data = dat, method = "BFGS", reltol = 0.5)
    
     --- value of length: 2 type: logical ---
     [1] FALSE FALSE
     --- function from context ---
     function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
     nGHK = NULL, ...)
     {
     if (!is.matrix(sigma)) {
     stop("argument 'sigma' must be a matrix")
     }
     else if (nrow(sigma) != ncol(sigma)) {
     stop("argument 'sigma' must be a square matrix")
     }
     else if (!all(is.numeric(sigma))) {
     stop("argument 'sigma must be a numeric matrix")
     }
     if (!all(is.numeric(lower))) {
     stop("argument 'lower' must be numeric")
     }
     else if (length(lower) == 1) {
     lower <- rep(lower, nrow(sigma))
     }
     else if (length(lower) != nrow(sigma)) {
     stop("argument 'lower' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     if (!all(is.numeric(upper))) {
     stop("argument 'upper' must be numeric")
     }
     else if (length(upper) == 1) {
     upper <- rep(upper, nrow(sigma))
     }
     else if (length(upper) != nrow(sigma)) {
     stop("argument 'upper' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     algOkay <- TRUE
     ghk <- FALSE
     if (!is.list(algorithm) && length(algorithm) != 1) {
     stop("argument 'algorithm' must be a single function",
     " or a single character string")
     }
     else if (is.function(algorithm)) {
     algResult <- do.call(algorithm, list())
     if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
     "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (is.character(algorithm)) {
     if (tolower(algorithm) == "ghk") {
     ghk <- TRUE
     }
     else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     if (!algOkay) {
     stop("argument 'algorithm' must be either one of the functions",
     " 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
     " \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
     }
     if (length(random.seed) != 1) {
     stop("argument 'random.seed' must be a single numerical value")
     }
     else if (!is.numeric(random.seed)) {
     stop("argument 'random.seed' must be numerical")
     }
     if (exists(".Random.seed")) {
     savedSeed <- .Random.seed
     }
     set.seed(random.seed)
     if (exists("savedSeed")) {
     on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
     }
     else {
     on.exit(rm(.Random.seed, envir = sys.frame()))
     }
     if (ghk) {
     if (is.null(nGHK)) {
     stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
     }
     else if (length(nGHK) != 1) {
     stop("argument 'nGHK' must be a single integer value")
     }
     else if (!is.numeric(nGHK)) {
     stop("argument 'nGHK' must be numeric")
     }
     else if (nGHK <= 0) {
     stop("argument 'nGHK' must be positive")
     }
     L <- try(t(chol(sigma)), silent = TRUE)
     if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
     }
     trunpt <- rep(NA, length(lower))
     above <- rep(NA, length(lower))
     for (i in 1:length(lower)) {
     if (lower[i] == -Inf) {
     trunpt[i] <- upper[i]
     above[i] <- 1
     }
     else if (upper[i] == Inf) {
     trunpt[i] <- lower[i]
     above[i] <- 0
     }
     else {
     stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
     " or the upper truncation point must be 'Inf'")
     }
     }
     sink(tempfile())
     on.exit(sink(), add = TRUE)
     result <- ghkvec(L = L, trunpt = trunpt, above = above,
     r = nGHK)
     result <- drop(result)
     }
     else {
     result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
     algorithm = algorithm, ...)
     }
     return(result)
     }
     <bytecode: 0x2a96c80>
     <environment: namespace:mvProbit>
     --- function search by body ---
     Function pmvnormWrap in namespace mvProbit has this body.
     ----------- END OF FAILURE REPORT --------------
     Error in if (class(L) == "try-error") { : the condition has length > 1
     Calls: mvProbit ... logLikFunc -> fnOrig -> mvProbitLogLikInternal -> pmvnormWrap
     Execution halted
    Running the tests in 'tests/mvProbitMargEff2.R' failed.
    Complete output:
     > # I thank Mohit Batham for providing this R script that demonstrated
     > # a bug in mvProbitMargEff() when called with 2 dependent variables
     > library( "mvProbit" )
     Loading required package: mvtnorm
     Loading required package: maxLik
     Loading required package: miscTools
    
     Please cite the 'maxLik' package as:
     Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
    
     If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
     https://r-forge.r-project.org/projects/maxlik/
     Loading required package: abind
     > nObs <- 100
     > set.seed( 123 )
     > xData <- data.frame(
     + const = rep( 1, nObs ),
     + x1 = as.numeric( rnorm( nObs ) > 0 ),
     + x2 = as.numeric( rnorm( nObs ) > 0 ),
     + x3 = rnorm( nObs ),
     + x4 = rnorm( nObs ))
     > beta <- c( 0.8, 1.2, -1.0, 1.4, -0.8,
     + -0.6, 1.0, 0.6, -1.2, -1.6)
     > sigma <- symMatrix( c( 1, 0.2, 1))
     > margEffCond <- try( mvProbitMargEff( ~ x1 + x2 + x3 + x4 , coef = beta,
     + sigma = sigma, data = xData, cond = TRUE ) )
     ----------- FAILURE REPORT --------------
     --- failure: the condition has length > 1 ---
     --- srcref ---
     :
     --- package (from environment) ---
     mvProbit
     --- call from context ---
     pmvnormWrap(upper = xBeta[i, ], sigma = coef$sigma, algorithm = algorithm,
     nGHK = nGHK, random.seed = random.seed, ...)
     --- call from argument ---
     if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
     }
     --- R stacktrace ---
     where 1: pmvnormWrap(upper = xBeta[i, ], sigma = coef$sigma, algorithm = algorithm,
     nGHK = nGHK, random.seed = random.seed, ...)
     where 2: mvProbitExpInternal(yMat = yMat, xMat = xMat, coef = betaLower,
     sigma = coef$sigma, cond = cond, algorithm = algorithm, nGHK = nGHK,
     random.seed = random.seed, ...)
     where 3: mvProbitMargEffInternal(yMat = yMat, xMat = xMat, coef = coef,
     sigma = sigma, cond = cond, algorithm = algorithm, nGHK = nGHK,
     eps = eps, dummyVars = dummyVars, random.seed = random.seed,
     ...)
     where 4: mvProbitMargEff(~x1 + x2 + x3 + x4, coef = beta, sigma = sigma,
     data = xData, cond = TRUE)
     where 5: doTryCatch(return(expr), name, parentenv, handler)
     where 6: tryCatchOne(expr, names, parentenv, handlers[[1L]])
     where 7: tryCatchList(expr, classes, parentenv, handlers)
     where 8: 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 9: try(mvProbitMargEff(~x1 + x2 + x3 + x4, coef = beta, sigma = sigma,
     data = xData, cond = TRUE))
    
     --- value of length: 2 type: logical ---
     [1] FALSE FALSE
     --- function from context ---
     function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
     nGHK = NULL, ...)
     {
     if (!is.matrix(sigma)) {
     stop("argument 'sigma' must be a matrix")
     }
     else if (nrow(sigma) != ncol(sigma)) {
     stop("argument 'sigma' must be a square matrix")
     }
     else if (!all(is.numeric(sigma))) {
     stop("argument 'sigma must be a numeric matrix")
     }
     if (!all(is.numeric(lower))) {
     stop("argument 'lower' must be numeric")
     }
     else if (length(lower) == 1) {
     lower <- rep(lower, nrow(sigma))
     }
     else if (length(lower) != nrow(sigma)) {
     stop("argument 'lower' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     if (!all(is.numeric(upper))) {
     stop("argument 'upper' must be numeric")
     }
     else if (length(upper) == 1) {
     upper <- rep(upper, nrow(sigma))
     }
     else if (length(upper) != nrow(sigma)) {
     stop("argument 'upper' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     algOkay <- TRUE
     ghk <- FALSE
     if (!is.list(algorithm) && length(algorithm) != 1) {
     stop("argument 'algorithm' must be a single function",
     " or a single character string")
     }
     else if (is.function(algorithm)) {
     algResult <- do.call(algorithm, list())
     if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
     "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (is.character(algorithm)) {
     if (tolower(algorithm) == "ghk") {
     ghk <- TRUE
     }
     else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     if (!algOkay) {
     stop("argument 'algorithm' must be either one of the functions",
     " 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
     " \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
     }
     if (length(random.seed) != 1) {
     stop("argument 'random.seed' must be a single numerical value")
     }
     else if (!is.numeric(random.seed)) {
     stop("argument 'random.seed' must be numerical")
     }
     if (exists(".Random.seed")) {
     savedSeed <- .Random.seed
     }
     set.seed(random.seed)
     if (exists("savedSeed")) {
     on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
     }
     else {
     on.exit(rm(.Random.seed, envir = sys.frame()))
     }
     if (ghk) {
     if (is.null(nGHK)) {
     stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
     }
     else if (length(nGHK) != 1) {
     stop("argument 'nGHK' must be a single integer value")
     }
     else if (!is.numeric(nGHK)) {
     stop("argument 'nGHK' must be numeric")
     }
     else if (nGHK <= 0) {
     stop("argument 'nGHK' must be positive")
     }
     L <- try(t(chol(sigma)), silent = TRUE)
     if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
     }
     trunpt <- rep(NA, length(lower))
     above <- rep(NA, length(lower))
     for (i in 1:length(lower)) {
     if (lower[i] == -Inf) {
     trunpt[i] <- upper[i]
     above[i] <- 1
     }
     else if (upper[i] == Inf) {
     trunpt[i] <- lower[i]
     above[i] <- 0
     }
     else {
     stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
     " or the upper truncation point must be 'Inf'")
     }
     }
     sink(tempfile())
     on.exit(sink(), add = TRUE)
     result <- ghkvec(L = L, trunpt = trunpt, above = above,
     r = nGHK)
     result <- drop(result)
     }
     else {
     result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
     algorithm = algorithm, ...)
     }
     return(result)
     }
     <bytecode: 0x282fba8>
     <environment: namespace:mvProbit>
     --- function search by body ---
     Function pmvnormWrap in namespace mvProbit has this body.
     ----------- END OF FAILURE REPORT --------------
     Error in if (class(L) == "try-error") { : the condition has length > 1
     > round( margEffCond, 3 )
     Error in round(margEffCond, 3) :
     non-numeric argument to mathematical function
     Execution halted
    Running the tests in 'tests/mvProbitTest.R' failed.
    Complete output:
     > library( "mvProbit" )
     Loading required package: mvtnorm
     Loading required package: maxLik
     Loading required package: miscTools
    
     Please cite the 'maxLik' package as:
     Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
    
     If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
     https://r-forge.r-project.org/projects/maxlik/
     Loading required package: abind
     > options( digits = 4 )
     >
     > ## generate a simulated data set
     > set.seed( 123 )
     > # number of observations
     > nObs <- 10
     >
     > # generate explanatory variables
     > xMat <- cbind(
     + const = rep( 1, nObs ),
     + x1 = as.numeric( rnorm( nObs ) > 0 ),
     + x2 = as.numeric( rnorm( nObs ) > 0 ),
     + x3 = rnorm( nObs ),
     + x4 = rnorm( nObs ) )
     >
     > # model coefficients
     > beta <- cbind( c( 0.8, 1.2, -1.0, 1.4, -0.8 ),
     + c( -0.6, 1.0, 0.6, -1.2, -1.6 ),
     + c( 0.5, -0.6, -0.7, 1.1, 1.2 ) )
     >
     > # covariance matrix of error terms
     > sigma <- miscTools::symMatrix( c( 1, 0.2, 0.4, 1, -0.1, 1 ) )
     >
     > # all parameters in a vector
     > allCoef <- c( c( beta ), sigma[ lower.tri( sigma ) ] )
     >
     > # generate dependent variables
     > yMatLin <- xMat %*% beta
     > yMat <- ( yMatLin + rmvnorm( nObs, sigma = sigma, pre0.9_9994 = TRUE ) ) > 0
     > colnames( yMat ) <- paste( "y", 1:3, sep = "" )
     > # (yMatLin > 0 )== yMat
     >
     > # unconditional expectations of dependent variables
     > yExp <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
     + sigma = sigma, data = as.data.frame( xMat ) )
     > round( yExp, 3 )
     V1 V2 V3
     1 0.021 0.725 0.194
     2 0.394 0.768 0.214
     3 0.125 0.788 0.196
     4 0.235 0.681 0.292
     5 0.680 0.435 0.579
     6 0.028 0.973 0.034
     7 0.958 0.186 0.784
     8 0.856 0.247 0.724
     9 0.061 0.968 0.034
     10 0.998 0.067 0.923
     > yExpA <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = allCoef,
     + data = as.data.frame( xMat ) )
     > all.equal( yExp, yExpA )
     [1] TRUE
     > yExp2 <- pnorm( yMatLin )
     > all.equal( yExp, as.data.frame( yExp2 ) )
     [1] TRUE
     >
     >
     > # conditional expectations of dependent variables
     > # (assuming that all other dependent variables are one)
     > yExpCond <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
     + sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
     + algorithm = GenzBretz() )
     > round( yExpCond, 3 )
     V1 V2 V3
     1 0.072 0.834 0.524
     2 0.660 0.779 0.314
     3 0.297 0.838 0.399
     4 0.448 0.727 0.461
     5 0.847 0.442 0.618
     6 0.139 0.985 0.162
     7 0.991 0.179 0.749
     8 0.950 0.244 0.710
     9 0.244 0.978 0.133
     10 1.000 0.065 0.892
     > yExpCondA <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = allCoef,
     + data = as.data.frame( xMat ), cond = TRUE,
     + algorithm = GenzBretz() )
     > all.equal( yExpCond, yExpCondA )
     [1] TRUE
     > yExpCond2 <- matrix( NA, nrow = nObs, ncol = ncol( yMat ) )
     > for( i in 1:nObs ) {
     + for( k in 1:ncol( yMat ) ) {
     + set.seed( 123 )
     + numerator <- pmvnorm( upper = yMatLin[ i, ], sigma = sigma )
     + set.seed( 123 )
     + denominator <- pmvnorm( upper = yMatLin[ i, -k ], sigma = sigma[ -k, -k ] )
     + yExpCond2[ i, k ] <- numerator / denominator
     + }
     + }
     > all.equal( yExpCond, as.data.frame( yExpCond2 ) )
     [1] TRUE
     > # now with explicitly specifying the algorithm
     > yExpCond3 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
     + sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
     + algorithm = GenzBretz )
     > all.equal( yExpCond, yExpCond3 )
     [1] TRUE
     > identical( yExpCond, yExpCond3 )
     [1] TRUE
     > # now with integrals obtained by the Miwa algorithm
     > yExpCond4 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
     + sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
     + algorithm = Miwa )
     > all.equal( yExpCond, yExpCond4, tol = 1e-3 )
     [1] TRUE
     > # now with integrals obtained by the Miwa algorithm, less precise
     > yExpCond5 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
     + sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
     + algorithm = Miwa( steps = 32 ) )
     > all.equal( yExpCond4, yExpCond5, tol = 1e-3 )
     [1] TRUE
     > # now with integrals obtained by the TVPACK algorithm
     > yExpCond6 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
     + sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
     + algorithm = TVPACK )
     > all.equal( yExpCond, yExpCond6, tol = 1e-3 )
     [1] TRUE
     > # now with integrals obtained by the TVPACK algorithm, less precise
     > yExpCond7 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
     + sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
     + algorithm = TVPACK( abseps = 0.5 ) )
     > all.equal( yExpCond6, yExpCond7, tol = 1e-3 )
     [1] "Component \"V1\": Mean relative difference: 0.03018"
     [2] "Component \"V2\": Mean relative difference: 0.04273"
     [3] "Component \"V3\": Mean relative difference: 0.04242"
     > # now with integrals obtained by the GHK algorithm
     > yExpCond8 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
     + sigma = sigma, data = as.data.frame( xMat ), cond = TRUE )
     ----------- FAILURE REPORT --------------
     --- failure: the condition has length > 1 ---
     --- srcref ---
     :
     --- package (from environment) ---
     mvProbit
     --- call from context ---
     pmvnormWrap(upper = xBeta[i, ], sigma = coef$sigma, algorithm = algorithm,
     nGHK = nGHK, random.seed = random.seed, ...)
     --- call from argument ---
     if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
     }
     --- R stacktrace ---
     where 1: pmvnormWrap(upper = xBeta[i, ], sigma = coef$sigma, algorithm = algorithm,
     nGHK = nGHK, random.seed = random.seed, ...)
     where 2: mvProbitExpInternal(yMat = yMat, xMat = xMat, coef = coef, sigma = sigma,
     cond = cond, algorithm = algorithm, nGHK = nGHK, random.seed = random.seed,
     ...)
     where 3: mvProbitExp(~x1 + x2 + x3 + x4, coef = c(beta), sigma = sigma,
     data = as.data.frame(xMat), cond = TRUE)
    
     --- value of length: 2 type: logical ---
     [1] FALSE FALSE
     --- function from context ---
     function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
     nGHK = NULL, ...)
     {
     if (!is.matrix(sigma)) {
     stop("argument 'sigma' must be a matrix")
     }
     else if (nrow(sigma) != ncol(sigma)) {
     stop("argument 'sigma' must be a square matrix")
     }
     else if (!all(is.numeric(sigma))) {
     stop("argument 'sigma must be a numeric matrix")
     }
     if (!all(is.numeric(lower))) {
     stop("argument 'lower' must be numeric")
     }
     else if (length(lower) == 1) {
     lower <- rep(lower, nrow(sigma))
     }
     else if (length(lower) != nrow(sigma)) {
     stop("argument 'lower' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     if (!all(is.numeric(upper))) {
     stop("argument 'upper' must be numeric")
     }
     else if (length(upper) == 1) {
     upper <- rep(upper, nrow(sigma))
     }
     else if (length(upper) != nrow(sigma)) {
     stop("argument 'upper' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     algOkay <- TRUE
     ghk <- FALSE
     if (!is.list(algorithm) && length(algorithm) != 1) {
     stop("argument 'algorithm' must be a single function",
     " or a single character string")
     }
     else if (is.function(algorithm)) {
     algResult <- do.call(algorithm, list())
     if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
     "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (is.character(algorithm)) {
     if (tolower(algorithm) == "ghk") {
     ghk <- TRUE
     }
     else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     if (!algOkay) {
     stop("argument 'algorithm' must be either one of the functions",
     " 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
     " \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
     }
     if (length(random.seed) != 1) {
     stop("argument 'random.seed' must be a single numerical value")
     }
     else if (!is.numeric(random.seed)) {
     stop("argument 'random.seed' must be numerical")
     }
     if (exists(".Random.seed")) {
     savedSeed <- .Random.seed
     }
     set.seed(random.seed)
     if (exists("savedSeed")) {
     on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
     }
     else {
     on.exit(rm(.Random.seed, envir = sys.frame()))
     }
     if (ghk) {
     if (is.null(nGHK)) {
     stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
     }
     else if (length(nGHK) != 1) {
     stop("argument 'nGHK' must be a single integer value")
     }
     else if (!is.numeric(nGHK)) {
     stop("argument 'nGHK' must be numeric")
     }
     else if (nGHK <= 0) {
     stop("argument 'nGHK' must be positive")
     }
     L <- try(t(chol(sigma)), silent = TRUE)
     if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
     }
     trunpt <- rep(NA, length(lower))
     above <- rep(NA, length(lower))
     for (i in 1:length(lower)) {
     if (lower[i] == -Inf) {
     trunpt[i] <- upper[i]
     above[i] <- 1
     }
     else if (upper[i] == Inf) {
     trunpt[i] <- lower[i]
     above[i] <- 0
     }
     else {
     stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
     " or the upper truncation point must be 'Inf'")
     }
     }
     sink(tempfile())
     on.exit(sink(), add = TRUE)
     result <- ghkvec(L = L, trunpt = trunpt, above = above,
     r = nGHK)
     result <- drop(result)
     }
     else {
     result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
     algorithm = algorithm, ...)
     }
     return(result)
     }
     <bytecode: 0xa78670>
     <environment: namespace:mvProbit>
     --- function search by body ---
     Function pmvnormWrap in namespace mvProbit has this body.
     ----------- END OF FAILURE REPORT --------------
     Error in if (class(L) == "try-error") { : the condition has length > 1
     Calls: mvProbitExp -> mvProbitExpInternal -> pmvnormWrap
     Execution halted
    Running the tests in 'tests/pmvnormWrapTest.R' failed.
    Complete output:
     > library( "mvProbit" )
     Loading required package: mvtnorm
     Loading required package: maxLik
     Loading required package: miscTools
    
     Please cite the 'maxLik' package as:
     Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
    
     If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
     https://r-forge.r-project.org/projects/maxlik/
     Loading required package: abind
     > library( "miscTools" )
     >
     > # covariance matrix
     > sigma <- miscTools::symMatrix( c( 1, 0.2, 0.4, 1, -0.1, 1 ) )
     >
     > ######## only upper ##########
     > upper <- c( -0.3, 0.7, -0.5 )
     > # Genz + Bretz (default)
     > pug <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
     + algorithm = GenzBretz(), random.seed = 123 )
     > print( pug )
     [1] 0.1360783
     attr(,"error")
     [1] 0.000118605
     attr(,"msg")
     [1] "Normal Completion"
     >
     > # Miwa (as function)
     > pum <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
     + algorithm = Miwa, random.seed = 123 )
     > print( pum )
     [1] 0.136069
     attr(,"error")
     [1] NA
     attr(,"msg")
     [1] "Normal Completion"
     > all.equal( pug, pum, check.attributes = FALSE, tol = 1e-4 )
     [1] TRUE
     >
     > # Miwa (as object returned from function Miwa())
     > pum1 <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
     + algorithm = Miwa(), random.seed = 123 )
     > all.equal( pum, pum1 )
     [1] TRUE
     >
     > # Miwa (as character string)
     > pum2 <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
     + algorithm = "Miwa", random.seed = 123 )
     > all.equal( pum, pum2 )
     [1] TRUE
     >
     > # TVPACK
     > put <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
     + algorithm = TVPACK, random.seed = 123 )
     > print( put )
     [1] 0.136069
     attr(,"error")
     [1] 1e-06
     attr(,"msg")
     [1] "Normal Completion"
     > all.equal( pug, put, check.attributes = FALSE, tol = 1e-4 )
     [1] TRUE
     > all.equal( pum, put, check.attributes = FALSE, tol = 1e-6 )
     [1] TRUE
     >
     > # GHK
     > pughk <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
     + algorithm = "ghk", random.seed = 123, nGHK = 1000 )
     ----------- FAILURE REPORT --------------
     --- failure: the condition has length > 1 ---
     --- srcref ---
     :
     --- package (from environment) ---
     mvProbit
     --- call from context ---
     mvProbit:::pmvnormWrap(upper = upper, sigma = sigma, algorithm = "ghk",
     random.seed = 123, nGHK = 1000)
     --- call from argument ---
     if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
     }
     --- R stacktrace ---
     where 1: mvProbit:::pmvnormWrap(upper = upper, sigma = sigma, algorithm = "ghk",
     random.seed = 123, nGHK = 1000)
    
     --- value of length: 2 type: logical ---
     [1] FALSE FALSE
     --- function from context ---
     function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
     nGHK = NULL, ...)
     {
     if (!is.matrix(sigma)) {
     stop("argument 'sigma' must be a matrix")
     }
     else if (nrow(sigma) != ncol(sigma)) {
     stop("argument 'sigma' must be a square matrix")
     }
     else if (!all(is.numeric(sigma))) {
     stop("argument 'sigma must be a numeric matrix")
     }
     if (!all(is.numeric(lower))) {
     stop("argument 'lower' must be numeric")
     }
     else if (length(lower) == 1) {
     lower <- rep(lower, nrow(sigma))
     }
     else if (length(lower) != nrow(sigma)) {
     stop("argument 'lower' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     if (!all(is.numeric(upper))) {
     stop("argument 'upper' must be numeric")
     }
     else if (length(upper) == 1) {
     upper <- rep(upper, nrow(sigma))
     }
     else if (length(upper) != nrow(sigma)) {
     stop("argument 'upper' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     algOkay <- TRUE
     ghk <- FALSE
     if (!is.list(algorithm) && length(algorithm) != 1) {
     stop("argument 'algorithm' must be a single function",
     " or a single character string")
     }
     else if (is.function(algorithm)) {
     algResult <- do.call(algorithm, list())
     if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
     "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (is.character(algorithm)) {
     if (tolower(algorithm) == "ghk") {
     ghk <- TRUE
     }
     else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     if (!algOkay) {
     stop("argument 'algorithm' must be either one of the functions",
     " 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
     " \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
     }
     if (length(random.seed) != 1) {
     stop("argument 'random.seed' must be a single numerical value")
     }
     else if (!is.numeric(random.seed)) {
     stop("argument 'random.seed' must be numerical")
     }
     if (exists(".Random.seed")) {
     savedSeed <- .Random.seed
     }
     set.seed(random.seed)
     if (exists("savedSeed")) {
     on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
     }
     else {
     on.exit(rm(.Random.seed, envir = sys.frame()))
     }
     if (ghk) {
     if (is.null(nGHK)) {
     stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
     }
     else if (length(nGHK) != 1) {
     stop("argument 'nGHK' must be a single integer value")
     }
     else if (!is.numeric(nGHK)) {
     stop("argument 'nGHK' must be numeric")
     }
     else if (nGHK <= 0) {
     stop("argument 'nGHK' must be positive")
     }
     L <- try(t(chol(sigma)), silent = TRUE)
     if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
     }
     trunpt <- rep(NA, length(lower))
     above <- rep(NA, length(lower))
     for (i in 1:length(lower)) {
     if (lower[i] == -Inf) {
     trunpt[i] <- upper[i]
     above[i] <- 1
     }
     else if (upper[i] == Inf) {
     trunpt[i] <- lower[i]
     above[i] <- 0
     }
     else {
     stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
     " or the upper truncation point must be 'Inf'")
     }
     }
     sink(tempfile())
     on.exit(sink(), add = TRUE)
     result <- ghkvec(L = L, trunpt = trunpt, above = above,
     r = nGHK)
     result <- drop(result)
     }
     else {
     result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
     algorithm = algorithm, ...)
     }
     return(result)
     }
     <bytecode: 0x46c91e0>
     <environment: namespace:mvProbit>
     --- function search by body ---
     Function pmvnormWrap in namespace mvProbit has this body.
     ----------- END OF FAILURE REPORT --------------
     Error in if (class(L) == "try-error") { : the condition has length > 1
     Calls: <Anonymous>
     Execution halted
Flavor: r-devel-linux-x86_64-debian-clang

Version: 0.1-8
Check: examples
Result: ERROR
    Running examples in ‘mvProbit-Ex.R’ failed
    The error most likely occurred in:
    
    > ### Name: mvProbit
    > ### Title: Estimation of Multivariate Probit Models
    > ### Aliases: mvProbit print.mvProbit
    > ### Keywords: models regression
    >
    > ### ** Examples
    >
    > ## generate a simulated data set
    > set.seed( 123 )
    > # number of observations
    > nObs <- 50
    >
    > # generate explanatory variables
    > xMat <- cbind(
    + const = rep( 1, nObs ),
    + x1 = as.numeric( rnorm( nObs ) > 0 ),
    + x2 = rnorm( nObs ) )
    >
    > # model coefficients
    > beta <- cbind( c( 0.8, 1.2, -0.8 ),
    + c( -0.6, 1.0, -1.6 ),
    + c( 0.5, -0.6, 1.2 ) )
    >
    > # covariance matrix of error terms
    > library( miscTools )
    > sigma <- symMatrix( c( 1, 0.2, 0.4, 1, -0.1, 1 ) )
    >
    > # generate dependent variables
    > yMatLin <- xMat %*% beta
    > yMat <- ( yMatLin + rmvnorm( nObs, sigma = sigma ) ) > 0
    > colnames( yMat ) <- paste( "y", 1:3, sep = "" )
    >
    > # estimation (BHHH optimizer and GHK algorithm)
    > estResult <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
    + data = as.data.frame( cbind( xMat, yMat ) ), iterlim = 1, nGHK = 50 )
    Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
     ----------- FAILURE REPORT --------------
     --- failure: the condition has length > 1 ---
     --- srcref ---
    :
     --- package (from environment) ---
    mvProbit
     --- call from context ---
    pmvnormWrap(upper = xBetaTmp, sigma = sigmaTmp, algorithm = algorithm,
     nGHK = nGHK, random.seed = randomSeed, ...)
     --- call from argument ---
    if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
    }
     --- R stacktrace ---
    where 1: pmvnormWrap(upper = xBetaTmp, sigma = sigmaTmp, algorithm = algorithm,
     nGHK = nGHK, random.seed = randomSeed, ...)
    where 2: mvProbitLogLikInternal(yMat = yMat, xMat = xMat, coef = param,
     sigma = NULL, algorithm = llAlgorithm, nGHK = llNGHK, returnGrad = llIntGrad,
     oneSidedGrad = llOneSidedGrad, eps = llEps, randomSeed = llRandom.seed,
     ...)
    where 3: fnOrig(theta, ...)
    where 4: fn(start1, fixed = fixed, sumObs = TRUE, returnHessian = returnHessian,
     ...)
    where 5: maxNRCompute(fn = function (theta, fnOrig, gradOrig = NULL, hessOrig = NULL,
     fixed, sumObs = FALSE, returnHessian = TRUE, ...)
    {
     nParam <- length(theta)
     f <- fnOrig(theta, ...)
     if (any(is.na(f))) {
     attr(f, "gradient") <- NA
     attr(f, "hessian") <- NA
     return(f)
     }
     gr <- attr(f, "gradient")
     if (is.null(gr)) {
     if (!is.null(gradOrig)) {
     gr <- gradOrig(theta, ...)
     }
     else {
     gr <- numericGradient(f = fnOrig, t0 = theta, fixed = fixed,
     ...)
     }
     }
     if (is.matrix(gr)) {
     if (ncol(gr) != length(theta)) {
     stop(paste0("if gradient is a matrix, it must have length(parameter) colums (currently ",
     length(theta), "), not ", ncol(gr)))
     }
     activeGr <- gr[, !fixed]
     }
     else {
     activeGr <- gr[!fixed]
     }
     if (any(is.na(activeGr))) {
     attr(f, "gradient") <- gr
     attr(f, "hessian") <- NA
     return(f)
     }
     if (observationGradient(gr, length(theta))) {
     gr <- as.matrix(gr)
     }
     if (is.null(dim(gr))) {
     gr[fixed] <- NA
     }
     else {
     gr[, fixed] <- NA
     }
     if (isTRUE(returnHessian)) {
     h <- attr(f, "hessian")
     if (is.null(h)) {
     if (!is.null(hessOrig)) {
     h <- as.matrix(hessOrig(theta, ...))
     }
     else {
     llFunc <- function(theta, ...) {
     return(sum(fnOrig(theta, ...)))
     }
     if (!is.null(attr(f, "gradient"))) {
     gradFunc <- function(theta, ...) {
     return(sumGradients(attr(fnOrig(theta, ...),
     "gradient"), nParam))
     }
     }
     else if (!is.null(gradOrig)) {
     gradFunc <- function(theta, ...) {
     return(sumGradients(gradOrig(theta, ...),
     nParam))
     }
     }
     else {
     gradFunc <- NULL
     }
     h <- numericHessian(f = llFunc, grad = gradFunc,
     t0 = theta, fixed = fixed, ...)
     }
     }
     if ((dim(h)[1] != nParam) | (dim(h)[2] != nParam)) {
     stop("Wrong hessian dimension. Needed ", nParam,
     "x", nParam, " but supplied ", dim(h)[1], "x",
     dim(h)[2])
     }
     else {
     h[fixed, ] <- NA
     h[, fixed] <- NA
     }
     }
     else if (tolower(returnHessian) == "bhhh") {
     h <- NULL
     if (is.null(dim(gr)) & any(is.na(gr[!fixed]))) {
     h <- NA
     }
     else if (is.matrix(gr)) {
     if (any(is.na(gr[, !fixed]))) {
     h <- NA
     }
     }
     if (is.null(h)) {
     checkBhhhGrad(g = gr, theta = theta, analytic = (!is.null(attr(f,
     "gradient")) || !is.null(gradOrig)), fixed = fixed)
     h <- -crossprod(gr)
     }
     attr(h, "type") = "BHHH"
     }
     else {
     h <- NULL
     }
     if (sumObs) {
     f <- sumKeepAttr(f)
     }
     if (sumObs) {
     gr <- sumGradients(gr, nParam)
     }
     if (!is.null(gradOrig) && !is.null(attr(f, "gradient"))) {
     attr(f, "gradBoth") <- TRUE
     }
     if (!is.null(hessOrig) && !is.null(attr(f, "hessian"))) {
     attr(f, "hessBoth") <- TRUE
     }
     attr(f, "gradient") <- gr
     attr(f, "hessian") <- h
     return(f)
    }, fnOrig = function (param, yMat, xMat, llAlgorithm, llNGHK,
     llIntGrad, llOneSidedGrad, llEps, llRandom.seed, ...)
    {
     logLikVal <- mvProbitLogLikInternal(yMat = yMat, xMat = xMat,
     coef = param, sigma = NULL, algorithm = llAlgorithm,
     nGHK = llNGHK, returnGrad = llIntGrad, oneSidedGrad = llOneSidedGrad,
     eps = llEps, randomSeed = llRandom.seed, ...)
     return(logLikVal)
    }, start = c(b_1_0 = 0.598378200054935, b_1_1 = 6.336646267785,
    b_1_2 = -1.01231147815159, b_2_0 = -1.03564793824529, b_2_1 = 1.26302210717712,
    b_2_2 = -1.70638777681512, b_3_0 = 0.410893287564513, b_3_1 = -0.373151392833251,
    b_3_2 = 2.92074017047352, R_1_2 = 0.172617843070694, R_1_3 = -0.0160393443285098,
    R_2_3 = -0.0704265674853818), finalHessian = "BHHH", bhhhHessian = TRUE,
     fixed = c(b_1_0 = FALSE, b_1_1 = FALSE, b_1_2 = FALSE, b_2_0 = FALSE,
     b_2_1 = FALSE, b_2_2 = FALSE, b_3_0 = FALSE, b_3_1 = FALSE,
     b_3_2 = FALSE, R_1_2 = FALSE, R_1_3 = FALSE, R_2_3 = FALSE
     ), control = new("MaxControl", tol = 1e-08, reltol = 1.49011611938477e-08,
     gradtol = 1e-06, steptol = 1e-10, lambdatol = 1e-06,
     qrtol = 1e-10, qac = "stephalving", marquardt_lambda0 = 0.01,
     marquardt_lambdaStep = 2, marquardt_maxLambda = 1e+12,
     nm_alpha = 1, nm_beta = 0.5, nm_gamma = 2, sann_cand = NULL,
     sann_temp = 10, sann_tmax = 10L, sann_randomSeed = 123L,
     iterlim = 1L, printLevel = 0L), yMat = c(0, 1, 1, 1,
     1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1,
     1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1,
     1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1,
     0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1,
     1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1,
     1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0,
     1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1,
     1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0), xMat = c(1, 1, 1,
     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0,
     1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0,
     1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0,
     1, 0, 0.253318513994755, -0.028546755348703, -0.0428704572913161,
     1.36860228401446, -0.225770985659268, 1.51647060442954, -1.54875280423022,
     0.584613749636069, 0.123854243844614, 0.215941568743973,
     0.379639482759882, -0.502323453109302, -0.33320738366942,
     -1.01857538310709, -1.07179122647558, 0.303528641404258,
     0.448209778629426, 0.0530042267305041, 0.922267467879738,
     2.05008468562714, -0.491031166056535, -2.30916887564081,
     1.00573852446226, -0.709200762582393, -0.688008616467358,
     1.0255713696967, -0.284773007051009, -1.22071771225454, 0.18130347974915,
     -0.138891362439045, 0.00576418589988693, 0.38528040112633,
     -0.370660031792409, 0.644376548518833, -0.220486561818751,
     0.331781963915697, 1.09683901314935, 0.435181490833803, -0.325931585531227,
     1.14880761845109, 0.993503855962119, 0.54839695950807, 0.238731735111441,
     -0.627906076039371, 1.36065244853001, -0.600259587147127,
     2.18733299301658, 1.53261062618519, -0.235700359100477, -1.02642090030678
     ), llAlgorithm = "GHK", llNGHK = 50, llIntGrad = TRUE, llOneSidedGrad = FALSE,
     llEps = 1e-06, llRandom.seed = 123)
    where 6: eval(as.call(cl))
    where 7: eval(as.call(cl))
    where 8: maxNR(fn = fn, grad = grad, hess = hess, start = start, finalHessian = finalHessian,
     bhhhHessian = TRUE, ...)
    where 9: maxRoutine(fn = logLik, grad = grad, hess = hess, start = start,
     constraints = constraints, ...)
    where 10: maxLik(logLik = logLik, start = start, method = method, finalHessian = finalHessian,
     yMat = yMat, xMat = xMat, llAlgorithm = algorithm, llNGHK = nGHK,
     llIntGrad = intGrad, llOneSidedGrad = oneSidedGrad, llEps = eps,
     llRandom.seed = random.seed, ...)
    where 11: mvProbit(cbind(y1, y2, y3) ~ x1 + x2, data = as.data.frame(cbind(xMat,
     yMat)), iterlim = 1, nGHK = 50)
    
     --- value of length: 2 type: logical ---
    [1] FALSE FALSE
     --- function from context ---
    function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
     nGHK = NULL, ...)
    {
     if (!is.matrix(sigma)) {
     stop("argument 'sigma' must be a matrix")
     }
     else if (nrow(sigma) != ncol(sigma)) {
     stop("argument 'sigma' must be a square matrix")
     }
     else if (!all(is.numeric(sigma))) {
     stop("argument 'sigma must be a numeric matrix")
     }
     if (!all(is.numeric(lower))) {
     stop("argument 'lower' must be numeric")
     }
     else if (length(lower) == 1) {
     lower <- rep(lower, nrow(sigma))
     }
     else if (length(lower) != nrow(sigma)) {
     stop("argument 'lower' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     if (!all(is.numeric(upper))) {
     stop("argument 'upper' must be numeric")
     }
     else if (length(upper) == 1) {
     upper <- rep(upper, nrow(sigma))
     }
     else if (length(upper) != nrow(sigma)) {
     stop("argument 'upper' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     algOkay <- TRUE
     ghk <- FALSE
     if (!is.list(algorithm) && length(algorithm) != 1) {
     stop("argument 'algorithm' must be a single function",
     " or a single character string")
     }
     else if (is.function(algorithm)) {
     algResult <- do.call(algorithm, list())
     if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
     "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (is.character(algorithm)) {
     if (tolower(algorithm) == "ghk") {
     ghk <- TRUE
     }
     else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     if (!algOkay) {
     stop("argument 'algorithm' must be either one of the functions",
     " 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
     " \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
     }
     if (length(random.seed) != 1) {
     stop("argument 'random.seed' must be a single numerical value")
     }
     else if (!is.numeric(random.seed)) {
     stop("argument 'random.seed' must be numerical")
     }
     if (exists(".Random.seed")) {
     savedSeed <- .Random.seed
     }
     set.seed(random.seed)
     if (exists("savedSeed")) {
     on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
     }
     else {
     on.exit(rm(.Random.seed, envir = sys.frame()))
     }
     if (ghk) {
     if (is.null(nGHK)) {
     stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
     }
     else if (length(nGHK) != 1) {
     stop("argument 'nGHK' must be a single integer value")
     }
     else if (!is.numeric(nGHK)) {
     stop("argument 'nGHK' must be numeric")
     }
     else if (nGHK <= 0) {
     stop("argument 'nGHK' must be positive")
     }
     L <- try(t(chol(sigma)), silent = TRUE)
     if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
     }
     trunpt <- rep(NA, length(lower))
     above <- rep(NA, length(lower))
     for (i in 1:length(lower)) {
     if (lower[i] == -Inf) {
     trunpt[i] <- upper[i]
     above[i] <- 1
     }
     else if (upper[i] == Inf) {
     trunpt[i] <- lower[i]
     above[i] <- 0
     }
     else {
     stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
     " or the upper truncation point must be 'Inf'")
     }
     }
     sink(tempfile())
     on.exit(sink(), add = TRUE)
     result <- ghkvec(L = L, trunpt = trunpt, above = above,
     r = nGHK)
     result <- drop(result)
     }
     else {
     result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
     algorithm = algorithm, ...)
     }
     return(result)
    }
    <bytecode: 0x3655878>
    <environment: namespace:mvProbit>
     --- function search by body ---
    Function pmvnormWrap in namespace mvProbit has this body.
     ----------- END OF FAILURE REPORT --------------
    Error in if (class(L) == "try-error") { : the condition has length > 1
    Calls: mvProbit ... fn -> fnOrig -> mvProbitLogLikInternal -> pmvnormWrap
    Execution halted
Flavor: r-devel-linux-x86_64-fedora-clang

Version: 0.1-8
Check: tests
Result: ERROR
     Running ‘mvProbitEst.R’ [34s/36s]
     Running ‘mvProbitMargEff2.R’
     Running ‘mvProbitTest.R’
     Running ‘mvProbitTest5.R’ [145s/148s]
     Comparing ‘mvProbitTest5.Rout’ to ‘mvProbitTest5.Rout.save’ ...1604c1604
    < 5 0.00 -0.01 0 0.00 0.00 0.00 0.00
    ---
    > 5 0.00 0.00 0 0.00 0.00 0.00 0.00
    1618,1620c1618,1620
    < 1 0 0.00 0 0 0 0 0
    < 5 0 -0.01 0 0 0 0 0
    < 10 0 0.00 0 0 0 0 0
    ---
    > 1 0 0 0 0 0 0 0
    > 5 0 0 0 0 0 0 0
    > 10 0 0 0 0 0 0 0
    1728c1728
    < 5 -0.01
    ---
    > 5 0.00
     Running ‘pmvnormWrapTest.R’
    Running the tests in ‘tests/mvProbitEst.R’ failed.
    Complete output:
     > library( "mvProbit" )
     Loading required package: mvtnorm
     Loading required package: maxLik
     Loading required package: miscTools
    
     Please cite the 'maxLik' package as:
     Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
    
     If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
     https://r-forge.r-project.org/projects/maxlik/
     Loading required package: abind
     > options( digits = 4 )
     >
     > ## generate a simulated data set
     > set.seed( 123 )
     > # number of observations
     > nObs <- 50
     >
     > # generate explanatory variables
     > xMat <- cbind(
     + const = rep( 1, nObs ),
     + x1 = as.numeric( rnorm( nObs ) > 0 ),
     + x2 = rnorm( nObs ) )
     >
     > # model coefficients
     > beta <- cbind( c( 0.8, 1.2, -0.8 ),
     + c( -0.6, 1.0, -1.6 ),
     + c( 0.5, -0.6, 1.2 ) )
     >
     > # covariance matrix of error terms
     > sigma <- miscTools::symMatrix( c( 1, 0.2, 0.4, 1, -0.1, 1 ) )
     >
     > # generate dependent variables
     > yMatLin <- xMat %*% beta
     > yMat <- ( yMatLin + rmvnorm( nObs, sigma = sigma, pre0.9_9994 = TRUE ) ) > 0
     > colnames( yMat ) <- paste( "y", 1:3, sep = "" )
     >
     > # create data frame
     > dat <- as.data.frame( cbind( xMat, yMat ) )
     >
     > # estimation with the BHHH algorithm, two-sided gradients
     > estResultBHHH <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + start = c( beta ), startSigma = sigma,
     + data = dat, tol = 0.5,
     + algorithm = GenzBretz() )
     > print( estResultBHHH )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, algorithm = GenzBretz(), tol = 0.5)
    
     Coefficients:
     b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
     1.02508 1.30536 -1.02809 -0.78248 1.14318 -1.66938 0.58229 -0.69667
     b_3_2 R_1_2 R_1_3 R_2_3
     1.23283 0.86322 0.36022 -0.08406
    
     > summary( estResultBHHH )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, algorithm = GenzBretz(), tol = 0.5)
    
     Coefficients:
     Estimate Std. error t value Pr(> t)
     b_1_0 1.02508 0.44548 2.301 0.0214 *
     b_1_1 1.30536 0.73511 1.776 0.0758 .
     b_1_2 -1.02809 0.68851 -1.493 0.1354
     b_2_0 -0.78248 0.38798 -2.017 0.0437 *
     b_2_1 1.14318 0.53026 2.156 0.0311 *
     b_2_2 -1.66938 0.81308 -2.053 0.0401 *
     b_3_0 0.58229 0.64107 0.908 0.3637
     b_3_1 -0.69667 0.89818 -0.776 0.4380
     b_3_2 1.23283 0.62452 1.974 0.0484 *
     R_1_2 0.86322 18.60071 0.046 0.9630
     R_1_3 0.36022 0.64712 0.557 0.5778
     R_2_3 -0.08406 0.34909 -0.241 0.8097
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     BHHH maximisation, 1 iterations
     Return code 2: successive function values within tolerance limit
     Log-likelihood: -56.91 on 12 Df
    
     > logLik( estResultBHHH )
     'log Lik.' -56.91 (df=12)
     > estResultBHHHA <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + start = c( beta, sigma[ lower.tri( sigma ) ] ),
     + data = dat, tol = 0.5,
     + algorithm = GenzBretz() )
     > all.equal( estResultBHHH, estResultBHHHA )
     [1] "Component \"call\": target, current do not match when deparsed"
     >
     > # estimation with the BHHH algorithm, one-sided gradients
     > estResultBHHH1 <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + start = c( beta ), startSigma = sigma,
     + data = dat, tol = 0.5,
     + algorithm = GenzBretz(), oneSidedGrad = TRUE )
     > print( estResultBHHH1 )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, algorithm = GenzBretz(), oneSidedGrad = TRUE,
     tol = 0.5)
    
     Coefficients:
     b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
     1.02508 1.30536 -1.02809 -0.78248 1.14318 -1.66938 0.58229 -0.69667
     b_3_2 R_1_2 R_1_3 R_2_3
     1.23283 0.86322 0.36022 -0.08406
    
     > summary( estResultBHHH1 )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, algorithm = GenzBretz(), oneSidedGrad = TRUE,
     tol = 0.5)
    
     Coefficients:
     Estimate Std. error t value Pr(> t)
     b_1_0 1.02508 0.44548 2.301 0.0214 *
     b_1_1 1.30536 0.73511 1.776 0.0758 .
     b_1_2 -1.02809 0.68851 -1.493 0.1354
     b_2_0 -0.78248 0.38797 -2.017 0.0437 *
     b_2_1 1.14318 0.53026 2.156 0.0311 *
     b_2_2 -1.66938 0.81307 -2.053 0.0401 *
     b_3_0 0.58229 0.64106 0.908 0.3637
     b_3_1 -0.69667 0.89817 -0.776 0.4380
     b_3_2 1.23283 0.62452 1.974 0.0484 *
     R_1_2 0.86322 18.60041 0.046 0.9630
     R_1_3 0.36022 0.64712 0.557 0.5778
     R_2_3 -0.08406 0.34909 -0.241 0.8097
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     BHHH maximisation, 1 iterations
     Return code 2: successive function values within tolerance limit
     Log-likelihood: -56.91 on 12 Df
    
     > logLik( estResultBHHH1 )
     'log Lik.' -56.91 (df=12)
     > all.equal( estResultBHHH, estResultBHHH1, tol = 1e-5 )
     [1] "Component \"call\": target, current do not match when deparsed"
     >
     > # estimation with the BFGS algorithm, two-sided gradients
     > estResultBFGS <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + start = c( beta ), startSigma = sigma, method = "BFGS",
     + data = dat,
     + reltol = 0.5, algorithm = GenzBretz() )
     > print( estResultBFGS )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, method = "BFGS", algorithm = GenzBretz(),
     reltol = 0.5)
    
     Coefficients:
     b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
     0.53607 1.16101 -1.14212 -0.44952 1.05496 -1.57624 0.40751 -0.85400
     b_3_2 R_1_2 R_1_3 R_2_3
     1.25958 0.29483 0.32822 -0.02382
    
     > summary( estResultBFGS )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, method = "BFGS", algorithm = GenzBretz(),
     reltol = 0.5)
    
     Coefficients:
     Estimate Std. error t value Pr(> t)
     b_1_0 0.53607 0.82709 0.648 0.5169
     b_1_1 1.16101 0.91515 1.269 0.2046
     b_1_2 -1.14212 0.93847 -1.217 0.2236
     b_2_0 -0.44952 0.79215 -0.567 0.5704
     b_2_1 1.05496 0.88398 1.193 0.2327
     b_2_2 -1.57624 0.73187 -2.154 0.0313 *
     b_3_0 0.40751 0.48897 0.833 0.4046
     b_3_1 -0.85400 0.84871 -1.006 0.3143
     b_3_2 1.25958 0.68040 1.851 0.0641 .
     R_1_2 0.29483 1.85817 0.159 0.8739
     R_1_3 0.32822 0.70923 0.463 0.6435
     R_2_3 -0.02382 0.34466 -0.069 0.9449
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     BFGS maximization, 4 iterations
     Return code 0: successful convergence
     Log-likelihood: -52.94 on 12 Df
    
     > logLik( estResultBFGS )
     'log Lik.' -52.94 (df=12)
     >
     > # estimation with the BFGS algorithm, one-sided gradients
     > estResultBFGS1 <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + start = c( beta ), startSigma = sigma, method = "BFGS",
     + data = dat,
     + reltol = 0.5, algorithm = GenzBretz(), oneSidedGrad = TRUE )
     > print( estResultBFGS1 )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, method = "BFGS", algorithm = GenzBretz(),
     oneSidedGrad = TRUE, reltol = 0.5)
    
     Coefficients:
     b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
     0.53607 1.16101 -1.14213 -0.44952 1.05496 -1.57624 0.40751 -0.85400
     b_3_2 R_1_2 R_1_3 R_2_3
     1.25958 0.29483 0.32822 -0.02382
    
     > summary( estResultBFGS1 )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, method = "BFGS", algorithm = GenzBretz(),
     oneSidedGrad = TRUE, reltol = 0.5)
    
     Coefficients:
     Estimate Std. error t value Pr(> t)
     b_1_0 0.53607 0.82709 0.648 0.5169
     b_1_1 1.16101 0.91515 1.269 0.2046
     b_1_2 -1.14213 0.93847 -1.217 0.2236
     b_2_0 -0.44952 0.79215 -0.567 0.5704
     b_2_1 1.05496 0.88398 1.193 0.2327
     b_2_2 -1.57624 0.73187 -2.154 0.0313 *
     b_3_0 0.40751 0.48897 0.833 0.4046
     b_3_1 -0.85400 0.84871 -1.006 0.3143
     b_3_2 1.25958 0.68040 1.851 0.0641 .
     R_1_2 0.29483 1.85817 0.159 0.8739
     R_1_3 0.32822 0.70923 0.463 0.6435
     R_2_3 -0.02382 0.34466 -0.069 0.9449
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     BFGS maximization, 4 iterations
     Return code 0: successful convergence
     Log-likelihood: -52.94 on 12 Df
    
     > logLik( estResultBFGS1 )
     'log Lik.' -52.94 (df=12)
     > all.equal( estResultBFGS, estResultBFGS1, tol = 1e-5 )
     [1] "Component \"call\": target, current do not match when deparsed"
     >
     > # estimation with the BFGS algorithm, one-sided gradients, no starting values
     > estResultBFGS1a <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + data = dat, method = "BFGS",
     + reltol = 0.5, algorithm = GenzBretz(), oneSidedGrad = TRUE )
     > print( estResultBFGS1a )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, method = "BFGS",
     algorithm = GenzBretz(), oneSidedGrad = TRUE, reltol = 0.5)
    
     Coefficients:
     b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
     0.79353 1.88367 -1.74922 -0.40036 0.91577 -1.52106 0.97790 -1.69433
     b_3_2 R_1_2 R_1_3 R_2_3
     1.57625 0.75270 0.23721 -0.03453
    
     > summary( estResultBFGS1a )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, method = "BFGS",
     algorithm = GenzBretz(), oneSidedGrad = TRUE, reltol = 0.5)
    
     Coefficients:
     Estimate Std. error t value Pr(> t)
     b_1_0 0.79353 0.90412 0.878 0.3801
     b_1_1 1.88367 0.90668 2.078 0.0378 *
     b_1_2 -1.74922 0.97400 -1.796 0.0725 .
     b_2_0 -0.40036 0.69056 -0.580 0.5621
     b_2_1 0.91577 0.97998 0.934 0.3501
     b_2_2 -1.52106 0.72006 -2.112 0.0347 *
     b_3_0 0.97790 0.64497 1.516 0.1295
     b_3_1 -1.69433 0.94788 -1.787 0.0739 .
     b_3_2 1.57625 0.73399 2.148 0.0318 *
     R_1_2 0.75270 4.32729 0.174 0.8619
     R_1_3 0.23721 0.90535 0.262 0.7933
     R_2_3 -0.03453 0.35889 -0.096 0.9234
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     BFGS maximization, 3 iterations
     Return code 0: successful convergence
     Log-likelihood: -49.36 on 12 Df
    
     > logLik( estResultBFGS1a )
     'log Lik.' -49.36 (df=12)
     >
     > # estimation with the BFGS algorithm, one-sided gradients, no starting values for beta
     > estResultBFGS1b <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + startSigma = sigma, data = dat,
     + method = "BFGS", reltol = 0.5, algorithm = GenzBretz(), oneSidedGrad = TRUE )
     > print( estResultBFGS1b )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, startSigma = sigma,
     method = "BFGS", algorithm = GenzBretz(), oneSidedGrad = TRUE,
     reltol = 0.5)
    
     Coefficients:
     b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1 b_3_2 R_1_2
     0.8110 1.9060 -1.7627 -0.3770 0.9441 -1.5332 0.9372 -1.7498 1.6024 0.6817
     R_1_3 R_2_3
     0.4044 0.1663
    
     > summary( estResultBFGS1b )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, startSigma = sigma,
     method = "BFGS", algorithm = GenzBretz(), oneSidedGrad = TRUE,
     reltol = 0.5)
    
     Coefficients:
     Estimate Std. error t value Pr(> t)
     b_1_0 0.8110 0.9297 0.872 0.3830
     b_1_1 1.9060 0.9160 2.081 0.0375 *
     b_1_2 -1.7627 1.0131 -1.740 0.0819 .
     b_2_0 -0.3770 0.7049 -0.535 0.5928
     b_2_1 0.9441 1.0345 0.913 0.3615
     b_2_2 -1.5332 0.7517 -2.040 0.0414 *
     b_3_0 0.9372 0.6291 1.490 0.1363
     b_3_1 -1.7498 0.9996 -1.751 0.0800 .
     b_3_2 1.6024 0.7848 2.042 0.0412 *
     R_1_2 0.6817 3.4094 0.200 0.8415
     R_1_3 0.4044 0.8024 0.504 0.6143
     R_2_3 0.1663 0.3389 0.491 0.6235
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     BFGS maximization, 3 iterations
     Return code 0: successful convergence
     Log-likelihood: -49.57 on 12 Df
    
     > logLik( estResultBFGS1b )
     'log Lik.' -49.57 (df=12)
     >
     > # estimation with the BFGS algorithm, one-sided gradients, no starting values for sigma
     > estResultBFGS1s <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + start = c( beta ), data = dat,
     + method = "BFGS", reltol = 0.5, algorithm = GenzBretz(), oneSidedGrad = TRUE )
     > print( estResultBFGS1s )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     method = "BFGS", algorithm = GenzBretz(), oneSidedGrad = TRUE,
     reltol = 0.5)
    
     Coefficients:
     b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
     0.55814 1.14958 -1.12747 -0.44426 1.07325 -1.57226 0.37255 -0.85470
     b_3_2 R_1_2 R_1_3 R_2_3
     1.23011 0.34981 -0.01177 0.09435
    
     > summary( estResultBFGS1s )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     method = "BFGS", algorithm = GenzBretz(), oneSidedGrad = TRUE,
     reltol = 0.5)
    
     Coefficients:
     Estimate Std. error t value Pr(> t)
     b_1_0 0.55814 0.86615 0.644 0.5193
     b_1_1 1.14958 0.90607 1.269 0.2045
     b_1_2 -1.12747 0.94191 -1.197 0.2313
     b_2_0 -0.44426 0.85526 -0.519 0.6035
     b_2_1 1.07325 0.99314 1.081 0.2798
     b_2_2 -1.57226 0.78557 -2.001 0.0453 *
     b_3_0 0.37255 0.56172 0.663 0.5072
     b_3_1 -0.85470 0.95214 -0.898 0.3694
     b_3_2 1.23011 0.71506 1.720 0.0854 .
     R_1_2 0.34981 2.02289 0.173 0.8627
     R_1_3 -0.01177 0.70776 -0.017 0.9867
     R_2_3 0.09435 0.33920 0.278 0.7809
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     BFGS maximization, 4 iterations
     Return code 0: successful convergence
     Log-likelihood: -52.78 on 12 Df
    
     > logLik( estResultBFGS1s )
     'log Lik.' -52.78 (df=12)
     >
     > # estimation with the BFGS algorithm, Miwa algorithm for obtaining integrals
     > estResultBFGSm <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + start = c( beta ), startSigma = sigma,
     + data = dat, method = "BFGS",
     + reltol = 0.5, algorithm = Miwa( steps = 64 ) )
     There were 13 warnings (use warnings() to see them)
     > print( estResultBFGSm )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, method = "BFGS", algorithm = Miwa(steps = 64),
     reltol = 0.5)
    
     Coefficients:
     b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
     0.53609 1.16102 -1.14213 -0.44952 1.05496 -1.57624 0.40753 -0.85400
     b_3_2 R_1_2 R_1_3 R_2_3
     1.25960 0.29491 0.32819 -0.02386
    
     > summary( estResultBFGSm )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, method = "BFGS", algorithm = Miwa(steps = 64),
     reltol = 0.5)
    
     Coefficients:
     Estimate Std. error t value Pr(> t)
     b_1_0 0.53609 0.82702 0.648 0.5168
     b_1_1 1.16102 0.91511 1.269 0.2045
     b_1_2 -1.14213 0.93837 -1.217 0.2235
     b_2_0 -0.44952 0.79215 -0.567 0.5704
     b_2_1 1.05496 0.88386 1.194 0.2326
     b_2_2 -1.57624 0.73173 -2.154 0.0312 *
     b_3_0 0.40753 0.48894 0.834 0.4046
     b_3_1 -0.85400 0.84860 -1.006 0.3142
     b_3_2 1.25960 0.68024 1.852 0.0641 .
     R_1_2 0.29491 1.85789 0.159 0.8739
     R_1_3 0.32819 0.70903 0.463 0.6435
     R_2_3 -0.02386 0.34468 -0.069 0.9448
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     BFGS maximization, 4 iterations
     Return code 0: successful convergence
     Log-likelihood: -52.94 on 12 Df
    
     > logLik( estResultBFGSm )
     'log Lik.' -52.94 (df=12)
     > all.equal( estResultBFGS, estResultBFGSm, tol = 1e-3 )
     [1] "Component \"call\": target, current do not match when deparsed"
     >
     > # estimation with the BFGS algorithm, GHK algorithm for obtaining integrals
     > estResultBFGSg <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + start = c( beta ), startSigma = sigma,
     + data = dat, method = "BFGS",
     + reltol = 0.5 )
     ----------- FAILURE REPORT --------------
     --- failure: the condition has length > 1 ---
     --- srcref ---
     :
     --- package (from environment) ---
     mvProbit
     --- call from context ---
     pmvnormWrap(upper = xBetaTmp, sigma = sigmaTmp, algorithm = algorithm,
     nGHK = nGHK, random.seed = randomSeed, ...)
     --- call from argument ---
     if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
     }
     --- R stacktrace ---
     where 1: pmvnormWrap(upper = xBetaTmp, sigma = sigmaTmp, algorithm = algorithm,
     nGHK = nGHK, random.seed = randomSeed, ...)
     where 2: mvProbitLogLikInternal(yMat = yMat, xMat = xMat, coef = param,
     sigma = NULL, algorithm = llAlgorithm, nGHK = llNGHK, returnGrad = llIntGrad,
     oneSidedGrad = llOneSidedGrad, eps = llEps, randomSeed = llRandom.seed,
     ...)
     where 3: fnOrig(theta, ...)
     where 4: logLikFunc(theta, fnOrig = function (param, yMat, xMat, llAlgorithm,
     llNGHK, llIntGrad, llOneSidedGrad, llEps, llRandom.seed,
     ...)
     {
     logLikVal <- mvProbitLogLikInternal(yMat = yMat, xMat = xMat,
     coef = param, sigma = NULL, algorithm = llAlgorithm,
     nGHK = llNGHK, returnGrad = llIntGrad, oneSidedGrad = llOneSidedGrad,
     eps = llEps, randomSeed = llRandom.seed, ...)
     return(logLikVal)
     }, gradOrig = ..2, hessOrig = ..3, yMat = c(1, 1, 1, 1, 1, 1,
     1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1,
     1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0,
     1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0,
     0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0,
     0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1,
     0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1,
     0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0), xMat = c(1,
     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1,
     0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1,
     1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0.253318513994755,
     -0.028546755348703, -0.0428704572913161, 1.36860228401446, -0.225770985659268,
     1.51647060442954, -1.54875280423022, 0.584613749636069, 0.123854243844614,
     0.215941568743973, 0.379639482759882, -0.502323453109302, -0.33320738366942,
     -1.01857538310709, -1.07179122647558, 0.303528641404258, 0.448209778629426,
     0.0530042267305041, 0.922267467879738, 2.05008468562714, -0.491031166056535,
     -2.30916887564081, 1.00573852446226, -0.709200762582393, -0.688008616467358,
     1.0255713696967, -0.284773007051009, -1.22071771225454, 0.18130347974915,
     -0.138891362439045, 0.00576418589988693, 0.38528040112633, -0.370660031792409,
     0.644376548518833, -0.220486561818751, 0.331781963915697, 1.09683901314935,
     0.435181490833803, -0.325931585531227, 1.14880761845109, 0.993503855962119,
     0.54839695950807, 0.238731735111441, -0.627906076039371, 1.36065244853001,
     -0.600259587147127, 2.18733299301658, 1.53261062618519, -0.235700359100477,
     -1.02642090030678), llAlgorithm = "GHK", llNGHK = 1000, llIntGrad = TRUE,
     llOneSidedGrad = FALSE, llEps = 1e-06, llRandom.seed = 123)
     where 5: eval(f, sys.frame(sys.parent()))
     where 6: eval(f, sys.frame(sys.parent()))
     where 7: callWithoutArgs(theta, fName = fName, args = names(formals(sumt)),
     ...)
     where 8: (function (theta, fName, ...)
     {
     return(callWithoutArgs(theta, fName = fName, args = names(formals(sumt)),
     ...))
     })(theta = c(b_1_0 = 0.8, b_1_1 = 1.2, b_1_2 = -0.8, b_2_0 = -0.6,
     b_2_1 = 1, b_2_2 = -1.6, b_3_0 = 0.5, b_3_1 = -0.6, b_3_2 = 1.2,
     R_1_2 = 0.2, R_1_3 = 0.4, R_2_3 = -0.1), fName = "logLikFunc",
     fnOrig = function (param, yMat, xMat, llAlgorithm, llNGHK,
     llIntGrad, llOneSidedGrad, llEps, llRandom.seed, ...)
     {
     logLikVal <- mvProbitLogLikInternal(yMat = yMat, xMat = xMat,
     coef = param, sigma = NULL, algorithm = llAlgorithm,
     nGHK = llNGHK, returnGrad = llIntGrad, oneSidedGrad = llOneSidedGrad,
     eps = llEps, randomSeed = llRandom.seed, ...)
     return(logLikVal)
     }, gradOrig = NULL, hessOrig = NULL, yMat = c(1, 1, 1, 1,
     1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0,
     1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0,
     0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1,
     1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1,
     0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1,
     0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0,
     0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1,
     1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0), xMat = c(1, 1, 1,
     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0,
     1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0,
     1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0,
     1, 0, 0.253318513994755, -0.028546755348703, -0.0428704572913161,
     1.36860228401446, -0.225770985659268, 1.51647060442954, -1.54875280423022,
     0.584613749636069, 0.123854243844614, 0.215941568743973,
     0.379639482759882, -0.502323453109302, -0.33320738366942,
     -1.01857538310709, -1.07179122647558, 0.303528641404258,
     0.448209778629426, 0.0530042267305041, 0.922267467879738,
     2.05008468562714, -0.491031166056535, -2.30916887564081,
     1.00573852446226, -0.709200762582393, -0.688008616467358,
     1.0255713696967, -0.284773007051009, -1.22071771225454, 0.18130347974915,
     -0.138891362439045, 0.00576418589988693, 0.38528040112633,
     -0.370660031792409, 0.644376548518833, -0.220486561818751,
     0.331781963915697, 1.09683901314935, 0.435181490833803, -0.325931585531227,
     1.14880761845109, 0.993503855962119, 0.54839695950807, 0.238731735111441,
     -0.627906076039371, 1.36065244853001, -0.600259587147127,
     2.18733299301658, 1.53261062618519, -0.235700359100477, -1.02642090030678
     ), llAlgorithm = "GHK", llNGHK = 1000, llIntGrad = TRUE,
     llOneSidedGrad = FALSE, llEps = 1e-06, llRandom.seed = 123)
     where 9: do.call(callWithoutSumt, argList)
     where 10: maxOptim(fn = fn, grad = grad, hess = hess, start = start, method = "BFGS",
     fixed = fixed, constraints = constraints, finalHessian = finalHessian,
     parscale = parscale, control = mControl, ...)
     where 11: maxRoutine(fn = logLik, grad = grad, hess = hess, start = start,
     constraints = constraints, ...)
     where 12: maxLik(logLik = logLik, start = start, method = method, finalHessian = finalHessian,
     yMat = yMat, xMat = xMat, llAlgorithm = algorithm, llNGHK = nGHK,
     llIntGrad = intGrad, llOneSidedGrad = oneSidedGrad, llEps = eps,
     llRandom.seed = random.seed, ...)
     where 13: mvProbit(cbind(y1, y2, y3) ~ x1 + x2, start = c(beta), startSigma = sigma,
     data = dat, method = "BFGS", reltol = 0.5)
    
     --- value of length: 2 type: logical ---
     [1] FALSE FALSE
     --- function from context ---
     function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
     nGHK = NULL, ...)
     {
     if (!is.matrix(sigma)) {
     stop("argument 'sigma' must be a matrix")
     }
     else if (nrow(sigma) != ncol(sigma)) {
     stop("argument 'sigma' must be a square matrix")
     }
     else if (!all(is.numeric(sigma))) {
     stop("argument 'sigma must be a numeric matrix")
     }
     if (!all(is.numeric(lower))) {
     stop("argument 'lower' must be numeric")
     }
     else if (length(lower) == 1) {
     lower <- rep(lower, nrow(sigma))
     }
     else if (length(lower) != nrow(sigma)) {
     stop("argument 'lower' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     if (!all(is.numeric(upper))) {
     stop("argument 'upper' must be numeric")
     }
     else if (length(upper) == 1) {
     upper <- rep(upper, nrow(sigma))
     }
     else if (length(upper) != nrow(sigma)) {
     stop("argument 'upper' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     algOkay <- TRUE
     ghk <- FALSE
     if (!is.list(algorithm) && length(algorithm) != 1) {
     stop("argument 'algorithm' must be a single function",
     " or a single character string")
     }
     else if (is.function(algorithm)) {
     algResult <- do.call(algorithm, list())
     if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
     "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (is.character(algorithm)) {
     if (tolower(algorithm) == "ghk") {
     ghk <- TRUE
     }
     else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     if (!algOkay) {
     stop("argument 'algorithm' must be either one of the functions",
     " 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
     " \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
     }
     if (length(random.seed) != 1) {
     stop("argument 'random.seed' must be a single numerical value")
     }
     else if (!is.numeric(random.seed)) {
     stop("argument 'random.seed' must be numerical")
     }
     if (exists(".Random.seed")) {
     savedSeed <- .Random.seed
     }
     set.seed(random.seed)
     if (exists("savedSeed")) {
     on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
     }
     else {
     on.exit(rm(.Random.seed, envir = sys.frame()))
     }
     if (ghk) {
     if (is.null(nGHK)) {
     stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
     }
     else if (length(nGHK) != 1) {
     stop("argument 'nGHK' must be a single integer value")
     }
     else if (!is.numeric(nGHK)) {
     stop("argument 'nGHK' must be numeric")
     }
     else if (nGHK <= 0) {
     stop("argument 'nGHK' must be positive")
     }
     L <- try(t(chol(sigma)), silent = TRUE)
     if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
     }
     trunpt <- rep(NA, length(lower))
     above <- rep(NA, length(lower))
     for (i in 1:length(lower)) {
     if (lower[i] == -Inf) {
     trunpt[i] <- upper[i]
     above[i] <- 1
     }
     else if (upper[i] == Inf) {
     trunpt[i] <- lower[i]
     above[i] <- 0
     }
     else {
     stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
     " or the upper truncation point must be 'Inf'")
     }
     }
     sink(tempfile())
     on.exit(sink(), add = TRUE)
     result <- ghkvec(L = L, trunpt = trunpt, above = above,
     r = nGHK)
     result <- drop(result)
     }
     else {
     result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
     algorithm = algorithm, ...)
     }
     return(result)
     }
     <bytecode: 0x1074890>
     <environment: namespace:mvProbit>
     --- function search by body ---
     Function pmvnormWrap in namespace mvProbit has this body.
     ----------- END OF FAILURE REPORT --------------
     Error in if (class(L) == "try-error") { : the condition has length > 1
     Calls: mvProbit ... logLikFunc -> fnOrig -> mvProbitLogLikInternal -> pmvnormWrap
     Execution halted
    Running the tests in ‘tests/mvProbitMargEff2.R’ failed.
    Complete output:
     > # I thank Mohit Batham for providing this R script that demonstrated
     > # a bug in mvProbitMargEff() when called with 2 dependent variables
     > library( "mvProbit" )
     Loading required package: mvtnorm
     Loading required package: maxLik
     Loading required package: miscTools
    
     Please cite the 'maxLik' package as:
     Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
    
     If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
     https://r-forge.r-project.org/projects/maxlik/
     Loading required package: abind
     > nObs <- 100
     > set.seed( 123 )
     > xData <- data.frame(
     + const = rep( 1, nObs ),
     + x1 = as.numeric( rnorm( nObs ) > 0 ),
     + x2 = as.numeric( rnorm( nObs ) > 0 ),
     + x3 = rnorm( nObs ),
     + x4 = rnorm( nObs ))
     > beta <- c( 0.8, 1.2, -1.0, 1.4, -0.8,
     + -0.6, 1.0, 0.6, -1.2, -1.6)
     > sigma <- symMatrix( c( 1, 0.2, 1))
     > margEffCond <- try( mvProbitMargEff( ~ x1 + x2 + x3 + x4 , coef = beta,
     + sigma = sigma, data = xData, cond = TRUE ) )
     ----------- FAILURE REPORT --------------
     --- failure: the condition has length > 1 ---
     --- srcref ---
     :
     --- package (from environment) ---
     mvProbit
     --- call from context ---
     pmvnormWrap(upper = xBeta[i, ], sigma = coef$sigma, algorithm = algorithm,
     nGHK = nGHK, random.seed = random.seed, ...)
     --- call from argument ---
     if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
     }
     --- R stacktrace ---
     where 1: pmvnormWrap(upper = xBeta[i, ], sigma = coef$sigma, algorithm = algorithm,
     nGHK = nGHK, random.seed = random.seed, ...)
     where 2: mvProbitExpInternal(yMat = yMat, xMat = xMat, coef = betaLower,
     sigma = coef$sigma, cond = cond, algorithm = algorithm, nGHK = nGHK,
     random.seed = random.seed, ...)
     where 3: mvProbitMargEffInternal(yMat = yMat, xMat = xMat, coef = coef,
     sigma = sigma, cond = cond, algorithm = algorithm, nGHK = nGHK,
     eps = eps, dummyVars = dummyVars, random.seed = random.seed,
     ...)
     where 4: mvProbitMargEff(~x1 + x2 + x3 + x4, coef = beta, sigma = sigma,
     data = xData, cond = TRUE)
     where 5: doTryCatch(return(expr), name, parentenv, handler)
     where 6: tryCatchOne(expr, names, parentenv, handlers[[1L]])
     where 7: tryCatchList(expr, classes, parentenv, handlers)
     where 8: 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 9: try(mvProbitMargEff(~x1 + x2 + x3 + x4, coef = beta, sigma = sigma,
     data = xData, cond = TRUE))
    
     --- value of length: 2 type: logical ---
     [1] FALSE FALSE
     --- function from context ---
     function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
     nGHK = NULL, ...)
     {
     if (!is.matrix(sigma)) {
     stop("argument 'sigma' must be a matrix")
     }
     else if (nrow(sigma) != ncol(sigma)) {
     stop("argument 'sigma' must be a square matrix")
     }
     else if (!all(is.numeric(sigma))) {
     stop("argument 'sigma must be a numeric matrix")
     }
     if (!all(is.numeric(lower))) {
     stop("argument 'lower' must be numeric")
     }
     else if (length(lower) == 1) {
     lower <- rep(lower, nrow(sigma))
     }
     else if (length(lower) != nrow(sigma)) {
     stop("argument 'lower' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     if (!all(is.numeric(upper))) {
     stop("argument 'upper' must be numeric")
     }
     else if (length(upper) == 1) {
     upper <- rep(upper, nrow(sigma))
     }
     else if (length(upper) != nrow(sigma)) {
     stop("argument 'upper' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     algOkay <- TRUE
     ghk <- FALSE
     if (!is.list(algorithm) && length(algorithm) != 1) {
     stop("argument 'algorithm' must be a single function",
     " or a single character string")
     }
     else if (is.function(algorithm)) {
     algResult <- do.call(algorithm, list())
     if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
     "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (is.character(algorithm)) {
     if (tolower(algorithm) == "ghk") {
     ghk <- TRUE
     }
     else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     if (!algOkay) {
     stop("argument 'algorithm' must be either one of the functions",
     " 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
     " \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
     }
     if (length(random.seed) != 1) {
     stop("argument 'random.seed' must be a single numerical value")
     }
     else if (!is.numeric(random.seed)) {
     stop("argument 'random.seed' must be numerical")
     }
     if (exists(".Random.seed")) {
     savedSeed <- .Random.seed
     }
     set.seed(random.seed)
     if (exists("savedSeed")) {
     on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
     }
     else {
     on.exit(rm(.Random.seed, envir = sys.frame()))
     }
     if (ghk) {
     if (is.null(nGHK)) {
     stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
     }
     else if (length(nGHK) != 1) {
     stop("argument 'nGHK' must be a single integer value")
     }
     else if (!is.numeric(nGHK)) {
     stop("argument 'nGHK' must be numeric")
     }
     else if (nGHK <= 0) {
     stop("argument 'nGHK' must be positive")
     }
     L <- try(t(chol(sigma)), silent = TRUE)
     if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
     }
     trunpt <- rep(NA, length(lower))
     above <- rep(NA, length(lower))
     for (i in 1:length(lower)) {
     if (lower[i] == -Inf) {
     trunpt[i] <- upper[i]
     above[i] <- 1
     }
     else if (upper[i] == Inf) {
     trunpt[i] <- lower[i]
     above[i] <- 0
     }
     else {
     stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
     " or the upper truncation point must be 'Inf'")
     }
     }
     sink(tempfile())
     on.exit(sink(), add = TRUE)
     result <- ghkvec(L = L, trunpt = trunpt, above = above,
     r = nGHK)
     result <- drop(result)
     }
     else {
     result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
     algorithm = algorithm, ...)
     }
     return(result)
     }
     <bytecode: 0x29ccd30>
     <environment: namespace:mvProbit>
     --- function search by body ---
     Function pmvnormWrap in namespace mvProbit has this body.
     ----------- END OF FAILURE REPORT --------------
     Error in if (class(L) == "try-error") { : the condition has length > 1
     > round( margEffCond, 3 )
     Error in round(margEffCond, 3) :
     non-numeric argument to mathematical function
     Execution halted
    Running the tests in ‘tests/mvProbitTest.R’ failed.
    Complete output:
     > library( "mvProbit" )
     Loading required package: mvtnorm
     Loading required package: maxLik
     Loading required package: miscTools
    
     Please cite the 'maxLik' package as:
     Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
    
     If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
     https://r-forge.r-project.org/projects/maxlik/
     Loading required package: abind
     > options( digits = 4 )
     >
     > ## generate a simulated data set
     > set.seed( 123 )
     > # number of observations
     > nObs <- 10
     >
     > # generate explanatory variables
     > xMat <- cbind(
     + const = rep( 1, nObs ),
     + x1 = as.numeric( rnorm( nObs ) > 0 ),
     + x2 = as.numeric( rnorm( nObs ) > 0 ),
     + x3 = rnorm( nObs ),
     + x4 = rnorm( nObs ) )
     >
     > # model coefficients
     > beta <- cbind( c( 0.8, 1.2, -1.0, 1.4, -0.8 ),
     + c( -0.6, 1.0, 0.6, -1.2, -1.6 ),
     + c( 0.5, -0.6, -0.7, 1.1, 1.2 ) )
     >
     > # covariance matrix of error terms
     > sigma <- miscTools::symMatrix( c( 1, 0.2, 0.4, 1, -0.1, 1 ) )
     >
     > # all parameters in a vector
     > allCoef <- c( c( beta ), sigma[ lower.tri( sigma ) ] )
     >
     > # generate dependent variables
     > yMatLin <- xMat %*% beta
     > yMat <- ( yMatLin + rmvnorm( nObs, sigma = sigma, pre0.9_9994 = TRUE ) ) > 0
     > colnames( yMat ) <- paste( "y", 1:3, sep = "" )
     > # (yMatLin > 0 )== yMat
     >
     > # unconditional expectations of dependent variables
     > yExp <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
     + sigma = sigma, data = as.data.frame( xMat ) )
     > round( yExp, 3 )
     V1 V2 V3
     1 0.021 0.725 0.194
     2 0.394 0.768 0.214
     3 0.125 0.788 0.196
     4 0.235 0.681 0.292
     5 0.680 0.435 0.579
     6 0.028 0.973 0.034
     7 0.958 0.186 0.784
     8 0.856 0.247 0.724
     9 0.061 0.968 0.034
     10 0.998 0.067 0.923
     > yExpA <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = allCoef,
     + data = as.data.frame( xMat ) )
     > all.equal( yExp, yExpA )
     [1] TRUE
     > yExp2 <- pnorm( yMatLin )
     > all.equal( yExp, as.data.frame( yExp2 ) )
     [1] TRUE
     >
     >
     > # conditional expectations of dependent variables
     > # (assuming that all other dependent variables are one)
     > yExpCond <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
     + sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
     + algorithm = GenzBretz() )
     > round( yExpCond, 3 )
     V1 V2 V3
     1 0.072 0.834 0.524
     2 0.660 0.779 0.314
     3 0.297 0.838 0.399
     4 0.448 0.727 0.461
     5 0.847 0.442 0.618
     6 0.139 0.985 0.162
     7 0.991 0.179 0.749
     8 0.950 0.244 0.710
     9 0.244 0.978 0.133
     10 1.000 0.065 0.892
     > yExpCondA <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = allCoef,
     + data = as.data.frame( xMat ), cond = TRUE,
     + algorithm = GenzBretz() )
     > all.equal( yExpCond, yExpCondA )
     [1] TRUE
     > yExpCond2 <- matrix( NA, nrow = nObs, ncol = ncol( yMat ) )
     > for( i in 1:nObs ) {
     + for( k in 1:ncol( yMat ) ) {
     + set.seed( 123 )
     + numerator <- pmvnorm( upper = yMatLin[ i, ], sigma = sigma )
     + set.seed( 123 )
     + denominator <- pmvnorm( upper = yMatLin[ i, -k ], sigma = sigma[ -k, -k ] )
     + yExpCond2[ i, k ] <- numerator / denominator
     + }
     + }
     > all.equal( yExpCond, as.data.frame( yExpCond2 ) )
     [1] TRUE
     > # now with explicitly specifying the algorithm
     > yExpCond3 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
     + sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
     + algorithm = GenzBretz )
     > all.equal( yExpCond, yExpCond3 )
     [1] TRUE
     > identical( yExpCond, yExpCond3 )
     [1] TRUE
     > # now with integrals obtained by the Miwa algorithm
     > yExpCond4 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
     + sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
     + algorithm = Miwa )
     > all.equal( yExpCond, yExpCond4, tol = 1e-3 )
     [1] TRUE
     > # now with integrals obtained by the Miwa algorithm, less precise
     > yExpCond5 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
     + sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
     + algorithm = Miwa( steps = 32 ) )
     > all.equal( yExpCond4, yExpCond5, tol = 1e-3 )
     [1] TRUE
     > # now with integrals obtained by the TVPACK algorithm
     > yExpCond6 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
     + sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
     + algorithm = TVPACK )
     > all.equal( yExpCond, yExpCond6, tol = 1e-3 )
     [1] TRUE
     > # now with integrals obtained by the TVPACK algorithm, less precise
     > yExpCond7 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
     + sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
     + algorithm = TVPACK( abseps = 0.5 ) )
     > all.equal( yExpCond6, yExpCond7, tol = 1e-3 )
     [1] "Component \"V1\": Mean relative difference: 0.03018"
     [2] "Component \"V2\": Mean relative difference: 0.04273"
     [3] "Component \"V3\": Mean relative difference: 0.04242"
     > # now with integrals obtained by the GHK algorithm
     > yExpCond8 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
     + sigma = sigma, data = as.data.frame( xMat ), cond = TRUE )
     ----------- FAILURE REPORT --------------
     --- failure: the condition has length > 1 ---
     --- srcref ---
     :
     --- package (from environment) ---
     mvProbit
     --- call from context ---
     pmvnormWrap(upper = xBeta[i, ], sigma = coef$sigma, algorithm = algorithm,
     nGHK = nGHK, random.seed = random.seed, ...)
     --- call from argument ---
     if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
     }
     --- R stacktrace ---
     where 1: pmvnormWrap(upper = xBeta[i, ], sigma = coef$sigma, algorithm = algorithm,
     nGHK = nGHK, random.seed = random.seed, ...)
     where 2: mvProbitExpInternal(yMat = yMat, xMat = xMat, coef = coef, sigma = sigma,
     cond = cond, algorithm = algorithm, nGHK = nGHK, random.seed = random.seed,
     ...)
     where 3: mvProbitExp(~x1 + x2 + x3 + x4, coef = c(beta), sigma = sigma,
     data = as.data.frame(xMat), cond = TRUE)
    
     --- value of length: 2 type: logical ---
     [1] FALSE FALSE
     --- function from context ---
     function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
     nGHK = NULL, ...)
     {
     if (!is.matrix(sigma)) {
     stop("argument 'sigma' must be a matrix")
     }
     else if (nrow(sigma) != ncol(sigma)) {
     stop("argument 'sigma' must be a square matrix")
     }
     else if (!all(is.numeric(sigma))) {
     stop("argument 'sigma must be a numeric matrix")
     }
     if (!all(is.numeric(lower))) {
     stop("argument 'lower' must be numeric")
     }
     else if (length(lower) == 1) {
     lower <- rep(lower, nrow(sigma))
     }
     else if (length(lower) != nrow(sigma)) {
     stop("argument 'lower' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     if (!all(is.numeric(upper))) {
     stop("argument 'upper' must be numeric")
     }
     else if (length(upper) == 1) {
     upper <- rep(upper, nrow(sigma))
     }
     else if (length(upper) != nrow(sigma)) {
     stop("argument 'upper' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     algOkay <- TRUE
     ghk <- FALSE
     if (!is.list(algorithm) && length(algorithm) != 1) {
     stop("argument 'algorithm' must be a single function",
     " or a single character string")
     }
     else if (is.function(algorithm)) {
     algResult <- do.call(algorithm, list())
     if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
     "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (is.character(algorithm)) {
     if (tolower(algorithm) == "ghk") {
     ghk <- TRUE
     }
     else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     if (!algOkay) {
     stop("argument 'algorithm' must be either one of the functions",
     " 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
     " \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
     }
     if (length(random.seed) != 1) {
     stop("argument 'random.seed' must be a single numerical value")
     }
     else if (!is.numeric(random.seed)) {
     stop("argument 'random.seed' must be numerical")
     }
     if (exists(".Random.seed")) {
     savedSeed <- .Random.seed
     }
     set.seed(random.seed)
     if (exists("savedSeed")) {
     on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
     }
     else {
     on.exit(rm(.Random.seed, envir = sys.frame()))
     }
     if (ghk) {
     if (is.null(nGHK)) {
     stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
     }
     else if (length(nGHK) != 1) {
     stop("argument 'nGHK' must be a single integer value")
     }
     else if (!is.numeric(nGHK)) {
     stop("argument 'nGHK' must be numeric")
     }
     else if (nGHK <= 0) {
     stop("argument 'nGHK' must be positive")
     }
     L <- try(t(chol(sigma)), silent = TRUE)
     if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
     }
     trunpt <- rep(NA, length(lower))
     above <- rep(NA, length(lower))
     for (i in 1:length(lower)) {
     if (lower[i] == -Inf) {
     trunpt[i] <- upper[i]
     above[i] <- 1
     }
     else if (upper[i] == Inf) {
     trunpt[i] <- lower[i]
     above[i] <- 0
     }
     else {
     stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
     " or the upper truncation point must be 'Inf'")
     }
     }
     sink(tempfile())
     on.exit(sink(), add = TRUE)
     result <- ghkvec(L = L, trunpt = trunpt, above = above,
     r = nGHK)
     result <- drop(result)
     }
     else {
     result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
     algorithm = algorithm, ...)
     }
     return(result)
     }
     <bytecode: 0x1b424a0>
     <environment: namespace:mvProbit>
     --- function search by body ---
     Function pmvnormWrap in namespace mvProbit has this body.
     ----------- END OF FAILURE REPORT --------------
     Error in if (class(L) == "try-error") { : the condition has length > 1
     Calls: mvProbitExp -> mvProbitExpInternal -> pmvnormWrap
     Execution halted
    Running the tests in ‘tests/pmvnormWrapTest.R’ failed.
    Complete output:
     > library( "mvProbit" )
     Loading required package: mvtnorm
     Loading required package: maxLik
     Loading required package: miscTools
    
     Please cite the 'maxLik' package as:
     Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
    
     If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
     https://r-forge.r-project.org/projects/maxlik/
     Loading required package: abind
     > library( "miscTools" )
     >
     > # covariance matrix
     > sigma <- miscTools::symMatrix( c( 1, 0.2, 0.4, 1, -0.1, 1 ) )
     >
     > ######## only upper ##########
     > upper <- c( -0.3, 0.7, -0.5 )
     > # Genz + Bretz (default)
     > pug <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
     + algorithm = GenzBretz(), random.seed = 123 )
     > print( pug )
     [1] 0.1360783
     attr(,"error")
     [1] 0.000118605
     attr(,"msg")
     [1] "Normal Completion"
     >
     > # Miwa (as function)
     > pum <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
     + algorithm = Miwa, random.seed = 123 )
     > print( pum )
     [1] 0.136069
     attr(,"error")
     [1] NA
     attr(,"msg")
     [1] "Normal Completion"
     > all.equal( pug, pum, check.attributes = FALSE, tol = 1e-4 )
     [1] TRUE
     >
     > # Miwa (as object returned from function Miwa())
     > pum1 <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
     + algorithm = Miwa(), random.seed = 123 )
     > all.equal( pum, pum1 )
     [1] TRUE
     >
     > # Miwa (as character string)
     > pum2 <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
     + algorithm = "Miwa", random.seed = 123 )
     > all.equal( pum, pum2 )
     [1] TRUE
     >
     > # TVPACK
     > put <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
     + algorithm = TVPACK, random.seed = 123 )
     > print( put )
     [1] 0.136069
     attr(,"error")
     [1] 1e-06
     attr(,"msg")
     [1] "Normal Completion"
     > all.equal( pug, put, check.attributes = FALSE, tol = 1e-4 )
     [1] TRUE
     > all.equal( pum, put, check.attributes = FALSE, tol = 1e-6 )
     [1] TRUE
     >
     > # GHK
     > pughk <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
     + algorithm = "ghk", random.seed = 123, nGHK = 1000 )
     ----------- FAILURE REPORT --------------
     --- failure: the condition has length > 1 ---
     --- srcref ---
     :
     --- package (from environment) ---
     mvProbit
     --- call from context ---
     mvProbit:::pmvnormWrap(upper = upper, sigma = sigma, algorithm = "ghk",
     random.seed = 123, nGHK = 1000)
     --- call from argument ---
     if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
     }
     --- R stacktrace ---
     where 1: mvProbit:::pmvnormWrap(upper = upper, sigma = sigma, algorithm = "ghk",
     random.seed = 123, nGHK = 1000)
    
     --- value of length: 2 type: logical ---
     [1] FALSE FALSE
     --- function from context ---
     function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
     nGHK = NULL, ...)
     {
     if (!is.matrix(sigma)) {
     stop("argument 'sigma' must be a matrix")
     }
     else if (nrow(sigma) != ncol(sigma)) {
     stop("argument 'sigma' must be a square matrix")
     }
     else if (!all(is.numeric(sigma))) {
     stop("argument 'sigma must be a numeric matrix")
     }
     if (!all(is.numeric(lower))) {
     stop("argument 'lower' must be numeric")
     }
     else if (length(lower) == 1) {
     lower <- rep(lower, nrow(sigma))
     }
     else if (length(lower) != nrow(sigma)) {
     stop("argument 'lower' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     if (!all(is.numeric(upper))) {
     stop("argument 'upper' must be numeric")
     }
     else if (length(upper) == 1) {
     upper <- rep(upper, nrow(sigma))
     }
     else if (length(upper) != nrow(sigma)) {
     stop("argument 'upper' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     algOkay <- TRUE
     ghk <- FALSE
     if (!is.list(algorithm) && length(algorithm) != 1) {
     stop("argument 'algorithm' must be a single function",
     " or a single character string")
     }
     else if (is.function(algorithm)) {
     algResult <- do.call(algorithm, list())
     if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
     "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (is.character(algorithm)) {
     if (tolower(algorithm) == "ghk") {
     ghk <- TRUE
     }
     else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     if (!algOkay) {
     stop("argument 'algorithm' must be either one of the functions",
     " 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
     " \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
     }
     if (length(random.seed) != 1) {
     stop("argument 'random.seed' must be a single numerical value")
     }
     else if (!is.numeric(random.seed)) {
     stop("argument 'random.seed' must be numerical")
     }
     if (exists(".Random.seed")) {
     savedSeed <- .Random.seed
     }
     set.seed(random.seed)
     if (exists("savedSeed")) {
     on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
     }
     else {
     on.exit(rm(.Random.seed, envir = sys.frame()))
     }
     if (ghk) {
     if (is.null(nGHK)) {
     stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
     }
     else if (length(nGHK) != 1) {
     stop("argument 'nGHK' must be a single integer value")
     }
     else if (!is.numeric(nGHK)) {
     stop("argument 'nGHK' must be numeric")
     }
     else if (nGHK <= 0) {
     stop("argument 'nGHK' must be positive")
     }
     L <- try(t(chol(sigma)), silent = TRUE)
     if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
     }
     trunpt <- rep(NA, length(lower))
     above <- rep(NA, length(lower))
     for (i in 1:length(lower)) {
     if (lower[i] == -Inf) {
     trunpt[i] <- upper[i]
     above[i] <- 1
     }
     else if (upper[i] == Inf) {
     trunpt[i] <- lower[i]
     above[i] <- 0
     }
     else {
     stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
     " or the upper truncation point must be 'Inf'")
     }
     }
     sink(tempfile())
     on.exit(sink(), add = TRUE)
     result <- ghkvec(L = L, trunpt = trunpt, above = above,
     r = nGHK)
     result <- drop(result)
     }
     else {
     result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
     algorithm = algorithm, ...)
     }
     return(result)
     }
     <bytecode: 0x3dcc488>
     <environment: namespace:mvProbit>
     --- function search by body ---
     Function pmvnormWrap in namespace mvProbit has this body.
     ----------- END OF FAILURE REPORT --------------
     Error in if (class(L) == "try-error") { : the condition has length > 1
     Calls: <Anonymous>
     Execution halted
Flavor: r-devel-linux-x86_64-fedora-clang

Version: 0.1-8
Check: examples
Result: ERROR
    Running examples in ‘mvProbit-Ex.R’ failed
    The error most likely occurred in:
    
    > ### Name: mvProbit
    > ### Title: Estimation of Multivariate Probit Models
    > ### Aliases: mvProbit print.mvProbit
    > ### Keywords: models regression
    >
    > ### ** Examples
    >
    > ## generate a simulated data set
    > set.seed( 123 )
    > # number of observations
    > nObs <- 50
    >
    > # generate explanatory variables
    > xMat <- cbind(
    + const = rep( 1, nObs ),
    + x1 = as.numeric( rnorm( nObs ) > 0 ),
    + x2 = rnorm( nObs ) )
    >
    > # model coefficients
    > beta <- cbind( c( 0.8, 1.2, -0.8 ),
    + c( -0.6, 1.0, -1.6 ),
    + c( 0.5, -0.6, 1.2 ) )
    >
    > # covariance matrix of error terms
    > library( miscTools )
    > sigma <- symMatrix( c( 1, 0.2, 0.4, 1, -0.1, 1 ) )
    >
    > # generate dependent variables
    > yMatLin <- xMat %*% beta
    > yMat <- ( yMatLin + rmvnorm( nObs, sigma = sigma ) ) > 0
    > colnames( yMat ) <- paste( "y", 1:3, sep = "" )
    >
    > # estimation (BHHH optimizer and GHK algorithm)
    > estResult <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
    + data = as.data.frame( cbind( xMat, yMat ) ), iterlim = 1, nGHK = 50 )
    Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
     ----------- FAILURE REPORT --------------
     --- failure: the condition has length > 1 ---
     --- srcref ---
    :
     --- package (from environment) ---
    mvProbit
     --- call from context ---
    pmvnormWrap(upper = xBetaTmp, sigma = sigmaTmp, algorithm = algorithm,
     nGHK = nGHK, random.seed = randomSeed, ...)
     --- call from argument ---
    if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
    }
     --- R stacktrace ---
    where 1: pmvnormWrap(upper = xBetaTmp, sigma = sigmaTmp, algorithm = algorithm,
     nGHK = nGHK, random.seed = randomSeed, ...)
    where 2: mvProbitLogLikInternal(yMat = yMat, xMat = xMat, coef = param,
     sigma = NULL, algorithm = llAlgorithm, nGHK = llNGHK, returnGrad = llIntGrad,
     oneSidedGrad = llOneSidedGrad, eps = llEps, randomSeed = llRandom.seed,
     ...)
    where 3: fnOrig(theta, ...)
    where 4: fn(start1, fixed = fixed, sumObs = TRUE, returnHessian = returnHessian,
     ...)
    where 5: maxNRCompute(fn = function (theta, fnOrig, gradOrig = NULL, hessOrig = NULL,
     fixed, sumObs = FALSE, returnHessian = TRUE, ...)
    {
     nParam <- length(theta)
     f <- fnOrig(theta, ...)
     if (any(is.na(f))) {
     attr(f, "gradient") <- NA
     attr(f, "hessian") <- NA
     return(f)
     }
     gr <- attr(f, "gradient")
     if (is.null(gr)) {
     if (!is.null(gradOrig)) {
     gr <- gradOrig(theta, ...)
     }
     else {
     gr <- numericGradient(f = fnOrig, t0 = theta, fixed = fixed,
     ...)
     }
     }
     if (is.matrix(gr)) {
     if (ncol(gr) != length(theta)) {
     stop(paste0("if gradient is a matrix, it must have length(parameter) colums (currently ",
     length(theta), "), not ", ncol(gr)))
     }
     activeGr <- gr[, !fixed]
     }
     else {
     activeGr <- gr[!fixed]
     }
     if (any(is.na(activeGr))) {
     attr(f, "gradient") <- gr
     attr(f, "hessian") <- NA
     return(f)
     }
     if (observationGradient(gr, length(theta))) {
     gr <- as.matrix(gr)
     }
     if (is.null(dim(gr))) {
     gr[fixed] <- NA
     }
     else {
     gr[, fixed] <- NA
     }
     if (isTRUE(returnHessian)) {
     h <- attr(f, "hessian")
     if (is.null(h)) {
     if (!is.null(hessOrig)) {
     h <- as.matrix(hessOrig(theta, ...))
     }
     else {
     llFunc <- function(theta, ...) {
     return(sum(fnOrig(theta, ...)))
     }
     if (!is.null(attr(f, "gradient"))) {
     gradFunc <- function(theta, ...) {
     return(sumGradients(attr(fnOrig(theta, ...),
     "gradient"), nParam))
     }
     }
     else if (!is.null(gradOrig)) {
     gradFunc <- function(theta, ...) {
     return(sumGradients(gradOrig(theta, ...),
     nParam))
     }
     }
     else {
     gradFunc <- NULL
     }
     h <- numericHessian(f = llFunc, grad = gradFunc,
     t0 = theta, fixed = fixed, ...)
     }
     }
     if ((dim(h)[1] != nParam) | (dim(h)[2] != nParam)) {
     stop("Wrong hessian dimension. Needed ", nParam,
     "x", nParam, " but supplied ", dim(h)[1], "x",
     dim(h)[2])
     }
     else {
     h[fixed, ] <- NA
     h[, fixed] <- NA
     }
     }
     else if (tolower(returnHessian) == "bhhh") {
     h <- NULL
     if (is.null(dim(gr)) & any(is.na(gr[!fixed]))) {
     h <- NA
     }
     else if (is.matrix(gr)) {
     if (any(is.na(gr[, !fixed]))) {
     h <- NA
     }
     }
     if (is.null(h)) {
     checkBhhhGrad(g = gr, theta = theta, analytic = (!is.null(attr(f,
     "gradient")) || !is.null(gradOrig)), fixed = fixed)
     h <- -crossprod(gr)
     }
     attr(h, "type") = "BHHH"
     }
     else {
     h <- NULL
     }
     if (sumObs) {
     f <- sumKeepAttr(f)
     }
     if (sumObs) {
     gr <- sumGradients(gr, nParam)
     }
     if (!is.null(gradOrig) && !is.null(attr(f, "gradient"))) {
     attr(f, "gradBoth") <- TRUE
     }
     if (!is.null(hessOrig) && !is.null(attr(f, "hessian"))) {
     attr(f, "hessBoth") <- TRUE
     }
     attr(f, "gradient") <- gr
     attr(f, "hessian") <- h
     return(f)
    }, fnOrig = function (param, yMat, xMat, llAlgorithm, llNGHK,
     llIntGrad, llOneSidedGrad, llEps, llRandom.seed, ...)
    {
     logLikVal <- mvProbitLogLikInternal(yMat = yMat, xMat = xMat,
     coef = param, sigma = NULL, algorithm = llAlgorithm,
     nGHK = llNGHK, returnGrad = llIntGrad, oneSidedGrad = llOneSidedGrad,
     eps = llEps, randomSeed = llRandom.seed, ...)
     return(logLikVal)
    }, start = c(b_1_0 = 0.598378200054935, b_1_1 = 6.336646267785,
    b_1_2 = -1.01231147815159, b_2_0 = -1.03564793824529, b_2_1 = 1.26302210717712,
    b_2_2 = -1.70638777681512, b_3_0 = 0.410893287564513, b_3_1 = -0.373151392833251,
    b_3_2 = 2.92074017047352, R_1_2 = 0.172617843070694, R_1_3 = -0.0160393443285098,
    R_2_3 = -0.0704265674853818), finalHessian = "BHHH", bhhhHessian = TRUE,
     fixed = c(b_1_0 = FALSE, b_1_1 = FALSE, b_1_2 = FALSE, b_2_0 = FALSE,
     b_2_1 = FALSE, b_2_2 = FALSE, b_3_0 = FALSE, b_3_1 = FALSE,
     b_3_2 = FALSE, R_1_2 = FALSE, R_1_3 = FALSE, R_2_3 = FALSE
     ), control = new("MaxControl", tol = 1e-08, reltol = 1.49011611938477e-08,
     gradtol = 1e-06, steptol = 1e-10, lambdatol = 1e-06,
     qrtol = 1e-10, qac = "stephalving", marquardt_lambda0 = 0.01,
     marquardt_lambdaStep = 2, marquardt_maxLambda = 1e+12,
     nm_alpha = 1, nm_beta = 0.5, nm_gamma = 2, sann_cand = NULL,
     sann_temp = 10, sann_tmax = 10L, sann_randomSeed = 123L,
     iterlim = 1L, printLevel = 0L), yMat = c(0, 1, 1, 1,
     1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1,
     1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1,
     1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1,
     0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1,
     1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1,
     1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0,
     1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1,
     1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0), xMat = c(1, 1, 1,
     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0,
     1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0,
     1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0,
     1, 0, 0.253318513994755, -0.028546755348703, -0.0428704572913161,
     1.36860228401446, -0.225770985659268, 1.51647060442954, -1.54875280423022,
     0.584613749636069, 0.123854243844614, 0.215941568743973,
     0.379639482759882, -0.502323453109302, -0.33320738366942,
     -1.01857538310709, -1.07179122647558, 0.303528641404258,
     0.448209778629426, 0.0530042267305041, 0.922267467879738,
     2.05008468562714, -0.491031166056535, -2.30916887564081,
     1.00573852446226, -0.709200762582393, -0.688008616467358,
     1.0255713696967, -0.284773007051009, -1.22071771225454, 0.18130347974915,
     -0.138891362439045, 0.00576418589988693, 0.38528040112633,
     -0.370660031792409, 0.644376548518833, -0.220486561818751,
     0.331781963915697, 1.09683901314935, 0.435181490833803, -0.325931585531227,
     1.14880761845109, 0.993503855962119, 0.54839695950807, 0.238731735111441,
     -0.627906076039371, 1.36065244853001, -0.600259587147127,
     2.18733299301658, 1.53261062618519, -0.235700359100477, -1.02642090030678
     ), llAlgorithm = "GHK", llNGHK = 50, llIntGrad = TRUE, llOneSidedGrad = FALSE,
     llEps = 1e-06, llRandom.seed = 123)
    where 6: eval(as.call(cl))
    where 7: eval(as.call(cl))
    where 8: maxNR(fn = fn, grad = grad, hess = hess, start = start, finalHessian = finalHessian,
     bhhhHessian = TRUE, ...)
    where 9: maxRoutine(fn = logLik, grad = grad, hess = hess, start = start,
     constraints = constraints, ...)
    where 10: maxLik(logLik = logLik, start = start, method = method, finalHessian = finalHessian,
     yMat = yMat, xMat = xMat, llAlgorithm = algorithm, llNGHK = nGHK,
     llIntGrad = intGrad, llOneSidedGrad = oneSidedGrad, llEps = eps,
     llRandom.seed = random.seed, ...)
    where 11: mvProbit(cbind(y1, y2, y3) ~ x1 + x2, data = as.data.frame(cbind(xMat,
     yMat)), iterlim = 1, nGHK = 50)
    
     --- value of length: 2 type: logical ---
    [1] FALSE FALSE
     --- function from context ---
    function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
     nGHK = NULL, ...)
    {
     if (!is.matrix(sigma)) {
     stop("argument 'sigma' must be a matrix")
     }
     else if (nrow(sigma) != ncol(sigma)) {
     stop("argument 'sigma' must be a square matrix")
     }
     else if (!all(is.numeric(sigma))) {
     stop("argument 'sigma must be a numeric matrix")
     }
     if (!all(is.numeric(lower))) {
     stop("argument 'lower' must be numeric")
     }
     else if (length(lower) == 1) {
     lower <- rep(lower, nrow(sigma))
     }
     else if (length(lower) != nrow(sigma)) {
     stop("argument 'lower' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     if (!all(is.numeric(upper))) {
     stop("argument 'upper' must be numeric")
     }
     else if (length(upper) == 1) {
     upper <- rep(upper, nrow(sigma))
     }
     else if (length(upper) != nrow(sigma)) {
     stop("argument 'upper' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     algOkay <- TRUE
     ghk <- FALSE
     if (!is.list(algorithm) && length(algorithm) != 1) {
     stop("argument 'algorithm' must be a single function",
     " or a single character string")
     }
     else if (is.function(algorithm)) {
     algResult <- do.call(algorithm, list())
     if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
     "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (is.character(algorithm)) {
     if (tolower(algorithm) == "ghk") {
     ghk <- TRUE
     }
     else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     if (!algOkay) {
     stop("argument 'algorithm' must be either one of the functions",
     " 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
     " \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
     }
     if (length(random.seed) != 1) {
     stop("argument 'random.seed' must be a single numerical value")
     }
     else if (!is.numeric(random.seed)) {
     stop("argument 'random.seed' must be numerical")
     }
     if (exists(".Random.seed")) {
     savedSeed <- .Random.seed
     }
     set.seed(random.seed)
     if (exists("savedSeed")) {
     on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
     }
     else {
     on.exit(rm(.Random.seed, envir = sys.frame()))
     }
     if (ghk) {
     if (is.null(nGHK)) {
     stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
     }
     else if (length(nGHK) != 1) {
     stop("argument 'nGHK' must be a single integer value")
     }
     else if (!is.numeric(nGHK)) {
     stop("argument 'nGHK' must be numeric")
     }
     else if (nGHK <= 0) {
     stop("argument 'nGHK' must be positive")
     }
     L <- try(t(chol(sigma)), silent = TRUE)
     if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
     }
     trunpt <- rep(NA, length(lower))
     above <- rep(NA, length(lower))
     for (i in 1:length(lower)) {
     if (lower[i] == -Inf) {
     trunpt[i] <- upper[i]
     above[i] <- 1
     }
     else if (upper[i] == Inf) {
     trunpt[i] <- lower[i]
     above[i] <- 0
     }
     else {
     stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
     " or the upper truncation point must be 'Inf'")
     }
     }
     sink(tempfile())
     on.exit(sink(), add = TRUE)
     result <- ghkvec(L = L, trunpt = trunpt, above = above,
     r = nGHK)
     result <- drop(result)
     }
     else {
     result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
     algorithm = algorithm, ...)
     }
     return(result)
    }
    <bytecode: 0x36e8a08>
    <environment: namespace:mvProbit>
     --- function search by body ---
    Function pmvnormWrap in namespace mvProbit has this body.
     ----------- END OF FAILURE REPORT --------------
    Error in if (class(L) == "try-error") { : the condition has length > 1
    Calls: mvProbit ... fn -> fnOrig -> mvProbitLogLikInternal -> pmvnormWrap
    Execution halted
Flavor: r-devel-linux-x86_64-fedora-gcc

Version: 0.1-8
Check: tests
Result: ERROR
     Running ‘mvProbitEst.R’ [33s/50s]
     Running ‘mvProbitMargEff2.R’
     Running ‘mvProbitTest.R’
     Running ‘mvProbitTest5.R’ [144s/204s]
     Comparing ‘mvProbitTest5.Rout’ to ‘mvProbitTest5.Rout.save’ ...1604c1604
    < 5 0.00 -0.01 0 0.00 0.00 0.00 0.00
    ---
    > 5 0.00 0.00 0 0.00 0.00 0.00 0.00
    1618,1620c1618,1620
    < 1 0 0.00 0 0 0 0 0
    < 5 0 -0.01 0 0 0 0 0
    < 10 0 0.00 0 0 0 0 0
    ---
    > 1 0 0 0 0 0 0 0
    > 5 0 0 0 0 0 0 0
    > 10 0 0 0 0 0 0 0
    1728c1728
    < 5 -0.01
    ---
    > 5 0.00
     Running ‘pmvnormWrapTest.R’
    Running the tests in ‘tests/mvProbitEst.R’ failed.
    Complete output:
     > library( "mvProbit" )
     Loading required package: mvtnorm
     Loading required package: maxLik
     Loading required package: miscTools
    
     Please cite the 'maxLik' package as:
     Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
    
     If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
     https://r-forge.r-project.org/projects/maxlik/
     Loading required package: abind
     > options( digits = 4 )
     >
     > ## generate a simulated data set
     > set.seed( 123 )
     > # number of observations
     > nObs <- 50
     >
     > # generate explanatory variables
     > xMat <- cbind(
     + const = rep( 1, nObs ),
     + x1 = as.numeric( rnorm( nObs ) > 0 ),
     + x2 = rnorm( nObs ) )
     >
     > # model coefficients
     > beta <- cbind( c( 0.8, 1.2, -0.8 ),
     + c( -0.6, 1.0, -1.6 ),
     + c( 0.5, -0.6, 1.2 ) )
     >
     > # covariance matrix of error terms
     > sigma <- miscTools::symMatrix( c( 1, 0.2, 0.4, 1, -0.1, 1 ) )
     >
     > # generate dependent variables
     > yMatLin <- xMat %*% beta
     > yMat <- ( yMatLin + rmvnorm( nObs, sigma = sigma, pre0.9_9994 = TRUE ) ) > 0
     > colnames( yMat ) <- paste( "y", 1:3, sep = "" )
     >
     > # create data frame
     > dat <- as.data.frame( cbind( xMat, yMat ) )
     >
     > # estimation with the BHHH algorithm, two-sided gradients
     > estResultBHHH <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + start = c( beta ), startSigma = sigma,
     + data = dat, tol = 0.5,
     + algorithm = GenzBretz() )
     > print( estResultBHHH )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, algorithm = GenzBretz(), tol = 0.5)
    
     Coefficients:
     b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
     1.02508 1.30536 -1.02809 -0.78248 1.14318 -1.66938 0.58229 -0.69667
     b_3_2 R_1_2 R_1_3 R_2_3
     1.23283 0.86322 0.36022 -0.08406
    
     > summary( estResultBHHH )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, algorithm = GenzBretz(), tol = 0.5)
    
     Coefficients:
     Estimate Std. error t value Pr(> t)
     b_1_0 1.02508 0.44548 2.301 0.0214 *
     b_1_1 1.30536 0.73511 1.776 0.0758 .
     b_1_2 -1.02809 0.68851 -1.493 0.1354
     b_2_0 -0.78248 0.38798 -2.017 0.0437 *
     b_2_1 1.14318 0.53026 2.156 0.0311 *
     b_2_2 -1.66938 0.81308 -2.053 0.0401 *
     b_3_0 0.58229 0.64107 0.908 0.3637
     b_3_1 -0.69667 0.89818 -0.776 0.4380
     b_3_2 1.23283 0.62452 1.974 0.0484 *
     R_1_2 0.86322 18.60071 0.046 0.9630
     R_1_3 0.36022 0.64712 0.557 0.5778
     R_2_3 -0.08406 0.34909 -0.241 0.8097
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     BHHH maximisation, 1 iterations
     Return code 2: successive function values within tolerance limit
     Log-likelihood: -56.91 on 12 Df
    
     > logLik( estResultBHHH )
     'log Lik.' -56.91 (df=12)
     > estResultBHHHA <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + start = c( beta, sigma[ lower.tri( sigma ) ] ),
     + data = dat, tol = 0.5,
     + algorithm = GenzBretz() )
     > all.equal( estResultBHHH, estResultBHHHA )
     [1] "Component \"call\": target, current do not match when deparsed"
     >
     > # estimation with the BHHH algorithm, one-sided gradients
     > estResultBHHH1 <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + start = c( beta ), startSigma = sigma,
     + data = dat, tol = 0.5,
     + algorithm = GenzBretz(), oneSidedGrad = TRUE )
     > print( estResultBHHH1 )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, algorithm = GenzBretz(), oneSidedGrad = TRUE,
     tol = 0.5)
    
     Coefficients:
     b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
     1.02508 1.30536 -1.02809 -0.78248 1.14318 -1.66938 0.58229 -0.69667
     b_3_2 R_1_2 R_1_3 R_2_3
     1.23283 0.86322 0.36022 -0.08406
    
     > summary( estResultBHHH1 )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, algorithm = GenzBretz(), oneSidedGrad = TRUE,
     tol = 0.5)
    
     Coefficients:
     Estimate Std. error t value Pr(> t)
     b_1_0 1.02508 0.44548 2.301 0.0214 *
     b_1_1 1.30536 0.73511 1.776 0.0758 .
     b_1_2 -1.02809 0.68851 -1.493 0.1354
     b_2_0 -0.78248 0.38797 -2.017 0.0437 *
     b_2_1 1.14318 0.53026 2.156 0.0311 *
     b_2_2 -1.66938 0.81307 -2.053 0.0401 *
     b_3_0 0.58229 0.64106 0.908 0.3637
     b_3_1 -0.69667 0.89817 -0.776 0.4380
     b_3_2 1.23283 0.62452 1.974 0.0484 *
     R_1_2 0.86322 18.60041 0.046 0.9630
     R_1_3 0.36022 0.64712 0.557 0.5778
     R_2_3 -0.08406 0.34909 -0.241 0.8097
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     BHHH maximisation, 1 iterations
     Return code 2: successive function values within tolerance limit
     Log-likelihood: -56.91 on 12 Df
    
     > logLik( estResultBHHH1 )
     'log Lik.' -56.91 (df=12)
     > all.equal( estResultBHHH, estResultBHHH1, tol = 1e-5 )
     [1] "Component \"call\": target, current do not match when deparsed"
     >
     > # estimation with the BFGS algorithm, two-sided gradients
     > estResultBFGS <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + start = c( beta ), startSigma = sigma, method = "BFGS",
     + data = dat,
     + reltol = 0.5, algorithm = GenzBretz() )
     > print( estResultBFGS )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, method = "BFGS", algorithm = GenzBretz(),
     reltol = 0.5)
    
     Coefficients:
     b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
     0.53607 1.16101 -1.14212 -0.44952 1.05496 -1.57624 0.40751 -0.85400
     b_3_2 R_1_2 R_1_3 R_2_3
     1.25958 0.29483 0.32822 -0.02382
    
     > summary( estResultBFGS )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, method = "BFGS", algorithm = GenzBretz(),
     reltol = 0.5)
    
     Coefficients:
     Estimate Std. error t value Pr(> t)
     b_1_0 0.53607 0.82709 0.648 0.5169
     b_1_1 1.16101 0.91515 1.269 0.2046
     b_1_2 -1.14212 0.93847 -1.217 0.2236
     b_2_0 -0.44952 0.79215 -0.567 0.5704
     b_2_1 1.05496 0.88398 1.193 0.2327
     b_2_2 -1.57624 0.73187 -2.154 0.0313 *
     b_3_0 0.40751 0.48897 0.833 0.4046
     b_3_1 -0.85400 0.84871 -1.006 0.3143
     b_3_2 1.25958 0.68040 1.851 0.0641 .
     R_1_2 0.29483 1.85817 0.159 0.8739
     R_1_3 0.32822 0.70923 0.463 0.6435
     R_2_3 -0.02382 0.34466 -0.069 0.9449
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     BFGS maximization, 4 iterations
     Return code 0: successful convergence
     Log-likelihood: -52.94 on 12 Df
    
     > logLik( estResultBFGS )
     'log Lik.' -52.94 (df=12)
     >
     > # estimation with the BFGS algorithm, one-sided gradients
     > estResultBFGS1 <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + start = c( beta ), startSigma = sigma, method = "BFGS",
     + data = dat,
     + reltol = 0.5, algorithm = GenzBretz(), oneSidedGrad = TRUE )
     > print( estResultBFGS1 )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, method = "BFGS", algorithm = GenzBretz(),
     oneSidedGrad = TRUE, reltol = 0.5)
    
     Coefficients:
     b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
     0.53607 1.16101 -1.14213 -0.44952 1.05496 -1.57624 0.40751 -0.85400
     b_3_2 R_1_2 R_1_3 R_2_3
     1.25958 0.29483 0.32822 -0.02382
    
     > summary( estResultBFGS1 )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, method = "BFGS", algorithm = GenzBretz(),
     oneSidedGrad = TRUE, reltol = 0.5)
    
     Coefficients:
     Estimate Std. error t value Pr(> t)
     b_1_0 0.53607 0.82709 0.648 0.5169
     b_1_1 1.16101 0.91515 1.269 0.2046
     b_1_2 -1.14213 0.93847 -1.217 0.2236
     b_2_0 -0.44952 0.79215 -0.567 0.5704
     b_2_1 1.05496 0.88398 1.193 0.2327
     b_2_2 -1.57624 0.73187 -2.154 0.0313 *
     b_3_0 0.40751 0.48897 0.833 0.4046
     b_3_1 -0.85400 0.84871 -1.006 0.3143
     b_3_2 1.25958 0.68040 1.851 0.0641 .
     R_1_2 0.29483 1.85817 0.159 0.8739
     R_1_3 0.32822 0.70923 0.463 0.6435
     R_2_3 -0.02382 0.34466 -0.069 0.9449
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     BFGS maximization, 4 iterations
     Return code 0: successful convergence
     Log-likelihood: -52.94 on 12 Df
    
     > logLik( estResultBFGS1 )
     'log Lik.' -52.94 (df=12)
     > all.equal( estResultBFGS, estResultBFGS1, tol = 1e-5 )
     [1] "Component \"call\": target, current do not match when deparsed"
     >
     > # estimation with the BFGS algorithm, one-sided gradients, no starting values
     > estResultBFGS1a <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + data = dat, method = "BFGS",
     + reltol = 0.5, algorithm = GenzBretz(), oneSidedGrad = TRUE )
     > print( estResultBFGS1a )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, method = "BFGS",
     algorithm = GenzBretz(), oneSidedGrad = TRUE, reltol = 0.5)
    
     Coefficients:
     b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
     0.79353 1.88367 -1.74922 -0.40036 0.91577 -1.52106 0.97790 -1.69433
     b_3_2 R_1_2 R_1_3 R_2_3
     1.57625 0.75270 0.23721 -0.03453
    
     > summary( estResultBFGS1a )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, method = "BFGS",
     algorithm = GenzBretz(), oneSidedGrad = TRUE, reltol = 0.5)
    
     Coefficients:
     Estimate Std. error t value Pr(> t)
     b_1_0 0.79353 0.90412 0.878 0.3801
     b_1_1 1.88367 0.90668 2.078 0.0378 *
     b_1_2 -1.74922 0.97400 -1.796 0.0725 .
     b_2_0 -0.40036 0.69056 -0.580 0.5621
     b_2_1 0.91577 0.97998 0.934 0.3501
     b_2_2 -1.52106 0.72006 -2.112 0.0347 *
     b_3_0 0.97790 0.64497 1.516 0.1295
     b_3_1 -1.69433 0.94788 -1.787 0.0739 .
     b_3_2 1.57625 0.73399 2.148 0.0318 *
     R_1_2 0.75270 4.32729 0.174 0.8619
     R_1_3 0.23721 0.90535 0.262 0.7933
     R_2_3 -0.03453 0.35889 -0.096 0.9234
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     BFGS maximization, 3 iterations
     Return code 0: successful convergence
     Log-likelihood: -49.36 on 12 Df
    
     > logLik( estResultBFGS1a )
     'log Lik.' -49.36 (df=12)
     >
     > # estimation with the BFGS algorithm, one-sided gradients, no starting values for beta
     > estResultBFGS1b <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + startSigma = sigma, data = dat,
     + method = "BFGS", reltol = 0.5, algorithm = GenzBretz(), oneSidedGrad = TRUE )
     > print( estResultBFGS1b )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, startSigma = sigma,
     method = "BFGS", algorithm = GenzBretz(), oneSidedGrad = TRUE,
     reltol = 0.5)
    
     Coefficients:
     b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1 b_3_2 R_1_2
     0.8110 1.9060 -1.7627 -0.3770 0.9441 -1.5332 0.9372 -1.7498 1.6024 0.6817
     R_1_3 R_2_3
     0.4044 0.1663
    
     > summary( estResultBFGS1b )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, startSigma = sigma,
     method = "BFGS", algorithm = GenzBretz(), oneSidedGrad = TRUE,
     reltol = 0.5)
    
     Coefficients:
     Estimate Std. error t value Pr(> t)
     b_1_0 0.8110 0.9297 0.872 0.3830
     b_1_1 1.9060 0.9160 2.081 0.0375 *
     b_1_2 -1.7627 1.0131 -1.740 0.0819 .
     b_2_0 -0.3770 0.7049 -0.535 0.5928
     b_2_1 0.9441 1.0345 0.913 0.3615
     b_2_2 -1.5332 0.7517 -2.040 0.0414 *
     b_3_0 0.9372 0.6291 1.490 0.1363
     b_3_1 -1.7498 0.9996 -1.751 0.0800 .
     b_3_2 1.6024 0.7848 2.042 0.0412 *
     R_1_2 0.6817 3.4094 0.200 0.8415
     R_1_3 0.4044 0.8024 0.504 0.6143
     R_2_3 0.1663 0.3389 0.491 0.6235
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     BFGS maximization, 3 iterations
     Return code 0: successful convergence
     Log-likelihood: -49.57 on 12 Df
    
     > logLik( estResultBFGS1b )
     'log Lik.' -49.57 (df=12)
     >
     > # estimation with the BFGS algorithm, one-sided gradients, no starting values for sigma
     > estResultBFGS1s <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + start = c( beta ), data = dat,
     + method = "BFGS", reltol = 0.5, algorithm = GenzBretz(), oneSidedGrad = TRUE )
     > print( estResultBFGS1s )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     method = "BFGS", algorithm = GenzBretz(), oneSidedGrad = TRUE,
     reltol = 0.5)
    
     Coefficients:
     b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
     0.55814 1.14958 -1.12747 -0.44426 1.07325 -1.57226 0.37255 -0.85470
     b_3_2 R_1_2 R_1_3 R_2_3
     1.23011 0.34981 -0.01177 0.09435
    
     > summary( estResultBFGS1s )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     method = "BFGS", algorithm = GenzBretz(), oneSidedGrad = TRUE,
     reltol = 0.5)
    
     Coefficients:
     Estimate Std. error t value Pr(> t)
     b_1_0 0.55814 0.86615 0.644 0.5193
     b_1_1 1.14958 0.90607 1.269 0.2045
     b_1_2 -1.12747 0.94191 -1.197 0.2313
     b_2_0 -0.44426 0.85526 -0.519 0.6035
     b_2_1 1.07325 0.99314 1.081 0.2798
     b_2_2 -1.57226 0.78557 -2.001 0.0453 *
     b_3_0 0.37255 0.56172 0.663 0.5072
     b_3_1 -0.85470 0.95214 -0.898 0.3694
     b_3_2 1.23011 0.71506 1.720 0.0854 .
     R_1_2 0.34981 2.02289 0.173 0.8627
     R_1_3 -0.01177 0.70776 -0.017 0.9867
     R_2_3 0.09435 0.33920 0.278 0.7809
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     BFGS maximization, 4 iterations
     Return code 0: successful convergence
     Log-likelihood: -52.78 on 12 Df
    
     > logLik( estResultBFGS1s )
     'log Lik.' -52.78 (df=12)
     >
     > # estimation with the BFGS algorithm, Miwa algorithm for obtaining integrals
     > estResultBFGSm <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + start = c( beta ), startSigma = sigma,
     + data = dat, method = "BFGS",
     + reltol = 0.5, algorithm = Miwa( steps = 64 ) )
     There were 13 warnings (use warnings() to see them)
     > print( estResultBFGSm )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, method = "BFGS", algorithm = Miwa(steps = 64),
     reltol = 0.5)
    
     Coefficients:
     b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
     0.53609 1.16102 -1.14213 -0.44952 1.05496 -1.57624 0.40753 -0.85400
     b_3_2 R_1_2 R_1_3 R_2_3
     1.25960 0.29491 0.32819 -0.02386
    
     > summary( estResultBFGSm )
    
     Call:
     mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
     startSigma = sigma, method = "BFGS", algorithm = Miwa(steps = 64),
     reltol = 0.5)
    
     Coefficients:
     Estimate Std. error t value Pr(> t)
     b_1_0 0.53609 0.82702 0.648 0.5168
     b_1_1 1.16102 0.91511 1.269 0.2045
     b_1_2 -1.14213 0.93837 -1.217 0.2235
     b_2_0 -0.44952 0.79215 -0.567 0.5704
     b_2_1 1.05496 0.88386 1.194 0.2326
     b_2_2 -1.57624 0.73173 -2.154 0.0312 *
     b_3_0 0.40753 0.48894 0.834 0.4046
     b_3_1 -0.85400 0.84860 -1.006 0.3142
     b_3_2 1.25960 0.68024 1.852 0.0641 .
     R_1_2 0.29491 1.85789 0.159 0.8739
     R_1_3 0.32819 0.70903 0.463 0.6435
     R_2_3 -0.02386 0.34468 -0.069 0.9448
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     BFGS maximization, 4 iterations
     Return code 0: successful convergence
     Log-likelihood: -52.94 on 12 Df
    
     > logLik( estResultBFGSm )
     'log Lik.' -52.94 (df=12)
     > all.equal( estResultBFGS, estResultBFGSm, tol = 1e-3 )
     [1] "Component \"call\": target, current do not match when deparsed"
     >
     > # estimation with the BFGS algorithm, GHK algorithm for obtaining integrals
     > estResultBFGSg <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
     + start = c( beta ), startSigma = sigma,
     + data = dat, method = "BFGS",
     + reltol = 0.5 )
     ----------- FAILURE REPORT --------------
     --- failure: the condition has length > 1 ---
     --- srcref ---
     :
     --- package (from environment) ---
     mvProbit
     --- call from context ---
     pmvnormWrap(upper = xBetaTmp, sigma = sigmaTmp, algorithm = algorithm,
     nGHK = nGHK, random.seed = randomSeed, ...)
     --- call from argument ---
     if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
     }
     --- R stacktrace ---
     where 1: pmvnormWrap(upper = xBetaTmp, sigma = sigmaTmp, algorithm = algorithm,
     nGHK = nGHK, random.seed = randomSeed, ...)
     where 2: mvProbitLogLikInternal(yMat = yMat, xMat = xMat, coef = param,
     sigma = NULL, algorithm = llAlgorithm, nGHK = llNGHK, returnGrad = llIntGrad,
     oneSidedGrad = llOneSidedGrad, eps = llEps, randomSeed = llRandom.seed,
     ...)
     where 3: fnOrig(theta, ...)
     where 4: logLikFunc(theta, fnOrig = function (param, yMat, xMat, llAlgorithm,
     llNGHK, llIntGrad, llOneSidedGrad, llEps, llRandom.seed,
     ...)
     {
     logLikVal <- mvProbitLogLikInternal(yMat = yMat, xMat = xMat,
     coef = param, sigma = NULL, algorithm = llAlgorithm,
     nGHK = llNGHK, returnGrad = llIntGrad, oneSidedGrad = llOneSidedGrad,
     eps = llEps, randomSeed = llRandom.seed, ...)
     return(logLikVal)
     }, gradOrig = ..2, hessOrig = ..3, yMat = c(1, 1, 1, 1, 1, 1,
     1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1,
     1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0,
     1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0,
     0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0,
     0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1,
     0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1,
     0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0), xMat = c(1,
     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1,
     0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1,
     1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0.253318513994755,
     -0.028546755348703, -0.0428704572913161, 1.36860228401446, -0.225770985659268,
     1.51647060442954, -1.54875280423022, 0.584613749636069, 0.123854243844614,
     0.215941568743973, 0.379639482759882, -0.502323453109302, -0.33320738366942,
     -1.01857538310709, -1.07179122647558, 0.303528641404258, 0.448209778629426,
     0.0530042267305041, 0.922267467879738, 2.05008468562714, -0.491031166056535,
     -2.30916887564081, 1.00573852446226, -0.709200762582393, -0.688008616467358,
     1.0255713696967, -0.284773007051009, -1.22071771225454, 0.18130347974915,
     -0.138891362439045, 0.00576418589988693, 0.38528040112633, -0.370660031792409,
     0.644376548518833, -0.220486561818751, 0.331781963915697, 1.09683901314935,
     0.435181490833803, -0.325931585531227, 1.14880761845109, 0.993503855962119,
     0.54839695950807, 0.238731735111441, -0.627906076039371, 1.36065244853001,
     -0.600259587147127, 2.18733299301658, 1.53261062618519, -0.235700359100477,
     -1.02642090030678), llAlgorithm = "GHK", llNGHK = 1000, llIntGrad = TRUE,
     llOneSidedGrad = FALSE, llEps = 1e-06, llRandom.seed = 123)
     where 5: eval(f, sys.frame(sys.parent()))
     where 6: eval(f, sys.frame(sys.parent()))
     where 7: callWithoutArgs(theta, fName = fName, args = names(formals(sumt)),
     ...)
     where 8: (function (theta, fName, ...)
     {
     return(callWithoutArgs(theta, fName = fName, args = names(formals(sumt)),
     ...))
     })(theta = c(b_1_0 = 0.8, b_1_1 = 1.2, b_1_2 = -0.8, b_2_0 = -0.6,
     b_2_1 = 1, b_2_2 = -1.6, b_3_0 = 0.5, b_3_1 = -0.6, b_3_2 = 1.2,
     R_1_2 = 0.2, R_1_3 = 0.4, R_2_3 = -0.1), fName = "logLikFunc",
     fnOrig = function (param, yMat, xMat, llAlgorithm, llNGHK,
     llIntGrad, llOneSidedGrad, llEps, llRandom.seed, ...)
     {
     logLikVal <- mvProbitLogLikInternal(yMat = yMat, xMat = xMat,
     coef = param, sigma = NULL, algorithm = llAlgorithm,
     nGHK = llNGHK, returnGrad = llIntGrad, oneSidedGrad = llOneSidedGrad,
     eps = llEps, randomSeed = llRandom.seed, ...)
     return(logLikVal)
     }, gradOrig = NULL, hessOrig = NULL, yMat = c(1, 1, 1, 1,
     1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0,
     1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0,
     0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1,
     1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1,
     0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1,
     0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0,
     0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1,
     1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0), xMat = c(1, 1, 1,
     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0,
     1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0,
     1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0,
     1, 0, 0.253318513994755, -0.028546755348703, -0.0428704572913161,
     1.36860228401446, -0.225770985659268, 1.51647060442954, -1.54875280423022,
     0.584613749636069, 0.123854243844614, 0.215941568743973,
     0.379639482759882, -0.502323453109302, -0.33320738366942,
     -1.01857538310709, -1.07179122647558, 0.303528641404258,
     0.448209778629426, 0.0530042267305041, 0.922267467879738,
     2.05008468562714, -0.491031166056535, -2.30916887564081,
     1.00573852446226, -0.709200762582393, -0.688008616467358,
     1.0255713696967, -0.284773007051009, -1.22071771225454, 0.18130347974915,
     -0.138891362439045, 0.00576418589988693, 0.38528040112633,
     -0.370660031792409, 0.644376548518833, -0.220486561818751,
     0.331781963915697, 1.09683901314935, 0.435181490833803, -0.325931585531227,
     1.14880761845109, 0.993503855962119, 0.54839695950807, 0.238731735111441,
     -0.627906076039371, 1.36065244853001, -0.600259587147127,
     2.18733299301658, 1.53261062618519, -0.235700359100477, -1.02642090030678
     ), llAlgorithm = "GHK", llNGHK = 1000, llIntGrad = TRUE,
     llOneSidedGrad = FALSE, llEps = 1e-06, llRandom.seed = 123)
     where 9: do.call(callWithoutSumt, argList)
     where 10: maxOptim(fn = fn, grad = grad, hess = hess, start = start, method = "BFGS",
     fixed = fixed, constraints = constraints, finalHessian = finalHessian,
     parscale = parscale, control = mControl, ...)
     where 11: maxRoutine(fn = logLik, grad = grad, hess = hess, start = start,
     constraints = constraints, ...)
     where 12: maxLik(logLik = logLik, start = start, method = method, finalHessian = finalHessian,
     yMat = yMat, xMat = xMat, llAlgorithm = algorithm, llNGHK = nGHK,
     llIntGrad = intGrad, llOneSidedGrad = oneSidedGrad, llEps = eps,
     llRandom.seed = random.seed, ...)
     where 13: mvProbit(cbind(y1, y2, y3) ~ x1 + x2, start = c(beta), startSigma = sigma,
     data = dat, method = "BFGS", reltol = 0.5)
    
     --- value of length: 2 type: logical ---
     [1] FALSE FALSE
     --- function from context ---
     function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
     nGHK = NULL, ...)
     {
     if (!is.matrix(sigma)) {
     stop("argument 'sigma' must be a matrix")
     }
     else if (nrow(sigma) != ncol(sigma)) {
     stop("argument 'sigma' must be a square matrix")
     }
     else if (!all(is.numeric(sigma))) {
     stop("argument 'sigma must be a numeric matrix")
     }
     if (!all(is.numeric(lower))) {
     stop("argument 'lower' must be numeric")
     }
     else if (length(lower) == 1) {
     lower <- rep(lower, nrow(sigma))
     }
     else if (length(lower) != nrow(sigma)) {
     stop("argument 'lower' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     if (!all(is.numeric(upper))) {
     stop("argument 'upper' must be numeric")
     }
     else if (length(upper) == 1) {
     upper <- rep(upper, nrow(sigma))
     }
     else if (length(upper) != nrow(sigma)) {
     stop("argument 'upper' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     algOkay <- TRUE
     ghk <- FALSE
     if (!is.list(algorithm) && length(algorithm) != 1) {
     stop("argument 'algorithm' must be a single function",
     " or a single character string")
     }
     else if (is.function(algorithm)) {
     algResult <- do.call(algorithm, list())
     if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
     "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (is.character(algorithm)) {
     if (tolower(algorithm) == "ghk") {
     ghk <- TRUE
     }
     else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     if (!algOkay) {
     stop("argument 'algorithm' must be either one of the functions",
     " 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
     " \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
     }
     if (length(random.seed) != 1) {
     stop("argument 'random.seed' must be a single numerical value")
     }
     else if (!is.numeric(random.seed)) {
     stop("argument 'random.seed' must be numerical")
     }
     if (exists(".Random.seed")) {
     savedSeed <- .Random.seed
     }
     set.seed(random.seed)
     if (exists("savedSeed")) {
     on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
     }
     else {
     on.exit(rm(.Random.seed, envir = sys.frame()))
     }
     if (ghk) {
     if (is.null(nGHK)) {
     stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
     }
     else if (length(nGHK) != 1) {
     stop("argument 'nGHK' must be a single integer value")
     }
     else if (!is.numeric(nGHK)) {
     stop("argument 'nGHK' must be numeric")
     }
     else if (nGHK <= 0) {
     stop("argument 'nGHK' must be positive")
     }
     L <- try(t(chol(sigma)), silent = TRUE)
     if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
     }
     trunpt <- rep(NA, length(lower))
     above <- rep(NA, length(lower))
     for (i in 1:length(lower)) {
     if (lower[i] == -Inf) {
     trunpt[i] <- upper[i]
     above[i] <- 1
     }
     else if (upper[i] == Inf) {
     trunpt[i] <- lower[i]
     above[i] <- 0
     }
     else {
     stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
     " or the upper truncation point must be 'Inf'")
     }
     }
     sink(tempfile())
     on.exit(sink(), add = TRUE)
     result <- ghkvec(L = L, trunpt = trunpt, above = above,
     r = nGHK)
     result <- drop(result)
     }
     else {
     result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
     algorithm = algorithm, ...)
     }
     return(result)
     }
     <bytecode: 0x2797898>
     <environment: namespace:mvProbit>
     --- function search by body ---
     Function pmvnormWrap in namespace mvProbit has this body.
     ----------- END OF FAILURE REPORT --------------
     Error in if (class(L) == "try-error") { : the condition has length > 1
     Calls: mvProbit ... logLikFunc -> fnOrig -> mvProbitLogLikInternal -> pmvnormWrap
     Execution halted
    Running the tests in ‘tests/mvProbitMargEff2.R’ failed.
    Complete output:
     > # I thank Mohit Batham for providing this R script that demonstrated
     > # a bug in mvProbitMargEff() when called with 2 dependent variables
     > library( "mvProbit" )
     Loading required package: mvtnorm
     Loading required package: maxLik
     Loading required package: miscTools
    
     Please cite the 'maxLik' package as:
     Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
    
     If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
     https://r-forge.r-project.org/projects/maxlik/
     Loading required package: abind
     > nObs <- 100
     > set.seed( 123 )
     > xData <- data.frame(
     + const = rep( 1, nObs ),
     + x1 = as.numeric( rnorm( nObs ) > 0 ),
     + x2 = as.numeric( rnorm( nObs ) > 0 ),
     + x3 = rnorm( nObs ),
     + x4 = rnorm( nObs ))
     > beta <- c( 0.8, 1.2, -1.0, 1.4, -0.8,
     + -0.6, 1.0, 0.6, -1.2, -1.6)
     > sigma <- symMatrix( c( 1, 0.2, 1))
     > margEffCond <- try( mvProbitMargEff( ~ x1 + x2 + x3 + x4 , coef = beta,
     + sigma = sigma, data = xData, cond = TRUE ) )
     ----------- FAILURE REPORT --------------
     --- failure: the condition has length > 1 ---
     --- srcref ---
     :
     --- package (from environment) ---
     mvProbit
     --- call from context ---
     pmvnormWrap(upper = xBeta[i, ], sigma = coef$sigma, algorithm = algorithm,
     nGHK = nGHK, random.seed = random.seed, ...)
     --- call from argument ---
     if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
     }
     --- R stacktrace ---
     where 1: pmvnormWrap(upper = xBeta[i, ], sigma = coef$sigma, algorithm = algorithm,
     nGHK = nGHK, random.seed = random.seed, ...)
     where 2: mvProbitExpInternal(yMat = yMat, xMat = xMat, coef = betaLower,
     sigma = coef$sigma, cond = cond, algorithm = algorithm, nGHK = nGHK,
     random.seed = random.seed, ...)
     where 3: mvProbitMargEffInternal(yMat = yMat, xMat = xMat, coef = coef,
     sigma = sigma, cond = cond, algorithm = algorithm, nGHK = nGHK,
     eps = eps, dummyVars = dummyVars, random.seed = random.seed,
     ...)
     where 4: mvProbitMargEff(~x1 + x2 + x3 + x4, coef = beta, sigma = sigma,
     data = xData, cond = TRUE)
     where 5: doTryCatch(return(expr), name, parentenv, handler)
     where 6: tryCatchOne(expr, names, parentenv, handlers[[1L]])
     where 7: tryCatchList(expr, classes, parentenv, handlers)
     where 8: 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 9: try(mvProbitMargEff(~x1 + x2 + x3 + x4, coef = beta, sigma = sigma,
     data = xData, cond = TRUE))
    
     --- value of length: 2 type: logical ---
     [1] FALSE FALSE
     --- function from context ---
     function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
     nGHK = NULL, ...)
     {
     if (!is.matrix(sigma)) {
     stop("argument 'sigma' must be a matrix")
     }
     else if (nrow(sigma) != ncol(sigma)) {
     stop("argument 'sigma' must be a square matrix")
     }
     else if (!all(is.numeric(sigma))) {
     stop("argument 'sigma must be a numeric matrix")
     }
     if (!all(is.numeric(lower))) {
     stop("argument 'lower' must be numeric")
     }
     else if (length(lower) == 1) {
     lower <- rep(lower, nrow(sigma))
     }
     else if (length(lower) != nrow(sigma)) {
     stop("argument 'lower' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     if (!all(is.numeric(upper))) {
     stop("argument 'upper' must be numeric")
     }
     else if (length(upper) == 1) {
     upper <- rep(upper, nrow(sigma))
     }
     else if (length(upper) != nrow(sigma)) {
     stop("argument 'upper' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     algOkay <- TRUE
     ghk <- FALSE
     if (!is.list(algorithm) && length(algorithm) != 1) {
     stop("argument 'algorithm' must be a single function",
     " or a single character string")
     }
     else if (is.function(algorithm)) {
     algResult <- do.call(algorithm, list())
     if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
     "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (is.character(algorithm)) {
     if (tolower(algorithm) == "ghk") {
     ghk <- TRUE
     }
     else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     if (!algOkay) {
     stop("argument 'algorithm' must be either one of the functions",
     " 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
     " \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
     }
     if (length(random.seed) != 1) {
     stop("argument 'random.seed' must be a single numerical value")
     }
     else if (!is.numeric(random.seed)) {
     stop("argument 'random.seed' must be numerical")
     }
     if (exists(".Random.seed")) {
     savedSeed <- .Random.seed
     }
     set.seed(random.seed)
     if (exists("savedSeed")) {
     on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
     }
     else {
     on.exit(rm(.Random.seed, envir = sys.frame()))
     }
     if (ghk) {
     if (is.null(nGHK)) {
     stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
     }
     else if (length(nGHK) != 1) {
     stop("argument 'nGHK' must be a single integer value")
     }
     else if (!is.numeric(nGHK)) {
     stop("argument 'nGHK' must be numeric")
     }
     else if (nGHK <= 0) {
     stop("argument 'nGHK' must be positive")
     }
     L <- try(t(chol(sigma)), silent = TRUE)
     if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
     }
     trunpt <- rep(NA, length(lower))
     above <- rep(NA, length(lower))
     for (i in 1:length(lower)) {
     if (lower[i] == -Inf) {
     trunpt[i] <- upper[i]
     above[i] <- 1
     }
     else if (upper[i] == Inf) {
     trunpt[i] <- lower[i]
     above[i] <- 0
     }
     else {
     stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
     " or the upper truncation point must be 'Inf'")
     }
     }
     sink(tempfile())
     on.exit(sink(), add = TRUE)
     result <- ghkvec(L = L, trunpt = trunpt, above = above,
     r = nGHK)
     result <- drop(result)
     }
     else {
     result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
     algorithm = algorithm, ...)
     }
     return(result)
     }
     <bytecode: 0x37a32c0>
     <environment: namespace:mvProbit>
     --- function search by body ---
     Function pmvnormWrap in namespace mvProbit has this body.
     ----------- END OF FAILURE REPORT --------------
     Error in if (class(L) == "try-error") { : the condition has length > 1
     > round( margEffCond, 3 )
     Error in round(margEffCond, 3) :
     non-numeric argument to mathematical function
     Execution halted
    Running the tests in ‘tests/mvProbitTest.R’ failed.
    Complete output:
     > library( "mvProbit" )
     Loading required package: mvtnorm
     Loading required package: maxLik
     Loading required package: miscTools
    
     Please cite the 'maxLik' package as:
     Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
    
     If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
     https://r-forge.r-project.org/projects/maxlik/
     Loading required package: abind
     > options( digits = 4 )
     >
     > ## generate a simulated data set
     > set.seed( 123 )
     > # number of observations
     > nObs <- 10
     >
     > # generate explanatory variables
     > xMat <- cbind(
     + const = rep( 1, nObs ),
     + x1 = as.numeric( rnorm( nObs ) > 0 ),
     + x2 = as.numeric( rnorm( nObs ) > 0 ),
     + x3 = rnorm( nObs ),
     + x4 = rnorm( nObs ) )
     >
     > # model coefficients
     > beta <- cbind( c( 0.8, 1.2, -1.0, 1.4, -0.8 ),
     + c( -0.6, 1.0, 0.6, -1.2, -1.6 ),
     + c( 0.5, -0.6, -0.7, 1.1, 1.2 ) )
     >
     > # covariance matrix of error terms
     > sigma <- miscTools::symMatrix( c( 1, 0.2, 0.4, 1, -0.1, 1 ) )
     >
     > # all parameters in a vector
     > allCoef <- c( c( beta ), sigma[ lower.tri( sigma ) ] )
     >
     > # generate dependent variables
     > yMatLin <- xMat %*% beta
     > yMat <- ( yMatLin + rmvnorm( nObs, sigma = sigma, pre0.9_9994 = TRUE ) ) > 0
     > colnames( yMat ) <- paste( "y", 1:3, sep = "" )
     > # (yMatLin > 0 )== yMat
     >
     > # unconditional expectations of dependent variables
     > yExp <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
     + sigma = sigma, data = as.data.frame( xMat ) )
     > round( yExp, 3 )
     V1 V2 V3
     1 0.021 0.725 0.194
     2 0.394 0.768 0.214
     3 0.125 0.788 0.196
     4 0.235 0.681 0.292
     5 0.680 0.435 0.579
     6 0.028 0.973 0.034
     7 0.958 0.186 0.784
     8 0.856 0.247 0.724
     9 0.061 0.968 0.034
     10 0.998 0.067 0.923
     > yExpA <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = allCoef,
     + data = as.data.frame( xMat ) )
     > all.equal( yExp, yExpA )
     [1] TRUE
     > yExp2 <- pnorm( yMatLin )
     > all.equal( yExp, as.data.frame( yExp2 ) )
     [1] TRUE
     >
     >
     > # conditional expectations of dependent variables
     > # (assuming that all other dependent variables are one)
     > yExpCond <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
     + sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
     + algorithm = GenzBretz() )
     > round( yExpCond, 3 )
     V1 V2 V3
     1 0.072 0.834 0.524
     2 0.660 0.779 0.314
     3 0.297 0.838 0.399
     4 0.448 0.727 0.461
     5 0.847 0.442 0.618
     6 0.139 0.985 0.162
     7 0.991 0.179 0.749
     8 0.950 0.244 0.710
     9 0.244 0.978 0.133
     10 1.000 0.065 0.892
     > yExpCondA <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = allCoef,
     + data = as.data.frame( xMat ), cond = TRUE,
     + algorithm = GenzBretz() )
     > all.equal( yExpCond, yExpCondA )
     [1] TRUE
     > yExpCond2 <- matrix( NA, nrow = nObs, ncol = ncol( yMat ) )
     > for( i in 1:nObs ) {
     + for( k in 1:ncol( yMat ) ) {
     + set.seed( 123 )
     + numerator <- pmvnorm( upper = yMatLin[ i, ], sigma = sigma )
     + set.seed( 123 )
     + denominator <- pmvnorm( upper = yMatLin[ i, -k ], sigma = sigma[ -k, -k ] )
     + yExpCond2[ i, k ] <- numerator / denominator
     + }
     + }
     > all.equal( yExpCond, as.data.frame( yExpCond2 ) )
     [1] TRUE
     > # now with explicitly specifying the algorithm
     > yExpCond3 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
     + sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
     + algorithm = GenzBretz )
     > all.equal( yExpCond, yExpCond3 )
     [1] TRUE
     > identical( yExpCond, yExpCond3 )
     [1] TRUE
     > # now with integrals obtained by the Miwa algorithm
     > yExpCond4 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
     + sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
     + algorithm = Miwa )
     > all.equal( yExpCond, yExpCond4, tol = 1e-3 )
     [1] TRUE
     > # now with integrals obtained by the Miwa algorithm, less precise
     > yExpCond5 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
     + sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
     + algorithm = Miwa( steps = 32 ) )
     > all.equal( yExpCond4, yExpCond5, tol = 1e-3 )
     [1] TRUE
     > # now with integrals obtained by the TVPACK algorithm
     > yExpCond6 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
     + sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
     + algorithm = TVPACK )
     > all.equal( yExpCond, yExpCond6, tol = 1e-3 )
     [1] TRUE
     > # now with integrals obtained by the TVPACK algorithm, less precise
     > yExpCond7 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
     + sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
     + algorithm = TVPACK( abseps = 0.5 ) )
     > all.equal( yExpCond6, yExpCond7, tol = 1e-3 )
     [1] "Component \"V1\": Mean relative difference: 0.03018"
     [2] "Component \"V2\": Mean relative difference: 0.04273"
     [3] "Component \"V3\": Mean relative difference: 0.04242"
     > # now with integrals obtained by the GHK algorithm
     > yExpCond8 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
     + sigma = sigma, data = as.data.frame( xMat ), cond = TRUE )
     ----------- FAILURE REPORT --------------
     --- failure: the condition has length > 1 ---
     --- srcref ---
     :
     --- package (from environment) ---
     mvProbit
     --- call from context ---
     pmvnormWrap(upper = xBeta[i, ], sigma = coef$sigma, algorithm = algorithm,
     nGHK = nGHK, random.seed = random.seed, ...)
     --- call from argument ---
     if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
     }
     --- R stacktrace ---
     where 1: pmvnormWrap(upper = xBeta[i, ], sigma = coef$sigma, algorithm = algorithm,
     nGHK = nGHK, random.seed = random.seed, ...)
     where 2: mvProbitExpInternal(yMat = yMat, xMat = xMat, coef = coef, sigma = sigma,
     cond = cond, algorithm = algorithm, nGHK = nGHK, random.seed = random.seed,
     ...)
     where 3: mvProbitExp(~x1 + x2 + x3 + x4, coef = c(beta), sigma = sigma,
     data = as.data.frame(xMat), cond = TRUE)
    
     --- value of length: 2 type: logical ---
     [1] FALSE FALSE
     --- function from context ---
     function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
     nGHK = NULL, ...)
     {
     if (!is.matrix(sigma)) {
     stop("argument 'sigma' must be a matrix")
     }
     else if (nrow(sigma) != ncol(sigma)) {
     stop("argument 'sigma' must be a square matrix")
     }
     else if (!all(is.numeric(sigma))) {
     stop("argument 'sigma must be a numeric matrix")
     }
     if (!all(is.numeric(lower))) {
     stop("argument 'lower' must be numeric")
     }
     else if (length(lower) == 1) {
     lower <- rep(lower, nrow(sigma))
     }
     else if (length(lower) != nrow(sigma)) {
     stop("argument 'lower' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     if (!all(is.numeric(upper))) {
     stop("argument 'upper' must be numeric")
     }
     else if (length(upper) == 1) {
     upper <- rep(upper, nrow(sigma))
     }
     else if (length(upper) != nrow(sigma)) {
     stop("argument 'upper' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     algOkay <- TRUE
     ghk <- FALSE
     if (!is.list(algorithm) && length(algorithm) != 1) {
     stop("argument 'algorithm' must be a single function",
     " or a single character string")
     }
     else if (is.function(algorithm)) {
     algResult <- do.call(algorithm, list())
     if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
     "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (is.character(algorithm)) {
     if (tolower(algorithm) == "ghk") {
     ghk <- TRUE
     }
     else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     if (!algOkay) {
     stop("argument 'algorithm' must be either one of the functions",
     " 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
     " \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
     }
     if (length(random.seed) != 1) {
     stop("argument 'random.seed' must be a single numerical value")
     }
     else if (!is.numeric(random.seed)) {
     stop("argument 'random.seed' must be numerical")
     }
     if (exists(".Random.seed")) {
     savedSeed <- .Random.seed
     }
     set.seed(random.seed)
     if (exists("savedSeed")) {
     on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
     }
     else {
     on.exit(rm(.Random.seed, envir = sys.frame()))
     }
     if (ghk) {
     if (is.null(nGHK)) {
     stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
     }
     else if (length(nGHK) != 1) {
     stop("argument 'nGHK' must be a single integer value")
     }
     else if (!is.numeric(nGHK)) {
     stop("argument 'nGHK' must be numeric")
     }
     else if (nGHK <= 0) {
     stop("argument 'nGHK' must be positive")
     }
     L <- try(t(chol(sigma)), silent = TRUE)
     if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
     }
     trunpt <- rep(NA, length(lower))
     above <- rep(NA, length(lower))
     for (i in 1:length(lower)) {
     if (lower[i] == -Inf) {
     trunpt[i] <- upper[i]
     above[i] <- 1
     }
     else if (upper[i] == Inf) {
     trunpt[i] <- lower[i]
     above[i] <- 0
     }
     else {
     stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
     " or the upper truncation point must be 'Inf'")
     }
     }
     sink(tempfile())
     on.exit(sink(), add = TRUE)
     result <- ghkvec(L = L, trunpt = trunpt, above = above,
     r = nGHK)
     result <- drop(result)
     }
     else {
     result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
     algorithm = algorithm, ...)
     }
     return(result)
     }
     <bytecode: 0x20ade30>
     <environment: namespace:mvProbit>
     --- function search by body ---
     Function pmvnormWrap in namespace mvProbit has this body.
     ----------- END OF FAILURE REPORT --------------
     Error in if (class(L) == "try-error") { : the condition has length > 1
     Calls: mvProbitExp -> mvProbitExpInternal -> pmvnormWrap
     Execution halted
    Running the tests in ‘tests/pmvnormWrapTest.R’ failed.
    Complete output:
     > library( "mvProbit" )
     Loading required package: mvtnorm
     Loading required package: maxLik
     Loading required package: miscTools
    
     Please cite the 'maxLik' package as:
     Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
    
     If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
     https://r-forge.r-project.org/projects/maxlik/
     Loading required package: abind
     > library( "miscTools" )
     >
     > # covariance matrix
     > sigma <- miscTools::symMatrix( c( 1, 0.2, 0.4, 1, -0.1, 1 ) )
     >
     > ######## only upper ##########
     > upper <- c( -0.3, 0.7, -0.5 )
     > # Genz + Bretz (default)
     > pug <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
     + algorithm = GenzBretz(), random.seed = 123 )
     > print( pug )
     [1] 0.1360783
     attr(,"error")
     [1] 0.000118605
     attr(,"msg")
     [1] "Normal Completion"
     >
     > # Miwa (as function)
     > pum <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
     + algorithm = Miwa, random.seed = 123 )
     > print( pum )
     [1] 0.136069
     attr(,"error")
     [1] NA
     attr(,"msg")
     [1] "Normal Completion"
     > all.equal( pug, pum, check.attributes = FALSE, tol = 1e-4 )
     [1] TRUE
     >
     > # Miwa (as object returned from function Miwa())
     > pum1 <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
     + algorithm = Miwa(), random.seed = 123 )
     > all.equal( pum, pum1 )
     [1] TRUE
     >
     > # Miwa (as character string)
     > pum2 <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
     + algorithm = "Miwa", random.seed = 123 )
     > all.equal( pum, pum2 )
     [1] TRUE
     >
     > # TVPACK
     > put <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
     + algorithm = TVPACK, random.seed = 123 )
     > print( put )
     [1] 0.136069
     attr(,"error")
     [1] 1e-06
     attr(,"msg")
     [1] "Normal Completion"
     > all.equal( pug, put, check.attributes = FALSE, tol = 1e-4 )
     [1] TRUE
     > all.equal( pum, put, check.attributes = FALSE, tol = 1e-6 )
     [1] TRUE
     >
     > # GHK
     > pughk <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
     + algorithm = "ghk", random.seed = 123, nGHK = 1000 )
     ----------- FAILURE REPORT --------------
     --- failure: the condition has length > 1 ---
     --- srcref ---
     :
     --- package (from environment) ---
     mvProbit
     --- call from context ---
     mvProbit:::pmvnormWrap(upper = upper, sigma = sigma, algorithm = "ghk",
     random.seed = 123, nGHK = 1000)
     --- call from argument ---
     if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
     }
     --- R stacktrace ---
     where 1: mvProbit:::pmvnormWrap(upper = upper, sigma = sigma, algorithm = "ghk",
     random.seed = 123, nGHK = 1000)
    
     --- value of length: 2 type: logical ---
     [1] FALSE FALSE
     --- function from context ---
     function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
     nGHK = NULL, ...)
     {
     if (!is.matrix(sigma)) {
     stop("argument 'sigma' must be a matrix")
     }
     else if (nrow(sigma) != ncol(sigma)) {
     stop("argument 'sigma' must be a square matrix")
     }
     else if (!all(is.numeric(sigma))) {
     stop("argument 'sigma must be a numeric matrix")
     }
     if (!all(is.numeric(lower))) {
     stop("argument 'lower' must be numeric")
     }
     else if (length(lower) == 1) {
     lower <- rep(lower, nrow(sigma))
     }
     else if (length(lower) != nrow(sigma)) {
     stop("argument 'lower' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     if (!all(is.numeric(upper))) {
     stop("argument 'upper' must be numeric")
     }
     else if (length(upper) == 1) {
     upper <- rep(upper, nrow(sigma))
     }
     else if (length(upper) != nrow(sigma)) {
     stop("argument 'upper' must either be a single numeric value",
     " or a vector with length equal to the number of rows/columns",
     " of argument 'sigma'")
     }
     algOkay <- TRUE
     ghk <- FALSE
     if (!is.list(algorithm) && length(algorithm) != 1) {
     stop("argument 'algorithm' must be a single function",
     " or a single character string")
     }
     else if (is.function(algorithm)) {
     algResult <- do.call(algorithm, list())
     if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
     "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (is.character(algorithm)) {
     if (tolower(algorithm) == "ghk") {
     ghk <- TRUE
     }
     else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     }
     else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
     algOkay <- FALSE
     }
     if (!algOkay) {
     stop("argument 'algorithm' must be either one of the functions",
     " 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
     " \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
     }
     if (length(random.seed) != 1) {
     stop("argument 'random.seed' must be a single numerical value")
     }
     else if (!is.numeric(random.seed)) {
     stop("argument 'random.seed' must be numerical")
     }
     if (exists(".Random.seed")) {
     savedSeed <- .Random.seed
     }
     set.seed(random.seed)
     if (exists("savedSeed")) {
     on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
     }
     else {
     on.exit(rm(.Random.seed, envir = sys.frame()))
     }
     if (ghk) {
     if (is.null(nGHK)) {
     stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
     }
     else if (length(nGHK) != 1) {
     stop("argument 'nGHK' must be a single integer value")
     }
     else if (!is.numeric(nGHK)) {
     stop("argument 'nGHK' must be numeric")
     }
     else if (nGHK <= 0) {
     stop("argument 'nGHK' must be positive")
     }
     L <- try(t(chol(sigma)), silent = TRUE)
     if (class(L) == "try-error") {
     warning("the correlation matrix is not positive definite")
     return(NA)
     }
     trunpt <- rep(NA, length(lower))
     above <- rep(NA, length(lower))
     for (i in 1:length(lower)) {
     if (lower[i] == -Inf) {
     trunpt[i] <- upper[i]
     above[i] <- 1
     }
     else if (upper[i] == Inf) {
     trunpt[i] <- lower[i]
     above[i] <- 0
     }
     else {
     stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
     " or the upper truncation point must be 'Inf'")
     }
     }
     sink(tempfile())
     on.exit(sink(), add = TRUE)
     result <- ghkvec(L = L, trunpt = trunpt, above = above,
     r = nGHK)
     result <- drop(result)
     }
     else {
     result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
     algorithm = algorithm, ...)
     }
     return(result)
     }
     <bytecode: 0x39b0a58>
     <environment: namespace:mvProbit>
     --- function search by body ---
     Function pmvnormWrap in namespace mvProbit has this body.
     ----------- END OF FAILURE REPORT --------------
     Error in if (class(L) == "try-error") { : the condition has length > 1
     Calls: <Anonymous>
     Execution halted
Flavor: r-devel-linux-x86_64-fedora-gcc