(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: