Hi all, I defined the function integrand as follows. library(mnormt) integrand <- function (rho, a, b,z, x, y){ Sigma <- matrix(c(1, rho, rho, 1), 2,2) mu <- rep(0,2) f <- pnorm(c((z-a*x)/b, (z-a*y)/b), mu, Sigma)*dnorm(c(0,0), mu, diag(2)) f } Then I would like to do the following: 1. Take the integral of the function integrand with respect to x and y. The integration limit is (-Inf, 1) for both x and y. 2. This integral is then a function of (rho, a, b, z). I then want to find the value of z that would make this integral to be equal to a constant, say, 0.1, for specific rho, a and b. Can anyone give some help on this? Thank you very much. Hanna [[alternative HTML version deleted]]