STM1001: Computer Lab 5B Solutions

Data Science Module

Topic 5B: Simulations in R

Example R code solutions for the Data Science Computer Lab 5, which uses the R package dplyr (Wickham et al. 2021) 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

set.seed(1)
n <- 5 # Specify sample size
x <- rnorm(n) # Randomly generate ns samples from the normal distribution
norm_mean <- mean(x) # Compute the mean of the samples, and store it in a new object

2.5

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, and store it 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, 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(-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, 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(-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

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

set.seed(1)
n <- 100
p = 0.3
x <- rbinom(n, 1, p)

hist(x, freq = FALSE, col = "chartreuse3", xlim = c(0, 1),
     main = paste("Histogram of random sample \n from Bernoulli distribution, n = ", n)) 
curve(dnorm(x, p, sd = sqrt(p*(1 - p))), add = TRUE, col = "blue", lwd = 2)

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 approximately normal (i.e. symmetric). This results adheres with the Central Limit Theorem.

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., 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.