Problem with Numerical derivatives (numDeriv) and mvtnorm

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: rvaradhan at jhmi.edu


----- Original Message -----
From: SL <sl465 at yahoo.fr>
Date: Saturday, November 21, 2009 2:42 pm
Subject: Re: [R] Problem with Numerical derivatives (numDeriv) and mvtnorm
To: r-help at 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 at r-project.org mailing list

PLEASE do read the posting guide 
and provide commented, minimal, self-contained, reproducible code.