##############################
#
# Script: Solutions 1 (assignment01a.R)
#
##############################



####################
# Problem 01.1

# P[Z>1] 
#   = 1-P[Z<1]
#   = 1-0.8433
#   = 0.1587
1-pnorm(1)

# P[|Z|>1] 
#   = P[Z>1]   + P[-Z>1]
#   = P[Z>1]   + P[Z<-1]
#   = 1-P[Z<1] + P[Z<-1]
#   = 1-0.8433 + 0.1567
#   = 0.3173
1-pnorm(1)+ pnorm(-1)

# P[|T|<1] 
#   = 1 - P[|T|>1]
#   = 1 - (P[T>1]   + P[-T>1])
#   = 1 - (P[T>1]   + P[T<-1])
#   = 1 - (1-P[T<1] + P[T<-1])
#   = 1 - (0.1870   + 0.1870)
#   = 0.6261
1 - (1-pt(1,df=4) + pt(-1,df=4))

# P[|X+1|<6]
#   = 1 - P[|X+1|>6]
#   = 1 - (P[X+1>6] + P[-X-1>6])
#   = 1 - (P[X>5]   + P[X<-7])
#   = 1 - (1-P[X<5] + P[X<-7])
#   = 1 - (0.5000   + 3.17e-5)
#   = 0.5000
1 - (1-pnorm(5, m=5,s=3) + pnorm(-7, m=5,s=3) )

# P[3 < R < 5]
#   = P[R<5] - P[R<3]
#   = 0.4562 - 0.1912
#   = 0.2650
pchisq(5,df=6) - pchisq(3,df=6)



##############################
# Problem 01.2

# Let us use what we did in class on Thursday:

# The function to claculate the power

powercalc <- function(m0,mA,s0,sA,alpha) {
	cv <- qnorm(1-alpha, mean=m0,sd=s0)		# Calculate the CV (of the null)
	power <- 1-pnorm(cv, mean=mA,sd=sA)		# Calculate 1-beta (of the alternative)
	return(power)
}

# The command is
powercalc(2,4,1,1,0.05)   # = 0.63876



##############################
# Problem 01.3


# The function to claculate the power

powercalc <- function(m0,mA,s0,sA,alpha) {
	cv <- qnorm(1-alpha, mean=m0,sd=s0)		# Calculate the CV (of the null)
	power <- 1-pnorm(cv, mean=mA,sd=sA)		# Calculate 1-beta (of the alternative)
	return(power)
}


# Null hypothesis parameters
m0 <- 0
s0 <- 4

# Alternative hypothesis parameters
mA <- 0
sA <- 4

# Set the standard alpha level
alpha <- 0.05

# Clear space for these two variables
xval <- numeric()
pwr  <- numeric()

# Loop to calculate the power at various values of the mean
for(i in 1:100) {
  mA <- i/5
  pwr[i] <- powercalc(m0,mA,s0,sA,alpha)
  xval[i] <- mA
}

# Plot our results
plot(xval,pwr, type="l",ylab="Power",las=1,ylim=c(0,1),yaxs="i",xaxs="i",main="Power Curve",xlab="Difference in means")
axis(2,at=pwr[1],labels=round(pwr[1],2),las=1,cex.axis=0.8)