############### compute power

n <- 26         
p.vec <-  seq(from=.4,by=.01,to=1)           # range of p's (Ha)
power <- rep(0,length(p.vec))

for (i in 1:length(p.vec)) {
	p <- 	p.vec[i] 
	s <- sqrt(p*(1-p)/n)/sqrt(1/(4*n))
	Ep <- (p-.5)/sqrt(1/(4*n))

	z <- (1.96 - Ep)/s
	power[i] <- 1-pnorm(z)
}
plot(p.vec,power,type="l",xlab="p",ylab="Power",cex.lab=1.2,cex.axis=1.2)


###### repeat for different values of n

plot(p.vec,power,type="n",xlab="p",ylab="Power",cex.lab=1.2,cex.axis=1.2)

p.vec <-  seq(from=.4,by=.01,to=1)
n.vec <- seq(from=20,by=50,to=250)
n.vec <- 2^(4:10)
for (j in 1:length(n.vec)){
	n <- n.vec[j]
	power <- rep(0,length(p.vec))

	for (i in 1:length(p.vec)) {
		p <- 	p.vec[i] 
		s <- sqrt(p*(1-p)/n)/sqrt(1/(4*n))
		Ep <- (p-.5)/sqrt(1/(4*n))


		z <- (1.96 - Ep)/s
		power[i] <- 1-pnorm(z)
	}
lines(p.vec,power,col=j)
}
legend(x=.8,y=.2,legend=as.factor(n.vec),col=1:length(n.vec),lty=1,title="Sample size")