# How to use Gamma distributions to estimate the number of failures?

I need to calculate the expected number of failures of a product within 6 years. The time until failure is said to be gamma distributed with $$\alpha=2$$ and $$\beta=0.5$$ for a mean time between failures of 1 year.

The product will be replaced after each failure, so I need to form a simulation to determine how many times I can expect to replace that product over the course of 6 years. How can I calculate the expected number of failures in 6 years with this information? I will be implementing this in R.

I don't want to sound like I'm demanding an answer. There's much more to the problem in question. I'm only looking for a starting point as I have not used the Gamma distribution for this purpose before and after searching online I couldn't find much help.

• "Estimate" suggests you are basing a calculation on data, but what data do you have? If you are just trying to compute with the Gamma distribution, then would "estimate" mean "approximate" or would it mean "calculate"? Finally, your last paragraph suggests you are familiar with distributions and using distribution functions, so what about the Gamma distribution strikes you as being any different from other distributions? – whuber Jun 11 at 17:28
• Thanks for replying. Poor wording on my part, I've edited my question to "calculate the expected number". I admittedly do not have much experience at all with the gamma distribution and I'm not familiar with its usefulness. I know that the poisson distribution is used to find the number of occurrences within a time period, which is why I'm confused as to why the problem specifies gamma. – Jacob Myer Jun 11 at 17:50
• The smartass answer is that if we are talking about one instance of the product, then when it fails, it fails. With an MTBF of 1 year, you could just say that you expect one failure in 6 years. More seriously, I assume the question is asking for something like the percent of the products that will fail in 6 years. It strikes me that R has a gamma RNG, so you just plug the parameters into the RNG. You can then simulate ... you'll have to fill in the rest. Please clarify if this is not what you need. – Weiwen Ng Jun 11 at 18:03
• Could you explain the meaning of "expected number of failures in 6 years"? Are you told that the Gamma distribution describes the failure times of products in a population, so that (for instance) 99.992% of all the products will have failed by six years; or are you putting a product into service, replacing it each time it fails, and calculating the expected number of those products that will have failed for you during a six year period (which would intuitively be near five to six, since the mean time between failures is one year)? – whuber Jun 11 at 21:15
• @whuber the second situation you described. The product will be replaced after each failure, so I need to form a simulation to determine how many times I can expect to replace that product over the course of 6 years. I am told the time until failure is gamma distributed with a mean of 1 year. – Jacob Myer Jun 11 at 21:41

Let's consider the problem generally so that we can delay the point at which details of Gamma distributions need to be known. I will proceed in three stages: the first develops a general formula, the second exploits properties of Gamma functions (but without doing much calculation), and the third applies these results to the data in the question.

### Relating service times to numbers of failures

At time $$0,$$ a product is put into service. It lasts a random amount of time $$T_1$$ governed by some distribution law $$F$$ (which obviously must assign zero probability to negative times!). At the moment $$0+T_1,$$ if the product is still needed, it is replaced by another. It lasts a random amount of time $$T_2$$--independent of $$T_1$$--governed by the same law $$F.$$ At the moment $$0+T_1+T_2,$$ if the product is still needed this second product is replaced by a third. This process continues until time $$\tau \gt 0,$$ representing how long some product needs to be in service. Let $$N(\tau)$$ be the total number of products put into service between times $$0$$ and $$\tau,$$ so that

$$0 + T_1 + T_2 + \cdots + T_{N(\tau)} \lt \tau.\tag{1}$$

$$N(\tau)$$ is a random variable. $$N(\tau)-1$$ is the number of product failures. What is its expected value, $$E[N(\tau)-1] = E[N(\tau)]-1$$?

Writing $$F_{N(\tau)}$$ for the cumulative distribution of $$N(\tau),$$ a convenient way to find this expectation begins with the formula

$$E[N(\tau)]-1 = -1+\sum_{k=0}^\infty \Pr(N(\tau) \gt k) = \sum_{k=1}^\infty (1 - F_{N(\tau)}(k)).$$

Now, according to definition $$(1)$$ above, the event "$$N(\tau)\le k$$" -- which means $$k$$ or fewer products serially lasted at least for a duration $$\tau$$ -- can also be understood as meaning the total lifetime of $$k$$ products -- whether indeed $$k$$ were actually needed -- was $$\tau$$ or greater. Thus

$$F_{N(\tau)}(k) = \Pr(N(\tau)\le k) = \Pr(0+T_1+T_2+\cdots+T_k \ge\tau).\tag{2}$$

Writing $$S_k = 0 + T_1 + T_2 + \cdots + T_k$$ and $$G_k$$ for its cumulative distribution function,

$$G_k(\tau) = \Pr(0 + T_1 + T_2 + \cdots + T_k \le\tau) = \Pr(S_k \le \tau),$$

it is apparent that the right hand side of $$(2)$$ is closely related to the CDF of $$S_k$$ evaluated at $$\tau,$$

$$\Pr(0+T_1+T_2+\cdots+T_k \ge\tau) = \Pr(S_k \ge \tau) = 1 - \Pr(S_k \lt \tau) = 1 - (G_k(\tau) - \Pr(S_k=\tau)).$$

For a continuous time model--where a product could fail at any moment but has only infinitesimal chance of failing at any given moment--the last term $$\Pr(S_k=\tau)$$ is zero. Together, the previous equations yield

$$E[N(\tau)-1] = \sum_{k=1}^\infty G_k(\tau).\tag{3}$$

This invites us to work out the CDFs $$G_k(\tau)$$ for $$k = 1, 2,$$ etc.

### Using Gamma distributions of service times

That's about as far as the general analysis can go. One reason why people hypothesize Gamma distributions for inter-arrival times is that they enjoy a very nice summation property:

The distribution of the sum of $$k$$ independent random variables with $$\Gamma(\alpha,\beta)$$ distributions has a $$\Gamma(k\alpha,\beta)$$ distribution.

In this notation, $$\alpha$$ is the "shape" parameter and $$\beta$$ merely establishes a scale--a unit of measurement--for the values. (There are two conventions: in one, which I will use, $$\beta$$ is the unit of measurement; in the other, $$1/\beta$$ is the unit. Obviously, switching between the two conventions is easy but when looking things up, be sure to check which convention is in use!)

Consequently, when $$F$$ is a $$\Gamma(\alpha,\beta)$$ distribution, then each $$G_k$$ is a $$\Gamma(k\alpha,\beta)$$ distribution. Because these distributions have density functions, their CDFs can be expressed as integrals:

$$G_k(\tau) = \frac{1}{\Gamma(k\alpha)}\int_0^\tau x^{k\alpha} e^{-x} \frac{\mathrm{d}x}{x}.$$

In this formula I have dropped all references to $$\beta$$ because we can re-introduce it at the end, knowing it's just the unit of measurement of $$\tau$$ (namely, time).

In R, formula $$(3)$$ is readily evaluated by summing values of pgamma. As an example:

#-- Specify the parameters of the problem
tau <- 6
alpha <- 2
beta <- 1/2

#-- Determine how far out to take the sum
u <- (rho + sqrt(6^2 + 4*tau/beta))/2
u <- pmax(1, round(u^2/alpha))

#-- Compute the sum
sum(pgamma(tau, alpha*seq_len(u), scale=beta))


In some cases, we can obtain simpler formulas. The summation in $$(3)$$ can be interchanged with the integration, giving

$$E[N(\tau)-1] = \int_0^\tau \left(\sum_{k=1}^\infty \frac{x^{k\alpha}}{\Gamma(k\alpha)}\right)\ e^{-x} \frac{\mathrm{d}x}{x}.\tag{4}$$

The summation has nice properties (it converges absolutely for all $$x$$) and is relatively straightforward to evaluate numerically.

### Calculations for the situation posed in the question

For small integral values of $$\alpha,$$ it's straightforward to obtain explicit formulas. In the case of the question, with $$\alpha=2,$$ use the facts $$0 = x^n-(-x)^n$$ when $$n$$ is even and $$2x^n=x^n-(-x)^n$$ when $$n$$ is odd to re-express the inner sum of $$(4)$$ as

$$\sum_{k=1}^\infty \frac{x^{k\alpha}}{\Gamma(k\alpha)} = \sum_{k=1}^\infty \frac{x^{2k}}{(2k-1)!} = \frac{x}{2}\sum_{j=0}^\infty \frac{x^{j} - (-x)^{j}}{j!} = \frac{x}{2}\left(e^x - e^{-x}\right).$$

(See my account of decimation for the theory behind this technique.)

Thus, $$(4)$$ becomes

$$E[N(\tau)-1] = \int_0^\tau \left(\frac{x}{2}\left(e^x - e^{-x}\right)\right)\ e^{-x} \frac{\mathrm{d}x}{x} = \frac{1}{2}\left(\tau - \frac{1}{2}\left(1 - e^{-2\tau}\right)\right).$$

Also in the question, $$\beta=1/2$$ years and the duration is $$6$$ years, so using $$\beta$$ as the time units we must write $$\tau = 6/\beta = 12$$ time units. Thus the answer is

$$E[N(12)-1] = \frac{1}{4}\left(23 + e^{-24}\right) \approx 5.750000000009438$$

(to double precision). As far as any simulation is concerned, then, expect the average number of products failing during six years to be near $$5.75.$$