# Parameter estimation for sum of independent exponential rv

I have a process which I suspect involves two exponentially distributed waiting times of different rates. That is, I am measuring $T = T_\lambda + T_\mu$ with $T_\lambda \sim Exp(\lambda)$ and $T_\mu \sim Exp(\mu)$. Writing down the distribution for $T$ is not problem:

$$p(T | \mu, \lambda) = \lambda \mu \int_0^t e^{-\mu \tau} e^{-\lambda (t - \tau)} d\tau = \frac{\lambda \mu}{\lambda - \mu} \left(e^{-\mu t} - e^{-\lambda t} \right)$$

So far so good. The problem is estimating the parameters, $\lambda$ and $\mu$.

If $\lambda = \mu$, then $p(T|\mu,\lambda) = \Gamma(T|\alpha = 2, \beta = \lambda)$, and the problem is easy, but I have no reason to believe this is the case.

It's also the case that $p(T|\mu,\lambda) = p(T|\lambda, \mu)$, so that the posterior distribution $p(\lambda, \mu | \text{Data})$ will be symmetric along $\mu = \lambda$. I could pick, say, $\mu = \lambda + \epsilon$, where $\epsilon > 0$ to remove this redundancy.

In any case, here's what I attempted:

1. Find the posterior distribution: Hope springs eternal, so I tried to actually do the integral of the likelihood times the prior. I'm using the improper priors $\lambda \sim \lambda^{-1}$ and $\epsilon \sim \epsilon^{-1}$. Maybe there's a nice conjugate prior somewhere out there, but I haven't found it.

$$Z(\{t_1, \dots, t_N\}) = \int_0^\infty\int_0^\infty \frac{\lambda^{N-1} (\lambda + \epsilon)^{N}}{\epsilon^{N+1}} e^{-\lambda \sum t_i} \prod_{i = 1}^N \left(1 - e^{-\epsilon t_i} \right) d\lambda d\epsilon$$

This is an absolute mess, due to that nasty product on the right. You can integrate out $\lambda$ first, and you get

$$\epsilon^{-(N+1)}\prod_{i = 1}^N \left(1 - e^{-\epsilon t_i} \right) \sum_{k=1}^N {N \choose k} \epsilon^k \int_0^\infty \lambda^{2N - k -1} e^{-\lambda \sum t_i} d \lambda = \prod_{i = 1}^N \left(1 - e^{-\epsilon t_i} \right) \sum_{k=1}^N {N \choose k} \epsilon^{k - N - 1} \frac{(2N - k - 1)!}{(\sum t_i)^{2N - k - 1} }$$

...okay. The product can be expanded out. It is an alternating sum:

$$\prod_{i = 1}^N \left(1 - e^{-\epsilon t_i} \right) = 1 - \sum_{j=1}^{2^N}z_j e^{a_j \epsilon}$$

where the set $\{a_j\}$ is the set of all $2^N$ partial sums of $\{t_i\}$, and $z_j = +1$ if $a_j$ is a sum of an even number of elements, or $-1$ otherwise. In any case, what this means is that the integral over $\epsilon$ should be

$$\sum_{k=0}^N {N \choose k} \frac{(2N - k - 1)!}{(\sum t_i)^{2N - k - 1} } \int_0^\infty \epsilon^{k-N-1} \left(1 - \sum_{j=1}^{2^N}z_j e^{a_j \epsilon} \right) d \epsilon$$

And here I get stuck. The exponent of $\epsilon$, $(k - N - 1) \in [-(N + 1), -1]$. It's a sum of diverging integrals. My expectation is that these divergences somehow "cancel out", since it's an alternating sum, and by performing the integral from $\delta$ to $\infty$, and then taking the limit $\delta \to 0$, things would work out. I have not done this because it seemed daunting and because the best case scenario is a messy analytical formula that involves $2^N$ partial sums.

1. MAP or MLE: The problem is the same with both. I'll do MLE as example:

$$\ln p(t|\lambda, \mu) = N\ln \lambda + N\ln \mu - N\ln (\lambda -\mu) + \sum_{i=0}^N\ln(e^{-t_i \mu} - e^{-t_i \lambda}).$$

Deriving with respect to $\lambda$ and setting to zero, we are left with:

$$\frac{N}{\lambda} - \frac{N}{\lambda - \mu} + \sum_{i=0}^N \frac{t_i e^{-t_i \lambda}}{e^{-t_i \mu} - e^{-t_i \lambda}} = 0$$

which I suspect does not have a closed solution for $\lambda$.

1. MCMC: My last resort. I tried to use PyMC to get samples from the posterior distribution. I have had trouble with this because the log likelihood involves the log difference of exponentials. This throws underflow errors, and results in nonsense. I know this is a known problem when calculating the log sum of exponentials, and there are functions like scipy.misc.logsumexp that get around this, but I am using the difference of exponentials, not the sum, and I don't know of any function that gets around that.
• I'd have reparameterized to $\mu = k\lambda$; not only would the algebra be a bit neater but that reparamaterization has a more useful interpretation (being the ratio of two rates) – Glen_b Feb 25 '17 at 5:43
• Have you checked out this paper eajournals.org/wp-content/uploads/… which appears solving exactly the problem you have? The keyword is "convolution". – Vladislavs Dovgalecs Feb 26 '17 at 7:03
• This answer might also help math.stackexchange.com/questions/474775/… – Vladislavs Dovgalecs Feb 26 '17 at 7:05
• It seems unlikely you can use improper priors on both $\lambda$ and $\mu$ and observe only $T_1+T_2$. – Xi'an Feb 27 '17 at 14:48
• A natural MCMC solution would include the missing $T_1$ and $T_2$ inside the simulation. – Xi'an Feb 27 '17 at 14:48