R help - Nov 2009 - Problem with Numerical derivatives (numDeriv) and mvtnorm

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

This is not a problem of numDerv.  The problem is that the function PP1 is
stochastic, i.e. it gives different values for the same argument p.  This is due
to the function pmvnorm(), which I presume uses some kind of Monte Carlo
sampling to compute the integral of a multivariate normal. Therefore, you cannot
evaluate the derivative of a stochastic function.

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: Stephane LUCHINI <stephane.luchini at univmed.fr>
Date: Friday, November 20, 2009 6:56 am
Subject: [R] Problem with Numerical derivatives (numDeriv) and mvtnorm
To: r-help at r-project.org

> 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 at r-project.org mailing list
>  
>  PLEASE do read the posting guide 
>  and provide commented, minimal, self-contained, reproducible code.

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