Estimating parameters of mixture of 2 exponential random variables (ideally in Python) Imagine a simulation experiment in which the output is n total numbers, where k of them are sampled from an exponential random variable with rate a and n-k are sampled from an exponential random variable with rate b. The constraints are that 0 < a ≤ b and 0 ≤ k ≤ n, but a, b, and k are all unknown. Further, because of details of the simulation experiment, when a << b, k ≈ 0, and when a = b, k ≈ n/2.
My goal is to estimate either a or b (don't care about k, and I don't need to estimate both a and b: just one of the two is fine). From speculation, it seems as though estimating just b might be the easiest path (when a << b, there is pretty much nothing to use to estimate a and plenty to estimate b, and when a = b, both there is still plenty to estimate b). I want to do it in Python ideally, but I am open to any free software. If someone could point me to an algorithm/statistical approach that handles this, I could always implement it myself, but an existing implementation would be nice.
My first approach was to optimize the following log-likelihood function:
$$
\log{L} = \sum_{i=1}^{n}{\max{\left( P(X=x_i|\lambda=a),P(X=x_i|\lambda=b) \right)}}
$$
This unfortunately didn't work very well. For some selections of the parameters, it's within an order of magnitude, but for others, it's absurdly off. Given my problem (with its constraints) and my goal of estimating the larger parameter of the two exponentials (without caring about the smaller parameter nor the number of points that came from either), any ideas?
 A: You can use an EM algorithm. Let $\kappa$ be the proportion of random variables from the distribution with rate parameter $a$. Initialize $a_0$, $b_0$ and $\kappa_0$. Iterate between the following steps for $T$ iterations:
E-Step:
For $i \in 1,...,n$ set:
$$
\nu_i^{t} = \frac{a_{t-1} e^{-a_{t-1} x_i}\kappa_{t-1}}{a_{t-1} e^{-a_{t-1} x_i}\kappa_{t-1} + b_{t-1}e^{-b_{t-1}x_i}(1-\kappa_{t-1})}
$$
M-Step:
Set:
$$
\kappa_t = \frac{1}{n}\sum_{i=1}^{n} \nu_i^t \qquad
a_t = \frac{\sum_{i=1}^{n} \nu_i^{t}}{\sum_{i=1}^{n}\nu_i^{t}x_i} \qquad
b_t = \frac{\sum_{i=1}^{n} (1-\nu_i^{t})}{\sum_{i=1}^{n}(1-\nu_i^t)x_i}
$$
You can probably improve your estimate of $\kappa$ by using your knowledge regarding the relationship between it and the rate parameters.
A: I tried the proposed EM algorithm, and it worked quite well!
I also tried the following, which worked quite well for a wide range of a and b (and a/b) that I tested:
First, for every single integer percentile (1st percentile, 2nd percentile, ..., 99th percentile), I compute the estimate of b using the quantile closed-form equation (where the i-th quantile is the (i *100)-th percentile) for an exponential distribution:
$$q_i = -\frac{\ln(1 − i)}{b}$$
so:
$$b = -\frac{\ln(1 − i)}{q_i}$$
The result is a list where each i-th element corresponds to the b estimate using the (i+1)-th percentile.
Then, I perform peak-calling on this list using the Python implementation of the Matlab peak-calling function. Then, I take the list of resulting peaks and return the minimum. It seems to work fairly well.
