STM1001: Computer Lab 5B Solutions

Data Science Stream

Topic 5B: Simulations in R

Example R code solutions for the Data Science Computer Lab 5, which uses the R packages dplyr (Wickham et al. 2021) and ggplot2 (Wickham et al. 2022), and penguin data from the palmerpenguins R package (Horst, Hill, and Gorman 2020), are presented below.

1 Sampling Distributions

1.1 Sampling

No answer required.

1.2 The Central Limit Theorem

No answer required.

2 Simulations in R

2.1

set.seed(1)
n <- 100
x <- rnorm(n)

2.2

hist(x, freq = FALSE, col = "chartreuse3", xlim = c(-3, 3),
     main = paste("Histogram of random sample \n from standard Normal distribution, n = ", n)) 
# Overlay standard normal curve
curve(dnorm(x), add = TRUE, col = "blue", lwd = 2)
# Add line denoting sample mean
abline(v = round(mean(x), 3), col = "red", lwd = 2, lty = 2)

2.3

No answer required.

2.4

No answer required.

2.5 for loops

No answer required.

2.6

#Specify sample size
ns <- 5

# Specify number of times to conduct process  
trials <- 10000

# create object in which to store sample means once they are generated
norm.means <- rep(0, trials) # Note length is equal to the number of trials

# Run for loop  
for(i in 1:trials){
    # Randomly generate ns samples from the normal distribution
    x <- rnorm(ns) 
    # Compute mean of samples, store in the norm.means object, in ith position
    norm.means[i] <- mean(x)}

2.7

hist(norm.means, freq = FALSE, breaks = 20, col = "red", 
     xlim = c(-3, 3),
     xlab = expression(bar(x)),
     main = paste("Histogram of means, n = ", ns)) 

curve(dnorm(x, 
      mean = 0, sd = 1/sqrt(ns)), 
      add = TRUE, col = "blue", lwd = 2) 

The simulated sample means are centred on the population mean of 0, but do display some variance.

2.8

#Specify sample size
ns <- 30
  
# Specify number of times to conduct process  
trials <- 10000

# create object in which to store sample means once they are generated
norm.means <- rep(0, trials) # Note length is equal to the number of trials

# Run for loop  
for(i in 1:trials){
    # Randomly generate ns samples from the normal distribution
    x <- rnorm(ns) 
    # Compute mean of samples, store in the norm.means object, in ith position
    norm.means[i] <- mean(x)}

hist(norm.means, freq = FALSE, breaks = 20, col = "red", 
     xlim = c(-3, 3),
     xlab = expression(bar(x)),
     main = paste("Histogram of means, n = ", ns)) 
curve(dnorm(x, 
      mean = 0, sd = 1/sqrt(ns)), 
      add = TRUE, col = "blue", lwd = 2) 

As we increase the sample size, the variance of the sample means decreases.

2.9

#Specify sample size
ns <- 60
  
# Specify number of times to conduct process  
trials <- 10000

# create object in which to store sample means once they are generated
norm.means <- rep(0, trials) # Note length is equal to the number of trials

# Run for loop  
for(i in 1:trials){
    # Randomly generate ns samples from the normal distribution
    x <- rnorm(ns) 
    # Compute mean of samples, store in the norm.means object, in ith position 
    norm.means[i] <- mean(x)}

hist(norm.means, freq = FALSE, breaks = 20, col = "red", 
     xlim = c(-3, 3),
     xlab = expression(bar(x)),
     main = paste("Histogram of means, n = ", ns)) 
curve(dnorm(x, 
      mean = 0, sd = 1/sqrt(ns)), 
      add = TRUE, col = "blue", lwd = 2) 

For a sample size of 60, the variance of the sample means is very low, providing us with an accurate estimate of the population mean.

3 CLT Simulation - Exponential Distribution

3.1

set.seed(1)
n <- 100
rate = 10
x <- rexp(n, rate = rate)

3.2

hist(x, freq = FALSE, col = "chartreuse3", xlim = c(0, 0.5),
     main = paste("Histogram of random sample \n from Exponential distribution, n = ", n))

curve(dnorm(x, 1/rate, 1/rate), add = TRUE, col = "blue", lwd = 2)

3.3 The Rate Parameter and Population Mean

No answer required.

3.4

#Specify sample size
ns <- 5
  
# Specify rate
rate = 10

# Specify number of times to conduct process  
trials <- 10000

# create object in which to store sample means once they are generated
norm.means <- rep(0, trials) # Note length is equal to the number of trials

# Run for loop  
for(i in 1:trials){
    # Randomly generate ns samples from the exponential distribution
    x <- rexp(ns, rate = rate)
    # Compute mean of samples, and store it in the norm.means object, in ith position 
    norm.means[i] <- mean(x) 
    }

hist(norm.means, freq = FALSE, breaks = 20, col = "red", 
     xlim = c(0, 0.4),
     xlab = expression(bar(x)),
     main = paste("Histogram of means, n = ", ns)) 

curve(dnorm(x, mean = 1/rate, sd = 1/rate/sqrt(ns)), add = TRUE,
        col = "blue", lwd = 2)

Here we use an x-axis range of 0 to 0.4.

3.5

The resulting data set is not as skewed as the underlying distribution, but there is still a noticeable skew (this is because of the low sample size).

3.6

#Specify sample size
ns <- 30
  
# Specify rate
rate = 10

# Specify number of times to conduct process  
trials <- 10000

# create object in which to store sample means once they are generated
norm.means <- rep(0, trials) # Note length is equal to the number of trials

# Run for loop  
for(i in 1:trials){
    # Randomly generate ns samples from the exponential distribution
    x <- rexp(ns, rate = rate)
    # Compute mean of samples, and store it in the norm.means object, in ith position 
    norm.means[i] <- mean(x) 
    }

hist(norm.means, freq = FALSE, breaks = 20, col = "red", 
     xlim = c(0, 0.2),
     xlab = expression(bar(x)),
     main = paste("Histogram of means, n = ", ns)) 

curve(dnorm(x, mean = 1/rate, sd = 1/rate/sqrt(ns)), add = TRUE,
        col = "blue", lwd = 2)

Here, we change the x axis range to be from 0 to 0.2. The sample means now display very little variance, and provide a much more accurate estimate of the population mean of 0.1. Note that the Normal density curve also fits the data well. These results are due to the sample size being bigger (30 instead of 5).

3.7

#Specify sample size
ns <- 60
  
# Specify rate
rate = 10

# Specify number of times to conduct process  
trials <- 10000

# create object in which to store sample means once they are generated
norm.means <- rep(0, trials) # Note length is equal to the number of trials

# Run for loop  
for(i in 1:trials){
    # Randomly generate ns samples from the exponential distribution
    x <- rexp(ns, rate = rate)
    # Compute mean of samples, and store it in the norm.means object, in ith position 
    norm.means[i] <- mean(x) 
    }

hist(norm.means, freq = FALSE, breaks = 20, col = "red", 
     xlim = c(0, 0.2),
     xlab = expression(bar(x)),
     main = paste("Histogram of means, n = ", ns)) 

curve(dnorm(x, mean = 1/rate, sd = 1/rate/sqrt(ns)), add = TRUE,
        col = "blue", lwd = 2)

For a sample size of 60, the resultant sample means are more accurate, and the corresponding normal density curve fits these distribution of these sample means extremely well (as expected via the CLT).

4 CLT Simulation - Bernoulli Distribution

4.1

No answer required.

4.2

#Specify sample size
ns <- 5

# Specify success rate
p = 0.3

# Specify number of times to conduct process  
trials <- 10000

# create object in which to store sample means once they are generated
norm.means <- rep(0, trials) # Note length is equal to the number of trials

# Run for loop  
for(i in 1:trials){
    # Randomly generate ns samples from the Bernoulli distribution
    x <- rbinom(ns, 1, p)
    # Compute mean of samples, and store it in the norm.means object, in ith position 
    norm.means[i] <- mean(x) 
    }

hist(norm.means, freq = FALSE, breaks = 20, col = "red", 
     xlim = c(0, 1),
     xlab = expression(bar(x)),
     main = paste("Histogram of means, n = ", ns)) 

# Set mu and sigma for the approximating normal distribution
mu <- p
sigma <- sqrt(p*(1 - p))
  
# Overlay the normal density
curve(dnorm(x, mu, sigma/sqrt(ns)), add = TRUE, col = "blue", lwd = 2)

4.3

The plot of the sample means is clearly not normally distributed.

4.4

#Specify sample size
ns <- 30
  
# Specify success rate
p = 0.3

# Specify number of times to conduct process  
trials <- 10000

# create object in which to store sample means once they are generated
norm.means <- rep(0, trials) # Note length is equal to the number of trials

# Run for loop  
for(i in 1:trials){
    # Randomly generate ns samples from the Bernoulli distribution
    x <- rbinom(ns, 1, p)
    # Compute mean of samples, and store it in the norm.means object, in ith position 
    norm.means[i] <- mean(x) 
    }

hist(norm.means, freq = FALSE, breaks = 20, col = "red", 
     xlim = c(0, 1),
     xlab = expression(bar(x)),
     main = paste("Histogram of means, n = ", ns)) 

# Set mu and sigma for the approximating normal distribution
mu <- p
sigma <- sqrt(p*(1 - p))
  
# Overlay the normal density
curve(dnorm(x, mu, sigma/sqrt(ns)), add = TRUE, col = "blue", lwd = 2)

All of a sudden, the distribution of the sample means begins to look much more normally distributed.

4.5

#Specify sample size
ns <- 60
  
# Specify success rate
p = 0.3

# Specify number of times to conduct process  
trials <- 10000

# create object in which to store sample means once they are generated
norm.means <- rep(0, trials) # Note length is equal to the number of trials

# Run for loop  
for(i in 1:trials){
    # Randomly generate ns samples from the Bernoulli distribution
    x <- rbinom(ns, 1, p)
    # Compute mean of samples, and store it in the norm.means object, in ith position 
    norm.means[i] <- mean(x) 
    }

hist(norm.means, freq = FALSE, breaks = 20, col = "red", 
     xlim = c(0, 1),
     xlab = expression(bar(x)),
     main = paste("Histogram of means, n = ", ns)) 

# Set mu and sigma for the approximating normal distribution
mu <- p
sigma <- sqrt(p*(1 - p))
  
# Overlay the normal density
curve(dnorm(x, mu, sigma/sqrt(ns)), add = TRUE, col = "blue", lwd = 2)

For a sample size of 60, the sample means are now looking as if they almost could have been generated from a normal distribution, which is amazing given they are obtained from data generated from a discrete distribution.

5 Extension - Average Diamond Prices

The R code below should have been run:

install.packages("ggplot2")
library(ggplot2)

5.1

hist(diamonds$price, col = "skyblue", xlab = "Diamond Price ($)",
     main = "Histogram of Diamond Prices ($)", freq = F)

5.2

Example R code is provided below:

#Specify sample size
ns <- 5
  
# Specify number of times to conduct process  
trials <- 100

# create object in which to store sample means once they are generated
norm.means <- rep(0, trials) # Note length is equal to the number of trials

# Run for loop  
for(i in 1:trials){
    # Randomly sample ns values from the diamonds data set
    x <- sample(diamonds$price, ns)
    # Compute mean of samples, and store it in the norm.means object, in ith position 
    norm.means[i] <- mean(x) 
    }

hist(norm.means, freq = FALSE, col = "skyblue",
     xlab = expression(bar(x)),
     main = paste("Histogram of means, n = ", ns)) 

We note that the distribution of the sample means is quite skewed.

5.3

#Specify sample size
ns <- 30
  
# Specify number of times to conduct process  
trials <- 100

# create object in which to store sample means once they are generated
norm.means2 <- rep(0, trials) # Note length is equal to the number of trials

# Run for loop  
for(i in 1:trials){
    # Randomly sample ns values from the diamonds data set
    x <- sample(diamonds$price, ns)
    # Compute mean of samples, and store it in the norm.means object, in ith position 
    norm.means2[i] <- mean(x) 
    }

hist(norm.means2, freq = FALSE, col = "skyblue",
     xlab = expression(bar(x)),
     main = paste("Histogram of means, n = ", ns)) 

By increasing our sample size to 30, the distribution of the sample means now appears more symmetric and closer to that of a normal distribution. This results adheres with the Central Limit Theorem. Further increases in the sample size would result in a sample mean distribution that even more closely matched that of a normal distribution.

Great job, that’s everything for today!

Hopefully, this computer lab has helped cement your understanding of the Central Limit Theorem.

References

Horst, Allison Marie, Alison Presmanes Hill, and Kristen B Gorman. 2020. Palmerpenguins: Palmer Archipelago (Antarctica) Penguin Data. https://doi.org/10.5281/zenodo.3960218.

Wickham, H., W. Chang, L. Henry, T. L. Pedersen, K. Takahashi, C. Wilke, K. Woo, H. Yutani, D. Dunnington, and RStudio. 2022. ggplot2: Create Elegant Data Visualisations Using the Grammar of Graphics. https://cran.r-project.org/web/packages/ggplot2/index.html.

Wickham, H., R. Francois, L. Henry, K. Muller, and RStudio. 2021. dplyr: A Grammar of Data Manipulation. https://dplyr.tidyverse.org, https://github.com/tidyverse/dplyr.

These notes have been prepared by Rupert Kuveke and Amanda Shaker. The copyright for the material in these notes resides with the authors named above, with the Department of Mathematical and Physical Sciences and with La Trobe University. Copyright in this work is vested in La Trobe University including all La Trobe University branding and naming. Unless otherwise stated, material within this work is licensed under a Creative Commons Attribution-Non Commercial-Non Derivatives License BY-NC-ND.