Here is a proof based on the discussion under David E Speyer's answer.
It uses the same method also to sharpen the proof of the lower bound. 

We will give a simple proof for the bound $$
 (n/e)^n\sqrt{2\pi n}
 \le
n! \le (n/e)^n(\sqrt{2\pi n}+1).
$$

We use the integral representation of $n!$:
\begin{align*}
    n!
    &= \int_0^\infty x^n e^{-x}
    %= \int_0^\infty e^{n \log(x) - x}
    = n\int_0^\infty e^{n \log(nx) - nx} dx
    %= n n^{n}\int_0^\infty e^{n (\log(x) - x)} dx
    = n (n/e)^{n}\int_{-1}^\infty e^{n (\log(1+x) - x)} dx.
\end{align*}
Note that $\log(1+x)-x$ is maximized at $x=0$.
Inspired by the series expansion $\log(1+x)-x=-x^2/2 + O(x^3)$ we make a change of variables,
$\log(1+x)-x=-z^2/2$, or:
$$
z = \text{sign}(x)\sqrt{2(x-\log(1+x))}.
$$
The differentials turn out to be pretty simple:
$$
d(-z^2/2)=
-zdz = d(\log(1+x)-x) = \frac{dx}{1+x} - dx = -\frac{x}{1+x}dx,
$$
which after adjusting the integration limits yield:
$$
    n!
    = (n/e)^{n} n\int_{-\infty}^\infty e^{-nz^2} (z+z/x) dz.
$$
We now need the bounds
\begin{align}
    1+2z/3
    &\le z+z/x \le \max\{1, 1+z\}
    \quad\text{for all }x\in\mathbb R.
    &&(0)
\end{align}
Since $\int_{-\infty}^\infty e^{-nz^2}dz=\sqrt{2\pi/n}$ and $\int_{0}^\infty e^{-nz^2} z\, dz=1/n$,
this immediately gives us
$$
\sqrt{2\pi n}
\le
n\int_{-\infty}^\infty e^{-nz^2} (z+z/x) dz
\le
\sqrt{2\pi n} + 1.
$$
Which is what we want.

It remains to show eq. (0),
which we will do using the following [standard logarithmic inequalities][1]:
\begin{align}
    \log(1+x) &\ge x(2+x)/(2+2x) \quad &\text{when } x<0
    &&(1)
    \\
    \log(1+x) &\ge x-x^2/2 \quad &\text{when } x>0
    &&(2)
    \\
    \log(1+x) &\le \frac{x(6+x)}{6+4x}
     \quad &\text{when } x>-1
    &&(3)
\end{align}
We note $z(1+1/x)$ is non-negative, since $1+1/x<0$ implies $\mathrm{sign}(x)=-1$.
Thus it suffices to show $z^2(1+1/x)^2\le 1$, which follows using eq. (1):
$$
z^2 = 2(x-\log(1+x))
\le \frac{x^2}{1+x}
\le \frac{x^2}{(1+x)^2}
= (1+1/x)^{-2},
$$
For $x>0$ we can use eq. (2) to bound
$
z = \sqrt{2(x-\log(1+x))}
\le \sqrt{2(x^2/2)} = x,
$
so $z/x\le 1$.
This proves the upper bound of (0).
For the lower bound, it suffices to show
$z^2(1/3+1/x)^2 \ge 1.$
Using eq. (2) we get
\begin{align*}
z^2(1/3+1/x)^2
%&= 2(x - \log(1+x))(1/3+1/x)^2
&\ge 2(x - x(6+x)/(6+4x))(1/3+1/x)^2
\\&= \frac{x^2}{9 + 6 x} + 1
\quad\ge 1 \quad \text{since } x\ge -1.
\end{align*}

---

The full Stirling approximation can similarly be derived from the series expansion:
$$
z(1+1/x) = 1 + \frac{2z}{3} + \frac{z^2}{12}
- \frac{2z^3}{135} + \frac{z^4}{864} + \dots
$$
This form is nice as it allows easy integration with respect to the  Gaussian pdf.
E.g.
$$
n\int_{-\infty}^{\infty}\frac{z^4}{864} e^{-n z^2/2} dz = \frac{\sqrt{2 \pi n}}{288n^{2}}
$$
matching the expected
$$
n! \sim \sqrt{2\pi n}\left(\frac{n}{e}\right)^n \left(1 +\frac{1}{12n}+\frac{1}{288n^2} - \frac{139}{51840n^3} -\frac{571}{2488320n^4}+ \cdots \right).
$$

I think one can even give a uniform bound of
$$
z(1+1/x) \le 1 + \frac{2z}{3} + \frac{z^2}{c},
$$
for some $c\approx 8$, strengthening the upper bound even more.
However, the logarithmic inequalities become quite tedious.

  [1]: https://en.wikipedia.org/wiki/List_of_logarithmic_identities#Inequalities