FIFO (M/D/1/N queue) Overflow probability I'm an engineer, trying to figure out the optimal FIFO depth for a particular application. The events I usually work with follow a Poisson distribution, and a fixed readout rate (thus the M/D/1/N naming).
In an example, events arrive with an expected rate of 1MHz, and they are to be read out with a fixed rate of 1Mevents/s. If the incoming events had a fixed timing distribution (1 event every 1us) there would be no need for event storage at all, but the Poissonian nature of the process dictates that it's possible to receive multiple events in any 1-us interval: it is necessary to store (buffer) these excess events (where the local event rate exceeds 1MHz) for later readout.
In this example, data is lost when the FIFO is full and I have a new event that needs to be written.
Until now, I've designed FIFOs by using simulations, and adjusting their size according to the data loss that I experience. But I wanted to be able to make educated guesses to the appropriate FIFO size beforehand.
I've found the following resources online, which however turned out to be of little help to me, given my poor knowledge of the topic:

*

*https://en.wikipedia.org/wiki/M/D/1_queue#Stationary_distribution

*https://www.physicsforums.com/threads/m-d-1-n-queueing-theory-help.600325/
Most often posts online focus on average waiting time or busy period, which I don't think apply to my problem.
What I thought would work
The probability of event loss can be described as (1) the probability that the FIFO is full, and (2) that I receive another event before the readout.
$P_{overflow} = P_{full} * P_{1\,event\,in\,\mu}$ where $\mu$ is the readout time (1us)
Now, $P_{1 event in \mu}$ is trivial, and given by the PDF of the Poisson distribution.
$P_{full}$... is another matter. I think it can be evaluated as:
$P_{full} = P_{full\,-\,1} * P_{1\,event\,in\,\mu}$. And, thus, recursively:
$P_{full} = \prod_i^{FIFO\,Depth} P_{1\,event\,in\,\mu} * P_{empty}$
Although I don't think this is actually right, and I don't know how to evaluate $P_{empty}$ so... I'm at loss.
Is this the right way to proceed? I'm guessing that the matter is far more complex than I anticipated. Which road should I pursue in studying it?
PS: I would eventually like to embed this in a small FIFO depth calculator in Python, so I can use SciPy primitives.
 A: I'm not sure your particular question has been considered in the literature, but we can easily derive the answer from the first principles. Assume that the readout rate is 1 entry per unit of time and the arrival intensity is $\lambda$ entries per unit of time (so my $\lambda$ is just the ratio of the arrival intensity to the readout rate in general, which is $1$ in the example you posted).
Let $m$ be the buffer length (I include the entry being read into the buffer, so your "no buffer" is my "buffer of length 1"). Let us consider the states of the buffer (the numbers of entries in it) at the moments of ending the readouts of an entry. The possible states are $0,1,\dots,m-1$ (the probability of an arrival exactly at the moment of finishing the readout is $0$, so we ignore it). The key observation is that the resulting process is a discrete Markov chain. We need now to understand its transition probabilities and the expected times for the jumps.
Let's start with the state $0$. It means that the reading device is currently just idle. The expected time until the next arrival is $1/\lambda$, after which it works for 1 unit of time until the next readout completion. So the total expected time of transition to the next state in this case is $1+\frac 1\lambda$. The probability $p_{0,k}$ to make a transition to the state $k$ for $k<m-1$ is exactly the probability to have $k$ arrivals within 1 unit of time after the first one, which is $\frac{\lambda^k}{k!}e^{-\lambda}$. The value $p_{0,m-1}$ can be obtained from making the sum of transition probabilities exactly $1$, but we'll not need it.
For any other state $k\ge 1$, the next readout starts right away, so the expected (actually, deterministic) time for the transition is $1$ in this case. The probability to jump to the state $k-1+\ell$ is that of exactly $\ell$ arrivals during the reading time, i.e., $\frac{\lambda^\ell}{\ell!}e^{-\lambda}$ as long as $k-1+\ell<m-1$. As before, $p_{k,m-1}$ complements the row sum to $1$.
Thus, the transition matrix is
$$
P=e^{-\lambda}\begin{bmatrix}
1 &\lambda &\lambda^2/2! & \lambda^3/3! &\dots & \lambda^{m-2}/(m-2)! & *
\\
1 &\lambda &\lambda^2/2! & \lambda^3/3! &\dots & \lambda^{m-2}/(m-2)! & *
\\
0 & 1 &\lambda &\lambda^2/2! &\dots & \lambda^{m-3}/(m-3)! & *
\\
0 & 0 & 1 &\lambda &\dots & \lambda^{m-4}/(m-4)! & *
\\
\vdots & \vdots &\vdots &\ddots &\ddots &\vdots &\vdots 
\\
0 & 0 &0 & 0 &\dots & 1 & * 
\end{bmatrix}
$$
(I denoted all entries in the last column by $*$; as we'll see, they won't be needed for the computation).
The loss rate is easy to count as the difference between the arrival rate $\lambda$ and the useful reading rate. To determine the latter notice that if $\pi_0,\dots,\pi_{m-1}$ is the steady state, then out of $t$ transitions on average $\pi_0t$ were from state $0$, so the typical (in any reasonable sense of this word, since the Law of Large Numbers is all powerful) time used for $t$ transitions (readouts) is $t(1+\pi_0\frac 1\lambda)$. Thus, the useful readout rate is $1$ over the factor in parentheses, i.e., $\frac{\lambda}{\lambda+\pi_0}$.
Thus, we have a neat formula
$$
\text{The loss rate (per unit time) }= \lambda\left(1-\frac 1{\lambda+\pi_0}\right)\,.
$$
and, accordingly,
$$
\text{The loss portion (per entry) }= 1-\frac 1{\lambda+\pi_0}\,.
$$
It remains to determine $\pi_0$ from the system $\pi P=\pi$. Note that it is always a degenerate system to which one needs to add one normalizing equation $\sum_k \pi_k=1$. That's why we do not need to bother with the last column in $P$.
From the first column, we get
$$
e^\lambda\pi_0=\pi_0+\pi_1
$$
so $\pi_1=(e^{\lambda}-1)\pi_0$, which makes it tempting to express each $p_k$ as $u_k\pi_0$ with $u_0=1, u_1=e^\lambda-1$. The equations on $u_k$ beyond the first one are
$$
e^\lambda u_{k-1}=u_k+\lambda u_{k-1}+\sum_{\ell=2}^{k-1}\frac{\lambda^\ell}{\ell!}u_{k-\ell}+\frac{\lambda^{k-1}}{(k-1)!}u_{0}
$$
which can be easily solved recurrently as
$$
u_k=(e^{\lambda}-\lambda)u_{k-1}-\sum_{\ell=2}^{k-1}\frac{\lambda^\ell}{\ell!}u_{k-\ell}-\frac{\lambda^{k-1}}{(k-1)!}u_{0}\,.
$$
For the buffer length $m$, we just have
$$
\pi_0=\left(\sum_{k=0}^{m-1}u_k\right)^{-1}\,.
$$
In general, that's it for an arbitrary arrival rate $\lambda$ (you are, probably, mainly interested in $\lambda\le 1$ but the computation works for all values and is quick and stable). However, your case $\lambda=1$ is rather special and deserves a separate treatment.
If we use generating function $f(z)=\sum_{k=0}^\infty u_kz^k$, then our infinite sequence of equations for $u_k$ can be written as
$$
[(u_0+u_1)+u_2z+u_3z^2+\dots]e^{-\lambda}[1+\lambda z+\tfrac{\lambda^2}{2!}z^2+\dots]
\\
=[u_0+u_1z+u_2z^2+\dots]\,,
$$
i.e., ($u_0=1$)
$$
\left(1+\frac{f(z)-1}z\right)e^{\lambda(z-1)}=f(z)\,
$$
whence
$$
f(z)=\frac{1-z}{1-ze^{\lambda(1-z)}}\,.
$$
Taking $\lambda=1$ and denoting $w=1-z$, we get
$$
ze^{1-z}=(1-w)e^w=(1-w)(1+w+w^2/2+w^3/6+\dots)
\\
=1-w^2/2-w^3/3-\dots\,,
$$
so
$$
f(z)=\frac{2}{w+2w^2/3-\dots}=\frac 2w-\frac 43+O(w)\,. 
$$
Since the partial sums of the coefficients of $f$ are the coefficients of $\frac{f(z)}{1-z}$, we conclude that up to a geometrically decaying sequence they are just the coefficients of
$$
\frac 2{(1-z)^2}-\frac 4{3(1-z)}\,,
$$
i.e., $u_0+\dots+u_{m-1}\approx 2[(m-1)+1]-\frac 43=2m-\frac 43$ and
correspondingly the loss rate is close to $\frac 1{2m-\frac 13}$ in this case. The direct computation shows that "close" means "very close for all practical purposes" for $m\ge 5$.
I hope, this helps a bit. Feel free to ask questions if anything is unclear.
