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 :
>
> 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
> Date: Saturday, November 21, 2009 8:15 pm
> Subject: Re: [R] Problem with Numerical derivatives (numDeriv) and mvtnorm
> To: SL
> 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=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
>>
>> PLEASE do read the posting guide
>> and provide commented, minimal, self-contained, reproducible code.
>
> __
> R-help@r-project.org mailing list
> https://st
