I have a bunch of machines. These machines fail from time to time, stay down for a while, recover and then commence running.
First, I wanted to model the failures as a Poisson point process. This assumes the machines are homogeneous in terms of their failure rates. In this model, the estimation becomes quite easy. You count the total failures across all the machines, $n$ and the total time the machines collectively ran, $t$. Then, the maximum likelihood estimate for the failure rate of the machines becomes:
$$\lambda = \frac{n}{t}$$
So far so good. But now I want to get rid of the assumption that the machines are homogeneous. So, there is some distribution from which their failure rates are drawn. A Gamma distribution is a natural choice for this distribution since its the conjugate prior. And this is covered in section 5.4.3 of the book, Introduction to probability models by Ross (10th ed). The overall point process becomes NBD or the negative binomial (the number of events in any interval of time is negative binomial instead of Poisson). Now, how do I take my data and fit the parameters of this Negative binomial point process?
I assume getting the total failures and total machine run time per machine will probably serve as a pre-processing step (just as the total failures and total machine time across all machines were used in the Poisson case). Not sure how to make progress from there.
Including relevant equations from the book:
The gamma distribution from which the rate, $\lambda$ of the Poisson is drawn has a PDF:
$$g(\lambda) = \theta e^{-\theta \lambda} \frac{(\theta \lambda)^{m-1}}{(m-1)!}$$
And this results in the PMF of the number of events:
$$P(N(t)=n) = {n+m-1 \choose n}\left(\frac{\theta}{t+\theta}\right)^{m} \left(\frac{t}{t+\theta}\right)^n$$
If for a machine you observe over a span of time k failures, the number of failures in a one-hour span can be modeled by $X\sim \mathscr {NegBinom}\left(t, \frac{t}{k+t}\right)$. As $t\rightarrow \infty$, $X$ will approach the distribution $\mathscr {Pois}\left(k\right)$. For example, if in a one hour span of time you observe $k=5$ failures from a particular machine, you can create the random variable $X\sim \mathscr{NB}\left(3600, \frac{3600}{3605}\right)$ to represent the number of failures observed in a one-hour period. For each of the machines, you will have a negative binomial random with a potentially different parameter of $p=\frac{3600}{k+3600}$, so you will have to figure out what the sum of negative binomial random variables with different probability of success turns out to be.