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
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.
