(On editing this, I realize that I should have computed the number of items---length(x)---in the var functions and then passed it into the sumsq function rather than compute it twice.)

We can see that the default var function in R uses the sample method:

Let's grab some numbers from rnorm and organize them as a matrix. Since we haven't specified a mean and standard deviation, they will get the default values μ = 0 and σ = 1.

Now we use the built-in function "apply", the second argument (=1) says to apply the indicated function to m by rows.

No question about it. Dividing by n underestimates the population variance, while n-1 looks good. It provides, as they say, an unbiased estimator. I did not appreciate that the sample means have a strong positive skew. That's interesting. Looks like some kind of gamma distribution.

It also works for the uniform distribution: