How exactly did statisticians agree to using (n-1) as the unbiased estimator for population variance without simulation? 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?
 A: Most proofs I have seen are simple enough that Gauss (however he did it) probably found it pretty easy to prove.
I've been looking for a derivation on CV that I could link you to (there are a number of links to proofs off-site, including at least one in answers here), but I haven't found one here on CV in a couple of searches, so for the sake of completeness, I'll give a simple one. Given its simplicity, it's easy to see how people would start to use what's usually called Bessel's correction.
This takes $E(X^2)=\text{Var}(X) + E(X)^2$ as assumed knowledge, and assumes that the first few basic variance properties are known.
\begin{eqnarray}
E[\sum_{i=1}^{n} (x_i-\bar x)^2] &=& E[\sum_{i=1}^{n} x_i^2-2\bar x\sum_{i=1}^{n} x_i+n\bar{x}^2]\\
&=& E[\sum_{i=1}^{n} x_i^2-n\bar{x}^2] \\
&=& n E[x_i^2]- n E[\bar{x}^2]\\
&=& n (\mu^2 + \sigma^2) - n(\mu^2+\sigma^2/n)\\
&=& (n-1) \sigma^2
\end{eqnarray}
A: 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.
A: According to Weisstein's World of Mathematics, it was first proved by Gauss in 1823. The reference is volume 4 of Gauss' Werke, which can be read at https://archive.org/details/werkecarlf04gausrich. The relevant pages seem to be 47-49. It seems that Gauss investigated the question and came up with a proof. I don't read Latin, but there is a German summary in the text. Pages 103-104 explain what he did (Edit: I added a rough translation):

Allein da man nicht berechtigt ist, die sichersten Werthe fuer die
wahren Werthe selbst zu halten, so ueberzeugt man sich leicht, dass
man durch dieses Verfahren allemal den wahrscheinlichsten und
mittleren Fehler zu klein finden muss, und daher die gegebenen
Resultaten eine groessere Genauigkeit beilegt, als sie wirklich
besitzen. [But since one is not entitled to treat the most probable values as though they were the actual values, one can easily convince oneself that one must always find that the most probable error and the average error are too small, and that therefore the given results possess a greater accuracy than they really have.]

from which it would seem that it was well-known that the sample variance is a biased estimate of the population variance. The article goes on to say that the difference between the two is usually ignored because it's not important if the sample size is big enough. Then it says:

Der Verfasser hat daher diesen Gegenstand eine besondere Untersuchung
unterworfen, die zu einem sehr Merkwuerdigen hoechst einfachen
Resultate gefuehrt hat. Man braucht nemlich den nach dem angezeigten
fahlerhaften Verfahren gefundenen mittleren Fehler, um ihn in die
richtigen zu verwandeln, nur mit
$$\sqrt{\frac{\pi-\rho}{\pi}}$$
zu multiplicieren, wo $\pi$ die Anzahl der beobachtungen (number of
observations) und $\rho$ die Anzahl der unbekannten groessen (number
of unknowns) bedeutet. [The author has therefore  made a special study of this object which has led to a very strange and extremely simple result. Namely, one needs only to multiply the average error found by the above erroneous process by (the given expression) to change it into the right one, where $\pi$ is the number of observations and $\rho$ is the number of unknown quantities.]

So if this is indeed the first time that the correction was found, then it seems that it was found by a clever calculation by Gauss, but people were already aware that some correction was required, so perhaps someone else could have found it empirically before this. Or possibly previous authors didn't care to derive the precise answer because they were working with fairly large data sets anyway.
Summary: manual, but people already knew that $n$ in the denominator wasn't quite right.
A: For me one piece of intuition is that
$$\begin{array}{c}
\mbox{The degree to which}\\
X_{i}\mbox{ varies from }\bar{X}
\end{array}+\begin{array}{c}
\mbox{The degree to which}\\
\bar{X}\mbox{ varies from }\mu
\end{array}=\begin{array}{c}
\mbox{The degree to which }\\
X_{i}\mbox{ varies from }\mu.
\end{array}$$
That is, 
$$\mathbf{E}\left[\left(X_{i}-\bar{X}\right)^{2}\right]+\mathbf{E}\left[\left(\bar{X}-\mu\right)^{2}\right]=\mathbf{E}\left[\left(X_{i}-\mu\right)^{2}\right].$$
Actually proving the above equation takes a bit of algebra (this algebra is very similar to @Glen_b's answer above). But assuming it is true, we can rearrange to get: 
$$\mathbf{E}\left[\left(X_{i}-\bar{X}\right)^{2}\right]=\underset{\sigma^{2}}{\underbrace{\mathbf{E}\left[\left(X_{i}-\mu\right)^{2}\right]}}-\underset{\frac{\sigma^{2}}{n}}{\underbrace{\mathbf{E}\left[\left(\bar{X}-\mu\right)^{2}\right]}}=\frac{n-1}{n}\sigma^2.$$
For me, another piece of intuition is that using $\bar{X}$ instead of $\mu$ introduces bias. And this bias is exactly equal to $\mathbf{E}\left[\left(\bar{X}-\mu\right)^{2}\right]=\frac{\sigma^2}{n}$.
A: Most of the answers have already elaborately explained it but apart from those there's one simple illustration that one could find helpful:
Suppose you are given that $n=4$ and the first three numbers are:
$8,4,6$,_
Now the fourth number can be anything since there are no constraints. Now consider the situation when you are given that $n=4$ and $\bar x=6$, then if the first three numbers are: $8,4,6$ then the fourth number has to be $6$. 
This is to say that if you know $n-1$ values and $\bar x$, then the $nth$ value has no freedom. Thus $n-1$ gives us an unbiased estimator.
