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

[Ravi Varadhan]

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.