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

[Ravi Varadhan]

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.