[R] question regarding qmvnorm

[Torsten Hothorn]

On Wed, 20 Apr 2011, Ravi Varadhan wrote:

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

thanks, Ravi, for the report -- indeed, `tail = "upper"' caused 
`approx_interval()' to act insane. Fixed in 0.9-99 on R-forge and soon on 
CRAN.

Best,

Torsten

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