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:
[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