`tetra1` <-
function(mu,r12,crit=1e-6, maxiter=20)
{
if ( missing(mu) ) {stop('You need to provide the mean vector')}
if (r12==0)
{
L <- c(0,0,0)
return(L)
# stop("The pairwise correlation 'r12' need to be > 0")
}
# find [a,b] that brackets the root
if (r12>0){a <- 0; b <- 1}
else {a <- -1; b <- 0}
u12 <- mu[1]*mu[2] + r12*sqrt( mu[1]*(1-mu[1]) * mu[2]*(1-mu[2]) )
x1 <- qnorm(mu[1])
x2 <- qnorm(mu[2])
tcc <- r12
#============================================================================#
# Repeated bisection to locate the root in [a,b].
#============================================================================#
niter <- 0
cdf <- pmvnorm( lower=c(-Inf, -Inf), upper=c(x1,x2), corr=toeplitz(c(1,tcc)) )
error <- cdf - u12
fail <- ( abs(error) >= crit )
while(fail && (niter < maxiter))
{
if (error > 0) {b <- tcc}
else {a <- tcc}
tcc <- a + 0.5*(b-a)
cdf <- pmvnorm( lower=c(-Inf, -Inf), upper=c(x1,x2),corr=toeplitz(c(1,tcc)) )
error <- cdf - u12
fail <- (abs(error) >= crit)
niter <- niter + 1;
}
L <- c(tcc, fail, niter)
return(L)
}