library(bbmle)
library(emdbook)
set.seed(1001)
a = c(20,10)
b=c(1,2)
k = 5
N <- 50
Nvec <- round(lseq(10,1000,length=20))
odd <- function(x) { x %% 2 == 1}
Nvec[odd(Nvec)] <- Nvec[odd(Nvec)]+1
nsim <- 200
power.N <- numeric(length(Nvec))
for (j in 1:length(Nvec)) {
  cat(j,"\n")
  N <- Nvec[j]
  pval <- numeric(nsim)
  for (i in 1:nsim) {
    cat("  ",i,"\n")
    g = rep(1:2,each=N/2)
    x= runif(N,min=0,max=5)
    y_det= a[g]*b/(b+x)
    y = rnbinom(N,mu=y_det,size=k)
    m0 = mle2(y~dnbinom(mu=a*b*x/(b+x),size=k),start=list(a=15,b=1,k=5))
    m1 = mle2(y~dnbinom(mu=a*b*x/(b+x),size=k),
      parameters=list(a~g,b~g),
      start=list(a=15,b=1,k=5))
    pval[i] <- anova(m0,m1)[2,"Pr(>Chisq)"]
  }
  power.N[j] <- sum(pval<0.05)/nsim
  save("power.N","Nvec",file="negbinhyp-batch2.RData")
}
save("power.N","Nvec",file="negbinhyp-batch2.RData")