Weird R behavior when simulating LLN for standard deviations

Question

I am simulating the law of large numbers and how it applies to standard deviation, as well. I wrote this code that works well but there is something that I am having a hard time understanding.

n <- 1000
dice.sd <- numeric(n)

sd(1:6) #1.870829

for (i in 1:n) {
  dice.sd[i] <- sd(sample(1:6, i, replace = TRUE))
}
plot(dice.sd)
abline(h=1.870829)
abline(h=1.7078)

As you can see, I made this loop to simulate LLN for standard deviations. According to the documentation, the sd() function uses n-1 for calculating the sample standard deviation, which should be about 1.87 for a die. However, when I run my simulation and graph the results, the standard deviation is converging to about 1.7078, which is the population standard deviation (using just n). Why is this the case? My loop originally was using the sample standard deviation, so why is it converging to the population standard deviation?

question from:https://stackoverflow.com/questions/65876392/weird-r-behavior-when-simulating-lln-for-standard-deviations

Answer 1

dice.sd[i] <- sd(sample(1:6, i, replace = TRUE))

This line creates a sample of size i, not a sample of size 6, which I assume is what you want. Your i goes up to 1000, which is large enough for the sample SD to be close to the population SD.

What you want is

dice.sd[i] <- sd(sample(1:6, 6, replace = TRUE))