Techniques/diagnostics for gaining confidence in normality assumptions and resulting confidence intervals

Question

I have data that is reasonably assumed to be iid samples from some distribution. Our goal is to put a confidence interval on the population mean and have something similar for the population variance. Notationally, we have IID $X_i, i = 1, ..., n$ with mean $\mu$, variance $\sigma^2$ unknown. Sample size $n$ varies from 200 to 20,000.

Plotting my data, it is trimodal, so definitely does not seem to be coming from a normal distribution.

Computing confidence intervals on the sample mean $\bar{X_n}$ is no problem. I'm just not sure how much to trust them.

Below is my planned diagnostics. Can you tell me if it make sense?

Here my reasoning is to see if many samples of $\bar{X_{n-1}}$ does indeed follow a normal to gain confidence that our $n$ is large enough that the central limit theorem approximation is valid. I can compute the sample mean $\bar{X_{n-1}}$ when I hold out one data point. That gives $n$ different samples of sample means $\bar{X_{n-1}}$. I can plot those and see if they follow a normal, or perform a Q-Q plot against a normal with standard deviation the computes standard error ($S_n/\sqrt{n}$) to see if they are close to normal. If so, then it is evidence that $\bar{X_n}$ will be normal and the confidence interval computed with variance $S_n^2/n$ is valid. Does this sound like a sound method?

Are there other diagnostics to gain confidence in the confidence intervals that are perhaps better?

---

related question: Can I use sub-sample means to verify that the sample mean is approximately normal?

Answer 1

Leaving one out is a creative solution, but unfortunately, it cannot help you to guess the form of the distribution of $\bar{X}$, because you will generate values that are equal to $(n\bar{X}-X_i)/(n-1)$, essentially a constant, $n\cdot \bar{X}$, minus the original sample, the whole thing scaled by $1/(n-1)$. So it will still look trimodal.

The idea to check artificially generated outputs of $\bar{X}$ is however good!

You can do something very closely related: simple bootstrapping, i.e. resampling your sample with replacement.

With such large sample sizes as you have, you should in most cases be on the safe side, unless you can spot that the data produce extreme outliers (looking at a boxplot can give an idea).

If you have R, you can experiment with following R code to get an impression:

# function to simulate from a nasty trimodal distribution
rtrimod <- function(n) rlnorm(n) + 10 * sample(0:2, size = n, replace=TRUE)
  
# sample size: 200

n <- 200

x <- rtrimod(n)
hist(x, breaks = 50)

# first we simulate the real sampling distribution

nsimul <- 10000
xbar <- numeric(nsimul) # reserve space
for (i in 1:nsimul) 
  xbar[i] <- mean(rtrimod(n))

hist(xbar, breaks = 50)
qqnorm(xbar, main = "qqnorm of sample means")

# techniques to try out
# leave one out:

xminus <- numeric(n) 
for (i in 1:n)
  xminus[i] <- mean(x[-i])

hist(xminus, breaks = 50) 
qqnorm(xminus) 

# now bootstrap

nboot <- 10000
xboot <- numeric(nboot) 

for (i in 1:nboot)
  xboot[i] <- mean(sample(x, replace = TRUE))

hist(xboot, breaks = 50) 
qqnorm(xboot, main = "qqnorm of bootstrap")

p.s. Pythonists: in R, qqnorm produces a Q-Q plot against a standard normal. This does not matter, since the qqnorm plot of any other large enough normal sample will look nicely linear, just not like the identity line, but with abscissa $\mu$ and slope $\sigma$.

p.p.s: Automatically translated code for python (not checked):

import numpy as np
import scipy.stats as stats

def rtrimod(n):
    return np.random.lognormal(size=n) + 10 * np.random.choice([0, 1, 2], size=n)

n = 200
x = rtrimod(n)
_ = plt.hist(x, bins=50)

nsimul = 10000
xbar = np.zeros(nsimul)
for i in range(nsimul):
    xbar[i] = np.mean(rtrimod(n))
_ = plt.hist(xbar, bins=50)
_ = stats.probplot(xbar, plot=plt)

xminus = np.zeros(n)
for i in range(n):
    xminus[i] = np.mean(np.delete(x, i))
_ = plt.hist(xminus, bins=50)
_ = stats.probplot(xminus, plot=plt)

nboot = 10000
xboot = np.zeros(nboot)
for i in range(nboot):
    xboot[i] = np.mean(np.random.choice(x, size=n, replace=True))
_ = plt.hist(xboot, bins=50)
_ = stats.probplot(xboot, plot=plt)


# or a qqplot: 
def qqplot(Y):   
    Y = sorted(Y)
    N = len(Y)
    x = [norm.ppf(i/(N+1)) for i in range(1, N+1)] 
    ## to compare to y = x: 
    line_values = np.arange(x[0], x[-1], .1)
    plt.plot(line_values, line_values, color = 'r')

    # plot our percentiles against each other: 
    plt.scatter(x, Y,  marker = "o")
    plt.ylabel("Sample's Emperical Quantiles" )
    plt.xlabel("Standard Normal's Quantiles")
    plt.show()

qqplot(xboot)