Es una mezcla de tres distribuciones, y se puede encontrar bastante fácilmente por fuerza bruta, si uno se permite pasar por alto algunos detalles importantes (por ejemplo, "es $\lambda < \mu$ ").
Dejemos que $n$ sea el número de clientes del sistema; $n=0$ significa que nadie está siendo atendido o esperando, $n=1$ significa que se está atendiendo a un cliente pero no hay nadie esperando, etc. Dejemos que $p_n$ sea la probabilidad en estado estacionario (que se supone de aquí en adelante) de que $n$ los clientes están en el sistema.
Evidentemente, si $n=0$ el tiempo $t$ hasta el siguiente inicio del servicio se distribuye $f(t|n=0) = \text{Exp}\{\lambda\}$ y si $n \geq 2$ el tiempo hasta el próximo inicio del servicio se distribuye $f(t|n \geq 2) = \text{Exp}\{\mu\}$ ya que el siguiente servicio se iniciará inmediatamente después de la finalización del actual. Si $n=1$ el tiempo hasta el siguiente inicio del servicio es el máximo del tiempo hasta la siguiente llegada y el tiempo hasta la siguiente finalización del servicio. Esta última distribución tiene la siguiente forma fácilmente derivada:
$f(t|n=1) = \lambda e^{-\lambda t} + \mu e^{-\mu t} - (\lambda + \mu)e^{-(\lambda+\mu)t}$
Las probabilidades de la mezcla corresponden a $p_0$ , $1-p_0-p_1$ y $p_1$ respectivamente. Las probabilidades $p_0$ , $p_1$ y $1-p_0-p_1$ se puede escribir como:
$$p_0 = \left[\sum_{k=0}^{\infty}\left({\lambda \over{\mu}}\right)^k\right]^{-1} = 1 - {\lambda\over{\mu}}$$ $$p_1 = \left({\lambda \over{\mu}}\right)p_0 = {\lambda\over{\mu}} - \left({\lambda\over{\mu}}\right)^2$$ $$1-p_0-p_1 = \left({\lambda\over{\mu}}\right)^2$$
Escribirlo todo, con algún reordenamiento de los términos, da:
$$f(t) = {\lambda(\mu-\lambda)\over{\mu}}\left(e^{-\lambda t} + e^{-\mu t}\right) + \left({\lambda\over{\mu}}\right)^2\left(\lambda e^{-\lambda t} + \mu e^{-\mu t} - (\lambda+\mu)e^{-(\lambda+\mu)t}\right)$$
Mi fuente de fórmulas de probabilidad fue Kleinrock, Sistemas de Colas .
Edición: La derivación de $f(t|n=1)$ es la siguiente, escrita como la derivación del máximo de dos variantes exponenciales independientes $x \sim \text{Exp}\{\lambda\}$ y $y \sim \text{Exp}\{\mu\}$ . Las FDA correspondientes son $F_X(x) = 1-\exp{\{-\lambda x\}}$ y $F_Y(y) = 1-\exp{\{-\mu y\}}$ .
Para ello, utilizaremos la "técnica de la función de distribución acumulativa". Obsérvese en primer lugar que la afirmación " $\max(x,y) \leq t$ " equivale a " $x \leq t$ y $y \leq t$ ". La probabilidad de que $\max(x,y) \leq t$ es sólo el producto de las probabilidades de que $x \leq t$ y $y \leq t$ (como $x$ y $y$ son independientes). Escribiendo esto se obtiene:
$$F_{\max(x,y)}(t) = \left(1-e^{-\lambda t}\right)\left(1-e^{-\mu t}\right) = 1 - e^{-\lambda t} - e^{-\mu t} + e^{-(\lambda+\mu) t}$$
y tomando la derivada con respecto a $t$ te lleva a la función de densidad.
