theoretical confidence interval depending on sample size [closed]

Question

I am using R and plain English to express my question. Let us say I have a "true"/made up population, which is normally distributed with a mean of 500000 and a standard deviation of 13000:

mean <- 500000
sd <- 13000
population <- rnorm(10000000, mean=mean, sd=sd)

I can then repeatedly sample from this population using different samples sizes ranging from 3 to 150. For each sample, I also calculate the 95% confidence interval using a t distribution. I think this is safe as at the beginning my sample size is pretty low (i.e. 3). This is not really exploited here but I think this would be correct?

results <- NULL

for (sample_size_ in 3:150) {
    
    for (sample_number_ in 1:10) {
        
        sample_ <- sample(population, sample_size_)
        ci_test <- t.test(sample_, conf.level=0.95)

        df <- data.frame(
            sample_size = sample_size_
            , low_ci = ci_test$conf.int[[1]] 
            , up_ci = ci_test$conf.int[[2]]
            , mean = mean(sample_)
        )

        results <- rbind(results, df)
        
    }
}

I can then plot the data frame, which nicely shows how the sample means gets closer to the population mean (red dotted line) as the sample size increases:

ggplot(results, aes(y = mean, x = sample_size)) +
    geom_point() +
    geom_hline(yintercept=mean, linetype="dotted", color = "red", size=3)

This is the output (the actual code uses a theme):

I guess this is straightforward but how can I determine these blue lines expressing the 95% theoretical uncertainty of the sample mean being inside the blue area?

Is there a formula and does it depend on the knowledge of that the true distribution is normal? I guess it should not, as the sample means will be normally distributed according to the CLT?

Thanks.

PS:

This is the current/final code:

mu_ <- 500000
sd_ <- 13000
z_value <- qnorm(.975)
max_sample_size = 500
repeats = 2

results <- NULL
upper_and_lower <- NULL

for (sample_size_ in 3:max_sample_size) {
    
    for (sample_number_ in 1:repeats) {
        
        sample_ <- rnorm(sample_size_, mean=mu_, sd=sd_)
        ci_test <- t.test(sample_, conf.level=0.95)

        df <- data.frame(
            sample_size = sample_size_
            , low_ci = ci_test$conf.int[[1]] 
            , up_ci = ci_test$conf.int[[2]]
            , mean = mean(sample_)
        )

        results <- rbind(results, df)
        
    }
    
    df <- data.frame(
        sample_size = sample_size_
        , upper = mu_ + (z_value * (sd_/sqrt(sample_size_)))
        , lower = mu_ - (z_value * (sd_/sqrt(sample_size_)))
    )
    upper_and_lower <- rbind(upper_and_lower, df)
}

ggplot() +
    geom_point(data=results, aes(y = mean, x = sample_size)) +
    geom_line(data=upper_and_lower, aes(y = upper, x = sample_size), color = "blue", size=2) +
    geom_line(data=upper_and_lower, aes(y = lower, x = sample_size), color = "blue", size=2)

Resulting graph with theoretical 95% uncertainty band:

Answer 1

Nice experiment. The blue lines will be at $\mu \pm z_{\alpha/2} \sigma/\sqrt{n}$ where $\alpha = 0.05$ and $\alpha \mapsto z_\alpha$ is the upper quantile function of the standard normal and $n$ is sample_size_. In this case, this is exact since you are generating from the normal distribution and the sample mean will have the distribution $N(\mu, \sigma^2/n)$. In general, it holds approximately based on the CLT for large $n$.

If you are trying to interpret the CI, that would be a bit different. Putting a CI around each of the points vertically, for each sample_size_, roughly 95% of those intervals will cross the dotted red line (you cover the true mean $\approx$ 95% of the time, hence the coverage probability is 0.95). Since your replication size is small (10), this is not that accurate. Try increasing sample_number_ to say 1000 to see it better. What you will observe is that although the length of the CIs will be smaller as you increase $n$, still $\approx$ 95% of them keep covering the red line no matter what $n$ is.

---

PS. I am assuming that instead of these two lines:

population <- rnorm(10000000, mean=mean, sd=sd)
sample_ <- sample(population, sample_size_)

you would do something like this:

sample_ <- rnorm(sample_size_, mean=mean, sd=sd)

That is, sample from the normal distribution, not a finite population drawn from the normal distribution as you are doing here, although the results will be close in your case (unless your sample size starts approaching 10000000).