[R] No convergence using ADAPT
Ravi Varadhan
rvaradhan at jhmi.edu
Mon Jul 9 17:27:06 CEST 2007
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.
More information about the R-help
mailing list