[R] Problem with Numerical derivatives (numDeriv) and mvtnorm

[Stephane LUCHINI]

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