模拟和可视化中心极限定理

模拟和可视化中心极限定理
Simulating and Visualising the Central Limit Theorem

原始链接: https://blog.foletta.net/post/2025-07-14-clt/

## 通过模拟探索中心极限定理这篇文章详细记录了作者深入理解中心极限定理（CLT）的个人历程——一个在计算机科学学位期间最初被回避的概念，但现在却觉得非常有趣。作者使用R中的模拟和可视化来超越理论理解，并亲眼看到中心极限定理的作用。中心极限定理指出，来自*任何*总体样本均值的分布，随着样本量的增加，将趋近于正态分布。为了验证这一点，作者从六种不同的分布（均匀分布、正态分布、二项分布、贝塔分布、指数分布和卡方分布）生成数据，并反复计算样本均值。初步模拟结果显示，样本量为6时，样本均值的分布*倾向于*正态分布，但偏态分布需要更大的样本量才能收敛。重要的是，使用正态分布来计算小样本的置信区间被证明是不正确的——需要使用t分布。将样本量增加到60可以提高准确性，但偏态分布仍然滞后。进一步探索，以2的幂次方增加样本量，展示了不同分布收敛到正态性的速度。均匀分布迅速正态化，而指数分布需要显著更大的样本量。这种实践方法提供了对中心极限定理的优势和局限性的实用理解，突出了它在各种总体分布中的行为。

I completed a Computer Science degree at uni, and bundled a lot of maths subjects in as electives: partial differential equations, vector calculus, discrete maths, linear algebra. For some reason however I always avoided statistics subjects. Maybe there’s a story to be told about a young person finding uncertainty uncomfortable, because twenty years later I find statistics, particularly the Bayesian flavour, really interesting.

One problem with a self-directed journey is that there’s foundational knowledge that has come to me in dribs and drabs, and one of the most foundational is the Central Limit Theorem (CLT). In this post I want to interrogate and explore the CLT using simulation and visualisation in an attempt to understand how it works in practice, not in theory. This is predominantly a process to help me better understand the CLT; you’re just here for the ride. Hopefully that ride can help you get where you need to go as well.

It’s been awhile since I’ve included any code in a post, so where it makes sense I’ll show the generating R code, with a liberal sprinking of comments so it’s hopefully not too inscrutable.

I don’t want this to be like an online recipe with pages of back story before you get to the meat and bones, but a brief summary of the CLT before we begin is unavoidable. In plain English the CLT can be described as such:

“If you take repeated samples of size n from a distribution and calculate the sample mean for each, as n gets approaches infinity, the distribution of sample means approaches a normal distribution.”

For the classic CLT there’s a couple of assumptions about the source distribution:

The sample is drawn independently (no autocorrelation like in a time series).
All the data points are drawn from the same distribution (“independent and identically distributed” or i.i.d).
The distribution has a finite mean and variance (e.g. no Cauchy or Pareto distributions).

There are other versions of the CLT in which some of these assumptions are relaxed, but we’ll focus on the ‘classic’ version.

Putting it math terms:

$$ \frac{\bar{X}_n - \mu}{\sigma/\sqrt{n}} \overset{d}\longrightarrow \mathcal{N}(0, 1) $$

We’re not the kind of people to just accept something because someone threw some fancy Greek symbols at us. We want to simulate it to give ourselves some confidence it works in practice.

Let’s create a tibble of ten-thousand random values from six different distributions, which we’ll call our ‘population’:

# 10,000 draws from six different distributions
population_data <- tibble(
    uniform = runif(10000, min = -20, max = 20),
    normal = rnorm(10000, mean = 0, sd = 4),
    binomial = rbinom(10000, size = 1, prob = .5),
    beta = rbeta(10000, shape1 = .9, shape2 = .5),
    exponential = rexp(10000, .4),
    chisquare = rchisq(10000, df = 2),
)

uniform	normal	binomial	beta	exponential	chisquare
Six 'Population' Distributions - Ten-Thousand Values - First Five Rows
-15.243338	-0.7022216	0	0.9473270	2.34764246	2.8571390
-9.788476	-3.2383413	1	0.0730441	6.00295709	0.9272733
-11.432831	-2.0977844	0	0.5160286	0.91033272	3.1607328
18.142428	2.5360707	1	0.9428468	3.08170008	1.9489974
-7.078033	-1.7224620	0	0.6134417	0.08687766	2.5317882

This ‘wide’ tibble is good for sampling from, but we’ll also transform it into a long version which will have other uses (note the ’_l’ suffix).

# Long version of random data
population_data_l <-
    population_data |>
    pivot_longer(cols = everything(), names_to = 'distribution')

distribution	value
Six Distributions - Post 'pivot_longer() - First Value of Each
beta	0.9473270
binomial	0.0000000
chisquare	2.8571390
exponential	2.3476425
normal	-0.7022216
uniform	-15.2433376

Here’s a histogram of each of the population distributions:

Let’s define a function take_random_sample_mean() which takes a sample from all of the population distributions and calculates the mean. If we use this function repeatedly, we should end up with a data set that demonstrates the central limit theorem.

Let’s take 20,000 sample means of size 60, bind it all together into a single data frame, and shape it into a long version.

# Define a function to take a random sample from our data
take_random_sample_mean <- function(data, sample_size) {
    slice_sample(.data = data, n = sample_size) |>
    summarise(across(everything(), list(sample_mean = mean, sample_sd = sd)))
}

# Draw 20,000 means of size 60 from our random data
sample_size <- 60
sample_means <- 
    map(1:20000, ~take_random_sample_mean(population_data, sample_size = sample_size)) |> 
    # Bind the sample means into a single tibble
    list_rbind() |> 
    # Move to a long version of the data
    pivot_longer(
        everything(),
        names_to = c('distribution', '.value'),
        names_pattern = '(\\w+)_(\\w+_\\w+)'
    )

Here’s the resulting histograms of the samples with only the x-axis free to change scale. The beta and binomial look pretty normal, and the others might be as well, but the differences in variance make it difficult to tell.

If you recall the formula at the start, to get to a standard normal we also need to subtract the population mean and divide by the population standard deviation over the square-root of n. Using the long version of the population data we can calculate the population mean and SD statistics:

population_data_stats <-
    population_data_l |>
    group_by(distribution) |>
    summarise(
        population_mean = mean(value),
        population_sd = sd(value),
        n = n()
    )

distribution	population_mean	population_sd	n
Random Means - Statistics
beta	0.64352237	0.3105243	10000
binomial	0.50260000	0.5000182	10000
chisquare	2.03125146	2.0189741	10000
exponential	2.48144359	2.5049145	10000
normal	-0.03205698	3.9813502	10000
uniform	-0.09148261	11.6405771	10000

Joining each of those statistics into the sample means by their distribution allows us to nicely standardise.

# CLT Calculation
clt <-
    sample_means |>
    # Join in the popultion mean, sd, and n
    left_join(population_data_stats, by = 'distribution') |>
    # Scale to the standard normal
    mutate(
        clt = (sample_mean - population_mean) / (population_sd / sqrt(sample_size) )
    )

That’s better: they now at least all look the same, except for the binomial which ends up having higher counts because of its discrete rather than continuous values. While you can kind of guess that they’re normal from a histogram, we can get a better sense quantile-quantile plot. The standard normal quantiles are on the x-axis, and our sample mean quantiles are on the y-axis.

They all track pretty closely to a standard normal, however the chi-square and exponential do tend to diverge slightly at the ends. We’ll dive a bit deeper into that later in the post.

Let’s summarise: we took twenty-thousand sample means of size sixty from six wildly different distributions, and we were able to see that the distributions of these sample means approximately followed a normal distribution. This is the essence of the central limit theorem, which allows us to use statistical methods for normal distributions on problems that involve population distributions that aren’t normal.

That’s all well and good if you’ve got the population mean and standard deviation, but in most cases you’re not going to have that. You’re also likely not going to have the resources to take twenty-thousand different samples.

Let’s take a a pretty classic example of interviewing six people about about something (weight, height, voting intentions, etc) and take the average: your sample mean. What the CLT can give you is an interval around your sample mean that would give you the a confidence interval (we’ll use the classic 95%) that the true mean of the population is somewhere in the interval.

So you use a bit of algebra and move some terms around in the CLT formula and you get this:

$$ \bar{X} \pm z_.025 \cdot \ \frac{s}{\sqrt{n}}
$$

where $\bar{X}$ is your sample mean, $z_.025$ is the critical value (aka qnorm() in R) and $s$ is the sample standard deviation. The sample standard deviation over the square root n is more commonly known as the standard error.

More simply, we take the sample mean and add/subtract the 97.5 quantile from a normal distribution times the standard error to get our 95% interval. Some may have noted a small mistake I’ve made here, which I’’ll leave in but will soon come to light.

We’re in frequentist territory here, so while it’s tempting to say “there’s a 95% probability that the true mean is in the interval”, we shouldn’t. Why? Because to a frequentist, the true mean $\mu$ is a fixed value: it’s either in the interval or it’s not, we can’t assign a probability to the population mean. What we should say is that if we were to repeat the sampling process many times, in the long run we should see the true mean within this confidence interval 95% of the time. We can simulate this to see if that is true.

We use a sample size of 6, taking the mean of these 6 samples from each of our population distributions 10,000 times.

# Ten-thousand random sample means of size six
small_sample_size <- 6 

repeated_samples <-
    map(1:10000, ~take_random_sample_mean(population_data, sample_size = small_sample_size)) |>
    # Bind the samples into a single tibble
    list_rbind() |> 
    # Wide to long tibble
    pivot_longer(
        everything(),
        names_to = c('distribution', '.value'),
        names_pattern = '(\\w+)_(\\w+_\\w+)'
    ) |>
    # Join with our population mean 
    left_join(population_data_stats, by = 'distribution')

For each of the samples from the distributions we calculate the 95% confidence interval. Because we know the true population mean, we can determine whether it is or is not within the interval. Then we calculate the percentage of samples for which the population mean fell within the 95% interval.

conf_intervals <-
    repeated_samples |>
    # Calculate CIs and whether CI contains the true population mean
    mutate(
        upper_ci = sample_mean + qnorm(0.975) * sample_sd / sqrt(small_sample_size),
        lower_ci = sample_mean - qnorm(0.975) * sample_sd / sqrt(small_sample_size),
        within_ci = population_mean < upper_ci & population_mean > lower_ci
    ) |> 
    # Determine percentage of CI that contain the true mean
    group_by(distribution) |> 
    summarise(percent_within_ci = mean(within_ci))

Distribution	% Within CI
normal	89.36%
uniform	88.95%
beta	87.87%
chisquare	83.32%
exponential	82.87%
binomial	77.90%

Uh oh! We’re way off our 95% here, what happened?

This is the mistake I was talking about earlier: we used a normal to calculate the CIs. We estimated the population standard deviation $\sigma$ using our sample standard deviation $s$. With our small sample size of 6, our CLT formula follows a t-distribution, not a normal.

We’ll re-run our confidence intervals, but this time we use qt(), the t-distribution quantile function instead of qnorm().

conf_intervals <-
    repeated_samples |> 
    mutate(
        upper_ci = sample_mean + qt(p = 0.975, df = small_sample_size - 1) * (sample_sd / sqrt(small_sample_size)),
        lower_ci = sample_mean - qt(p = 0.975, df = small_sample_size - 1) * (sample_sd / sqrt(small_sample_size)),
        within_ci = population_mean < upper_ci & population_mean > lower_ci
    ) |> 
    group_by(distribution) |> 
    summarise(percent_within_ci = mean(within_ci))

Distribution	% Within CI
binomial	96.83%
normal	95.05%
uniform	93.88%
beta	92.59%
chisquare	88.91%
exponential	88.47%

That looks a bit better! The binomial is higher than 95% due to its discrete nature (can’t smoothly hit all possible values), the normal and uniform distributions are close to our 95% value, but our heavily skewed beta, exponential a chi-square still aren’t up to scratch. With those heavily skewed population distributions, we need more samples before the central limit theorem ‘kicks in’.

Let’s re-run using the dataset from the start of the post, which used s sample size of 60.

conf_intervals <-
    # Using our original sample size of 60
    sample_means |> 
    left_join(population_data_stats, by = 'distribution') |>
    mutate(
        upper_ci = sample_mean + qnorm(0.975) * (sample_sd / sqrt(sample_size)),
        lower_ci = sample_mean - qnorm(0.975) * (sample_sd / sqrt(sample_size)),
        within_ci = population_mean < upper_ci & population_mean > lower_ci
    ) |> 
    group_by(distribution) |> 
    summarise(percent_within_ci = mean(within_ci))

Distribution	% Within CI
binomial	94.86%
normal	94.52%
beta	94.40%
uniform	94.36%
chisquare	93.60%
exponential	92.97%

Better, but even with a sample size of 60, the skewed distirbutions are still slightly under the mark.

This begs the question: how do the sample means of these skewed distributions behave as the sample size increased? What we can do is repeatedly sample, but increase the sample size each time. We’ll go up in powers of two from 1 to 1024 samples.

increasing_sample_size <-
    # Sample sizes increasing in powers of 2 (1, 2, 4, 8, ...)
    map(
        2^(0:10),
        # Anonymous function
        \(y) {
            # 1000 sample means
            map(1:1000, ~take_random_sample_mean(population_data, sample_size = y)) |>
            # Bind them all together
            list_rbind() |> 
            # Wide to long per distribution
            pivot_longer(
                everything(),
                names_to = c('distribution', '.value'),
                names_pattern = '(\\w+)_(\\w+_\\w+)'
            ) |>
            # Add in our population means and SDs
            left_join(population_data_stats, by = 'distribution') |>
            # Add sample size as a column
            mutate(sample_size = y) |> 
            # Normalise to -> N(0,1)
            mutate(clt = (sample_mean - population_mean) / (population_sd / sqrt(sample_size) )) 
        }
    ) |>
    list_rbind()

With this data we can create an animation of a Q-Q plot for the uniform versus exponential distributions, showing how the distribution of sample means changes as the sample size increases. You’ll see the distribution approaching the standard normal distribution.

ou’ll see the distribution of sample means from the uniform distribution approaches a normal much faster than the expontential. It’s very subjective, but I think the uniform stsrts looking reasonably good at a sample size of 8. The exponential however takes much longer to converge to a normal.

The central limit theorem is something I’ve read about many times and had a reasonable grasp on, but I was always wondering about how it behaved with different population distributions. Running simulations and visualising how it behaved in different scenarios has given me (and hopefully yourselves) a much clearer view on how it works, and probably more importantly where it doesn’t work well.

模拟和可视化中心极限定理 Simulating and Visualising the Central Limit Theorem

模拟和可视化中心极限定理
Simulating and Visualising the Central Limit Theorem