The formula for computing variance has $(n-1)$ in the denominator:

$s^2 = \frac{\sum_{i=1}^N (x_i - \bar{x})^2}{n-1}$

I've always wondered why. However, reading and watching a few good videos about "why" it is, it seems, $(n-1)$ is a good unbiased estimator of the population variance. Whereas $n$ underestimates and $(n-2)$ overestimates the population variance.

What I'm curious to know, is that in the era of no computers how exactly was this choice made?  Is there an actual mathematical proof proving this or was this purely empirical and statisticians made A LOT of calculations by hand to come up with the "best explanation" at the time? 

Just how did statisticians come up with this formula in the early 19th century with the aid of computers? Manual or there is more to it than meets the eye?

The correction is called Bessel's correction and it has a mathematical proof. Personally, I was taught it the easy way: using $n-1$ is how you correct the bias of $E[\frac{1}{n}\sum_1^n(x_i - \bar x)^2]$ (see here).

You can also explain the correction based on the concept of degrees of freedom, simulation isn't strictly needed.
