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 <rvarad...@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: rvarad...@jhmi.edu > > > ----- Original Message ----- > From: Ravi Varadhan <rvarad...@jhmi.edu> > Date: Saturday, November 21, 2009 8:15 pm > Subject: Re: [R] Problem with Numerical derivatives (numDeriv) and mvtnorm > To: SL <sl...@yahoo.fr> > Cc: r-help@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: rvarad...@jhmi.edu >> >> >> ----- Original Message ----- >> From: SL <sl...@yahoo.fr> >> 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 >> >> 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-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.