Go back to your calculus text and review the definition of derivative:
f'(x) = lim h -> 0 [f(x+h) - f(x)] / h
when f(x) and f(x + h) are random variables, the above limit does not exist.
In fact, f'(x) is also a random variable.
Now, if you want the derivative you have to use a multivariate integration
algorithm that yields a deterministic value. The function `sadmvn' in the
package "mnormt" can do this:
require(mnormt)
PP2 <- 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]
pp <- sadmvn(lower=rep(-Inf, 4), upper=c(thetac,thetae,thetab,thetao),
mean=rep(0,4), varcov=ssigma, maxpt=100000)
return(pp)
}
xx <- -0.6675762
P2(xx)
require(numDeriv)
grad(x=xx, func=PP2)
I hope this helps,
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: SL <sl...@yahoo.fr>
Date: Saturday, November 21, 2009 2:42 pm
Subject: Re: [R] Problem with Numerical derivatives (numDeriv) and mvtnorm
To: r-help@r-project.org
> 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
>
> ______________________________________________
> 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.
