STM1001: Computer Lab 5B

Data Science Stream

Topic 5B: Simulations in R

Welcome to the fifth computer lab for the Data Science stream of STM1001.

In Topic 4, we introduced sampling distributions and the concept of the Central Limit Theorem.

In this computer lab we will use R to conduct the simulations discussed in Chapter 3 of Topic 4.

• It might be helpful to keep this content open in a separate tab while you work through the lab material, in case you would like to quickly double-check your understanding of any concepts.

By the end of this lab, you should be able to conduct simple simulations in R. Let’s get started!

1 Sampling Distributions

1.1 Sampling

🏡 Recall that a sample is a random selection of units from a chosen population.

• For example, if the STM1001 cohort was our chosen population, then a random selection of STM1001 students would be a sample from this population.

Because it is typically unfeasible to collect data on an entire population of interest, we instead conduct analyses using a sample. We hope that the inferences we make using the sample data can be extrapolated to provide information about the population, and we can use statistical techniques to ascertain the accuracy of these inferences.

One such statistical technique we can utilise is the Central Limit Theorem.

1.2 The Central Limit Theorem

🏡 Recall the definition of the Central Limit Theorem (CLT):

Let \(X_1, \ldots, X_n\) be a random sample from a distribution with finite mean \(\mu\) and finite variance \(\sigma^2\).

For \(\overline{X}\) denoting the sample mean, if \(n\) is sufficiently large then \[\overline{X}\stackrel{\tiny \text{approx.}}\sim N\left(\mu,\frac{\sigma^2}{n}\right) , \] where \(\stackrel{\tiny \text{approx.}}\sim\) denotes ‘approximately distributed as’.

By the CLT, it follows that as we increase our sample size, the sample mean \(\overline{X}\) will get increasingly closer (i.e. converge) to the population mean \(\mu\). As \(n\) increases, the variance of the sample mean \(\overline{X}\) will also decrease.

In other words, while samples of individuals can exhibit significant variation, samples of means tend to have less variability (as each mean is calculated from a sample of individuals, reducing the impact of outliers).

Let’s visualise this process in R using some simulations.

2 Simulations in R

Carrying out simulations in R is actually quite straightforward. We will walk through the steps involved one at a time.

2.1

💻 Suppose that we are interested in the random variable \(X \sim N(0, 1)\) - in other words, here \(X\) follows the standard normal distribution.

We can generate a sample of 100 observations from this distribution using the following R code:

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

Run this code now.

Note: The default mean and standard deviation values used by rnorm are 0 and 1 respectively, so we do not need to specify them in this instance.

2.2

💻 We can plot our generated data using the hist function, as shown below.

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)

Note that we have also used the curve and dnorm functions here to overlay the standard normal distribution probability density curve.

2.3

💻 Our main interest here is the sample mean, not the variability of the sample values.

If we check the sample mean of our data set x, we find that this is equal to approximately 0.109 (shown by the red dashed line in the plot above). Since this is the sample mean for just one sample of 100 observations, we can’t be too confident in using this estimate as an accurate estimate of the population mean (which we know was specified to be 0).

Rather than relying on just the one sample of 100 observations, we will now simulate a number of random samples from \(X\), in order to obtain multiple sample mean estimates, which we can use in turn to obtain more accurate estimate of the population mean.

2.4

💻 To begin, suppose we sample just 5 observations of \(X\), rather than 100, and then compute the sample mean of these 5 observations.

Using the code in 2.1, it follows that this process is quite straightforward:

set.seed(1)
n <- 5 # Specify sample size
x <- rnorm(n) # Generate n samples from the standard normal distribution
norm_mean <- mean(x) # Store the mean of the samples in a new object

🎧 Online students

💬 Run your object norm_mean and post your result in the chat.

2.5 for loops

💻 Since our underlying data is generated from a standard normal distribution, it follows that the distribution of the sample mean will also be normal (we can think of this result as a precursor to the Central Limit Theorem).

To accurately visualise the distribution of the sample mean, we will need to repeat the process outlined in 2.4 many times. This is where our simulations can really shine.

Suppose, rather than conducting the process outlined in 2.4 once, we instead conduct it 10,000 times. It would take a while to write this code out line by line in R (over 9000 lines potentially!).

Instead, to save time, we will use a for loop. A for loop is a function that will repeat a specified set of operations a specified number of times.

The Code chunk below provides a simple example of the R syntax for a for loop.

for(i in 1:10){
  x <- rnorm(i)
  mean(x)}

How do we interpret this code? Well,we have specified that for some input i, where i can take the integer values 1 up to 10, the for loop will simulate i random values from the standard Normal distribution, store them in the object x, and then compute the mean of x.

Thus, our for loop runs like this:

• The for loop starts with the first specified value, 1, and simulates 1 random value from the standard Normal distribution.
• The mean of x is now computed. This is the last operation in the for loop.
• Since the specified operations inside the loop have been finished for that value of i, the loop starts again, with i=2 now.
• This process continues, until the final specified value of i, namely i=10, has been used.
• At this point, the for loop is complete.

Note: We do not have to use the letter i here, it is simply convention. You can use a different letter, but remember, make sure you have not used your chosen letter as a name for a previous object in which you have stored data.

2.6

💻 The Code chunk below expands upon the for loop example discussed above in 2.5. Take a look at this code, read the comments carefully, and make sure you understand what the code does (note some values are missing, as denoted by the ...).

#Specify sample size
ns <- ...
# Specify number of times to conduct process  
trials <- ...
# 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)}

Once you feel confident in your understanding of this code, replace the ... missing details for ns and trials, using the following details:

• You want to take 10,000 samples of randomly generated values from the standard Normal distribution.
• Each sample should have a sample size of 5.

Once you have filled in the details, run your code.

Hint: If you have given this a decent shot, but are stuck, you can check the Code chunk below:

#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

💻 To visualise our results, we can use the following code:

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

We can also overlay a normal density curve, that, since our underlying data is normal, will reflect the distribution of the sample mean.

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

Run these code chunks now, and check your results. Are our simulated sample means close to the population mean of 0?

2.8

💻Repeat steps 2.6 and 2.7, using a sample of 30 instead of 5. What do you notice about the distribution of the simulated sample means?

🎧 Online students

💬 Leave a comment in the chat about the changes you observe in the distribution, changing from \(n=5\) to \(n=30\).

2.9

💻 Finally, repeat steps 2.6 and 2.7, using a sample of 60. What do you notice about the distribution of the simulated sample means now?

🎧 Online students

💬 Leave a comment in the chat about the changes you observe in the distribution, changing from \(n=30\) to \(n=60\).

3 CLT Simulation - Exponential Distribution

💻 To demonstrate the versatility of the Central Limit Theorem, let’s carry out some simulations involving data generated from the Exponential distribution, which (unlike the Normal distribution) is highly skewed (i.e. asymmetric).

3.1

💻 Recall from Topic 3 that the Exponential distribution is defined by one parameter, known as the rate parameter. The rate parameter is denoted by the Greek letter \(\lambda\).

If a random variable \(X\) follows an Exponential distribution with rate parameter \(\lambda\), we write \(X\sim Exp(\lambda)\).

Suppose now that our random variable \(X\) follows an exponential distribution with rate parameter \(\lambda = 10\). In other words, suppose \(X \sim \text{Exp}(10)\).

We can generate a sample of 100 observations from this distribution using the following R code:

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

Note that we use rexp to randomly generate values from the exponential distribution, instead of using rnorm (which is for the normal distribution).

Run this code now.

3.2

💻 Our generated data from 3.1 is highly skewed, and clearly different to data generated from the symmetric normal distribution. This disparity can be seen in Figure 3.1 below, which shows a histogram of the data, with a normal distribution probability density curve overlaid in blue (similar to in 2.2), based on appropriate mean and variance values.

Figure 3.1: Data generated from an Exponential distribution, with a normal curve overlaid.

Clearly, we cannot fit a normal curve well to data generated from an Exponential distribution. However, as we work through this question, we will find that when we model the distribution of the sample means of data generated from an Exponential distribution, this distribution of sample means can be fitted well by a normal curve, especially as our sample size increases. This (somewhat unintuitive) result is due to the Central Limit Theorem!

3.3 The Rate Parameter and Population Mean

💻 Instead of considering the individual data points generated from the Exponential distribution, let’s shift our focus to the sample mean of these data points.

For data generated from an Exponential distribution, the population mean \(\mu\) is in theory equal to 1 divided by the rate parameter \(\lambda\). I.e., for \(X\sim Exp(\lambda)\), the population mean of \(X\) equals \(\frac{1}{\lambda}\).

If we check the sample mean of x, we find that this is equal to 0.1030676 - this sample mean is quite a good estimate of the population mean (\(\mu = 1/10 = 0.1\)), given the theory.

3.4

💻 Let us now conduct the simulation process used in 2.6, for our data sampled from \(X \sim \text{Exp}(10)\).

Using the code in 2.6 as a guide, simulate 10,000 sample means of \(X\), using samples of size 5.

Plot your results, using the code in 2.7 as a guide.

Note: You will have to change the x-axis range for the histogram (via the xlim argument). Try using a range of 0 to 0.4.

Note 2: Use the code in the Code chunk below to overlay the appropriate Normal distribution probability density curve:

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

3.5

💻 What do you notice about your resultant plot?

3.6

💻 Repeats 3.4 for \(X \sim Exp(10)\), using a sample of 30 instead of 5. What do you notice about the distribution of the simulated sample means?

🎧 Online students

💬 Leave a comment in the chat about the changes you observe in the distribution, changing from \(n=5\) to \(n=30\).

3.7

💻 Finally, repeat 3.4 for \(X \sim Exp(10)\), using a sample of 60. What do you notice about the distribution of the simulated sample means now?

🎧 Online students

💬 Leave a comment in the chat about the changes you observe in the distribution, changing from \(n=30\) to \(n=60\).

4 CLT Simulation - Bernoulli Distribution

💻 As discussed in Section 3.3 of Topic 4, the Central Limit Theorem even applies to discrete random variables.

Suppose that our random variable \(X\) now follows a Bernoulli distribution with success parameter \(p = 0.3\). In other words, suppose \(X \sim \text{BERN}(0.3)\).

4.1

💻 If we generate a sample of 100 observations from this distribution (see the R code below) and then visualise our generated data, we see a roughly 70/30 split of observations between \(x=0\) and \(x=1\). Clearly, (just as in the Exponential distribution case), the normal distribution probability density curve, with appropriate parameter specifications, (shown in blue), is an extremely poor fit for this data.

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

Note that we use rbinom to randomly generate values from the Bernoulli distribution, instead of using rnorm (which is for the normal distribution).

4.2

💻 However, despite our data being generated from a discrete distribution, we can still use the Central Limit Theorem to obtain an accurate estimate of the population mean by considering the distribution of \(\bar{X}\)!

Let us now conduct the simulation process used in 2.6, for our Bernoulli-distributed data.

Using the code in 2.6 as a guide, simulate 10000 sample means of \(X\), using samples of size 5. Then, plot your results, using the code in 2.7 as a guide.

Note: You will have to change the x-axis range for the histogram (via the xlim argument).

Note 2: Use the code in the Code chunk below to overlay a normal distribution probability density curve with the appropriate parameter specifications:

# 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

💻 What do you notice about your resultant plot?

4.4

💻 Repeat 3.4 for \(X \sim BERN(0.3)\), using a sample of 30 instead of 5. What do you notice about the distribution of the simulated sample means?

🎧 Online students

💬 Leave a comment in the chat about the changes you observe in the distribution, changing from \(n=5\) to \(n=30\).

4.5

💻 Finally, repeat 3.4 for \(X \sim BERN(0.3)\), using a sample of 60. What do you notice about the distribution of the simulated sample means now?

🎧 Online students

💬 Leave a comment in the chat about the changes you observe in the distribution, changing from \(n=30\) to \(n=60\).

5 Extension: Average Diamond Prices

💻 Let’s consider how we can apply the Central Limit Theorem to real data.

For this question, we will consider data on approximately 54,000 round cut diamonds. This data is stored in the diamonds data set in the ggplot2 R package. We will treat this data as the population data.

Install and load the ggplot2 R package now.

Hint: If you need a refresher on installing and loading packages in R, check the Code chunk below:

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

5.1

💻 If we plot the price of round cut diamonds using a histogram, we can see that the data is clearly not normally distributed.

Make sure to run this code before continuing.

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

5.2

💻 Suppose that you are interested in estimating the population mean price for round cut diamonds, but do not have access to the full diamonds data set.

Instead, you can only sample prices for 5 round cut diamonds at a time, represented by the code below:

x <- sample(diamonds$price, 5)

Using this code and the code in 2.6 as a guide, obtain 100 sample means from 100 simulated samples of 5 random round cut diamond sale prices. Plot a histogram of these sample means. What do you observe?

🎧 Online students

💬 Volunteer to share your screen to show your histogram and discuss your results, if you got to this stage before the end of the lab.

5.3

💻 Repeat 5.2, but this time suppose that you are able to take samples of 30 diamonds at a time, instead of 5. What changes do you observe?

Great job, that’s everything for today!

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

References

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.