[R] Computing multivariate normal probabilities. Was: Re: Problem with Numerical derivatives (numDeriv) and mvtnorm

[stephane Luchini]

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.
>