Stephane LUCHINI
2009-Nov-20 11:55 UTC
[R] Problem with Numerical derivatives (numDeriv) and mvtnorm
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
Ravi Varadhan
2009-Nov-21 14:43 UTC
[R] Problem with Numerical derivatives (numDeriv) and mvtnorm
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: rvaradhan at jhmi.edu ----- Original Message ----- From: Stephane LUCHINI <stephane.luchini at univmed.fr> Date: Friday, November 20, 2009 6:56 am Subject: [R] Problem with Numerical derivatives (numDeriv) and mvtnorm To: r-help at 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 at r-project.org mailing list > > PLEASE do read the posting guide > and provide commented, minimal, self-contained, reproducible code.
Thanks for you comment. There is certainly some "Monte Carlo sampling" involved in mvtnorm but why derivatives could not be computed? In theory, the derivatives exist (eg. bivariate probit). Moreover, when used with optim, there are some numerical derivatives computed... does it mean that mvtnorm cannot be used in an optimisation problem? I think it hard to believe. One possibility would be to use the analytical derivatives and then a do-it-yourself integration but i was looking for something a bit more comprehensive. The mvtnorm package uses a specific way to compute pmvnorm and I'm far to do a good enough job so that derivatives can compare with what mvtnorm can do. Stef