[R] question about a program

[Joris Meys]

Most of the computation time is in the functions qvnorm. You can win a
little bit by optimizing the code, but the gain is relatively small.
You can also decrease the interval used to evaluate qvnorm to win some
speed there. As you look for the upper tail, no need to evaluate the
function in negative numbers. Eg :

pair_kFWER2=function(m,alpha,rho,k,pvaluessort){
library(mvtnorm)
cc_z <- numeric(m)
Var <- matrix(c(1,rho,rho,1), nrow=2, ncol=2, byrow=T)

lpart <- 1:k      # first part
hpart <- (k+1):m  # second part

# make first part of the cc_z
cc_z[lpart] <- qmvnorm((k*(k-1))/(m*(m-1))*alpha, tail="upper",
interval=c(0,4),sigma=Var)$quantile

# make second part of the cc_z
cc_z[hpart] <- sapply(hpart,function(i){
    qmvnorm((k*(k-1))/((m-i+k)*(m-i+k-1))*alpha,interval=c(0,4),
tail="upper", sigma=Var)$quantile
})

crit_cons <- 1-pnorm(cc_z)

# does the test whether values of crit_cons[i] are smaller than
# values pvaluessort[i] and returns a vector
# which says at which positions does not occur
# tail takes the last position. I guess that's what you wanted
out <- tail(which(!(crit_cons < pvaluessort)),1)
return(out)
}

> system.time(replicate(100,pair_kFWER(m,alpha,rho,k,pvaluessort)))
   user  system elapsed
   5.97    0.04    6.04

> system.time(replicate(100,pair_kFWER2(m,alpha,rho,k,pvaluessort)))
   user  system elapsed
   4.43    0.00    4.44

Check whether the function works correctly. It gives the same result
as the other one in my tests. Only in the case where your function
returns 0, the modified returns integer(0)

Cheers
Joris


On Wed, Jun 23, 2010 at 2:21 PM, li li <hannah.hlx at gmail.com> wrote:
> Dear all,
>   I have the following program for a multiple comparison procedure.
> There are two functions for the two steps. First step is to calculate the
> critical values,
> while the second step is the actual procedure [see below: program with two
> functions].
>   This  work fine. However, However I want to put them into one function
> for the convenience
> of later use [see below: program with one function]. Some how the big
> function works extremely
> slow.  For example I chose
> m <- 10
> rho <- 0.1
> k <- 2
> alpha <- 0.05
> pvaluessort <- sort(1-pnorm(rnorm(10))
>
> Is there anything that I did wrong?
>  Thank you!
>                    Hannah
>
>
> ######Program with two functions############
> ## first step
> library(mvtnorm)
> cc_f <- function(m, rho, k, alpha) {
>
> 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)
>                             }
> ## second step
> pair_kFWER=function(m,crit_cons,pvaluessort){
> k=0
> while((k <m)&&(crit_cons[m-k] < pvaluessort[m-k])){
> k=k+1
> }
> return(m-k)
> }
> ###########################################
> ##############Program with one function ##############
>
> pair_kFWER=function(m,alpha,rho,k,pvaluessort){
> library(mvtnorm)
> 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}
>               }
> crit_cons <- 1-pnorm(cc_z)
>
> k=0
> while((k <m)&&(crit_cons[m-k] < pvaluessort[m-k])){
> k=k+1
> }
> return(m-k)
> }
> #####################################################
>
>        [[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.
>



-- 
Joris Meys
Statistical consultant

Ghent University
Faculty of Bioscience Engineering
Department of Applied mathematics, biometrics and process control

tel : +32 9 264 59 87
Joris.Meys at Ugent.be
-------------------------------
Disclaimer : http://helpdesk.ugent.be/e-maildisclaimer.php