Mixture distribution power

gistfile1.txt

# mixture distribution power sims
rm(list=ls())

# from http://stats.stackexchange.com/questions/27088/shapiro-francia-test-error
sf.testBIG=function (x) 
  
{
  DNAME <- deparse(substitute(x))
  x <- sort(x[complete.cases(x)])
  n <- length(x)
  if ((n < 5 || n >= 1000000))   ###            <-----here is my edit
    stop("sample size must be between 5 and 10,000")
  y <- qnorm(ppoints(n, a = 3/8))
  W <- cor(x, y)^2
  u <- log(n)
  v <- log(u)
  mu <- -1.2725 + 1.0521 * (v - u)
  sig <- 1.0308 - 0.26758 * (v + 2/u)
  z <- (log(1 - W) - mu)/sig
  pval <- pnorm(z, lower.tail = FALSE)
  RVAL <- list(statistic = c(W = W), p.value = pval, method = 
                 "Shapiro-Francia normality test", 
               data.name = DNAME)
  class(RVAL) <- "htest"
  return(RVAL)
}

nsims = 100
ns = c(50,100,500,1000,5000,10000,50000,100000,500000,1000000)
ms = c(.8,.85,.9,.95)

d <- data.frame()
for (n in ns) {
  for (m in ms) {
    for (s in 1:nsims) {
      d <- rbind(d, 
                 data.frame(n = n,
                            m = m, 
                            s = s,
                            p = sf.testBIG(c(rnorm(n*m),
                                               rnorm(n*(1-m), 1, 1))
                                             )$p.value < .05))
    }
  }
}

ms <- d %>% group_by(n, m) %>%
  summarise(p = mean(p))

quartz()

qplot(log10(n), p, col=factor(m), group=factor(m), geom="line", 
      data=ms) + 
  xlab("N") + ylab("Power") +
  theme_classic()