Deriving Stirling's approximation using the Poisson distribution

I discovered a neat way to derive Stirling’s approximation for the factorial in a succinct manner. I came across this neat derivation in Information Theory, Inference and Learning Algorithms by David J.C. MacKay (though I’ve also just noticed that it’s listed on the Wikipedia page too). This is an addendum of sorts to my Factorial Overflow post. Consider an indexed family of Poisson random variables $N_{\lambda }\sim \operatorname{Pn}(\lambda), \lambda \in \mathbb N$ with probability mass functions $$ f_{N_\lambda}(n) = \frac{e^{-\lambda} \lambda^n}{n!}, \quad n \in \{0, 1,\dots \} $$ Now, the Poisson distribution has a divisibility property which means that the sum of i.i.d Poisson random variables is itself a Poisson random variable. So we can decompose $N_\lambda$ as $N_{\lambda}=\sum_{i=1}^\lambda N_i$ where $N_i \sim \operatorname{Pn}(1)$ i.i.d and $\mathbb{E}[N_i]=\operatorname{Var}[N_i]=1$. Thus, by the Central Limit Theorem, as $\lambda \to \infty$: $$…