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
Problem with Numerical derivatives (numDeriv) and mvtnorm
6 messages · SL, Ravi Varadhan, stephane Luchini +1 more
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.
Hi Torsten, It would be useful to "warn" the users that the multivariate normal probability calculated by "pmvnorm" using the GenzBretz algorithm is "random", i.e. the result can vary between repeated executions of the function. This would prevent inappropriate use of pmvnorm such as computing derivatives of it (see this email thread). It seems that the other algorithm "Miwa" is deterministic, but not sure how reliable it is (I had some trouble with it). It would also be useful in the help page to provide a link to two other functions for evaluating multivariate normal probabilities: mnormt::sadmvn mprobit::mvnapp In particular, the `mvnapp' function of Harry Joe in "mprobit" package seems to be very interesting as it provides very accurate results using asymptotic expansions. Best, 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: Ravi Varadhan <rvaradhan at jhmi.edu> Date: Saturday, November 21, 2009 8:15 pm Subject: Re: [R] Problem with Numerical derivatives (numDeriv) and mvtnorm To: SL <sl465 at yahoo.fr> Cc: r-help at r-project.org
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. ______________________________________________ R-help at r-project.org mailing list PLEASE do read the posting guide and provide commented, minimal, self-contained, reproducible code.
I'm now making some trials with sadmvn which provides results similar to pmvnorm for optimization but I know compute my OPG estimator of the covariance matrix with sadmvn (by the way Ravi, when I was refering to "exist in theory" I was refering to the theory not to the computation - would an appropriate "random" computation of partial derivative work?). Interestingly, mprobit also provides derivatives, exactly what I need. Unfortunatly it fails to install on mac os X! (I don't want to install windows in my system and my linux server is off for the moment). Stephane 2009/11/22 Ravi Varadhan <rvaradhan at jhmi.edu>:
Hi Torsten, It would be useful to "warn" the users that the multivariate normal probability calculated by "pmvnorm" using the GenzBretz algorithm is "random", i.e. the result can vary between repeated executions of the function. ?This would prevent inappropriate use of pmvnorm such as computing derivatives of it (see this email thread). It seems that the other algorithm "Miwa" is deterministic, but not sure how reliable it is (I had some trouble with it). It would also be useful in the help page to provide a link to two other functions for evaluating multivariate normal probabilities: mnormt::sadmvn mprobit::mvnapp In particular, the `mvnapp' function of Harry Joe in "mprobit" package seems to be very interesting as it provides very accurate results using asymptotic expansions. Best, 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: Ravi Varadhan <rvaradhan at jhmi.edu>
Date: Saturday, November 21, 2009 8:15 pm
Subject: Re: [R] Problem with Numerical derivatives (numDeriv) and mvtnorm
To: SL <sl465 at yahoo.fr>
Cc: r-help at r-project.org
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.
______________________________________________
R-help at r-project.org mailing list
PLEASE do read the posting guide
and provide commented, minimal, self-contained, reproducible code.
______________________________________________
R-help at 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.
On Sun, 22 Nov 2009, Ravi Varadhan wrote:
Hi Torsten,
Hi Ravi,
It would be useful to "warn" the users that the multivariate normal probability calculated by "pmvnorm" using the GenzBretz algorithm is "random", i.e. the result can vary between repeated executions of the function.
only if a different seed is used.
This would prevent inappropriate use of pmvnorm such as computing derivatives of it (see this email thread).
?pmvt has "Randomized quasi-Monte Carlo methods are used for the computations." and appropriate references. In addition, the new book by Alan Genz and Frank Bretz covers all technical details in depth, so the procedures are well documented. Anyway, I'll add a statement to ?pmvnorm. Best wishes, Torsten
Hi Torsten, Thanks for you comment. If you have some free time to spare, partial derivatives with respect to bounds and correlation coefficients would be great for pmvnorm! In complex problems, optim is not very good at estimating the hessian numerically and first order derivatives help to build an OPG estimator, which is not very good as compared to an analytical hessian but still much better than the numerical hessian provided by optim i have found the problems I study. Best, Stephane 2009/11/23 Torsten Hothorn <Torsten.Hothorn at stat.uni-muenchen.de>:
On Sun, 22 Nov 2009, Ravi Varadhan wrote:
Hi Torsten,
Hi Ravi,
It would be useful to "warn" the users that the multivariate normal probability calculated by "pmvnorm" using the GenzBretz algorithm is "random", i.e. the result can vary between repeated executions of the function.
only if a different seed is used.
This would prevent inappropriate use of pmvnorm such as computing derivatives of it (see this email thread).
?pmvt has "Randomized quasi-Monte Carlo methods are used for the computations." and appropriate references. In addition, the new book by Alan Genz and Frank Bretz covers all technical details in depth, so the procedures are well documented. Anyway, I'll add a statement to ?pmvnorm. Best wishes, Torsten
______________________________________________ R-help at 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.