0
$\begingroup$

I have a real world scenario that involves a process behaving as follows: there are two kinds of machines and when these machines fail, their recoveries follow exponential distributions. Let's say the first kind has rate $\lambda$ of recovering and the second has rate $\mu$. Also, there is some censor level in the data collection process where we don't see the data if it $>c$. If the data is un-censored, we know what kind it is. If it is censored, we lose this information. It could have been any of the two kinds. It's easy to generate data from such a process. We just generate some samples from an exponential with rate $\lambda$ and collect the ones below $c$ in a vector (call it $t$) and the ones above $c$ in another vector (call it $x$). Then we generate some samples from an exponential with rate $\mu$ and collect the ones below $c$ into a vector $s$ and append the ones above $c$ into $x$. Now we have three vectors, $t$, $s$ and $x$. We should be able to devise a censored mixture model and frame a maximum likelihood scheme that allows us to retrieve the values of $\lambda$, $\mu$ and $p$ (the proportion of samples generated from $exp(\lambda)$) from the data, $t$, $s$ and $x$.

Here is my attempt at this. I think it's 66% working but there is one thorn in the side.

The likelihood function will be given by:

$$L(\lambda, \mu, p) = \left(\prod\limits_{s=s_1}^{s_S} \mu e^{-\mu s}\right) \left(\prod\limits_{t=t_1}^{t_T} \lambda e^{-\lambda t}\right)\left(\prod\limits_{x=x_1}^{x_X} (p e^{-\mu x} + (1-p)e^{-\lambda x})\right)$$

Where $S$, $T$ and $X$ are the lengths of vectors $s$, $t$ and $x$. Also, we consistently censor the data at $c$. So, the entire vector $x$ is just a bunch of $c$'s. This makes the expression above:

$$L(\lambda, \mu, p) = \left(\prod\limits_{s=s_1}^{s_S} \mu e^{-\mu s}\right) \left(\prod\limits_{t=t_1}^{t_T} \lambda e^{-\lambda t}\right)\left( (p e^{-\mu c} + (1-p)e^{-\lambda c})\right)^c$$

Taking logarithm on both sides, we get the log-likelihood function:

$$l(\mu, \lambda, p) = \sum\limits_{s}\mu e^{-\mu s} + \sum\limits_{t}\lambda e^{-\lambda t} + c\log(p e^{-\mu c} + (1-p)e^{-\lambda c})$$

The derivatives w.r.t. $\mu$ and $\lambda$ are kind of making sense. However, taking derivative with respect to $p$ we get:

$$\frac{\partial l}{\partial p} = c\frac{e^{-\mu c} - e^{-\lambda c}}{(p e^{-\mu c} + (1-p)e^{-\lambda c})}$$

If we want to set this equal to $0$, we would need:

$$e^{-\lambda c} = e^{-\mu c}$$

However, this is impossible. What might I be missing here?

See code that generates data and calculates gradients with respect to parameters here: https://github.com/ryu577/survival/blob/master/distributions/expmixture.py

$\endgroup$
2
$\begingroup$

Interesting problem! The first issue is how we decide to model this situation, we can consider our observations to be drawn from bivariate distribution $(T,C)$ with a continuous and a discrete component. That is, we have $n$ iid observations $(T_i, C_i)$ defined in the following way,

$$T_i = B_i X_i + (1-B_i)Y_i $$

$$C_i = \begin{cases}B_i, \ \ T_i < c \\ 2, \ \ T_i \geq 2 \end{cases} $$

Where $X_i \sim \text{Exp} (\lambda)$, $Y_i \sim \text{Exp} (\mu)$ and $B_i \sim \text{Bernoulli}(p)$ and they are all independent in the relevant ways. $X_i, Y_i, B_i$ are also all hidden from us, $X_i$ is the output from the first machine and is only given to us when $B_i = 0$ and $Y_i$ is the output from the second machine and is only given to us when $B_i = 1$. But when the output is greater than $c$, $C_i$ censors which machine it came from and so it outputs $2$ as we observe, otherwise we'd observe whatever $B_i$ is.

At this point instead of trying to find the density, we can use an Expectation Maximisation (EM) algorithm to estimate the paramters since this would be a classic application of such since our observations depend on hidden variables.

If we want to find the density instead we need to obtain the joint density. However it is likely that it will not be amenable to analytic solution for the MLE, I calculated the joint density by calculating,

$$f_{T,C}(t,j) = \frac{\text{d}}{\text{d}x} \mathbb{P}[ T_i \leq x, \ C_i = j]$$

We find our density to be,

$$f_{T,C}(t, j) = \begin{cases}\mu (1-p) e^{-\mu x} \cdot \mathbf{1}_{x \leq c}, \ \ j = 0 \\ \lambda p e^{-\lambda x} \cdot \mathbf{1}_{x \leq c}, \ \ j = 1 \\ (p \lambda e^{-\lambda x} + (1-p)\mu e^{-\mu x}) \cdot \mathbf{1}_{x > c}, \ \ j = 2 \end{cases} $$

This seems like it would be impossible to find an analytic solution, so instead we should use an EM algorithm to estimate the parameters of interest. I'll update this response with an explanation soon.

$\endgroup$
  • $\begingroup$ Thanks, looking forward to your approach with EM. In the meantime, I added a formulation which seems to work for $\lambda$ and $\mu$ but doesn't make sense for $p$. Added the edit to the question. $\endgroup$ – Rohit Pandey Feb 17 '19 at 3:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.