n <- 100
mu_true <- 20
sigma2_true <- 10
# Data generation:
set.seed(42)
y <- rnorm(n, mu_true, sqrt(sigma2_true))
alpha <- .05
mean(y) + c(-1,1)*qt(1-alpha/2, n-1)*sqrt(var(y)/n)
2.5 % 97.5 %
(Intercept) 19.44941 20.75624
################################################
alpha <- 1-c(.8, .85, .9, .95, .99)
CI_alpha <- matrix(NA, ncol=2, nrow=length(alpha))
for(i in 1:length(alpha)){
CI_alpha[i,] <- mean(y) + c(-1,1)*qt(1-alpha[i]/2, n-1)*sqrt(var(y)/n)
}
plot(CI_alpha[,1], (1-alpha), pch=19, xlim=c(18.5,21.5), cex=.5, xlab="mu")
points(CI_alpha[,2], (1-alpha), pch=19, cex=.5)
segments(x0=CI_alpha[,1], y0=(1-alpha), x1=CI_alpha[,2], y1=(1-alpha))
abline(v=mu_true, col=2, lwd=2)
abline(v=mean(y), lty=2)
y_10 <- rnorm(10, mu_true, sqrt(sigma2_true))
y_500 <- rnorm(500, mu_true, sqrt(sigma2_true))
CI_alpha_10 <- matrix(NA, ncol=2, nrow=length(alpha))
CI_alpha_500 <- matrix(NA, ncol=2, nrow=length(alpha))
for(i in 1:length(alpha)){
CI_alpha_500[i,] <- mean(y_500) + c(-1,1)*qt(1-alpha[i]/2, n-1)*sqrt(var(y_500)/500)
CI_alpha_10[i,] <- mean(y_10) + c(-1,1)*qt(1-alpha[i]/2, n-1)*sqrt(var(y_10)/10)
}
par(mfrow=c(1,3))
plot(CI_alpha_10[,1], (1-alpha), pch=19, xlim=c(18.5,21.5), cex=.5, xlab="mu", main="n=10")
points(CI_alpha_10[,2], (1-alpha), pch=19, cex=.5)
segments(x0=CI_alpha_10[,1], y0=(1-alpha), x1=CI_alpha_10[,2], y1=(1-alpha))
abline(v=mu_true, col=2)
plot(CI_alpha[,1], (1-alpha), pch=19, xlim=c(18.5,21.5), cex=.5, xlab="mu", main="n=100")
points(CI_alpha[,2], (1-alpha), pch=19, cex=.5)
segments(x0=CI_alpha[,1], y0=(1-alpha), x1=CI_alpha[,2], y1=(1-alpha))
abline(v=mu_true, col=2)
plot(CI_alpha_500[,1], (1-alpha), pch=19, xlim=c(18.5,21.5), cex=.5, xlab="mu", main="n=500")
points(CI_alpha_500[,2], (1-alpha), pch=19, cex=.5)
segments(x0=CI_alpha_500[,1], y0=(1-alpha), x1=CI_alpha_500[,2], y1=(1-alpha))
abline(v=mu_true, col=2)
par(mfrow=c(1,1))
################################################
nsim <- 5000
n <- 100
alpha <- 0.05
coverage <- numeric(nsim)
CI <- matrix(NA, ncol=2, nrow=nsim)
set.seed(396)
for(i in 1:nsim){
y <- rnorm(n, mu_true, sqrt(sigma2_true))
CI[i,] <- mean(y) + c(-1,1)*qt(1-alpha/2, n-1)*sqrt(var(y)/n)
coverage[i] <- CI[i,1] <= mu_true & CI[i,2] >= mu_true
}
table(coverage)/nsim
coverage
0 1
0.0492 0.9508
plot(CI[1:100,1], 1:100, pch=19, xlim=range(CI), cex=.5, col=as.factor(coverage+1), xlab="mu")
points(CI[1:100,2], 1:100, pch=19, cex=.5, col=as.factor(coverage+1))
#lines(x=CI[1:2,1], y=CI[1:2,2])
segments(x0=CI[1:100,1], y0=1:100, x1=CI[1:100,2], y1=1:100, col=as.factor(coverage+1))
abline(v=mu_true)
