In this article, we are going to derive the Poisson distribution from the binomial distribution. If you don't know what the binomial distribution is or how to derive it from "scratch", I wrote an introduction to it in a previous article. As for the Poisson distribution, it is, obviously, also a discrete probability distribution whose shtick is quantifying the likelihood of $k$ events happening in some interval provided you know the average number of such events happening in said interval and that each event is independent of the previous ones. In other words, their timing doesn't affect the timing of the subsequent ones.
One dark example of where the timing of a single event could be considered to affect the timing of future events are passenger aircraft crashes, even though it is often given as an example problem when discussing the Poisson distribution. However, in reality, such events cause global repercussions that affect airlines, plane manufacturers and regulatory bodies in considerable and unpredictable ways. Fleets of aircrafts get grounded, corrections get issued and the issue gets resolved, all of which affects the likelihood of future crashes. If one is pressed for an aviation related problem, a better example that actually obeys the Poisson distribution would be the number of engine shutdowns on an aircraft over a large number of flight hours or cycles.
The boilerplate example is estimating the likelihood of a certain number of floods given statistics accumulated over, say, 10 years. Deviating from the mean on either side will, obviously, adversely affect the likelihood of that number of floods happening during a period of 10 years. All we have to do now is find a way to quantify this concept.
Sending $n$ to $\infty$
As mentioned, we start with the binomial distribution. It tells us the likelihood of having $k$ successes out of $n$ trials. The Poisson distribution tells us the likelihood of $k$ events happening given the average of such events happening during some interval (which is implicit). We know the expected value or the mean of the binomial distribution is equal to $nq$ where n is the number of trials and $q$ is the probability of success of a single trial. The expected value makes sense; the more trials we conduct and the higher the probability of success, the more successes we ought to expect. Thus, another way of expressing the binomial distribution by aliasing $\mu$ as $\lambda$ is:
$$ p_q(k, n) = {n \choose k} q^k(1-q)^{n-k} = {n \choose k} {\left(\frac{\lambda}{n}\right)}^k\left(1-\frac{\lambda}{n}\right)^{n-k} $$
We don't know what the probability $q$ of any one single trial is, but we know an "overall trend" $\lambda$ which we effectively take to be determined over an $\infty$ number of trials, hence to express the Poisson distribution, we allow the number of trials of $n$ go to infinity:
$$ p(k) = \lim_{n\to\infty} p_q(k, n) = \lim_{n\to\infty} {n \choose k} {\left(\frac{\lambda}{n}\right)}^k\left(1-\frac{\lambda}{n}\right)^{n-k} = \lim_{n\to\infty} \frac{n!}{(n-k)!k!} {\left(\frac{\lambda}{n}\right)}^k\left(1-\frac{\lambda}{n}\right)^{n-k} $$
Taking the limit here is not difficult, as we'll show in detail. What creates difficulty for people is the fact that $\lambda = nq$ instigates an indeterminate form if substituted in. However, we make the legitimate assumption that as $n$ goes to infinity, $q$ tends to zero. In other words, we assume $\lambda$ remains constant.
We start off slow, by moving the constant factors outside:
$$ p(k) = \frac{\lambda^k}{k!} \lim_{n\to\infty} \frac{n!}{(n-k)!} {\left(\frac{1}{n}\right)}^k\left(1-\frac{\lambda}{n}\right)^{n-k} $$
Notice that the equation can be rewritten as follows:
$$ p(k) = \frac{\lambda^k}{k!} \lim_{n\to\infty} \left(1-\frac{\lambda}{n}\right)^{n-k} \prod_{i=0}^k { \frac{n - i}{n} } = \frac{\lambda^k}{k!} \lim_{n\to\infty} \left(1-\frac{\lambda}{n}\right)^{n-k} \prod_{i=0}^{k-1} {\left( 1 - \frac{i}{n} \right)}$$
If you're having trouble seeing where that product came from, here it is expanded:
$$ \frac{n!}{(n-k)!} \frac{1}{n^k} = \frac{n}{n} \cdot \frac{n - 1}{n} \cdot \frac{n - 2}{n} \cdots \frac{n - k + 1}{n} = \prod_{i=0}^{k-1} { \frac{n - i}{n} } $$
It is easy to see that the product approaches $1$ when we allow $n$ to tend to infinity, hence:
$$ \lim_{n\to\infty} \prod_{i=0}^{k-1} { \left( 1 - \frac{i}{n} \right) } = 1 \implies p(k) = \frac{\lambda^k}{k!} \lim_{n\to\infty} { \left(1-\frac{\lambda}{n}\right)^{n-k} }$$
We rewrite the expression further chasing simplification:
$$ p(k) = \frac{\lambda^k}{k!} \lim_{n\to\infty} \left(1-\frac{\lambda}{n}\right)^{n} \left(1-\frac{\lambda}{n}\right)^{-k}$$
It's immediately obvious that:
$$ \lim_{n\to\infty} \left(1-\frac{\lambda}{n}\right)^{-k} = 1 \implies p(k) = \frac{\lambda^k}{k!} \lim_{n\to\infty} \left(1-\frac{\lambda}{n}\right)^{n}$$
The remaining factor might look familiar to the way we define $e$ (taking the limit of $(1 + 1/x)^x$ as $x$ goes to infinity). It's not quite there yet, but it stands to reason we can transform it in a way that it can be expressed as a power of $e$:
$$ \lim_{n\to\infty} \left(1-\frac{\lambda}{n}\right)^{n} = \lim_{m\to\infty} \left(\left(1+\frac{1}{m}\right)^{m}\right)^{-\lambda} = e^{-\lambda} \iff m = -\frac{n}{\lambda} $$
Which resolves our distribution to its final form:
$$\Large { p(k) = \frac{\lambda^k}{k!} e^{-\lambda} }$$
While it might not surprise you that the the expected value of this probability distribution is $\lambda$, its variance is also $\lambda$. The moment $e$ crops up in your model, its properties are bound to be interesting.