Running ‘mvProbitEst.R’ [21s/30s]
  Running ‘mvProbitMargEff2.R’ [1s/2s]
  Running ‘mvProbitTest.R’ [2s/3s]
  Running ‘mvProbitTest5.R’ [88s/119s]
  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’ [1s/2s]
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: 0x55e93e9a6860>
  <environment: namespace:mvProbit>
   --- function search by body ---
  Function pmvnormWrap in namespace mvProbit has this body.
   ----------- END OF FAILURE REPORT -------------- 
  Fatal error: the condition has length > 1
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: 0x5558b75737e8>
  <environment: namespace:mvProbit>
   --- function search by body ---
  Function pmvnormWrap in namespace mvProbit has this body.
   ----------- END OF FAILURE REPORT -------------- 
  Fatal error: the condition has length > 1
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: 0x563c50e23728>
  <environment: namespace:mvProbit>
   --- function search by body ---
  Function pmvnormWrap in namespace mvProbit has this body.
   ----------- END OF FAILURE REPORT -------------- 
  Fatal error: the condition has length > 1
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: 0x56111f7fd1f8>
  <environment: namespace:mvProbit>
   --- function search by body ---
  Function pmvnormWrap in namespace mvProbit has this body.
   ----------- END OF FAILURE REPORT -------------- 
  Fatal error: the condition has length > 1