Combinatorial\Probabilistic Proof of Stirling's Approximation Stirling's approximation is the following well-known asymptotic result:
$$n! \approx \left(\frac{n}{e}\right)^n \sqrt{2 \pi n}$$
This result has several analytical proofs, for example via Laplace's method, the trapezoidal rule (and Euler–Maclaurin formula for an asymptotic expansion) or Hayman's method (essentially Cauchy's formula). There are also other (usually ad-hoc) proofs.
The result is applied often in combinatorics and probability, especially in the study of random walks. I want a result which is the other way around - a combinatorial\probabilistic proof for Stirling's approximation. I'm not sure if this is possible, but to convince you that it might be I'll give some partial results.
First, write $n! =  \left(\frac{n}{e}\right)^n f(n)$. We'll obtain bounds related to $f(n)$:


*

*Consider $\{S_n\}_{n\ge 0}$, the random walk on $\mathbb{Z}^{d}$, with $S_0=\vec{0}$. It is recurrent iff $\sum_{n\ge 1} P(S_{2n} = \vec{0})$ diverges. This series is $\sum_{n\ge 1} (2d)^{-2n} \sum_{a_1+\cdots+a_d=n} \binom{2n}{a_1,a_1,a_2,a_2,\cdots,a_d,a_d}$. For $d=1$, this becomes $\sum_{n\ge 1} \frac{\binom{2n}{n}}{4^n} = \sum_{n\ge 1} \frac{f(2n)}{f^2(n)}$. There's is Stirling-free proof for the recurrence of the 1-dimensional random walk, using a criterion due to Nash-Williams.

*Similarly, for $d=2$, there's a rotation trick that shows $\sum_{n\ge 1} P(S_{2n} = \vec{0})=\sum_{n\ge 1} \left( \frac{\binom{2n}{n}}{4^n} \right)^2= \sum_{n\ge 1} \left( \frac{f(2n)}{f^2(n)} \right)^2$. Again, the Nash-Williams criterion proves the recurrence of the 2-dimensional random walk without Stirling, and hence the divergence of said series.

*Back to the 1-dimensional case. Using the reflection principle, the probability that the random walk returns to the origin during the first $2n$ steps is $1-\frac{\binom{2n}{n}}{4^n} =  1-\frac{f(2n)}{f^2(n)}$, and this tends to 1 as the walk is recurrent.

*Still with the 1-dimensional case. One can calculate $E|\frac{S_n}{\sqrt{n}}|$ explicitly - it is $2^{1-n} \sqrt{n}\binom{n-1}{\lfloor \frac{n-1}{2} \rfloor}$. On the other hand, the CLT together with the notion of uniform integrability justifies the following limit: $\lim_{n\to\infty}  E|\frac{S_n}{\sqrt{n}}| = E|N(0,1)| = \frac{2}{\sqrt{2\pi}}$. When $n$ is odd this implies $\sqrt{n} \frac{f(n-1)}{f^2(\frac{n-1}{2})} \to \frac{2}{\sqrt{2\pi}}$.


So if we assume $f(n) \approx C \cdot n^{\alpha}$, the first results give bounds on $\alpha$ ($\alpha \le 1, \alpha \le 0.5, \alpha >0$) and the last one gives the true value of $\alpha$ and $C$. The fact that the 3-dimensional random walk is transient implies $\alpha \ge \frac{1}{3}$, but I don't know how to prove it without tools that prove Stirling's approximation.
Questions:
Can we remove\relax the assumption $f(n) \approx C \cdot n^{\alpha}$ somehow? Can we use less analysis than we used in the 3rd result? Are there other combinatorial approaches to the problem that give partial results (or even the full one)? Can the fact that the 3-dimensional random walk is transient be proved Stirling-free (and preferably combinatorially)?
 A: The Wikipedia page List of probabilistic proofs of non-probabilistic theorems has a reference to the paper:

Blyth, Colin R.; Pathak, Pramod K.
   A Note on Easy Proofs of Stirling's Theorem.
  Amer. Math. Monthly 93 (1986), no. 5, 376–379. MR 1540867 DOI 10.2307/2323600

in which several simple proofs of Stirling's approximation are given, using the central limit theorem on Gamma or Poisson random variables.
A: Let me finish Qiaochu's answer. 
Observe that the ratio $\mathbb P ( S_n=n+k) / \mathbb P(S_n=n)$ is close to $1$: It is 
$$\prod_{j=1}^k \frac{n}{n+j}$$.
Using $$ e^{-j/n} \leq \frac{n}{n+j} \leq 1$$
The probability ratio is bounded between $e^{- k^2/n}$ and $1$.
Hence by averaging over an interval:
$$ e^{-\epsilon^2} \mathbb P( S_n=n) \leq \frac{\mathbb P \left( n\leq S_n< n+ \epsilon \sqrt{n} \right)  }{\epsilon \sqrt{n}}  \leq \mathbb P( S_n =n )$$
Applying the estimate from the central limit theorem:
$$\mathbb P \left( n\leq S_n< n+ \epsilon \sqrt{n} \right) \cong \frac{\epsilon}{\sqrt{2\pi}}+ O\left(\epsilon^3\right)$$
we get the estimate
$$\frac{ n^n e^{-n}}{n!} = \mathbb  P( S_n=n) =\frac{1}{\sqrt{2\pi n}} + O \left( \epsilon ^2\right)$$
Letting $\epsilon$ tend slowly to $0$, we win.
A: Meanwhile, I'll answer the other question about random walks: Yes, it is quite easy to prove random walks in $\geq 3$ dimensions are transient without Stirling's approximation. Indeed, this is the standard proof I know, and was used in the second hit when I googled for it here. See Section 2 here. In summary, let $p_n$ be the probability that the random walk returns to $0$ at time $n$. One shows that the walk is transient if and only if $\sum p_n < \infty$.
Then one writes
$$p_n = \frac{1}{(2 \pi)^d} \int_{-\pi}^{\pi} \cdots \int_{-\pi}^{\pi} \left( \frac{1}{d} \sum \cos \theta_j \right)^n d \theta_1 \cdots d \theta_n$$
and approximates the integral by
$$\frac{1}{(2 \pi)^d} \int_{- \infty}^{\infty} \cdots \int_{- \infty}^{\infty} \exp\left(- \frac{n}{2} \sum \theta_i^2 \right) d \theta_1 \cdots d \theta_n$$
to show that $p_n \sim \frac{c}{n^{d/2}}$, so $\sum p_n$ converges for $d>3$ and diverges for $d=1$, $d=2$.
A: Here's the beginning of the proof using CLT and Poisson random variables, although I don't know how to justify the last step. 
Let $X_i$ be iid Poisson random variables with parameter $1$. Then $S_n = \sum_{i=1}^n X_i$ is Poisson with parameter $n$, and hence
$$\mathbb{P} \left( S_n = n \right) = \frac{n^n e^{-n}}{n!}.$$
On the other hand, by the central limit theorem, $\frac{S_n - n}{\sqrt{n}}$ is asymptotically normal with variance $1$, and hence for $\varepsilon > 0$, 
$$\mathbb{P} \left( \left| \frac{S_n - n}{\sqrt{n}} \right| \le \varepsilon \right) \approx \int_{-\varepsilon}^{\varepsilon} \frac{1}{\sqrt{2\pi}} e^{- \frac{x^2}{2}} \, dx = \frac{2 \varepsilon}{\sqrt{2\pi}} + O(\varepsilon^3).$$
At this point I'd really like to set $\varepsilon = \frac{1}{2 \sqrt{n}}$, but I don't think I can justify this; the LHS doesn't change if I instead use $\varepsilon = \frac{1}{3 \sqrt{n}}$, for example, but the RHS does. That suggests that the CLT isn't a strong enough result to allow me to pick $\varepsilon$ to depend on $n$, so it seems I need to finish with some limiting argument as $\varepsilon \to 0$ instead. Help! 
A: Here is a completely elementary way.
You can start by summing,
$$
\sum_{j \leq n} \log j
$$
using Euler-Maclaurin. This will reveal that 
$$
n! \sim C \sqrt{n} \Big ( \frac{n}{e} \Big )^{n}
$$
and so it remains to determine the value of the constant $C > 0$. 
One way to go about this is to consider a sum of $n$ i.i.d random
variables $X_i$ which are $0$ or $1$ with equal probability $\tfrac 12$.
Then, since $X_j$ and $1 - X_j$ have the same distribution, 
$$
\mathbb{P} \Big ( \sum_{j = 1}^{n} X_j \geq \frac{n}{2} \Big ) =
\mathbb{P} \Big ( \sum_{j = 1}^{n} (1 - X_j) \geq \frac{n}{2} \Big )
$$
which implies that 
$$
\mathbb{P} \Big ( \sum_{j = 1}^{n} X_j \geq \frac{n}{2} \Big ) 
= \mathbb{P} \Big ( \sum_{j=1}^{n} X_j < \frac{n}{2} \Big )
$$
and therefore (since the sum of the LHS and RHS is 1),
$$
\mathbb{P} \Big ( \sum_{j = 1}^{n} X_j \geq \frac{n}{2} \Big ) = \frac{1}{2}.
$$
On the other hand, we also know that
$$
\mathbb{P} \Big ( \sum_{j = 1}^n X_j \geq \frac{n}{2} \Big )= \frac{1}{2^n}
\sum_{j \geq \frac{n}{2}} \binom{n}{j}
$$
Now using our ``weak Stirling approximation'' we find that, 
when $j = n/2 + \kappa \sqrt{n/2}$ then,
$$
\binom{n}{j} \sim \frac{2^n}{C \sqrt{n}} \cdot e^{-\kappa^2 / 2}
$$
Inserting this into the previous equation and summing we get
$$
\mathbb{P} \Big ( \sum_{j = 1}^n X_j \geq \frac{n}{2} \Big )
\sim \frac{1}{C} \int_{0}^{\infty} e^{-x^2/2} dx
$$
as $n \rightarrow \infty$. Therefore combining this with an earlier equation
it follows that
$$
\int_{0}^{\infty} e^{-x^2/2} dx = \frac{C}{2}
$$
We know by the classical trick of squaring and using polar co-ordinates that
$$
\int_{0}^{\infty} e^{-x^2/2} dx = \frac{\sqrt{2\pi}}{2}
$$
and therefore $C = \sqrt{2\pi}$. This determines finally the constant in Stirling's approximation.
P.S: Note that one can also completely eliminate the probabilistic language in this approach, and just use binomial coefficients. 
