############ Set home directory to Mypackages in R ... ie this can be       #######################
############ done with File > Change dir...                                 #######################
############ Issue the command  rm(list=ls(all=TRUE)) at the r prompt       #######################
############ to clean everything up.                                        #######################
############ Creates package in Mypackages by pasting all of what follows   #######################
############ at the R prompt.  Please note that you will have to edit the   #######################
############ DESCRIPTION file and the .rd files next.  See instructions     #######################
############ below and the "Writing R Extensions" for more guidance.        #######################

##########################################################################################
#############       Alan T. Arnholt         ##############################################
#############       4/7/2005                ##############################################

nn <- function(mu0,mu1,sigma,Alpha,Beta,alternative = c("two.sided", "less", "greater"))
{
#
# my axis used later...
#
my.axis <- function(side, at, labels,...)
{
for(i in seq(along=at)) axis(side=side, at=at[i], labels=labels[i],...)
}
#
alternative <- match.arg(alternative)
#
#   Checkers
#
    if (!missing(mu0) && (length(mu0) != 1 || is.na(mu0)))
        stop("mu0 must be a single number")
    if (!missing(mu1) && (length(mu1) != 1 || is.na(mu1)))
        stop("mu1 must be a single number")
    if (!missing(sigma) && (length(sigma) != 1 || is.na(sigma)))
        stop("sigma must be a single number")
    if (!missing(Alpha) && (length(Alpha) != 1 || !is.finite(Alpha) ||
        Alpha < 0 || Alpha > 1))
        stop("Alpha must be a single number between 0 and 1")
    if (!missing(Beta) && (length(Beta) != 1 || !is.finite(Beta) ||
        Beta < 0 || Beta > 1))
        stop("Beta must be a single number between 0 and 1")
    if ( (alternative=="greater") && (mu0 > mu1) )
        stop("mu0 must be less than mu1 to test greater than alternative")
    if ( (alternative=="less") && (mu1 > mu0) )
        stop("mu0 must be greater than mu1 to test less than alternative")
    if( mu0==mu1)
        warning("Note: beta can be no smaller than 1-alpha regardless of n")


###########################################################################################
###     Start

if (alternative == "two.sided" && mu0==mu1)
{
#
n <- 1
cvu <- qnorm(1-Alpha/2,mu0,sigma/sqrt(n))
cvl <- qnorm(Alpha/2,mu0,sigma/sqrt(n))
Power <- (1-pnorm(cvu,mu1,sigma/sqrt(n))) + (pnorm(cvl,mu1,sigma/sqrt(n)))
Beta <- round(1-Power,2)
x <- seq(min(mu0,mu1)-3.4*sigma/sqrt(n),max(mu1,mu0)+3.4*sigma/sqrt(n),length=1000)
y0 <- dnorm(x,mu0,sigma/sqrt(n))
y1 <- dnorm(x,mu1,sigma/sqrt(n))
cvl <- qnorm(Alpha/2,mu0,sigma/sqrt(n))
cvu <- qnorm(1-Alpha/2,mu0,sigma/sqrt(n))
plot(x,y0,type="l",xlab="",ylab="",main="",col="blue",lwd=2, xaxt="n", yaxt="n")
title(main=list(expression(paste("For ", H[a],":",mu!=mu[0],", if ",mu[0]==mu[1], ", ",beta==1-alpha," for all ",n,"."  )  )))
legend( max(x)-.17*(max(x)-min(x)), dnorm(mu0,mu0,sigma/sqrt(n)), c("P(Type I error)","", "P(Type II error)"), col = c("lightblue",NA,"pink"), text.col=c("black",NA,"black"),pch=c(15,NA,15), bg='gray95',cex=.75,pt.cex=2.5,lwd=2)
size <- 1
text(min(x), .95*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(alpha==.(Alpha)),pos=4,col="black",cex=size)
text(min(x), .85*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(beta==.(Beta)),pos=4,col="black",cex=size)
text(min(x), .75*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(mu[0]==.(mu0)),pos=4,col="black",cex=size)
text(min(x), .65*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(mu[1]==.(mu1)),pos=4,col="black",cex=size)
text(min(x), .55*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(sigma==.(sigma)),pos=4,col="black",cex=size)
# text(min(x), .45*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(n==.(n)),pos=4,col="black",cex=size)
my.axis(1,at=c(mu0,cvl,cvu),labels=c(expression(mu[0]),expression(cv[alpha/2]),expression(cv[1-alpha/2])), cex.axis=1.2 )
lines(x,y1,col="red",lwd=2)
xa <- c(cvu, seq(cvu,max(x),length=500),max(x))
ya <- c(0,dnorm(seq(cvu,max(x),length=500),mu0,sigma/sqrt(n)),0)
polygon(xa,ya,col="lightblue")
xc <- c(cvl, seq(cvl,min(x),length=500),min(x))
yc <- c(0,dnorm(seq(cvl,min(x),length=500),mu0,sigma/sqrt(n)),0)
polygon(xc,yc,col="lightblue")
xb <- c(cvl,seq(cvl,cvu,length=500), cvu)
yb <- c(0,dnorm(seq(cvl,cvu,length=500),mu1,sigma/sqrt(n)),0)
polygon(xb,yb,col="pink")
lines(x,y1,col="red",lwd=3)
lines(x,y0,col="blue",lwd=3)
lines(c(mu0,mu0),c(0,dnorm(mu0,mu0,sigma/sqrt(n))),col="blue",lty=2)
lines(c(mu1,mu1),c(0,dnorm(mu1,mu1,sigma/sqrt(n))),col="red",lty=2)
lines(c(min(x),max(x)),c(0,0),lwd=3)
box()
}
else if (alternative=="two.sided")
{
n <- 0
Power <- 0
while(Power < 1-Beta)
{
n <- n + 1
cvu <- qnorm(1-Alpha/2,mu0,sigma/sqrt(n))
cvl <- qnorm(Alpha/2,mu0,sigma/sqrt(n))
Power <- (1-pnorm(cvu,mu1,sigma/sqrt(n))) + (pnorm(cvl,mu1,sigma/sqrt(n)))
}
Beta <- round(1-Power,2)
x <- seq(min(mu0,mu1)-3.4*sigma/sqrt(n),max(mu1,mu0)+3.4*sigma/sqrt(n),length=1000)
y0 <- dnorm(x,mu0,sigma/sqrt(n))
y1 <- dnorm(x,mu1,sigma/sqrt(n))
cvl <- qnorm(Alpha/2,mu0,sigma/sqrt(n))
cvu <- qnorm(1-Alpha/2,mu0,sigma/sqrt(n))
plot(x,y0,type="l",xlab="",ylab="",main="",col="blue",lwd=2, xaxt="n", yaxt="n")
title(main=list(expression(paste("Testing the Alternative Hypothesis  "  , H[a],":",mu!=mu[0])  )))
legend( max(x)-.17*(max(x)-min(x)), dnorm(mu0,mu0,sigma/sqrt(n)), c("P(Type I error)","", "P(Type II error)"), col = c("lightblue",NA,"pink"), text.col=c("black",NA,"black"),pch=c(15,NA,15), bg='gray95',cex=.75,pt.cex=2.5,lwd=2)
size <- 1
text(min(x), .95*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(alpha==.(Alpha)),pos=4,col="black",cex=size)
text(min(x), .85*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(beta==.(Beta)),pos=4,col="black",cex=size)
text(min(x), .75*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(mu[0]==.(mu0)),pos=4,col="black",cex=size)
text(min(x), .65*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(mu[1]==.(mu1)),pos=4,col="black",cex=size)
text(min(x), .55*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(sigma==.(sigma)),pos=4,col="black",cex=size)
text(min(x), .45*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(n==.(n)),pos=4,col="black",cex=size)
my.axis(1,at=c(mu0,cvl,cvu,mu1),labels=c(expression(mu[0]),expression(cv[alpha/2]),expression(cv[1-alpha/2]),expression(mu[1])), cex.axis=1.2 )
lines(x,y1,col="red",lwd=2)
xa <- c(cvu, seq(cvu,max(x),length=500),max(x))
ya <- c(0,dnorm(seq(cvu,max(x),length=500),mu0,sigma/sqrt(n)),0)
polygon(xa,ya,col="lightblue")
xc <- c(cvl, seq(cvl,min(x),length=500),min(x))
yc <- c(0,dnorm(seq(cvl,min(x),length=500),mu0,sigma/sqrt(n)),0)
polygon(xc,yc,col="lightblue")
xb <- c(cvl,seq(cvl,cvu,length=500), cvu)
yb <- c(0,dnorm(seq(cvl,cvu,length=500),mu1,sigma/sqrt(n)),0)
polygon(xb,yb,col="pink")
lines(x,y1,col="red",lwd=3)
lines(x,y0,col="blue",lwd=3)
lines(c(mu0,mu0),c(0,dnorm(mu0,mu0,sigma/sqrt(n))),col="blue",lty=2)
lines(c(mu1,mu1),c(0,dnorm(mu1,mu1,sigma/sqrt(n))),col="red",lty=2)
lines(c(min(x),max(x)),c(0,0),lwd=3)
box()
}

else
{
n <- round((( (sigma*(qnorm(1-Alpha)-qnorm(Beta))) / (mu1-mu0))^2+1/2),0)
}

if (alternative == "greater" && (mu0==mu1) )
{
n <- 1
x <- seq(mu0-3.4*sigma/sqrt(n),mu1+3.4*sigma/sqrt(n),length=1000)
y0 <- dnorm(x,mu0,sigma/sqrt(n))
y1 <- dnorm(x,mu1,sigma/sqrt(n))
cv <- qnorm(1-Alpha,mu0,sigma/sqrt(n))
Beta <- round(pnorm(cv,mu1,sigma/sqrt(n)),3)
plot(x,y0,type="l",xlab="",ylab="",main="",col="blue",lwd=2, xaxt="n", yaxt="n")
title(main=list(expression(paste("For ", H[a],":",mu > mu[0],", if ",mu[0]==mu[1], ", ",beta==1-alpha," for all ",n,"."  )  )))
legend( max(x)-.17*(max(x)-min(x)), dnorm(mu0,mu0,sigma/sqrt(n)), c("P(Type I error)","", "P(Type II error)"), col = c("lightblue",NA,"pink"), text.col=c("black",NA,"black"),pch=c(15,NA,15), bg='gray95',cex=.75,pt.cex=2.5,lwd=2)
size <- 1
text(min(x), .95*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(alpha==.(Alpha)),pos=4,col="black",cex=size)
text(min(x), .85*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(beta==.(Beta)),pos=4,col="black",cex=size)
text(min(x), .75*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(mu[0]==.(mu0)),pos=4,col="black",cex=size)
text(min(x), .65*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(mu[1]==.(mu1)),pos=4,col="black",cex=size)
text(min(x), .55*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(sigma==.(sigma)),pos=4,col="black",cex=size)
# text(min(x), .45*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(n==.(n)),pos=4,col="black",cex=size)
my.axis(1,at=c(mu0,cv),labels=c(expression(mu[0]),expression(cv[1-alpha])), cex.axis=1.2 )
lines(x,y1,col="red",lwd=2)
xa <- c(cv, seq(cv,max(x),length=500),max(x))
ya <- c(0,dnorm(seq(cv,max(x),length=500),mu0,sigma/sqrt(n)),0)
polygon(xa,ya,col="lightblue")
xb <- c(min(x),seq(min(x),cv,length=500), cv)
yb <- c(0,dnorm(seq(min(x),cv,length=500),mu1,sigma/sqrt(n)),0)
polygon(xb,yb,col="pink")
lines(x,y1,col="red",lwd=3)
lines(x,y0,col="blue",lwd=3)
lines(c(mu0,mu0),c(0,dnorm(mu0,mu0,sigma/sqrt(n))),col="blue",lty=2)
lines(c(mu1,mu1),c(0,dnorm(mu1,mu1,sigma/sqrt(n))),col="red",lty=2)
lines(c(min(x),max(x)),c(0,0),lwd=3)
box()
}
else if (alternative=="greater")
# do normal stuff
{
x <- seq(mu0-3.4*sigma/sqrt(n),mu1+3.4*sigma/sqrt(n),length=1000)
y0 <- dnorm(x,mu0,sigma/sqrt(n))
y1 <- dnorm(x,mu1,sigma/sqrt(n))
cv <- qnorm(1-Alpha,mu0,sigma/sqrt(n))
Beta <- round(pnorm(cv,mu1,sigma/sqrt(n)),3)
plot(x,y0,type="l",xlab="",ylab="",main="",col="blue",lwd=2, xaxt="n", yaxt="n")
title(main=list(expression(paste("Testing the Alternative Hypothesis  "  , H[a],":",mu>mu[0])  )))
legend( max(x)-.17*(max(x)-min(x)), dnorm(mu0,mu0,sigma/sqrt(n)), c("P(Type I error)","", "P(Type II error)"), col = c("lightblue",NA,"pink"), text.col=c("black",NA,"black"),pch=c(15,NA,15), bg='gray95',cex=.75,pt.cex=2.5,lwd=2)
size <- 1
text(min(x), .95*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(alpha==.(Alpha)),pos=4,col="black",cex=size)
text(min(x), .85*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(beta==.(Beta)),pos=4,col="black",cex=size)
text(min(x), .75*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(mu[0]==.(mu0)),pos=4,col="black",cex=size)
text(min(x), .65*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(mu[1]==.(mu1)),pos=4,col="black",cex=size)
text(min(x), .55*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(sigma==.(sigma)),pos=4,col="black",cex=size)
text(min(x), .45*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(n==.(n)),pos=4,col="black",cex=size)
my.axis(1,at=c(mu0,cv,mu1),labels=c(expression(mu[0]),expression(cv[1-alpha]),expression(mu[1])), cex.axis=1.2 )
lines(x,y1,col="red",lwd=2)
xa <- c(cv, seq(cv,max(x),length=500),max(x))
ya <- c(0,dnorm(seq(cv,max(x),length=500),mu0,sigma/sqrt(n)),0)
polygon(xa,ya,col="lightblue")
xb <- c(min(x),seq(min(x),cv,length=500), cv)
yb <- c(0,dnorm(seq(min(x),cv,length=500),mu1,sigma/sqrt(n)),0)
polygon(xb,yb,col="pink")
lines(x,y1,col="red",lwd=3)
lines(x,y0,col="blue",lwd=3)
lines(c(mu0,mu0),c(0,dnorm(mu0,mu0,sigma/sqrt(n))),col="blue",lty=2)
lines(c(mu1,mu1),c(0,dnorm(mu1,mu1,sigma/sqrt(n))),col="red",lty=2)
lines(c(min(x),max(x)),c(0,0),lwd=3)
box()
}

if (alternative =="less" && (mu0==mu1))
{
n <- 1
x <- seq(mu1-3.4*sigma/sqrt(n),mu0+3.4*sigma/sqrt(n),length=1000)
y0 <- dnorm(x,mu0,sigma/sqrt(n))
y1 <- dnorm(x,mu1,sigma/sqrt(n))
cv <- qnorm(Alpha,mu0,sigma/sqrt(n))
Beta <- 1-Alpha
plot(x,y0,type="l",xlab="",ylab="",main="",col="blue",lwd=2, xaxt="n", yaxt="n")
title(main=list(expression(paste("For ", H[a],":",mu < mu[0],", if ",mu[0]==mu[1], ", ",beta==1-alpha," for all ",n,"."  )  )))
legend( max(x)-.17*(max(x)-min(x)), dnorm(mu0,mu0,sigma/sqrt(n)), c("P(Type I error)","", "P(Type II error)"), col = c("lightblue",NA,"pink"), text.col=c("black",NA,"black"),pch=c(15,NA,15), bg='gray95',cex=.75,pt.cex=2.5,lwd=2)
size <- 1
text(min(x), .95*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(alpha==.(Alpha)),pos=4,col="black",cex=size)
text(min(x), .85*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(beta==.(Beta)),pos=4,col="black",cex=size)
text(min(x), .75*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(mu[0]==.(mu0)),pos=4,col="black",cex=size)
text(min(x), .65*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(mu[1]==.(mu1)),pos=4,col="black",cex=size)
text(min(x), .55*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(sigma==.(sigma)),pos=4,col="black",cex=size)
# text(min(x), .45*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(n==.(n)),pos=4,col="black",cex=size)
my.axis(1,at=c(cv,mu0),labels=c(expression(cv[alpha]),expression(mu[0])), cex.axis=1.2 )
lines(x,y1,col="red",lwd=2)
xa <- c(cv, seq(cv,min(x),length=500),min(x))
ya <- c(0,dnorm(seq(cv,min(x),length=500),mu0,sigma/sqrt(n)),0)
polygon(xa,ya,col="lightblue")
xb <- c(cv,seq(cv, max(x),length=500), max(x))
yb <- c(0,dnorm(seq(cv, max(x),length=500),mu1,sigma/sqrt(n)),0)
polygon(xb,yb,col="pink")
lines(x,y1,col="red",lwd=3)
lines(x,y0,col="blue",lwd=3)
lines(c(mu0,mu0),c(0,dnorm(mu0,mu0,sigma/sqrt(n))),col="blue",lty=2)
lines(c(mu1,mu1),c(0,dnorm(mu1,mu1,sigma/sqrt(n))),col="red",lty=2)
lines(c(min(x),max(x)),c(0,0),lwd=3)
box()
}
else if (alternative=="less")
{
x <- seq(mu1-3.4*sigma/sqrt(n),mu0+3.4*sigma/sqrt(n),length=1000)
y0 <- dnorm(x,mu0,sigma/sqrt(n))
y1 <- dnorm(x,mu1,sigma/sqrt(n))
cv <- qnorm(Alpha,mu0,sigma/sqrt(n))
plot(x,y0,type="l",xlab="",ylab="",main="",col="blue",lwd=2, xaxt="n", yaxt="n")
title(main=list(expression(paste("Testing the Alternative Hypothesis  "  , H[a],":",mu<mu[0])  )))
legend( max(x)-.17*(max(x)-min(x)), dnorm(mu0,mu0,sigma/sqrt(n)), c("P(Type I error)","", "P(Type II error)"), col = c("lightblue",NA,"pink"), text.col=c("black",NA,"black"),pch=c(15,NA,15), bg='gray95',cex=.75,pt.cex=2.5,lwd=2)
size <- 1
text(min(x), .95*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(alpha==.(Alpha)),pos=4,col="black",cex=size)
text(min(x), .85*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(beta==.(Beta)),pos=4,col="black",cex=size)
text(min(x), .75*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(mu[0]==.(mu0)),pos=4,col="black",cex=size)
text(min(x), .65*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(mu[1]==.(mu1)),pos=4,col="black",cex=size)
text(min(x), .55*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(sigma==.(sigma)),pos=4,col="black",cex=size)
text(min(x), .45*dnorm(mu0,mu0,sigma/sqrt(n)),bquote(n==.(n)),pos=4,col="black",cex=size)

my.axis(1,at=c(mu1,cv,mu0),labels=c(expression(mu[1]),expression(cv[alpha]),expression(mu[0])), cex.axis=1.2 )
lines(x,y1,col="red",lwd=2)
xa <- c(cv, seq(cv,min(x),length=500),min(x))
ya <- c(0,dnorm(seq(cv,min(x),length=500),mu0,sigma/sqrt(n)),0)
polygon(xa,ya,col="lightblue")
xb <- c(cv,seq(cv, max(x),length=500), max(x))
yb <- c(0,dnorm(seq(cv, max(x),length=500),mu1,sigma/sqrt(n)),0)
polygon(xb,yb,col="pink")
lines(x,y1,col="red",lwd=3)
lines(x,y0,col="blue",lwd=3)
lines(c(mu0,mu0),c(0,dnorm(mu0,mu0,sigma/sqrt(n))),col="blue",lty=2)
lines(c(mu1,mu1),c(0,dnorm(mu1,mu1,sigma/sqrt(n))),col="red",lty=2)
lines(c(min(x),max(x)),c(0,0),lwd=3)
box()
}

if(mu0==mu1)
{
cat("\n")
cat("When mu0=mu1, the probability of committing a type II error (Beta) can", "\n")
cat("be no smaller than 1-alpha =", Beta,"regardless of n.",  "\n")
cat("\n\n")
}
else
{
cat("\n")
cat("The required sample size (n) to test the user specified hypothesis for the ", "\n")
cat("probability of type I error (alpha) to equal", Alpha, "and the probability ", "\n")
cat("of a type II error (beta) to be no more than", Beta, "is", n,".", "\n")
cat("\n\n")
}
}

Alan <- c(12,14,16,18,20)
Susie <- c(11,13,15,17,19)

####### Package Now ####

KDATA <- objects()
package.skeleton(list=KDATA, name="EXAMPLE")

###  Manually remove KDATA at the r prompt with > rm(KDATA)
###  Edit the DESCRIPTION file located in the EXAMPLE folder under Mypackages
###  See "Writing R extensions" for more details on what needs to be edited or have a look
###  at the DESCRIPTION file of packages others have written.
###  Edit the .rd files in the MAN directory beneath the EXAMPLE directory.  Again,
###  see "Writing R extensions" for more details.
###  Use commands below from command line in DOS prompt (In the Mypackages folder/directory)
###  These commands build a binary library...read the manuals if you want a source file or 
###  if you need more help.
###  >RCMD build --binary EXAMPLE   # Type in the Mypackages folder/directory!
###  Use packages > Install package(s) from local zip files to install package.