Re: [R] Computing multivariate normal probabilities. Was: Re: Problem with Numerical derivatives (numDeriv) and mvtnorm
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=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 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] mac os X: mprobit fails to install
Hi all, any chance that someone got through the installation problem of mprobit on mac os X? 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.
Re: [R] mac os X: mprobit fails to install
Thanks. I have 10.5.8 with R 2.10.0 now (i still had 2.9 in my previous messages). I also have gcc-4.2 installed but no Xcode package. It still fails to install - can it be the Xcode package? Where can I find it - I don't have my install CDs with me and will not get them soon? Stephane 2009/11/22 David Winsemius dwinsem...@comcast.net: There were quite a few implicit declaration warning messages when I followed Phil's advice, but I do seem to get a complete build on a Mac 10.5.8 running 64 bit R 2.10.0. Have you installed the Xcode package? The gcc-4.2? -- David. On Nov 22, 2009, at 3:15 PM, SL wrote: I have tried your command but without success. Any idea? Here is my log: Macbook:$ PKG_CFLAGS=-I/usr/include/sys R CMD INSTALL mprobit_0.9-2.tar.gz * Installing to library ‘/Users/stephaneluchini/Library/R/2.9/library’ * Installing *source* package ‘mprobit’ ... ** libs ** arch - i386 sh: make: command not found ERREUR : compilation failed pour le package ‘mprobit’ * Removing ‘/Users/stephaneluchini/Library/R/2.9/library/mprobit’ 2009/11/22 Phil Spector spec...@stat.berkeley.edu: Stephane - The check log indicated that malloc.h couldn't be found. Since that header file is located in /usr/include/sys on Macs, you could do the following: 1. Download mprobit_0.9-2.tar.gz from your local CRAN mirror. 2. At a terminal, type PKG_CFLAGS=-I/usr/include/sys R CMD INSTALL mprobit_0.9-2.tar.gz They'll be some warning messages, but the package should get built. - Phil Spector Statistical Computing Facility Department of Statistics UC Berkeley spec...@stat.berkeley.edu On Sun, 22 Nov 2009, stephane Luchini wrote: Hi all, any chance that someone got through the installation problem of mprobit on mac os X? 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. __ 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. David Winsemius, MD Heritage Laboratories West Hartford, CT __ 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.
[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.
[R] Partial derivatives of the multivariate cumulative distribution
I'm currently using the mvtnorm package to model unobserved heterogeneity in a structural model and using optim to estimate the model. I have got good clues that convergence is not really a problem but the hessian matrix estimate is very bad. To overcome this problem, I'm constructing an OPG estimator of the information matrix and I was wondering if there were an easy way to obtain partial derivatives of say for instance: P1 - pmvnorm(lower=c(-Inf,-Inf,-Inf,-Inf),upper=c(theta1,theta2,theta3,theta4),corr=ssigma) with respect to the mean parameters theta1, theta2, theta3, theta4 and the non-diagonal parameters in sigma, hence $\partial P_1 / \partial \theta_1$, etc... I can deal with numerical or analytical partial derivatives - a gradient would be fine since all observations share the same partial derivative. 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.
