R/mudiff.modwoc.equalvar.R

In SampleSizeMeans: Sample Size Calculations for Normal Means

`mudiff.modwoc.equalvar` <-
function(len, alpha, beta, n01, n02, level = 0.95, worst.level = 0.95, equal=TRUE)
{
  l.prior <- 2*qt((1+level)/2,2*alpha)*sqrt(beta/alpha*(1/n01+1/n02))
  
  if (l.prior <= len)
  {
    0 # prior knowledge is sufficient
  }
  else
  {    
    step <- 2^3
    found.upper.bound <- FALSE
    found.lower.bound <- FALSE
    direction <- 1
    n1 <- 16*beta/len^2*qnorm((1+level)/2)^2/qchisq(1-worst.level,2*alpha)
    n1 <- ceiling(max(1,n1-n01))

    while(step >= 2){
      step <- ifelse(found.upper.bound & found.lower.bound, step/2,step*2)
      n1 <- n1 + direction * step
      if(n1 <= 0){
        found.lower.bound <- TRUE
        n1 <- 0
      }
      n2 <- ifelse(!equal,n1+n01-n02,n1)
      n2 <- max(0,n2)
      ndot <- n1+n2

      a <- qf(worst.level, ndot, 2 * alpha)
      denom <- (1 + ndot/2/alpha * a)*8*beta
      num <- 1/(1/(n1+n01)+1/(n2+n02))*len*len*(ndot + 2 * alpha)
      b <- pf(num/denom, 1, ndot + 2 * alpha)

      if(b >= level) {
        found.upper.bound <- TRUE
        direction <- -1
      }
      else {
        found.lower.bound <- TRUE
        direction <- 1
      }
    }
    n1[direction == 1] <- n1 + 1
    n2 <- ifelse(!equal,n1+n01-n02,n1)
    n2 <- max(0,n2)
    c(n1,n2)
  }
}