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. 
 A: 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%.
