Question

1. Let W1, . . . , Wn be i.i.d. N(µ, σ²). A 100(1 − α)% confidence interval for σ² is

( ([n − 1]*S² ) / (χ²_{1−α/2,n−1} ) , ([n − 1]S²) / ( χ²_α/2,n−1) ,

where χ²_u,k denotes the 100×uth percentile of the χ² distribution with k degrees of freedom, S² = (n − 1)⁻¹ * the sum of (Wi − Wbar )² (from 1 to n) is the sample variance and Wbar = n⁻¹ * the sum of Wi (from i = 1 to n) the sample mean.

(a) Suppose µ = 68, σ = 4, n = 12, and α = 0.05. Estimate the coverage probability of this random confidence interval for σ² using a Monte Carlo simulation with 1e4 replications. That is, estimate the probability that the confidence interval captures σ² . Denote this estimated probability by pbar. Comment on how 1 − α compares to ±1/ √ 1e4. You may use the qchisq function (hint: the χ² distribution is not symmetric) but you may not use the built in rnorm function.

(b) Estimate the coverage probability for the confidence interval when W1, . . . , Wn are i.i.d. Exp(1). Notice that in this case σ² = var(W1) = 1. Estimate the coverage probability pbar for each (n, α) ∈ {10, 50, 100} × {0.05, 0.1}. In each case, compare 1 − α to pbar ± 1/ √ 1e4. Again, you may not use the built in rexp function. Comment on the results.

Answer 1

Here W are IID normal with . Interval estimate of a population variance:

is sample size, is sample variance.

Here

a) The coverage probability is estimated below using Simulation R code below.

m <- 10000
n <- 12
sigma <- 4
mu <- 68
alpha <- 0.05
#X <- array (dim = m)
n1 <- 0
for (i in 1:m)
{
W <- rnorm(n, mean = mu, sd = sigma)
s <- sd(W)
if(s^2*(n-1)/qchisq(alpha/2,n-1) > sigma^2 & s^2*(n-1)/qchisq(1-alpha/2,n-1) < sigma^2)
{
n1 = n1 + 1
}
  
}
pbar <- n1/m
pbar

The estimated coverage probability is which is close to theoretical .

b) The cocverage probability when is calulated below.

m <- 10000
n <- 12
sigma <- 1
alpha <- 0.05
#X <- array (dim = m)
n1 <- 0
for (i in 1:m)
{
W <- rexp(n, 1)
s <- sd(W)
if(s^2*(n-1)/qchisq(alpha/2,n-1) > sigma^2 & s^2*(n-1)/qchisq(1-alpha/2,n-1) < sigma^2)
{
n1 = n1 + 1
}
  
}
pbar <- n1/m
pbar