I am trying calculate a probability using numerical integration. The first program I ran spit out an answer in a very short time. The program is below: ## START PROGRAM trial <- function(input) { pmvnorm(lower = c(0,0), upper = c(2, 2), mean = input, sigma = matrix(c(.1, 0, 0, .1), nrow = 2, ncol = 2, byrow = FALSE)) } require(mvtnorm) require(adapt) bottomB <- -5*sqrt(.1) topB <- 2 + 5*sqrt(.1) areaB <- (topB - bottomB)^2 unscaled.Po.in.a <- adapt(2, lo = c(bottomB, bottomB), up = c(topB, topB), minpts = 1000, eps = 1e-4, functn = trial) (1/areaB)*unscaled.Po.in.a$value ## FINISH PROGRAM I tried to run the program again changing a.) sigma in the trial function, b.) upper in the trial function, and c.) the bounds of integration; that is, bottomB and topB. The new program is below: ## START PROGRAM trial <- function(input) { pmvnorm(lower = c(0,0), upper = c(10, 10), mean = input, sigma = matrix(c(.01, 0, 0, .01), nrow = 2, ncol = 2, byrow = FALSE)) } require(mvtnorm) require(adapt) bottomB <- -5*sqrt(.01) topB <- 10 + 5*sqrt(.01) areaB <- (topB - bottomB)^2 unscaled.Po.in.a <- adapt(2, lo = c(bottomB, bottomB), up = c(topB, topB), minpts = 1000, eps = 1e-4, functn = trial) (1/areaB)*unscaled.Po.in.a$value ## FINISH PROGRAM Now, the program just runs and runs (48 hours at last count!). By playing around with the program, I have deduced the program is highly sensitive to changing the upper option in the trial function. For example, using a vector like (4, 4) causes no problems and the program quickly yields an answer. I have a couple of other programs where I can easily obtain a simulation-based answer, but I would ultimately like to know what's up with this program before I give up on it so I can learn a thing or two. Does anyone have any clues or tricks to get around this problem? My guess is that it will simply be very difficult (impossible?) to obtain this type of relative error (eps = 1e-4) and I will have no choice but to pursue the simulation approach. Thanks for any responses (philip.turk at nau.edu)! -- Phil Philip Turk Assistant Professor of Statistics Northern Arizona University Department of Mathematics and Statistics PO Box 5717 Flagstaff, AZ 86011 Phone: 928-523-6884 Fax: 928-523-5847 E-mail: Philip.Turk at nau.edu Web Site: http://jan.ucc.nau.edu/~stapjt-p/
What versions of "adapt" and R are you using? The current package was built with R-2.5.1. I tried your program with R-2.5.0, and got the answer 0.1501053 in just a few seconds. At 03:20 PM 7/7/2007, Philip wrote:>I am trying calculate a probability using numerical integration. The first >program I ran spit out an answer in a very short time. The program is below: > >## START PROGRAM > >trial <- function(input) > >{ >pmvnorm(lower = c(0,0), upper = c(2, 2), mean = input, sigma = >matrix(c(.1, 0, >0, .1), nrow = 2, ncol = 2, byrow = FALSE)) >} > >require(mvtnorm) >require(adapt) > >bottomB <- -5*sqrt(.1) >topB <- 2 + 5*sqrt(.1) >areaB <- (topB - bottomB)^2 > >unscaled.Po.in.a <- adapt(2, lo = c(bottomB, bottomB), up = c(topB, topB), >minpts = 1000, eps = 1e-4, functn = trial) > >(1/areaB)*unscaled.Po.in.a$value > >## FINISH PROGRAM > >I tried to run the program again changing a.) sigma in the trial >function, b.) >upper in the trial function, and c.) the bounds of integration; that is, >bottomB and topB. The new program is below: > >## START PROGRAM > >trial <- function(input) > >{ >pmvnorm(lower = c(0,0), upper = c(10, 10), mean = input, sigma = >matrix(c(.01, >0, 0, .01), nrow = 2, ncol = 2, byrow = FALSE)) >} > >require(mvtnorm) >require(adapt) > >bottomB <- -5*sqrt(.01) >topB <- 10 + 5*sqrt(.01) >areaB <- (topB - bottomB)^2 > >unscaled.Po.in.a <- adapt(2, lo = c(bottomB, bottomB), up = c(topB, topB), >minpts = 1000, eps = 1e-4, functn = trial) > >(1/areaB)*unscaled.Po.in.a$value > >## FINISH PROGRAM > >Now, the program just runs and runs (48 hours at last count!). By playing >around with the program, I have deduced the program is highly sensitive to >changing the upper option in the trial function. For example, using a vector >like (4, 4) causes no problems and the program quickly yields an answer. I >have a couple of other programs where I can easily obtain a simulation-based >answer, but I would ultimately like to know what's up with this >program before >I give up on it so I can learn a thing or two. Does anyone have any clues or >tricks to get around this problem? My guess is that it will simply be very >difficult (impossible?) to obtain this type of relative error (eps = >1e-4) and >I will have no choice but to pursue the simulation approach. > >Thanks for any responses (philip.turk at nau.edu)! > >-- Phil > >Philip Turk >Assistant Professor of Statistics >Northern Arizona University >Department of Mathematics and Statistics >PO Box 5717 >Flagstaff, AZ 86011 >Phone: 928-523-6884 >Fax: 928-523-5847 >E-mail: Philip.Turk at nau.edu >Web Site: http://jan.ucc.nau.edu/~stapjt-p/ > >______________________________________________ >R-help at stat.math.ethz.ch mailing list >https://stat.ethz.ch/mailman/listinfo/r-help >PLEASE do read the posting guide http://www.R-project.org/posting-guide.html >and provide commented, minimal, self-contained, reproducible code.===============================================================Robert A. LaBudde, PhD, PAS, Dpl. ACAFS e-mail: ral at lcfltd.com Least Cost Formulations, Ltd. URL: http://lcfltd.com/ 824 Timberlake Drive Tel: 757-467-0954 Virginia Beach, VA 23464-3239 Fax: 757-467-2947 "Vere scire est per causas scire"
Since the covariance is zero (i.e. you have independent normals), you can simplify your problem so that you just need to perform one-dimensional integration. Here is how you can do it: trial2 <- function(input) { #pmvnorm(lower = c(0,0), upper = c(10, 10), mean = input, sigma matrix(c(.01, 0, 0, .01), nrow = 2, ncol = 2, byrow = FALSE)) pnorm(q=10, mean = input, sd = sqrt(.01)) - pnorm(q=0, mean = input, sd sqrt(.01)) } bottomB <- -5*sqrt(.01) topB <- 10 + 5*sqrt(.01) areaB <- (topB - bottomB)^2 new.ans <- 1/areaB * (integrate(f=trial2, lo = bottomB, up = topB)$val)^2> new.ans[1] 0.8264463 Hope this is helpful, Ravi. ---------------------------------------------------------------------------- ------- Ravi Varadhan, Ph.D. Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University Ph: (410) 502-2619 Fax: (410) 614-9625 Email: rvaradhan at jhmi.edu Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html ---------------------------------------------------------------------------- -------- -----Original Message----- From: r-help-bounces at stat.math.ethz.ch [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Philip Turk Sent: Saturday, July 07, 2007 3:20 PM To: r-help at stat.math.ethz.ch Subject: [R] No convergence using ADAPT I am trying calculate a probability using numerical integration. The first program I ran spit out an answer in a very short time. The program is below: ## START PROGRAM trial <- function(input) { pmvnorm(lower = c(0,0), upper = c(2, 2), mean = input, sigma = matrix(c(.1, 0, 0, .1), nrow = 2, ncol = 2, byrow = FALSE)) } require(mvtnorm) require(adapt) bottomB <- -5*sqrt(.1) topB <- 2 + 5*sqrt(.1) areaB <- (topB - bottomB)^2 unscaled.Po.in.a <- adapt(2, lo = c(bottomB, bottomB), up = c(topB, topB), minpts = 1000, eps = 1e-4, functn = trial) (1/areaB)*unscaled.Po.in.a$value ## FINISH PROGRAM I tried to run the program again changing a.) sigma in the trial function, b.) upper in the trial function, and c.) the bounds of integration; that is, bottomB and topB. The new program is below: ## START PROGRAM trial <- function(input) { pmvnorm(lower = c(0,0), upper = c(10, 10), mean = input, sigma matrix(c(.01, 0, 0, .01), nrow = 2, ncol = 2, byrow = FALSE)) } require(mvtnorm) require(adapt) bottomB <- -5*sqrt(.01) topB <- 10 + 5*sqrt(.01) areaB <- (topB - bottomB)^2 unscaled.Po.in.a <- adapt(2, lo = c(bottomB, bottomB), up = c(topB, topB), minpts = 1000, eps = 1e-4, functn = trial) (1/areaB)*unscaled.Po.in.a$value ## FINISH PROGRAM Now, the program just runs and runs (48 hours at last count!). By playing around with the program, I have deduced the program is highly sensitive to changing the upper option in the trial function. For example, using a vector like (4, 4) causes no problems and the program quickly yields an answer. I have a couple of other programs where I can easily obtain a simulation-based answer, but I would ultimately like to know what's up with this program before I give up on it so I can learn a thing or two. Does anyone have any clues or tricks to get around this problem? My guess is that it will simply be very difficult (impossible?) to obtain this type of relative error (eps = 1e-4) and I will have no choice but to pursue the simulation approach. Thanks for any responses (philip.turk at nau.edu)! -- Phil Philip Turk Assistant Professor of Statistics Northern Arizona University Department of Mathematics and Statistics PO Box 5717 Flagstaff, AZ 86011 Phone: 928-523-6884 Fax: 928-523-5847 E-mail: Philip.Turk at nau.edu Web Site: http://jan.ucc.nau.edu/~stapjt-p/ ______________________________________________ R-help at stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.