2
$\begingroup$

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?

$\endgroup$

1 Answer 1

1
$\begingroup$

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)
$\endgroup$
14
  • $\begingroup$ awesome explanation. I'll try it out. $\endgroup$ Jun 5, 2023 at 14:46
  • $\begingroup$ Why are you suggesting bootstrapping with replacement, as opposed to sampling without replacement? One could consider my original suggestion (leave one out) as sampling without replacement. I'm guessing this help with the "near constant" potential problems you suggested. $\endgroup$ Jun 5, 2023 at 14:48
  • 1
    $\begingroup$ Tx :-) - Did you try it out? your suggestion sounds very reasonable, much more systematic. But it is in a way "too systematic" - it just reproduces exactly what you see in your sample. CLT concerns all possible samples. Bootstrap tries to reflect all possible samples as far as possible. This implies in particular, that the items are sampled independently, which is not the case when you leave one out, but it is the case, when you sample with replacement. (To be precise: conditionally independent, given the observed sample). $\endgroup$
    – Ute
    Jun 5, 2023 at 14:53
  • $\begingroup$ Ok, just tried something similar. It basically shrinks the original distribution in the x-direction but reproduces it--still tri modal. Makes sense that these $\bar{X}_{n-1}(i) $ are not IID, as all pairs share $n-1$ variables. $\endgroup$ Jun 19, 2023 at 22:03
  • 1
    $\begingroup$ size has to be equal to the original sample size. It is a feature that some values are multiply included. $\endgroup$
    – Ute
    Jun 20, 2023 at 0:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.