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]]
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.
On 21/04/11 09:59, li li wrote:> 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(1) In the first instance you should probably contact the maintainer of the mvtnorm package: require(utils) maintainer("mvtnorm") (2) The code you enclosed indicates that you don't understand anything about functions. You have m, rho, k, and alpha given as *arguments* to your function, and yet you specify them explicitly within the body of the function. This makes no sense at all. If you are going to do something it pays to *understand* what you are doing rather than just slapping down some code at random and hoping that it will work. Delete the first four (non-blank!) lines of your code and then call cc_f(10,0.1,2,0.5) This does indeed throw an error --- from uniroot(). As I say, contact the package maintainer. cheers, Rolf Turner> 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) > > }