title: “Sampling Distributions and Central Limit Theorem Analysis” output: html_document —
First, we load the necessary packages and dataset. We treat the
NCbirths dataset as our population and visualize the
distribution of birth weights (BirthWeightGm).
# Load required libraries
library(Stat2Data)
library(dplyr)
library(ggplot2)
# Load dataset
data("NCbirths")
# Plot histogram of the population variable
ggplot(NCbirths, aes(x = BirthWeightGm)) +
geom_histogram(binwidth = 250, fill = "skyblue", color = "black") +
labs(
title = "Population Distribution of Birth Weights",
x = "Birth Weight (grams)",
y = "Count"
) +
theme_minimal()
# Calculate and store population parameters
mu <- mean(NCbirths$BirthWeightGm, na.rm = TRUE)
sigma <- sd(NCbirths$BirthWeightGm, na.rm = TRUE)
mu
## [1] 3295.619
sigma
## [1] 632.977
Interpretation: The population distribution of birth weights is unimodal and slightly left-skewed. The population mean (\(\mu\)) is approximately 3295.62 grams, and the population standard deviation (\(\sigma\)) is approximately 632.98 grams.
Next, we draw 1,000 random samples of size \(n = 30\) with replacement from the population and plot the sampling distribution of the sample means.
# Set sample size and random seed
n <- 30
set.seed(123)
# Draw 1,000 samples of size n=30
sample_means <- replicate(1000, {
sample_data <- sample(NCbirths$BirthWeightGm[!is.na(NCbirths$BirthWeightGm)], size = n, replace = TRUE)
mean(sample_data)
})
# Data frame for plotting
sim_data <- data.frame(sample_means = sample_means)
# Plot histogram of sample means
ggplot(sim_data, aes(x = sample_means)) +
geom_histogram(binwidth = 25, fill = "lightgreen", color = "black") +
labs(
title = "Sampling Distribution of Sample Means (n = 30)",
x = "Sample Mean Birth Weight (grams)",
y = "Count"
) +
theme_minimal()
# Summary statistics of simulated sample means
mean_sample_means <- mean(sample_means)
sd_sample_means <- sd(sample_means)
mean_sample_means
## [1] 3302.995
sd_sample_means
## [1] 114.0927
Interpretation: The sampling distribution of sample means is unimodal and symmetric, forming a bell curve. Compared to the population distribution, the center remains essentially the same (\(\text{mean} \approx\) 3303 grams vs. \(\mu \approx\) 3295.62 grams), but the spread is significantly narrower (\(\text{SD} \approx\) 114.09 grams vs. \(\sigma \approx\) 632.98 grams).
According to the Central Limit Theorem (CLT), the standard error is \(\text{SE} = \frac{\sigma}{\sqrt{n}}\). We compare the theoretical normal model with our simulation.
# Compute Standard Error
SE <- sigma / sqrt(n)
SE
## [1] 115.5653
# Density histogram with superimposed CLT normal curve
ggplot(sim_data, aes(x = sample_means)) +
geom_histogram(aes(y = after_stat(density)), binwidth = 25, fill = "lightgreen", color = "black") +
stat_function(
fun = dnorm,
args = list(mean = mu, sd = SE),
color = "red",
linewidth = 1
) +
labs(
title = "Simulated Sample Means with CLT Normal Overlay",
x = "Sample Mean Birth Weight (grams)",
y = "Density"
) +
theme_minimal()
Interpretation: The theoretical CLT model matches the simulated sampling distribution extremely well. The theoretical mean (\(\mu =\) 3295.62) is virtually identical to the simulated mean (3303), and the theoretical standard error (\(\text{SE} =\) 115.57) closely matches the simulated standard deviation (114.09). The red normal curve overlays the green histogram seamlessly.
pnorm()We compute the theoretical probability that a random sample of \(n = 30\) babies has a mean birth weight
greater than \(3,600\text{ grams}\)
using pnorm(), and compare it to the proportion from our
simulation.
# Theoretical probability using CLT
prob_clt <- pnorm(3600, mean = mu, sd = SE, lower.tail = FALSE)
prob_clt
## [1] 0.004221207
# Empirical proportion from simulation
prop_sim <- mean(sample_means > 3600)
prop_sim
## [1] 0.003
Interpretation: The theoretical probability from the CLT model is approximately 0.0042, and the empirical proportion from the simulation is 0.003. These values are nearly identical because 1,000 random samples provide a very close approximation to the theoretical normal curve, with minor differences due to random sampling variability in the simulation.
qnorm()We find the cutoff value above which \(90\%\) of sample means fall (the 10th
percentile) using qnorm(), and compare it with the
empirical 10th percentile from our simulation.
# Theoretical cutoff using CLT (10th percentile)
cutoff_clt <- qnorm(0.10, mean = mu, sd = SE)
cutoff_clt
## [1] 3147.516
# Empirical 10th percentile from simulation
cutoff_sim <- quantile(sample_means, 0.10)
cutoff_sim
## 10%
## 3150.536
Interpretation: The theoretical 10th percentile from
qnorm() is 3147.5 grams, and the empirical 10th percentile
from the simulation is 3150.5 grams. In context, \(90\%\) of random samples of 30 babies will
have a mean birth weight above approximately \(3,334.1\text{ grams}\).