###############################################
## Author: Joshua M. Tebbs
## Date: 15 July 2021
## Update: 14 September 2021
## STAT 513 course notes: R Code
## Chapter 10
###############################################

# Figure 10.1
# Page 3
y = seq(38,62,0.01)
pdf = dnorm(y,50,sqrt(10))
plot(y,pdf,type="l",xlab=expression(bar(y)),ylab="Sampling distribution",cex.lab=1.25)
abline(h=0)
x = seq(55,62,0.01)
y = dnorm(x,50,sqrt(10))
polygon(c(55,x,62),c(0,y,0),col="lightblue")
points(x=55,y=0,pch=19,cex=1.25)


# Figure 10.2
# Page 5
y = seq(40,100,0.01)
pdf = 10*dnorm(y,50,10)*(pnorm(y,50,10))^9
plot(y,pdf,type="l",xlab=expression(y[(10)]),ylab="Sampling distribution",cex.lab=1.25)
abline(h=0)
x = seq(80,100,0.01)
y = 10*dnorm(x,50,10)*(pnorm(x,50,10))^9
polygon(c(80,x,100),c(0,y,0),col="lightblue")
points(x=80,y=0,pch=19,cex=1.25)

# Monte Carlo simulation
B = 100000 # MC samples
n = 10 # sample size
x = rnorm(B*n,50,10) # mean = 50 (i.e., H0 true)
M = matrix(x,B,n) # matrix of B rows and n columns
t = M[1,]*0 # initialise
for(i in 1:B) {
	t[i] <- max(M[i,]) # calculate the max for each sample
	} 
length(t[t>80])/B # proportion of samples where max > 80


# Figure 10.3
# Page 8
# Left
y = seq(0,400,0.1)
pdf = dgamma(y,15,1/10)
plot(y,pdf,type="l",xlab="t",ylab="Sampling distribution",cex.lab=1.25)
abline(h=0)
k = round(qgamma(0.05,15,1/10),1) # critical value
x = seq(0,k,0.1)
y = dgamma(x,15,1/10)
polygon(c(0,x,k),c(0,y,0),col="lightblue")
points(x=k,y=0,pch=19,cex=1.5)
# Right
y = seq(0,150,0.1)
pdf = dgamma(y,15,1/5)
plot(y,pdf,type="l",xlab="t",ylab="Sampling distribution",cex.lab=1.25)
abline(h=0)
x = seq(k,150,0.1)
y = dgamma(x,15,1/5)
polygon(c(k,x,150),c(0,y,0),col="lightblue")
points(x=k,y=0,pch=19,cex=1.5)


# Figure 10.4
# Page 10
# Left
y = seq(50,120,1)
prob = dpois(y,84)
plot(y,prob,type="h",xlab="t",ylab="Sampling distribution",ylim=c(0,max(prob)),cex.lab=1.25)
abline(h=0)
k = qpois(0.99,84) # critical value
for (i in k:120){
	lines(c(i,i),c(0,dpois(i,84)),lty=1,col="lightblue",lwd=2)
	}
points(x=qpois(0.99,84),y=0,pch=19,cex=1.5)

# Accident data
accidents = c(rep(0,32),rep(1,26),rep(2,12),rep(3,7),rep(4,4),rep(5,2),rep(6,1))
sum(accidents)


# Figure 10.5
# Page 12
x = seq(-5,5,0.001)
pdf = dt(x,10)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",ylim=c(0,0.4))
abline(h=0)
x = seq(qt(0.95,10),5,0.001)
y = dt(x,10)
polygon(c(qt(0.95,10),x,5),c(0,y,0),col="lightblue")
points(x=qt(0.95,10),y=0,pch=19,cex=1.5)
text(-0.03,0.09,expression(1-alpha),cex=1.25)
text(3.3,0.04,expression(alpha),cex=1.25)


# Figure 10.6
# Page 14
x = seq(0,30,0.001)
pdf = dchisq(x,10)
plot(x,pdf,type="l",lty=1,xlab="",ylab="",xaxt="n",yaxt="n",bty="n",ylim=c(0,0.1))
abline(h=0)
x = seq(0,qchisq(0.025,10),0.001)
y = dchisq(x,10)
polygon(c(0,x,qchisq(0.025,10)),c(0,y,0),col="lightblue")
points(x=qchisq(0.025,10),y=0,pch=19,cex=1.5)
x = seq(qchisq(0.975,10),30,0.001)
y = dchisq(x,10)
polygon(c(qchisq(0.975,10),x,30),c(0,y,0),col="lightblue")
points(x=qchisq(0.975,10),y=0,pch=19,cex=1.5)
text(9.8,0.02,expression(1-alpha),cex=1.25)
text(0.5,0.01,expression(alpha/2),cex=1.25)
text(25,0.0075,expression(alpha/2),cex=1.25)

# Filament data
bulbs = c(4.43,5.93,3.74,5.82,5.90,2.90,2.64,6.49,5.31,8.49,
	1.01,1.07,1.41,3.42,1.22,4.01,0.57,1.47,2.81,8.52,
	0.52,4.77,0.85,2.21,6.85,3.43,1.87,5.15,2.02,10.58)
sum(bulbs^2)


# Figure 10.7
# Page 16
x = seq(0,30,0.001)
pdf = dchisq(x,10)
plot(x,pdf,type="l",lty=1,xlab="",ylab="",xaxt="n",yaxt="n",bty="n",ylim=c(0,0.1))
abline(h=0)
x = seq(0,qchisq(0.025,10),0.001)
y = dchisq(x,10)
polygon(c(0,x,qchisq(0.025,10)),c(0,y,0),col="lightblue")
points(x=qchisq(0.025,10),y=0,pch=19,cex=1.5)
x = seq(qchisq(0.975,10),30,0.001)
y = dchisq(x,10)
polygon(c(qchisq(0.975,10),x,30),c(0,y,0),col="lightblue")
points(x=qchisq(0.975,10),y=0,pch=19,cex=1.5)
text(9.8,0.02,expression(1-alpha),cex=1.25)
text(0.5,0.01,expression(alpha/2),cex=1.25)
text(25,0.0075,expression(alpha/2),cex=1.25)


# Figure 10.8
# Page 17
# Seizure data 
control = c(49.2,44.5,45.8,95.8,30.1,36.5,82.3,87.9,105.0,95.2,97.5,105.0,58.1,86.6,58.3,
	72.8,116.7,45.1,70.3,77.4)
treatment = c(97.1,73.4,68.5,91.8,106.6,21.57,27.79,10.77,11.81,62.1,94.9,142.5,53,125,79.5,
	29.5,78.4,127.5,90.6,57.1,111.6,77.7,142,82.9,101.5)
boxplot(control-mean(control), treatment-mean(treatment),
 xlab="",names=c("Control","Treatment"),ylab="Serum ALP activity level (IU/L)",ylim=c(-100,100),col="lightblue")
test.stat = (sum((control-mean(control))^2)/20)/(sum((treatment-mean(treatment))^2)/25)
test.stat


# Figure 10.9
# Page 18
# Left
y = seq(0,1,0.001)
pdf = dbeta(y,8,1)
plot(y,pdf,type="l",xlab=expression(y[(n)]),ylab="Sampling distribution",cex.lab=1.25)
abline(h=0)
k = qbeta(0.95,8,1) # critical value
x = seq(k,1,0.001)
y = dbeta(x,8,1)
polygon(c(k,x,1),c(0,y,0),col="lightblue")
points(x=k,y=0,pch=19,cex=1.5)
# Right
n = 11
y = seq(0,1.25,0.001)
pdf = n*y^(n-1)/1.25^n
plot(y,pdf,type="l",xlab=expression(y[(n)]),ylab="Sampling distribution",xlim=c(0,1.25),cex.lab=1.25)
abline(h=0)
abline(v=1.25,lty=2)
k = qbeta(0.95,11,1) # critical value
x = seq(0,k,0.001)
y = n*x^(n-1)/1.25^n
polygon(c(0,x,k),c(0,y,0),col="lightblue")
points(x=k,y=0,pch=19,cex=1.5)


# Figure 10.10
# Page 20
x = seq(-5,5,0.001)
pdf = dt(x,10)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",ylim=c(0,0.4))
abline(h=0)
x = seq(-5,qt(0.025,10),0.001)
y = dt(x,10)
polygon(c(-5,x,qt(0.025,10)),c(0,y,0),col="lightblue")
points(x=qt(0.025,10),y=0,pch=19,cex=1.5)
x = seq(qt(0.975,10),5,0.001)
y = dt(x,10)
polygon(c(qt(0.975,10),x,5),c(0,y,0),col="lightblue")
points(x=qt(0.975,10),y=0,pch=19,cex=1.5)
text(-0.025,0.1,expression(1-alpha),cex=1.25)
text(-3.5,0.05,expression(alpha/2),cex=1.25)
text(3.5,0.05,expression(alpha/2),cex=1.25)


# Figure 10.11
# Page 21
# Left
x = seq(-5,5,0.001)
pdf = dt(x,10)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",ylim=c(0,0.4))
abline(h=0)
x = seq(qt(0.95,10),5,0.001)
y = dt(x,10)
polygon(c(qt(0.95,10),x,5),c(0,y,0),col="lightblue")
points(x=qt(0.95,10),y=0,pch=19,cex=1.5)
text(-0.03,0.09,expression(1-alpha),cex=2)
text(3.3,0.04,expression(alpha),cex=2)
# Right
x = seq(-5,5,0.001)
pdf = dt(x,10)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",ylim=c(0,0.4))
abline(h=0)
x = seq(-5,qt(0.05,10),0.001)
y = dt(x,10)
polygon(c(-5,x,qt(0.05,10)),c(0,y,0),col="lightblue")
points(x=qt(0.05,10),y=0,pch=19,cex=1.5)
text(-0.03,0.09,expression(1-alpha),cex=2)
text(-3.3,0.04,expression(alpha),cex=2)


# Figure 10.12
# Page 23
# Left
x = seq(-5,5,0.001)
pdf = dt(x,10)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",ylim=c(0,0.4))
abline(h=0)
x = seq(qt(0.95,10),5,0.001)
y = dt(x,10)
polygon(c(qt(0.95,10),x,5),c(0,y,0),col="lightblue")
points(x=qt(0.95,10),y=0,pch=19,cex=1.5)
text(-0.03,0.09,expression(1-alpha),cex=1.75)
text(3.3,0.04,expression(alpha),cex=1.75)
axis(side=1,at=c(0,qt(0.95,10)),labels=c(expression(theta[0]),expression(k)),cex.axis=1.75)
# Right
x = seq(-5,5,0.001)
pdf = dt(x,10)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",ylim=c(0,0.4))
abline(h=0)
x = seq(-5,qt(0.05,10),0.001)
y = dt(x,10)
polygon(c(-5,x,qt(0.05,10)),c(0,y,0),col="lightblue")
points(x=qt(0.05,10),y=0,pch=19,cex=1.5)
text(-0.03,0.09,expression(1-beta),cex=1.75)
text(-3.3,0.04,expression(beta),cex=1.75)
axis(side=1,at=c(qt(0.05,10),0),labels=c(expression(k),expression(theta[a])),cex.axis=1.75)


# Figure 10.13
# Page 27
x = seq(-5,5,0.001)
pdf = dt(x,10)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",ylim=c(0,0.4))
abline(h=0)
x = seq(-5,qt(0.025,10),0.001)
y = dt(x,10)
polygon(c(-5,x,qt(0.025,10)),c(0,y,0),col="lightblue")
points(x=qt(0.025,10),y=0,pch=19,cex=1.5)
x = seq(qt(0.975,10),5,0.001)
y = dt(x,10)
polygon(c(qt(0.975,10),x,5),c(0,y,0),col="lightblue")
points(x=qt(0.975,10),y=0,pch=19,cex=1.5)
text(-0.025,0.1,expression(1-alpha),cex=1.25)
text(-3.5,0.05,expression(alpha/2),cex=1.25)
text(3.5,0.05,expression(alpha/2),cex=1.25)


# Figure 10.14
# Page 28
x = seq(0,30,0.001)
pdf = dchisq(x,10)
plot(x,pdf,type="l",lty=1,xlab="",ylab="",xaxt="n",yaxt="n",bty="n",ylim=c(0,0.1))
abline(h=0)
x = seq(0,qchisq(0.025,10),0.001)
y = dchisq(x,10)
polygon(c(0,x,qchisq(0.025,10)),c(0,y,0),col="lightblue")
points(x=qchisq(0.025,10),y=0,pch=19,cex=1.5)
x = seq(qchisq(0.975,10),30,0.001)
y = dchisq(x,10)
polygon(c(qchisq(0.975,10),x,30),c(0,y,0),col="lightblue")
points(x=qchisq(0.975,10),y=0,pch=19,cex=1.5)
text(9.8,0.02,expression(1-alpha),cex=1.25)
text(0.5,0.01,expression(alpha/2),cex=1.25)
text(25,0.0075,expression(alpha/2),cex=1.25)


# Figure 10.15
# Page 31
x = seq(0,30,0.001)
pdf = dchisq(x,10)
plot(x,pdf,type="l",lty=1,xlab="",ylab="",xaxt="n",yaxt="n",bty="n",ylim=c(0,0.1))
abline(h=0)
x = seq(0,qchisq(0.025,10),0.001)
y = dchisq(x,10)
polygon(c(0,x,qchisq(0.025,10)),c(0,y,0),col="lightblue")
points(x=qchisq(0.025,10),y=0,pch=19,cex=1.5)
x = seq(qchisq(0.975,10),30,0.001)
y = dchisq(x,10)
polygon(c(qchisq(0.975,10),x,30),c(0,y,0),col="lightblue")
points(x=qchisq(0.975,10),y=0,pch=19,cex=1.5)
text(9.8,0.02,expression(1-alpha),cex=1.25)
text(0.5,0.01,expression(alpha/2),cex=1.25)
text(25,0.0075,expression(alpha/2),cex=1.25)


# Figure 10.16
# Page 32
# Left
boxplot(pre_men,post_men,
 	xlab="",names=c("Pre-menarche","Post-menarche"),ylab="Age at menarche (in years)",ylim=c(12,18),col="lightblue")
t.test(pre_men,post_men,conf.level=0.95,var.equal=TRUE,alternative="greater")
# Right
x = seq(-10,10,0.001)
pdf = dt(x,150)
plot(x,pdf,type="l",lty=1,xlab="",ylab="",ylim=c(0,0.4))
abline(h=0)
points(x=7.0583,y=0,pch=19,cex=1.5)

# Swimmer data 
pre_men = c(15.0, 17.1, 14.6, 15.2, 14.9, 14.4, 14.7, 15.3, 13.6, 15.1, 16.2, 15.9, 13.8, 15.0, 15.4, 14.9,
14.2, 16.5, 13.2, 16.8, 15.3, 14.7, 13.9, 16.1, 15.4, 14.6, 15.2, 14.8, 13.7, 16.3, 15.1, 14.5,
16.4, 13.6, 14.8, 15.5, 13.9, 15.9, 16.0, 14.6, 14.0, 15.1, 14.8, 15.0, 14.8, 15.3, 15.7, 14.3,
13.9, 15.6, 15.4, 14.6, 15.2, 14.8, 13.7, 16.3, 15.1, 14.5, 13.6, 15.1, 16.2, 15.9, 13.8, 15.0,
15.4, 14.9, 16.2, 15.9, 13.8, 15.0, 15.4, 14.9, 14.2, 16.5, 13.4, 16.5, 14.8, 15.1, 14.9, 13.7,
16.2, 15.8, 15.4, 14.7, 14.3, 15.2, 14.6, 13.7, 14.9, 15.8, 15.1, 14.6, 13.8, 16.0, 15.0, 14.6)

post_men = c(14.0, 16.1, 13.4, 14.6, 13.7, 13.2, 13.7, 14.3, 12.9, 14.1, 15.1, 14.8, 12.8, 14.2, 14.1, 13.6,
14.2, 15.8, 12.7, 15.6, 14.1, 13.0, 12.9, 15.1, 15.0, 13.6, 14.2, 13.8, 12.7, 15.3, 14.1, 13.5,
15.3, 12.6, 13.8, 14.4, 12.9, 14.6, 15.0, 13.8, 13.0, 14.1, 13.8, 14.2, 13.6, 14.1, 14.5, 13.1, 
12.8, 14.3, 14.2, 13.5, 14.1, 13.6, 12.4, 15.1)


# Figure 10.17
# Page 34
hist(p.value,xlab="p-value",ylab="",main="",col="lightblue",freq=FALSE,yaxp=c(0,1,1),ylim=c(0,1.25))
abline(h=1,lty=2)

# Monte Carlo simulation
n = 10 # sample size
B = 10000 # number of Monte Carlo samples
data = matrix(rnorm(n*B,0,1),nrow=B,ncol=n)
p.value = data[,1]*0 # initialise
for (i in 1:B){
	p.value[i] = t.test(data[i,],conf.level=0.95)$p.value
	}


# Figure 10.18
# Page 36
theta = seq(2,10,0.001)
theta.0=6
n=10
alpha=0.05
sigma.0=2
par(mar = c(4, 4.4, 4, 2)) # adjust plotting region to avoid cutting off labels
plot(theta,
	1-pnorm((theta.0+qnorm(1-alpha/2)*(sigma.0/sqrt(n))-theta)/(sigma.0/sqrt(n))) + 
	pnorm((theta.0-qnorm(1-alpha/2)*(sigma.0/sqrt(n))-theta)/(sigma.0/sqrt(n))),
		type="l",lty=1,xlab=expression(theta),ylab=expression(K(theta)),
		xaxt="n",yaxt="n",cex.lab=1.25,ylim=c(0,1))
abline(h=alpha,lty=2)
axis(side=1,at=c(6),labels=c(expression(theta[0])),cex.axis=1.25)
axis(side=2,at=c(0,alpha,1),labels=c(0,expression(alpha),1),cex.axis=1.25,las=2)


# Figure 10.19
# Page 37
# Left
mu = seq(40,70,0.001)
mu.0=50
n=5
alpha=0.05
sigma.0=10
plot(mu,1-pnorm((qnorm(1-alpha)*(sigma.0/sqrt(n))+mu.0-mu)/(sigma.0/sqrt(n))),
	type="l",lty=1,xlab=expression(mu),ylab=expression(K(mu)),xlim=c(40,70),ylim=c(0,1),cex.lab=1.5)
abline(h=alpha,lty=2)
# Right
plot(mu,1-pnorm((qnorm(1-alpha)*(sigma.0/sqrt(n))+mu.0-mu)/(sigma.0/sqrt(n))),
	type="l",lty=1,xlab=expression(mu),ylab=expression(K(mu)),xlim=c(40,70),ylim=c(0,1),cex.lab=1.5)
abline(h=alpha,lty=2)
lines(mu,1-pnorm((qnorm(1-alpha)*(sigma.0/sqrt(15))+mu.0-mu)/(sigma.0/sqrt(15))),lty=4)
lines(mu,1-pnorm((qnorm(1-alpha)*(sigma.0/sqrt(30))+mu.0-mu)/(sigma.0/sqrt(30))),lty=8)
# Add legend
legend(40,1,lty=c(1,4,8), 
	c(expression(paste(n,"=5")),expression(paste(n,"=15")),expression(paste(n,"=30"))))


# Figure 10.20
# Page 38
alpha = 0.0115
theta = seq(0.75,2,0.001)
plot(theta,1-ppois(105,84*theta),type="l",lty=1,xlab=expression(theta),ylab=expression(K(theta)),xlim=c(0.75,2),ylim=c(0,1),cex.lab=1.25)
abline(h=alpha,lty=2)


# Figure 10.21
# Page 41
y = seq(0,1,0.01)
plot(y,dbeta(y,1,2),type="l",lty=1,xlab="y",ylab="Probability density function",xaxt="n",yaxp=c(0,4,4),,ylim=c(0,4),cex.lab=1.25)
lines(y,dbeta(y,1,4),lty=4)
abline(h=0)
abline(v=0,lty=2)
axis(side=1,at=c(0,1/2,1),labels=c("0","1/2","1"))
legend(0.5,3,lty=c(1,4),c(expression(paste(theta,"=2")),expression(paste(theta,"=4"))))


# Figure 10.22
# Page 44
x = seq(0,30,0.001)
pdf = dchisq(x,10)
plot(x,pdf,type="l",lty=1,xlab="",ylab="",xaxt="n",yaxt="n",bty="n",ylim=c(0,0.1))
abline(h=0)
x = seq(0,qchisq(0.05,10),0.001)
y = dchisq(x,10)
polygon(c(0,x,qchisq(0.05,10)),c(0,y,0),col="lightblue")
points(x=qchisq(0.05,10),y=0,pch=19,cex=1.5)
text(9.8,0.02,expression(1-alpha),cex=1.25)
text(0.5,0.01,expression(alpha),cex=1.25)


# Figure 10.23
# Page 47
theta = seq(5,10,0.001)
theta.0=6
n=15
alpha=0.05
sigma.0=2
# Left
plot(theta,1-pnorm((theta.0+qnorm(1-alpha)*(sigma.0/sqrt(n))-theta)/(sigma.0/sqrt(n))),
	type="l",lty=1,xlab=expression(theta),ylab="Power function",
	xaxt="n",yaxt="n",cex.lab=1.5)
lines(theta,1-pnorm((qnorm(1-alpha)*(sigma.0/sqrt(6))+theta.0-theta)/(sigma.0/sqrt(6))),lty=4) #n=6
abline(h=alpha,lty=2)
axis(side=1,at=c(6),labels=c(expression(theta[0])),cex.axis=1.5)
axis(side=2,at=c(0,alpha,1),labels=c(0,expression(alpha),1),cex.axis=1.5,las=2)
text(7,0.82,expression(K(theta)),cex=1.5)
text(8,0.62,expression(paste(K,"*",(theta),sep="")),cex=1.5)
# Right
theta.0=9
plot(theta,pnorm((theta.0+qnorm(alpha)*(sigma.0/sqrt(n))-theta)/(sigma.0/sqrt(n))),
	type="l",lty=1,xlab=expression(theta),ylab="Power function",
	xaxt="n",yaxt="n",cex.lab=1.5)
lines(theta,pnorm((qnorm(alpha)*(sigma.0/sqrt(6))+theta.0-theta)/(sigma.0/sqrt(6))),lty=4) #n=6
abline(h=alpha,lty=2)
axis(side=1,at=c(9),labels=c(expression(theta[0])),cex.axis=1.5)
axis(side=2,at=c(0,alpha,1),labels=c(0,expression(alpha),1),cex.axis=1.5,las=2)
text(8.05,0.82,expression(K(theta)),cex=1.5)
text(7,0.62,expression(paste(K,"*",(theta),sep="")),cex=1.5)


# Figure 10.24
# Page 51
par(mar = c(4, 4.4, 4, 2)) # adjust plotting region to avoid cutting off labels
alpha = 0.0443
p = seq(0.4,1,0.001)
plot(p,1-pbinom(58,100,p),
	type="l",lty=1,xlab=expression(p),ylab=expression(K(p)),xlim=c(0.4,1),ylim=c(0,1),cex.lab=1.25)
abline(h=alpha,lty=2)


# Figure 10.25
# Page 53
x = seq(0,30,0.001)
pdf = dchisq(x,10)
plot(x,pdf,type="l",lty=1,xlab="",ylab="",xaxt="n",yaxt="n",bty="n",ylim=c(0,0.1))
abline(h=0)
x = seq(0,qchisq(0.05,10),0.001)
y = dchisq(x,10)
polygon(c(0,x,qchisq(0.05,10)),c(0,y,0),col="lightblue")
points(x=qchisq(0.05,10),y=0,pch=19,cex=1.5)
text(9.8,0.02,expression(1-alpha),cex=1.25)
text(0.5,0.01,expression(alpha),cex=1.25)


# Figure 10.26
# Page 54
par(mar = c(4, 4.7, 4, 2)) # adjust plotting region to avoid cutting off labels
sigma.sq = seq(0,12,0.001)
sigma.sq.0=9
n=20
alpha=0.05
plot(sigma.sq,pchisq((sigma.sq.0*qchisq(alpha,n)/sigma.sq),n),
	type="l",lty=1,xlab=expression(sigma^2),ylab=expression(K(sigma^2)),cex.lab=1.25)
abline(h=alpha,lty=2)


# Figure 10.27
# Page 59
y = seq(0,80,0.1)
pdf = dgamma(y,4,1/6)
plot(y,pdf,type="l",xlab="t",ylab="g(t)",xaxt="n",yaxt="n",cex.lab=1.25)
abline(h=0)
abline(h=0.01,lty=2)
axis(side=1,at=c(6,41),labels=c(expression(k[1]),expression(k[2])),cex.axis=1.25)
axis(side=2,at=c(0.01),labels=c(expression(k)),cex.axis=1.25,las=2)
points(x=6,y=0,pch=19,cex=1.25)
points(x=41,y=0,pch=19,cex=1.25)


# Figure 10.28
# Page 60
par(mar = c(4, 4.4, 4, 2)) # adjust plotting region to avoid cutting off labels
theta = seq(0,50,0.01)
theta.0=20
n=30
alpha=0.05
plot(theta,pchisq(theta.0*qchisq(alpha/2,2*n)/theta,2*n)+1-pchisq(theta.0*qchisq(1-alpha/2,2*n)/theta,2*n),
	type="l",lty=1,xlab=expression(theta),ylab=expression(K(theta)),xlim=c(0,50),ylim=c(0,1),cex.lab=1.5)
abline(h=alpha,lty=2)