### Hypothesis Testing Examples
### Plotting the ROC curve.

#######################################################
#######################################################
# Example 1: Normal assume theta0 < theta1

# How to make a normal table in R/S-Plus

#######################################################
# Suppose 

alpha = 0.05

theta0 = 25; sigma0 = 3; n = 30;

xx = seq(theta0-5*sigma0/sqrt(n),theta0+5*sigma0/sqrt(n),.01)

yy = dnorm(xx, mean = theta0, sd = sigma0/sqrt(n))

plot(xx,yy,type="l")  # plot the null distribution

#######################################################
# Note that the qnorm function can be used to locate the
# critical values for hypothesis testing.  This is much 
# better than the table look up method presented in
# HTExamples1.r

# critical value using qnorm()

cv1 = qnorm(alpha,mean = theta0, sd = sigma0/sqrt(n), lower.tail = FALSE)

rr = cv1
rr

# compare the null with an alternative distribution

theta1 = 26

theta.min = min(theta0,theta1)
theta.max = max(theta0,theta1)

xx = seq(theta.min-5*sigma0/sqrt(n),theta.max+5*sigma0/sqrt(n),.01) 
	# theta values

yy = dnorm(xx, mean = theta0, sd = sigma0/sqrt(n))  # null
yy1 = dnorm(xx, mean = theta1, sd = sigma0/sqrt(n)) # alternative
plot(xx,yy,type="l")  # plot the null distribution
par(new=T)
plot(xx,yy1,type="l")  # plot the alternative distribution

pow1 = 1 - pnorm(rr, theta1, sigma0/sqrt(n))  # compute the power for theta1
pow1

# make the power curve

pow2 = 1 - pnorm(rr, xx, sigma0/sqrt(n))  # power over the theta values

plot(xx, pow2, type='l', main="power curve for alpha = 0.05")

# make an ROC curve, this curve plots the detection rate (power) for different
# values of the false positive rate (alpha) for a give alternative value of 
# the parameter.

theta1 = 26

inc = 0.01  # increment
alpha = seq((0+inc),(1-inc),inc)  # false postive values

cv2 = qnorm(alpha, mean = theta0, sd = sigma0/sqrt(n), lower.tail = FALSE)
rr = cv2  

pow3 = 1 - pnorm(rr, theta1, sigma0/sqrt(n))

plot(alpha,pow3,type='l',main="ROC for a value in the alternative",
   xlab="false postive rate (alpha)",ylab="detection rate (power)")

#######################################################
#######################################################
# Example 2: Exponential assume lambda0 > lambda1

#######################################################
# Suppose 

alpha = 0.05

lambda0 = 1.3; n = 3;

m = n/lambda0; v = n/lambda0^n;

xx = seq(0,m+5*v,.01)

xx = seq(0,5,0.1)

yy = dgamma(xx, shape = n, rate = lambda0)

plot(xx,yy,type="l")


#######################################################
# critical value using qgamma()

cv1 = qgamma(alpha, shape = n, rate = lambda0, lower.tail = FALSE)

rr = cv1
rr

# power

lambda1 = 0.52

pow1 = 1 - pgamma(rr, shape = n, rate = lambda1)  # compute the power for theta1
pow1

#######################################################
#######################################################
# Example 3: Poisson assume lambda0 > lambda1

#######################################################
# Suppose 

alpha = 0.05

lambda0 = 1; n = 10;

xx = seq(0,3*n*lambda0,1)

yy = dpois(xx, lambda = n*lambda0)

plot(xx,yy,type="h")  # plot the null distribution

#######################################################
# critical value using qpois()

cv1 = qpois(alpha, n*lambda0, lower.tail = TRUE)

# for a conservative test -1, since the q function is left continuous.

rr = cv1 - 1
rr

ppois(rr, n*lambda0)  # check conservative

# power

lambda1 = 0.52

pow1 = ppois(rr, n*lambda1)  # compute the power for lamda1
pow1

#######################################################
#######################################################
# Example 4: Bernoulli assume p0 < p1

#######################################################
# Suppose 

alpha = 0.05

p0 = .5; n = 20;

xx = seq(0,n,1)

yy = dbinom(xx, n, p0)

plot(xx,yy,type="h")  # plot the null distribution

#######################################################
# critical value using qbinom()

cv1 = qbinom(alpha, n, p0, lower.tail = FALSE)

# for a conservative test -1 

rr = cv1 
rr

1 - pbinom(rr, n, p0)  # check conservative

# power

p1 = 0.73

pow1 = 1 - pbinom(rr, n, p1)  # compute the power for p1
pow1