This is not a problem of numDerv. The problem is that the function PP1 is stochastic, i.e. it gives different values for the same argument p. This is due to the function pmvnorm(), which I presume uses some kind of Monte Carlo sampling to compute the integral of a multivariate normal. Therefore, you cannot evaluate the derivative of a stochastic function.
Ravi. ____________________________________________________________________ Ravi Varadhan, Ph.D. Assistant Professor, Division of Geriatric Medicine and Gerontology School of Medicine Johns Hopkins University Ph. (410) 502-2619 email: rvarad...@jhmi.edu ----- Original Message ----- From: Stephane LUCHINI <stephane.luch...@univmed.fr> Date: Friday, November 20, 2009 6:56 am Subject: [R] Problem with Numerical derivatives (numDeriv) and mvtnorm To: r-help@r-project.org > I'm trying to obtain numerical derivative of a probability computed > with mvtnorm with respect to its parameters using grad() and > jacobian() from NumDeriv. > > To simplify the matter, here is an example: > > PP1 <- function(p){ > thetac <- p > thetae <- 0.323340333 > thetab <- -0.280970036 > thetao <- 0.770768082 > ssigma <- diag(4) > ssigma[1,2] <- 0.229502120 > ssigma[1,3] <- 0.677949335 > ssigma[1,4] <- 0.552907745 > ssigma[2,3] <- 0.784263100 > ssigma[2,4] <- 0.374065025 > ssigma[3,4] <- 0.799238700 > ssigma[2,1] <- ssigma[1,2] > ssigma[3,1] <- ssigma[1,3] > ssigma[4,1] <- ssigma[1,4] > ssigma[3,2] <- ssigma[2,3] > ssigma[4,2] <- ssigma[2,4] > ssigma[4,3] <- ssigma[3,4] > pp1 <- > > pmvnorm(lower=c(-Inf,-Inf,-Inf,-Inf),upper=c(thetac,thetae,thetab,thetao),corr=ssigma) > return(pp1)} > > xx <- -0.6675762 > > If I compute several times the probability PP1, i obtain some > slightly > different numbers but I suppose this is ok: > > > PP1(xx) > [1] 0.1697787 > attr(,"error") > [1] 0.000840748 > attr(,"msg") > [1] "Normal Completion" > > PP1(xx) > [1] 0.1697337 > attr(,"error") > [1] 0.0009363715 > attr(,"msg") > [1] "Normal Completion" > > PP1(xx) > [1] 0.1691539 > attr(,"error") > [1] 0.0006682577 > attr(,"msg") > [1] "Normal Completion" > > When I now turn to the numerical derivatives of the probability, I > obtain large discrepencies if I repeat the instruction "grad": > > > grad(PP1,xx) > [1] -42.58016 > > grad(PP1,xx) > [1] 9.297055 > > grad(PP1,xx) > [1] -6.736725 > > grad(PP1,xx) > [1] -20.71176 > > grad(PP1,xx) > [1] 18.51968 > > The "problem" is the same if I turn to the jacobian function (when I > > want to compute all partial derivatives in one shot) > > Someone can help? > > Stephane > > ______________________________________________ > R-help@r-project.org mailing list > > PLEASE do read the posting guide > and provide commented, minimal, self-contained, reproducible code. ______________________________________________ R-help@r-project.org 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.