Coverage probability for Wald confidence interval with small sample size

a) My understanding is that the sample variance (i.e. the squared deviation divided by n - 1 instead of n) constitutes a mean-unbiased estimator of population variance. However, I've run multiple simulations and it appears that sample variance is a median-biased estimator of population variance, especially for small sample sizes. Is this right?

For the simulations I've run, the average sample variance is quite close to the population variance, but the median sample variance is lower than the population variance, which means that in more cases than not the sample variance underestimates the population variance. For some simulations this is a difference of two or even three percent.

Here's an example in R:

set.seed(33)

distribution <- rnorm(10000, 4, 8)

set.seed(33)

data.frame(x <- sapply(1:10000, function(i){ var(sample(distribution, 32)) > (var(distribution) * (10000 - 1)/10000) })) %>% filter(x == "TRUE") %>% dim()

In this case, the sample variance is larger than the variance 46.72 percent of the time. I mention this median-bias because it seems that this could effect confidence intervals.

b) For the simulations I've run, coverage probability for small sample sizes, even when sampling from a normal distribution, is smaller than nominal coverage probability (as calculated by a Wald confidence interval), often by around 1 percent or more.

Here's an example:

set.seed(33)

x <- lapply(1:100000, function(i){ sample(distribution, 32, replace = TRUE) })

bind_cols(data.frame(unlist(lapply(x, mean))), data.frame(unlist(lapply(x, >sd)))) %>% rename(average = unlist.lapply.x..mean.., stan_dev = unlist.lapply.x..sd..) >%>% mutate(high = average + 1.96 * stan_dev/sqrt(32), low = average - 1.96 * >stan_dev/sqrt(32)) %>% filter(high > mean(distribution)) %>% filter(low < mean(distribution)) %>% dim()

Here, the coverage probability is only 94.167 percent.

I understand that sample standard deviation (sample variance squared) is a (slightly) mean-biased (?) estimator of population standard deviation. Is the coverage probability above related to this or to the median-bias of sample variance.

I recognize that there are significant coverage problems with the Wald confidence interval for the binomial distribution (see https://projecteuclid.org/euclid.ss/1009213286), Poisson distribution, etc. I didn't realize that this was the case even for the normal distribution.

Any help in understanding the above would be much appreciated. If I've simply made a coding error, please do point this out.

Otherwise, could someone please suggest a better confidence interval than the Wald for normal and other continuous distributions with a small sample size and/or refer me to any relevant literature?

Much appreciated.

EDITED: For clarity and brevity.

First, for normal data $$S = \frac{1}{n-1}\sum_{i=1}^n (X_i - \bar X)^2$$ is (mean) biased for $$\sigma,$$ especially so for small $$n.$$ The unbiasing constant depends on $$n.$$ For your sample size $$n = 32$$ my simulation below in R, illustrates this. With a million iterations it is reasonable to expect about 2-place accuracy.

set.seed(1123)
s = replicate(10^6, sd(rnorm(32,4,8)))
mean(s)
[1] 7.93691      # biased: < 8
mean(s)/(sqrt(2/31)*gamma(32/2)/gamma(31/2))
[1] 8.001164     # unbiased: aprx 8


The sample variance $$S^2$$ has median below $$\sigma^2.$$ This is not surprising because $$\frac{(n-1)S^2}{\sigma^2} \sim \mathsf{Chisq}(n-1),$$ a right-skewed distribution.

median(s^2)
[1] 62.64106     # < 64
mean(s^2 < 64)
[1] 0.533473     # > 1/2


I am not sure what formula you are using for your 95% confidence interval. The commonly used CI for normal $$\mu$$ is based on $$S$$ as defined above and Student's t distribution with $$n - 1$$ degrees of freedom: $$\bar X \pm t^*S/\sqrt{n},$$ where $$t^*$$ cuts probability 2.5 from the upper tail of $$\mathsf{T}(\nu = n-1).$$ A sample size $$n = 32$$ is not quite large enough to substitute the normal quantile 1.96 for $$t^*.$$ R code continues from above:

set.seed(1123)
a = replicate(10^6, mean(rnorm(32, 4, 8)))
mean(a);  median(a)
[1] 3.999486         # aprx 4
[1] 3.998304         # aprx 4

qt = qt(.975, 31);  qt
[1]  2.039513        # not 1.96
LCL = a - qt*s/sqrt(32)
UCL = a + qt*s/sqrt(32)
mean((LCL <= 4) &  (UCL >= 4))
[1] 0.94955          # aprx 0.95


I am not sure I follow the syntax of your uncommented R code, but it seems you may be using 1.96 and a version of the sample variance with denominator $$n$$ instead of $$n-1.$$ If that is the case, then it should be no surprise that your coverage probability is not exactly 95%.

• Regarding, coverage probability, apologies for the unclear code. The function "sd" in R takes the variance with denominator n - 1, but you're right that I took the normal quantile and if this was replaced with the t-statistic, coverage probability would be very close to 95 percent. Thanks! – num_39 Nov 23 '19 at 18:55
• Sample variance with $n-1$ in denominator is (mean) unbiased estimator of population variance, but its square root is biased for population SD. Unbiasedness does not 'survive' nonlinear operations. // Many Wald CIs are asymptotic: work well if $n$ is sufficiently large that estimator is almost unbiased and very nearly normal. Try your example with $n = 1000$ instead of $n = 32$ and I suppose you'll get coverage probability nearer 95%. – BruceET Nov 23 '19 at 20:06