Re: [R] 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=10)
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
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.
Re: [R] Problem with Numerical derivatives (numDeriv) and mvtnorm
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 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.
Re: [R] Problem with Numerical derivatives (numDeriv) and mvtnorm
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: rvarad...@jhmi.edu
- Original Message -
From: Stephane LUCHINI
Date: Friday, November 20, 2009 6:56 am
Subject: [R] Problem with Numerical derivatives (numDeriv) and mvtnorm
To: r-help@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@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.
[R] 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
__
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.
