From: Ravi Varadhan <rvaradhan_at_jhmi.edu>

Date: Wed, 20 Apr 2011 18:32:49 -0400

From: r-help-bounces_at_r-project.org [r-help-bounces_at_r-project.org] On Behalf Of li li [hannah.hlx_at_gmail.com] Sent: Wednesday, April 20, 2011 5:59 PM

To: r-help

Subject: [R] question regarding qmvnorm

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_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. Received on Wed 20 Apr 2011 - 22:34:54 GMT

Date: Wed, 20 Apr 2011 18:32:49 -0400

If you had told us what the error message was, my job would have been easier. But, at least you provied the code, so it was not hard for me to see where the problem was. There is a problem with the strategy used by `qmvnorm' to locate the initial interval in which the quantile is supposed to lie. In particular, the problem is with the approx_interval() function.

In your case, the interval computed by this function does not contain the root. Hence, uniroot() fails. The solution is to provide the interval that contains the roots. I am cc'ing Torsten so that he is aware of the problem.

The following works:

m <- 10

rho <- 0.1

k <- 2

alpha <- 0.05

cc_z <- rep(NA, length=m)

var <- matrix(c(1,rho,rho,1), nrow=2, ncol=2, byrow=T)
for (i in 1:m){

if (i <= k) {cc_z[i] <- qmvnorm((k*(k-1))/(m*(m-1))*alpha, tail="upper",
sigma=var, interval=c(0, 5))$quantile} else
{cc_z[i] <- qmvnorm((k*(k-1))/((m-i+k)*(m-i+k-1))*alpha, tail
="upper", sigma=var, interval=c(0,5))$quantile}
}

cc_z

*> cc_z
*

[1] 1.9438197 1.9438197 1.8910860 1.8303681 1.7590806 1.6730814 1.5652220
[8] 1.4219558 1.2116459 0.8267921

*>
*

Ravi.

From: r-help-bounces_at_r-project.org [r-help-bounces_at_r-project.org] On Behalf Of li li [hannah.hlx_at_gmail.com] Sent: Wednesday, April 20, 2011 5:59 PM

To: r-help

Subject: [R] question regarding qmvnorm

Dear all,

I wrote the following function previously. It worked fine with the old
mvtnorm package.

Somehow with the updated package, I got a error message when trying to use
the function.

I need some help. It is sort of urgent. Can anyone please take a look. The
function is the following.

Thank you very much!

Hannah

library(mvtnorm)

cc_f <- function(m, rho, k, alpha) {

m <- 10

rho <- 0.1

k <- 2

alpha <- 0.05

cc_z <- numeric(m)

var <- matrix(c(1,rho,rho,1), nrow=2, ncol=2, byrow=T)

for (i in 1:m){

if (i <= k) {cc_z[i] <- qmvnorm((k*(k-1))/(m*(m-1))*alpha, tail="upper", sigma=var)$quantile} else

{cc_z[i] <- qmvnorm((k*(k-1))/((m-i+k)*(m-i+k-1))*alpha, tail ="upper", sigma=var)$quantile}

}

cc <- 1-pnorm(cc_z)

return(cc)

} [[alternative HTML version deleted]] ______________________________________________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.

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. Received on Wed 20 Apr 2011 - 22:34:54 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Thu 21 Apr 2011 - 10:30:32 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*