Power Curve for Binomial Test



bin.pwr <- function(n,p0=1/2,RR=NULL,alpha=.05)
{
# plot power vs p for H0: p=p0, default p0=1/2
# RR if present is the chosen rejection region
# alpha is the target significance level

if(is.null(RR)){
    X = 0:n
    db = dbinom(X,n,p0)
    sdb = sort(db)
    cb = cumsum(sdb)
    m = max((1:(n+1))[cb <= alpha])
    R = X[db <= sdb[m]]
}
else{
    R=RR
}
Alpha <- sum(dbinom(R,n,p0))
p <- seq(.01,.99,.01)
P <- p
for(i in 1:length(p))
   {
   P[i] <- sum(dbinom(R,n,p[i]))
   }
plot(p,P,ylab="Power",type="l")
title(paste("Power curve for n=",n,"   alpha =",round(Alpha,4)))
}