nbr_students = 10
nbr_same = 2
birthdays <- sample(1:366, size = nbr_students, replace = TRUE)
max(table(birthdays)) >= nbr_same

## [1] FALSE

We are not interested in the birthdays themselves, just whether at least two are the same. table(birthdays) records the number of times each calendar day was a birthday.

The code results in either FALSE (equivalently 0) if there isn’t a shared birthday, or TRUE (equivalently 1) if there is.

nbr_students = 10
nbr_same = 2
numTrials <- 1000
numHits = 0

for (t in 1:numTrials) {
  birthdays <- sample(1:366, size =  nbr_students, replace = TRUE)
  numHits <- numHits + (max(table(birthdays)) >= nbr_same)
}
numHits/numTrials

## [1] 0.125

How do your simulation results compare to the actual probabilities for a shared birthday?

• 11.66% for a group of 10,

• 41.06% for a group of 20, and

• 89.05% for a group of 40 students.

FYI: The probability that in a group of N students at least 2 share a birthday is: \[ 1 - \frac{366}{366}\cdot\frac{366 - 1}{366} \cdots \frac{366 - N}{366} \] No need to worry about this formula beyond this: It follows by counting possible events, and the \(1-\ldots\) structure comes from considering the complementary case, the one in which no students share a birthday.

Estimated probabilities (sample size 10000):

EngWal <- read.csv("EnglandWalesBirthdays.csv")

nbr_students = 50
nbr_same = 4
numTrials <- 10000
numHits = 0

for (t in 1:numTrials) {
  birthdays <- sample(1:366, size =  nbr_students, replace = TRUE, prob = EngWal$X2016 )
  numHits <- numHits + (max(table(birthdays)) >= nbr_same)
}
numHits/numTrials

## [1] 0.004

hist(pop, breaks = 40, col = "lightblue", 
     main="Histogram of entire population", xlab = "Temperature")

qqnorm(pop)
qqline(pop, col = "deeppink")

SRS <- sample(pop, size = 100, replace = FALSE)
par(mfrow=c(1,2))
hist(SRS, col="lightblue", xlim=range(pop), freq=FALSE)
hist(pop,col="lightblue", main="Histogram of temperatures", xlim=range(pop), freq=FALSE)
par(mfrow=c(1,1))

• Generally, the sample will have similar characteristics as the population, like skew to the left. But it really depends on your particular sample!

• With a sample size of only 100 out of over 400,000 data points, the mean and standard deviation of the sample are nevertheless close to \(\mu\) and \(\sigma\).

mean(SRS)

## [1] 16.764

sd(SRS)

## [1] 8.783513

We use a for-loop to draw 1000 random samples of size 100 and store the resulting means.

trials    <- 1000
estimates <- numeric(trials)

for (i in 1:trials) {
  SRS <- sample(pop, size = 100, replace = FALSE)
  estimates[i] <- mean(SRS)
}

hist(estimates, col="lightblue", breaks=20, xlim=range(estimates),
     main="Histogram of sample means")
abline(v = mu, lwd = 2, col = "deeppink")

qqnorm(estimates)
qqline(estimates, col="deeppink")

The histogram of the sample means looks like a normal distribution, centred around the true population mean.

mean(estimates)

## [1] 16.26916

mu

## [1] 16.29877

trials <- 2000
estimates_more_trials <-numeric(trials)

for (i in 1:trials) {
  SRS <- sample(pop, size = 100, replace = FALSE)
  estimates_more_trials[i] <- mean(SRS)
}

par(mfrow=c(1,2))
hist(estimates_more_trials, breaks = 40, col = "lightblue",
     main="Sample means from more trials", xlim=range(estimates), freq = FALSE)
hist(estimates, breaks = 40, col = "lightblue",
     main="Sample means from fewer trials", xlim=range(estimates), freq = FALSE)
par(mfrow=c(1,1))

trials <- 1000
estimates_lss <-numeric(trials)

for (i in 1:trials) {
  SRS <- sample(pop, size = 500, replace = FALSE)
  estimates_lss[i] <- mean(SRS)
}

par(mfrow=c(1,2))
hist(estimates_lss, breaks=20, col="lightblue", xlim=range(estimates), freq=FALSE,
     main="Sample means larger sampling size")
hist(estimates, breaks=40, col="lightblue", xlim=range(estimates), freq=FALSE,
     main="Sample means smaller sampling size")
par(mfrow=c(1,1))

• If we repeatedly choose samples from the population, the sample mean will take different values in different samples.

• A simple idea: if we repeat sampling many times and the mean does not change substantially (i.e. if we get similar answers each time), our estimates are fairly reliable.

• A more advanced idea: The distribution of sample means tends towards a normal distribution. So we can use the 68-95-99 rule ta quantify how reliable we expect our estimates to be. (We will return to this idea in later classes…)

This sounds great! But is it practical?

• In order to find the standard deviation of the sampling distribution, we drew 1000 random samples of size 200.

• This is not something we would do in real life. No one is going to conduct 1000 polls of 200 people each just to quantify the uncertainty of one poll!

• We have used simulations to illustrate the relationship between samples and populations, but usually we are working with one sample.

sd(pop)

## [1] 9.43757

For samples of size \(n=100\) we found

sd(estimates)

## [1] 0.9747276

which differs by a factor of \(\sqrt{100}\).

The approximation \(\sigma/\sqrt{n}\) is called the standard error of the mean.