Central-limit theorem via sample size or sampling magnitude?

Question

I have a small application of the central limit theorem:

"Take a sample of 40 exponentially distributed random numbers and calculate their mean. Repeat 1000 times, record each measured mean and calculate the mean of those 1000 means."

Here's my Rcode:

nsim   = 1000
nsamp  = 40
lambda = 0.2

clt <- data.frame(m=numeric(nsim), s=numeric(nsim))

for(s in 1:nsim)
{
    clt$v[s] <- mean(rexp(nsamp, lambda))
    clt$s[s] <- sd(rexp(nsamp, lambda))
}   

cat(mean(clt$v), "\n")
    cat(mean(clt$s), "\n")

I can observe that with increasing number of simulations (i.e. increasing the repetitions of the sampling event), the mean of means and mean of variance more closely agree with the theoretical values (1/lambda).

What I wonder is, does the Central-limit theorem only depend on the number of sampling events (1000) or is it just as much dependent on sufficiently large sample size (i.e. 40 in this example)?

Answer 1

The CLT is generally considered a theoretical distribution. We don't actually generate 1000 (or 10,000 or 100,000) samples from it. Rather, we imagine what would happen if we were to replicate this experiment 1000 (or, technically, infinite) times.

The R code you post caps it at 1000, probably because by 1000 samples it will demonstrate the properties of the CLT. But, again, the CLT really says what would happen if you repeated it an infinite number of times.

So, to your question:

does the Central-limit theorem only depend on the number of sampling events (1000) or is it just as much dependent on sufficiently large sample size (i.e. 40 in this example)?

It depends on both (kind of). The number of sampling events actually is infinite (though the CLT properties hold with smaller distributions), and as the sample size increases, the more normal (and smaller the variance) the distribution of means becomes.

Answer 2

The classical central limit theorem says that, under certain conditions, some random variable $Z$ that depends on a parameter $n$ will approach an $N(0, 1)$ distribution as $n\rightarrow\infty$. If you sample from the distribution $Z$ has, say, $R$-many times, then the Glivenko-Cantelli theorem says that the empirical distribution of those draws will approach the true distribution of $Z$ as $R\rightarrow\infty$.

If you want to simulate this, you have to pay attention to the differences between $N$ and $R$. In particular, $R$ does not appear in the central limit theorem.

Let's look at the CLT function I wrote.

clt_exp <- function(N, R, rate = 1){
  
  # N: Samples drawn from the exponential distribution
  # R: Times to sample from the "Z" distribution
  # rate: Parameter of the exponential distribution
  
  # Define the population mean and standard deviation
  #
  the_mean <- 1/rate
  the_sd <- 1/rate
  
  z <- rep(NA, R)
  
  # Take R-many draws from the distribution that is 
  # asymptotically N(0, 1)
  #
  for (i in 1:R){
    
    # Simulate N-many draws from an exponential distribution
    #
    x <- rexp(N, rate)
    
    # Transform the mean of x as specified by the CLT
    #
    z[i] <- (mean(x) - the_mean)/the_sd * sqrt(N)
  }
  
  return(z)
  
}

What this basically defines is something like rcltexp, where the R input is the number of samples to draw like it is in rexp or rnorm, and then N and rate are parameters of the distribution. When we run a function like rexp, we don't expect rexp(3) or even rexp(10) to look particularly exponential. However, we do expect rexp(99999) to look rather exponential. This is Glivenko-Cantelli in action.

Next, let's vary N and R to see what happens. Let's start with N = 4, so that the mean $\bar x$ is calculated as the mean of four exponentials. Throughout, we will keep the rate = 1.

library(ggplot2)
set.seed(2024)

clt_exp <- function(N, R, rate = 1){
  
  # N: Samples drawn from the exponential distribution
  # R: Times to sample from the "Z" distribution
  # rate: Parameter of the exponential distribution
  
  # Define the population mean and standard deviation
  #
  the_mean <- 1/rate
  the_sd <- 1/rate
  
  z <- rep(NA, R)
  
  # Take R-many draws from the distribution that is 
  # asymptotically N(0, 1)
  #
  for (i in 1:R){
    
    # Simulate N-many draws from an exponential distribution
    #
    x <- rexp(N, rate)
    
    # Transform the mean of x as specified by the CLT
    #
    z[i] <- (mean(x) - the_mean)/the_sd * sqrt(N)
  }
  
  return(z)
  
}

z_4_10 <- clt_exp(4, 10)
d_4_10 <- data.frame(
  z = z_4_10,
  ECDF = ecdf(z_4_10)(z_4_10),
  N = "N = 4",
  R = "R = 10"
)
z_4_100 <- clt_exp(4, 100)
d_4_100 <- data.frame(
  z = z_4_100,
  ECDF = ecdf(z_4_100)(z_4_100),
  N = "N = 4",
  R = "R = 100"
)
z_4_1000 <- clt_exp(4, 1000)
d_4_1000 <- data.frame(
  z = z_4_1000,
  ECDF = ecdf(z_4_1000)(z_4_1000),
  N = "N = 4",
  R = "R = 1000"
)
z_4_10000 <- clt_exp(4, 10000)
d_4_10000 <- data.frame(
  z = z_4_10000,
  ECDF = ecdf(z_4_10000)(z_4_10000),
  N = "N = 4",
  R = "R = 10000"
)
d_4 <- rbind(
  d_4_10,
  d_4_100,
  d_4_1000,
  d_4_10000
)
ggplot(d_4, aes(x = z, y = ECDF)) +
  geom_point() +
  geom_line() +
  facet_grid(rows = vars(R))

With just ten observations, we do not get much resolution about the distribution. With R = 100, we begin to see what the CDF looks like, and R = 1000 and certainly R = 10000 give good resolution. However, with N = 4, the distribution is fairly obviously non-normal. However, as Glivenko-Cantelli says should happen, whatever the distribution is, as R gets large, the observations beging to match that distribution.

Next, let's up N to 20.

z_20_10 <- clt_exp(20, 10)
d_20_10 <- data.frame(
  z = z_20_10,
  ECDF = ecdf(z_20_10)(z_20_10),
  N = "N = 20",
  R = "R = 10"
)


z_20_100 <- clt_exp(20, 100)
d_20_100 <- data.frame(
  z = z_20_100,
  ECDF = ecdf(z_20_100)(z_20_100),
  N = "N = 20",
  R = "R = 100"

)

z_20_1000 <- clt_exp(20, 1000)
d_20_1000 <- data.frame(
  z = z_20_1000,
  ECDF = ecdf(z_20_1000)(z_20_1000),
  N = "N = 20",
  R = "R = 1000"
)

z_20_10000 <- clt_exp(20, 10000)
d_20_10000 <- data.frame(
  z = z_20_10000,
  ECDF = ecdf(z_20_10000)(z_20_10000),
  N = "N = 20",
  R = "R = 10000"
)

d_20 <- rbind(
  d_20_10,
  d_20_100,
  d_20_1000,
  d_20_10000
)
ggplot(d_20, aes(x = z, y = ECDF)) +
  geom_point() +
  geom_line() +
  facet_grid(rows = vars(R))

With N = 20, the distribution looks more normal, not perfect, but much more normal. Let's see if a much larger N will give something that looks really Gaussian. Let's try N = 5000 with the same R values.

z_5000_10 <- clt_exp(5000, 10)
d_5000_10 <- data.frame(
  z = z_5000_10,
  ECDF = ecdf(z_5000_10)(z_5000_10),
  N = "N = 5000",
  R = "R = 10"
)


z_5000_100 <- clt_exp(5000, 100)
d_5000_100 <- data.frame(
  z = z_5000_100,
  ECDF = ecdf(z_5000_100)(z_5000_100),
  N = "N = 5000",
  R = "R = 100"
)


z_5000_1000 <- clt_exp(5000, 1000)
d_5000_1000 <- data.frame(
  z = z_5000_1000,
  ECDF = ecdf(z_5000_1000)(z_5000_1000),
  N = "N = 5000",
  R = "R = 1000"
)


z_5000_10000 <- clt_exp(5000, 10000)
d_5000_10000 <- data.frame(
  z = z_5000_10000,
  ECDF = ecdf(z_5000_10000)(z_5000_10000),
  N = "N = 5000",
  R = "R = 10000"
)


d_5000 <- rbind(
  d_5000_10,
  d_5000_100,
  d_5000_1000,
  d_5000_10000
)
ggplot(d_5000, aes(x = z, y = ECDF)) +
  geom_point() +
  geom_line() +
  facet_grid(rows = vars(R))

Well doesn't that look Gaussian?

Overall, if you want to simulate the central limit theorem, make the R a huge number and vary the N. The N is a parameter of the $Z = \sqrt{N}\dfrac{\bar X_N -\mu}{\sigma}$. The R is just how many times you draw from the distribution of $Z$. Let's look at an example for a fixed, large R of 10000.

set.seed(2024)
s <- seq(-7, 7, 0.01)
z_1_10000 <- clt_exp(1, 10000)
d_1 <- data.frame(
  z = s,
  ECDF = ecdf(z_1_10000)(s),
  N = 1
)
#
z_2_10000 <- clt_exp(2, 10000)
d_2 <- data.frame(
  z = s,
  ECDF = ecdf(z_2_10000)(s),
  N = 2
)
#
z_5_10000 <- clt_exp(5, 10000)
d_5 <- data.frame(
  z = s,
  ECDF = ecdf(z_5_10000)(s),
  N = 5
)
#
z_9_10000 <- clt_exp(9, 10000)
d_9 <- data.frame(
  z = s,
  ECDF = ecdf(z_9_10000)(s),
  N = 9
)
#
z_99_10000 <- clt_exp(99, 10000)
d_99 <- data.frame(
  z = s,
  ECDF = ecdf(z_99_10000)(s),
  N = 99
)
#
z_999_10000 <- clt_exp(999, 10000)
d_999 <- data.frame(
  z = s,
  ECDF = ecdf(z_999_10000)(s),
  N = 999
)
#
d_n01 <- data.frame(
  z = s,
  ECDF = pnorm(s, 0, 1),
  N = "N(0, 1)"
)
#
d <- rbind(
  d_1,
  d_2,
  d_5,
  d_9,
  d_99,
  d_999,
  d_n01
)
ggplot(d, aes(x = z, y = ECDF, col = N)) +
  geom_line() +
  facet_grid(rows = vars(N)) +
  theme(legend.position = "bottom")

Look at that nice convergence toward the $N(0, 1)$.

Compare that to if R = 50.

set.seed(2024)
s <- seq(-7, 7, 0.01)
z_1_50 <- clt_exp(1, 50)
d_1 <- data.frame(
  z = s,
  ECDF = ecdf(z_1_50)(s),
  N = 1
)
#
z_2_50 <- clt_exp(2, 50)
d_2 <- data.frame(
  z = s,
  ECDF = ecdf(z_2_50)(s),
  N = 2
)
#
z_5_50 <- clt_exp(5, 50)
d_5 <- data.frame(
  z = s,
  ECDF = ecdf(z_5_50)(s),
  N = 5
)
#
z_9_50 <- clt_exp(9, 50)
d_9 <- data.frame(
  z = s,
  ECDF = ecdf(z_9_50)(s),
  N = 9
)
#
z_99_50 <- clt_exp(99, 50)
d_99 <- data.frame(
  z = s,
  ECDF = ecdf(z_99_50)(s),
  N = 99
)
#
z_999_50 <- clt_exp(999, 50)
d_999 <- data.frame(
  z = s,
  ECDF = ecdf(z_999_50)(s),
  N = 999
)
#
d_n01 <- data.frame(
  z = s,
  ECDF = pnorm(s, 0, 1),
  N = "N(0, 1)"
)
#
d <- rbind(
  d_1,
  d_2,
  d_5,
  d_9,
  d_99,
  d_999,
  d_n01
)
ggplot(d, aes(x = z, y = ECDF, col = N)) +
  geom_line() +
  facet_grid(rows = vars(N)) +
  theme(legend.position = "bottom")

The story about the convergence is the same, just with rougher graphs.

Answer 3

The central limit theorem states that you should take samples of sufficiently large size (40 in your example), but it does not provide a limit on the minimum number of such samples, i.e., the statement holds regardless of the number of samples. See point #1 in this answer.

Answer 4

It does not matter for law of large numbers

I observe that with increasing number of simulations (i.e. increasing the repetitions of the sampling event), the mean of means and mean of variance more closely agree with the theoretical values (1/lambda).

Your application and observation of the behaviour of the mean relates more to the law of large numbers than the central limit theorem.

For the mean of means, it is irrelevant how you create the average. If you average all your samples or first average groups does not matter, because the end result is the same.

Example

Say your sample is $1,2,3,4$ then you could group it and compute average of group averages, but the result is the same (if the group's sizes are equal).

$$\frac{(2+3)/2 + (4+1)/2}{2} = \frac{2+3+4+1}{4}$$

---

Also note that the mean of the estimated standard deviations may approach some value for larger nsim, but that doesn't need to relate to the simulations following approximately a normal distribution.

The behaviour that this value approaches some constant is due to the law of large numbers.

It matters for central limit theorem

The CLT relates to the distribution of a sample average. When you standardize this based on the population mean and deviation, then the final result approaches a Normal distribution under the right circumstances (independence, finite variance, etc.).

If your question would have been about the central limit theorem or about approximations with a normal distribution then the size of the samples matters (not the number of repetitions, although this may have an effect when you make a histogram).

There is already a question about this: Why does increasing the sample size of coin flips not improve the normal curve approximation?

It deals with the issue of the two different sizes (the sample size and the repetitions size), as well as with problems of drawing a histogram.