stor555_N-P-hypothesistest

a=seq(-10,0,by=0.1)
aa = as.data.frame(a)
alpha=seq(0.01, 0.91, by=0.1)

diff_fun = function(a) pnorm(a+qnorm(1-alpha/2))-pnorm(a+qnorm(1-alpha)) - pnorm(a+qnorm(alpha/2))
powerNPfun = function(a) 1-pnorm(a+qnorm(1-alpha))
powerIfun = function(a) 1-pnorm(a+qnorm(1-alpha/2))+pnorm(a+qnorm(alpha/2))


D =apply(aa,1,diff_fun)
powerNP=apply(aa,1,powerNPfun)
powerI=apply(aa,1,powerIfun)


par(mfrow=c(1,2))
plot(a,powerNP[1,],type='l',col=1)
for(i in 2:length(alpha)) points(a, powerNP[i,],type='l',col=i)
legend("bottomleft",legend=paste("alpha=",as.character(alpha[1:10])), pch="--", col=1:10)
title("N-P")

plot(a,powerI[1,],type='l',col=1)
for(i in 2:length(alpha)) points(a, powerI[i,],type='l',col=i)
legend("bottomleft",legend=paste("alpha=",as.character(alpha[1:10])), pch="--", col=1:10)
title("Intuitive")

par(mfrow=c(1,1))
plot(a, D[1,],type='l',col=1,ylim=c(0,0.2),ylab="Power(N-P) - Power(Intuitive)",xlab="a = (theta0-theta1)*sqrt(n)/sigma")
for(i in 2:length(alpha)) points(a, D[i,],type='l',col=i)
legend("topleft",legend=paste("alpha=",as.character(alpha[1:10])), pch="--", col=1:10)
title("Power(NP) - Power(Intivtive)")

par(mfrow=c(1,2))
persp(alpha,a,D,col=c(6,7),theta = 60,xlim=c(0,1),ylim=c(-10,0))
persp(alpha,a,D,col=c(6,7),theta = 240,xlim=c(0,1),ylim=c(-10,0))