On deriving the Poisson distribution from the binomial distribution My own personal understanding of modeling the number of independent uniformly distributed events occurring in a fixed time interval, is that we can first approximate this with a binomial distribution on the whole time interval.
$$p(k;n,p)=\binom{n}{k}p^k(1-p)^{n-k}$$
Where $n$ is the total number of possible events and $p$ is the probability that any one event occurs during that interval of time, where we have the expected number of events $\lambda = np$.
Now the problem with a single binomial model on the whole time interval (besides the fact that all events probably won't have the same probability of occurrence) is that any event may occur more than once at different points in time (a mailer may send you multiple pieces of mail at different times during the day).  Thus we form a better approximation by partitioning our fixed time interval into $N$ subintervals and putting a binomial model
$$p(k;n,\frac{p}{N})=\binom{n}{k}\Big(\frac{p}{N}\Big)^k(1-\frac{p}{N})^{n-k}$$
on each of these subintervals.  Where simply dividing $p$ by $N$ is justified since the events are uniformly distributed over the time interval.
Now we add up the distributions on these $N$ subintervals to obtain $p(k;Nn,\frac{p}{N})=p(k;Nn,\frac{\lambda}{Nn})$, and taking $N\rightarrow\infty$, we obtain the Poisson distribution with mean $Nn\frac{\lambda}{Nn}=\lambda$.
The above derivation seems to me to be far more coherent than the one given by the sources I've looked at, such as wikipedia, which all make some vague argument about how very small intervals are likely to contain at most one event (plausible but not rigorous), and where events are treated as entirely fungible, and so on each of the $N$ subintervals we have the distribution $B(k;1,\frac{\lambda}{N})$ which we sum over the subintervals to obtain $B(k;N,\frac{\lambda}{N})$ and then take $N\rightarrow\infty$ to obtain the Poisson distribution with mean $\lambda$.
My question is, is there some reason I'm not seeing for why all the standard sources use the latter argument, rather than the one I've presented?
 A: I believe the reason for the standard derivation as opposed to your own is really quite dissatisfying, one less parameter (no $n$).  It seems like your most substantial frustration with the standard derivation is the "vague argument about how very small intervals are likely to contain at most one event (plausible but not rigorous)."  Let's try to make this rigorous with a few reasonable assumptions.  Also, please note that you've only really alleviated this issue in the sense that you have weakened this assumption to "tiny intervals contain fewer than $n+1$ events."  
Let's start with a simple scenario where we have $m$ events which may or may not take place on the time interval $[0,1]$.  Let's assume event $i$ has probability $p_i$ that it occurs, and P(event $i$ occurs at $x$ | given event $i$ occurs) follows a uniform distribution on the unit interval.  The probability that event $i$ occurs within $\epsilon$ event $j$ is
$$
q_{i,j} <= p_i \cdot p_j * \int_0^1 \int_{y-\epsilon}^{y+\epsilon}1 dxdy = p_i\cdot p_j \cdot \epsilon
$$
(The inequality is just due to the endpoints.  It should be $\epsilon - \epsilon^2$).  Of course the probability that any two events happen within $\epsilon$ of each other is therefore bounded by $\epsilon \cdot \sum_i\sum_{i<j} p_i\cdot p_j$.  Since we have finitely many events the right hand side is bounded, so as $\epsilon \rightarrow 0$, this quantity goes to 0.  
Even if we consider an infinite number of possible events, we must be maintaining the assumption that the expected number of events is equal to 
$$
\lambda = \sum_i p_i < \infty.
$$
Note that the right hand side of the first equation is bounded above by $\epsilon \lambda^2$.  Lastly, let's loosen and/or justify our assumption for the uniform distribution for the timing of an event, conditional on its occurrence.  I believe that someone with more time on their hands could easily demonstrate that we could replace this assumption by the assumption that this conditional distribution has a well defined PDF on $[0,1]$, i.e. there are no point masses.  But more to the nature of the problem, we truly don't care when these events happen within the interval of interest.  We are interested in modelling the number of events.  So making a strong assumption about the nature of that distribution, such as uniformity, is OK in my book.
I hope this was helpful.
