set.seed(123)

par(mfrow=c(1,3))

CIlower_90<-numeric(10000)
CIupper_90<-numeric(10000)
pvalue1_90<-numeric(10000)
pvalue2_90<-numeric(10000)

for(j in 1:10000) {
  (sample_1<-rnorm(100,10,2))
  (s2_1<-var(sample_1))
  (chi_1<-(99*s2_1)/4)
  (pvalue1_90[j]<-pchisq(chi_1, 0.95, df=99, lower.tail=F))
  (pvalue2_90[j]<-pchisq(chi_1, 0.05, df=99, lower.tail=T))
  (chi_u_1<-qchisq(0.95, df=99))
  (chi_l_1<-qchisq(0.05, df=99))
  (CIlower_90[j]<-(99*s2_1)/chi_u_1)
  (CIupper_90[j]<-(99*s2_1)/chi_l_1)
}

#CIlower_90
#CIupper_90
#pvalue1_90
#pvalue2_90

reject1_90<-pvalue1_90 <= 0.05 | pvalue2_90 <= 0.05
table(reject1_90)

color<-rep(gray(.7),100)
color[reject1_90[1:100]]<-"blue"

plot(0, xlim=c(1,7), ylim=c(1,100), ylab="Sample No.", xlab="", main="90% Confidence Interval for sigma-square")
abline(v=4, lty=2)
for(j in 1:10000) {
  lines(c(CIlower_90[j], CIupper_90[j]), c(j,j), col=color[j], lwd=1)
}

set.seed(123)

CIlower_95<-numeric(10000)
CIupper_95<-numeric(10000)
pvalue1_95<-numeric(10000)
pvalue2_95<-numeric(10000)

for(j in 1:10000) {
  (sample_2<-rnorm(100,10,2))
  (s2_2<-var(sample_2))
  (chi_2<-(99*s2_2)/4)
  (pvalue1_95[j]<-pchisq(chi_2, 0.975, df=99, lower.tail=F))
  (pvalue2_95[j]<-pchisq(chi_2, 0.025, df=99, lower.tail=T))
  (chi_u_2<-qchisq(0.975, df=99))
  (chi_l_2<-qchisq(0.025, df=99))
  (CIlower_95[j]<-(99*s2_2)/chi_u_2)
  (CIupper_95[j]<-(99*s2_2)/chi_l_2)
}

#CIlower_95
#CIupper_95
#pvalue1_95
#pvalue2_95

reject1_95<-pvalue1_95 <= 0.025 | pvalue2_95 <= 0.025
table(reject1_95)

color<-rep(gray(.7),100)
color[reject1_95[1:100]]<-"black"

plot(0, xlim=c(1,7), ylim=c(1,100), ylab="Sample No.", xlab="", main="95% Confidence Interval for sigma-square")
abline(v=4, lty=2)
for(j in 1:10000) {
  lines(c(CIlower_95[j], CIupper_95[j]), c(j,j), col=color[j], lwd=1)
}

set.seed(123)

CIlower_99<-numeric(10000)
CIupper_99<-numeric(10000)
pvalue1_99<-numeric(10000)
pvalue2_99<-numeric(10000)

for(j in 1:10000) {
  (sample_3<-rnorm(100,10,2))
  (s2_3<-var(sample_3))
  (chi_3<-(99*s2_3)/4)
  (pvalue1_99[j]<-pchisq(chi_3, 0.995, df=99, lower.tail=F))
  (pvalue2_99[j]<-pchisq(chi_3, 0.005, df=99, lower.tail=T))
  (chi_u_3<-qchisq(0.995, df=99))
  (chi_l_3<-qchisq(0.005, df=99))
  (CIlower_99[j]<-(99*s2_3)/chi_u_3)
  (CIupper_99[j]<-(99*s2_3)/chi_l_3)
}

#CIlower_99
#CIupper_99
#pvalue1_99
#pvalue2_99

reject1_99<-pvalue1_99 <= 0.005 | pvalue2_99 <= 0.005
table(reject1_99)

color<-rep(gray(.7),100)
color[reject1_99[1:100]]<-"dark red"

plot(0, xlim=c(1,7), ylim=c(1,100), ylab="Sample No.", xlab="", main="99% Confidence Interval for sigma-square")
abline(v=4, lty=2)
for(j in 1:10000) {
  lines(c(CIlower_99[j], CIupper_99[j]), c(j,j), col=color[j], lwd=1)
}