Maximum likelihood estimators based on discrete component lifetimes of a k-out-of-n system

This paper deals with parametric inference about the independent and identically distributed discrete lifetimes of components of a k-out-of-n system. We consider the maximum likelihood estimation assuming that the available data consists of component failure times observed up to and including the moment of the breakdown of the system. First, we provide general conditions for the almost sure existence of a strongly consistent sequence of maximum likelihood estimators (MLE’s). Then, we focus on three typical discrete failure distributions—the Poisson, binomial and negative binomial distributions—and prove that in these cases the MLE’s are unique, provided they exist, and that they are strongly consistent. Finally, we complete our results by Monte Carlo simulation study. Interestingly, the inference considered in the paper can be viewed as equivalent to one based on Type-II right censored discrete data. Therefore, our results can as well be applied to the case when Type-II right censored sample from a discrete distribution is observed.


Introduction
An important class of systems studied in reliability theory is one containing so called k-out-of-n systems. Such systems consist of n elements and work as long as at least k of the elements function. As technical structures having some redundancy, they find various applications in engineering when highly reliable products are needed, for example, they are used in design of servers in internet service or in design of automotive and aeronautic engines. Consequently, they have attracted substantial interest-a vast literature on k-out-of-n systems is available. One stream of this literature concerns inference about component lifetimes based on failure times in a k-out-of-n system or in a sample of k-out-of-n systems. Classic works of Halperin (1952) and Bhattacharyya (1985) describe asymptotic properties of MLE's based on failure times of components of a k-out-of-n system. Generalizations of this results to the case when some of the failure times are censored can be found, among others, in Kong and Fei (1996) and Lin and Balakrishnan (2011). There are also works developing estimation methods for the distribution of components in a system based on a sample of lifetimes of systems; see, for example, Ng et al. (2012), Navarro et al. (2012), Hermanns and Cramer (2018) and the references therein.
All the above-mentioned results, however, concern only the case when the component lifetimes have absolutely continuous distributions. Yet, in some applications the continuity assumption is not adequate. This is the case, for instance, when the system performs a task repetitively and its components have certain probabilities of breakdown upon each cycle or when the component lifetimes represent the numbers of turn-on and switch-off up to failures. While reliability properties of k-out-of-n systems consisting of components with discrete lifetimes have been studied over the years; see Weiss (1962), Young (1970), Tank and Eryilmaz (2015), Dembińska (2018), Dembińska and Goroncy (2020) and Dembińska et al. (2019), to the best of our knowledge results concerning inference about discrete lifetimes of components based on failure times in k-out-of-n systems are not known.
The aim of this paper is to fill in this gap in the literature. We focus on maximum likelihood estimation of an unknown parameter of discrete distribution of component lifetimes of a k-out-of-n system. The estimation is based on failure times of components observed up to and including the system breakdown. In Sect. 2, we set our notation, describe the inference problem under consideration and point out that this problem can be viewed as equivalent to inference from a Type-II right censored sample. Next, in Sect. 3, we present a theorem asserting that under some mild regularity conditions the MLE's of interest exist almost surely for all sufficiently large n and are strongly consistent. The proof of this theorem is postponed to the "Appendix". In Sect. 4, we choose three typical discrete failure distributions-Poisson, binomial and negative binomial-to be the distributions of lifetimes of the components and show that then the MLE's are unique if they exist, their values can be obtained easily by numerical methods and obtained MLE's are strongly consistent. In Sect. 5, we perform Monte Carlo simulation study to investigate finite-sample properties of MLE's discussed in Sect. 4. Section 6 contains an illustrative example based on real failure data while in Sect. 7 we give concluding remarks and problems for future investigations.

Maximum likelihood point estimation
Let F = {F(θ, ·), θ ∈ Θ} be a family of discrete cumulative distribution functions (cdf's), where θ ∈ R is the parameter of interest. Consider a k-out-of-n system which consists of n two-state (i.e., working or failed) components. We assume that the lifetimes of the components, T 1 , T 2 , . . . , T n , are independent and identically distributed (iid) random variables (rv's) with the common cdf F(θ, ·) ∈ F, so that F(θ, t) = P θ (T 1 ≤ t). Next, we denote f (θ, t) = P θ (T 1 = t), i.e., f (θ, ·) is the probability mass function (pmf) corresponding to F(θ, ·), and F(θ, t) = 1 − F(θ, t). Moreover, for simplicity of notation we require that for any θ ∈ Θ the support of F(θ, ·), denoted by supp F(θ, ·), is of the form {0, 1, . . . , M}, where M ≤ ∞. Yet, it is easily seen that the results of Sect. 3 hold more generally in the case when Our aim is to use the maximum-likelihood approach to estimate the unknown parameter θ from the failure data collected up to and including a breakdown of a k-out-of-n system. Let T 1: n ≤ T 2: n ≤ · · · ≤ T n: n denote the order statistics corresponding to T 1 , T 2 , . . . , T n . A k-out-of-n system works as long as at least k of its components work. It fails when the (n − k + 1)th component failure occurs. Thus, the lifetime of k-outof-n system is the (n − k + 1)th smallest of the component lifetimes, i.e., T n−k+1: n . However, in the case of discretely operating elements if k = 1 then at the moment of the system failure we do not necessarily have exactly n − k + 1 inoperative elementsdue to possible ties between component failures with non-zero probability the number of inoperative elements can be larger than n −k +1; see Davies and Dembińska (2019) for details. Therefore, collecting data up to and including a breakdown of a k-out-of-n system we can register not only the values of T 1: n , T 2: n , . . . , T n−k+1: n but also the value of S-the number of failed components at the moment of failure of the system. This means that we observe S, T 1: n , T 2: n , . . . , T n−k+1: n , or equivalently, S, T 1: n , T 2: n , . . . , T S: n .
It is worth pointing out that the results presented in this paper, even though formulated in terms of inference from failure times of components of a k-out-of-n system up to and including its breakdown, can as well be applied to inference based on Type-II right censored discrete data. Indeed, during an experiment in which Type-II right censoring is applied n items with iid lifetimes T 1 , T 2 , . . . , T n are placed on a test. Due to budget or time limitations or on account of ethical decisions in biomedical problems, the experiment is terminated at the moment of the r th failure, where r < n is fixed in advance. If the lifetimes T i , i = 1, . . . , n, are discrete rv's, then with non-zero probability it may happen that at the moment of the r th failure more than r items are broken. Clearly, in order not to lose any information it is reasonable to make the inference not only from the values of T 1: n , T 2: n , . . . , T r : n but to include also the value of S-the number of failed items at the time of the r th failure. Therefore, the problem is equivalent to inference based on S, T 1: n , T 2: n , . . . , T r : n , and with r = n − k + 1 it is exactly the same problem as inference from failure times of components of a k-out-ofn system up to and including its breakdown. To the best of our knowledge maximum likelihood inference for discrete distributions based on censored data has not been studied before in the literature.

Strong consistency
The standard theorems of asymptotic theory of MLE's constructed from iid observations do not apply to our problem in which we make inference from dependent and non-identically distributed rv's S, T 1: n , . . . , T n−k+1: n . However, as will be shown later on, the basic machinery of proving that under some regularity conditions MLE's from iid observations are strongly consistent can be modified to derive the following analogous result for MLE's obtained from failure data of a k-out-of-n system.
Theorem 1 Assume that the family F = {F(λ), λ ∈ Θ} satisfies the following three conditions: stands for the largest integer not exceeding x, then there exists a sequence (θ n , n ≥ 1) such that, with P θ -probability 1, -for all sufficiently large n,θ n is a solution to the likelihood equation (3); -θ n → θ as n → ∞.
Proof See the "Appendix".
Theorem 1 can be used in practice, because for a given family F = {F(θ, ·), θ ∈ Θ} we can check if its assumptions are satisfied without knowing the value of the true parameter θ . In particular, this theorem will allow us to deduce that the MLE's obtained in the next section are strongly consistent.

MLE's for some specific families of distributions
In this section, we will consider the Poisson Poiss(θ ), θ > 0, binomial b(w, θ ), θ ∈ (0, 1), and negative binomial nb(w, θ ), θ ∈ (0, 1), distributions as possible component lifetime distributions of a k-out-of-n system. These three distributions, besides the geometric one, are listed by Barlow and Proschan (1996) as typical discrete failure distributions widely used in reliability engineering. We will prove that in the case of all these discrete distributions if the MLE of the parameter θ based on observed values of S, T 1: n , . . . , T S: n exists then it is unique. Hence, Theorem 1 will guarantee that in these cases the MLE of θ exists almost surely for sufficiently large n and is strongly consistent.
It is worth pointing out that, since the geometric distribution is a special case of the negative binomial distribution, results proved here for the negative binomial lifetimes of components hold in particular for geometrically distributed lifetimes. Yet the geometric case is easier-then a closed-form formula for the MLE of θ can be obtained and hence not only asymptotic but also exact distributional properties of this estimator can be given. For this purpose, the geometric case will be considered in details in a separate paper.
To prove results of this section, we will make use of the following two lemmas. The first one is taken from Pólya and Szegő (1998, p. 41).

Lemma 1 Let the radius of convergence of the power series
, let the number of its zeros in the interval 0 < x < ρ be Z and let the number of changes of sign in the sequence of its coefficients be C. Then Z ≤ C.
The second lemma concerns linear combinations of Bernstein polynomials and was first proved by Schoenberg (1959). Recall that Bernstein polynomials of degree w are defined by

Lemma 2 The number of zeros of a given nonzero linear combination of Bernstein polynomials B
, does not exceed the number of sign changes of the sequence β = (β 0 , . . . , β n ). The first and the last signs of the sum are identical to the signs of the first and the last nonzero element of β, respectively.

Poisson distribution
Let the component lifetimes where θ > 0 is the parameter to estimate. Then by (1) and (4) the observed likelihood function of S, T 1: n , . . . , T S: n can be written as From (6) it is easily seen that if s = n, then the function L(θ ), θ > 0, has a global maximum at δ/n if δ > 0, and does not attain a global maximum if δ = 0. Hence, the MLE does not exist when T 1 = T 2 = · · · = T n = 0. We see at once that the probability of such an event, e −nθ , decreases to 0 as n → ∞. It is worth also noting that if s = n and δ > 0, that is if we observe the following event {T n−k+1: n = T n−k+2: n = · · · = T n: n > 0}, then the MLE is just equal to the sample mean.
It remains to consider the case when s ∈ {n − k + 1, . . . , n − 1}. For this purpose, note that (6) can be rewritten as From (7) it is clear that d dθ log L(θ ) has the same sign as h(θ ). But h(θ ) can be represented as where and t z 1 +···+z i − j < 0 for i = 1, . . . , m and j ≥ t s + 2, which is due to the fact that 0 ≤ t z 1 < t z 1 +z 2 < · · · < t z 1 +···+z m = t s . Consequently, the number of sign changes in the sequence (α s,t 1 ,...,t s ( j), j ≥ 0) equals one. The radius of convergence of the power series in (8) is ρ = ∞. Therefore, Lemma 1 guarantees that the number of zeros of h(θ ) in the interval (0, ∞) is at most one. But from (7) we see that (10) Hence, h(θ ) has exactly one zero in (0, ∞), which by (7) shows that the likelihood equation has exactly one solution in (0, ∞). Moreover, the function L(θ ) is first increasing, then decreasing, which implies it attains its global maximum and the observed MLE of θ , being the solution to (11), is unique. From (9) and (10) we know that the observed MLE of θ belongs to the finite interval δ s , δ+(n−s)(t s +1) s and therefore can be found easily through numerical methods.
Thus, we have proved the following theorem Theorem 2 From the Poisson distribution with pmf given in (5), suppose we have observed failure times of components of a k-out-of-n system up to and including the breakdown of the system S = s, T 1: n = t 1 , . . . , T s: n = t s .

Binomial distribution
Now suppose that the component lifetimes T i , i = 1, . . . , n, of a k-out-of-n system have the binomial distribution b(w, θ ) with the following pmf where w ∈ {1, 2, . . .} is known and θ ∈ (0, 1) is the parameter to estimate. With the notation (4), the observed likelihood function (1) is given by where C 2 does not depend on θ . An easy computation shows that, for θ ∈ (0, 1), If s = n then we see from (14) that the function L(θ ), θ ∈ (0, 1), has a global maximum at δ/(wn) if δ > 0, and does not attain a global maximum if δ = 0. Hence, similarly to the Poisson case, the MLE does not exist when T 1 = T 2 = · · · = T n = 0 and it is easily seen that the probability of non-existence, (1 − θ) nw , approaches 0 as n → ∞. Moreover, if T n−k+1: n = T n−k+2: n = · · · = T n: n > 0, then the MLE is equal to the sample mean divided by w.
The case of s ∈ {n − k + 1, . . . , n − 1} requires more effort. Note that (14) can be rewritten as where Clearly g(θ ) has the same sign as d dθ log L(θ ). Moreover, it is easy to check that g(θ ) can be represented as the following linear combination of Bernstein polynomials The coefficient β s,t 1 ,...,t s (t s ) is positive since t s < w when s < n. Now for j = t s + 1, . . . , w − 1 we check that β s,t 1 ,...,t s ( j) < 0, which is equivalent to the inequality But t z 1 < t z 1 +z 2 < · · · < t s < w, which shows that for j = t s + 1, . . . , w − 1 the expressions in the braces in (16) are negative and hence (16) holds. Finally, we verify that β s,t 1 ,...,t s (w) < 0. This corresponds to the inequality which is true because of the relation (17). Summarizing we have proved that Lemma 2 now ensures that the sign of the derivative (15) is first positive and then negative on (0, 1). From this we conclude that the observed likelihood function (13) is first increasing and then decreasing there. Hence, it has a global maximum in (0, 1) which is attained at the point being the only solution to the observed likelihood equation. Therefore, the observed MLE of θ exists and is unique.
Thus, we have the following analogue of Theorem 2. (12), suppose we have observed failure times of components of a k-out-of-n system up to and including the breakdown of the system S = s, T 1: n = t 1 , . . . , T s: n = t s .

Negative binomial distribution
Consider a k-out-of-n system composed of n components whose lifetimes T i , i = 1, . . . , n, have the negative binomial distribution nb(w, θ ) with a pmf where w ∈ {1, 2, . . .} is known and θ ∈ (0, 1) is the parameter to estimate. Then the observed likelihood function (1) takes on the form where C 3 does not depend on θ and δ is given in (4). The observed likelihood equation (2) becomes (21) If s = n and δ = 0, then the function L(θ ) given in (19) is decreasing and consequently the observed MLE of θ does not exist. It is obvious that the probability of non-existence, P(T 1 = T 2 = · · · = T n = 0) = (1 − θ) nw , decreases to 0 as n → ∞. Otherwise, that is when s < n or δ > 0, we have and, since L(θ ), θ ∈ (0, 1), is continuous and positive, it has a global maximum. In the case when s = n and δ > 0 we easily see from (20) that this global maximum is attained at θ = δ/ (nw + δ) and hence that the observed MLE of θ is equal tô wheret is the observed sample mean. It remains to consider the case when s ∈ {n − k + 1, . . . , n − 1}. For this purpose, note that the left-hand side of (21) can be represented as the following power series where Indeed, if s < n then δ + (n − s)(t s + 1) > 0, which implies γ s,t 1 ,...,t s (t s + 1) > 0. Moreover, if j ≥ t s + 2 then j > t s > · · · > t z 1 +z 2 > t z 1 and consequently which shows that γ s,t 1 ,...,t s ( j) < 0 for j ≥ t s + 2.
Since the radius of convergence of the power series in (22) is ρ = 1, from Lemma 1 and (23) we obtain that the left-hand side of (21) [or equivalently of (20)] considered as a function of θ has at most one zero in the interval (0, 1). But from the previous discussion, we know that the function L(θ ), θ ∈ (0, 1), has a global maximum. Therefore, the left-hand side of (20) has exactly one zero in (0, 1) and this zero is a point at which the likelihood function L(θ ) attains its global maximum. The observed MLE of θ is unique and can be obtained easily by numerical methods.
Thus, we have proved the following result.
Theorem 4 From the negative binomial distribution with pmf given in (18), suppose we have observed failure times of components of a k-out-of-n system up to and including the breakdown of the system S = s, T 1: n = t 1 , . . . , T s: n = t s .

Monte Carlo simulation study
From Sect. 4 we know that in the case of Poisson Poiss(θ ), binomial b(w, θ ) and negative binomial nb(w, θ ) distributions the maximum likelihood estimators of θ based on failure times of components of a k-out-of-n system observed up to and including the breakdown of the system are strongly consistent as n → ∞ and k = [ pn], where p ∈ (0, 1) is fixed. The aim of this section is to investigate finite-sample properties of these estimators via Monte Carlo simulation study. For this purpose we assume Poisson Poiss(θ = 1), binomial b(w = 4, θ = 0.5) and negative binomial nb(w = 5, θ = 0.15) component lifetimes. The parameters of these distributions were chosen so that the corresponding variances are equal (in the case of the Poisson and binomial distributions) or approximately equal (in the case of the negative binomial distribution) to one. The almost equal variances allow to make comparisons between the three considered cases. Next, for each of the chosen distributions and for some selected values of n and k we generate N = 1000 times the failure times of components of a k-out-of-n system observed up to and including the system breakdown obtaining the data of the form s (i) , t (i) 1 ≤ · · · ≤ t (i) n−k+1 = · · · = t (i) s (i) , i = 1, . . . , N . For each i = 1, . . . , N we then computeθ (i) ML , the maximum likelihood estimator of θ , using numerical methods if necessary. More precisely, to solve the corresponding likelihood equation we use the method of finding the unique root of a continuous function in a finite interval, such as the bisection method. Finally, we compute the mean and standard deviation ofθ (i) ML , i = 1, . . . , N . These values can be treated as the simulated expectation and standard deviation ofθ ML . The obtained results are presented in Tables 1, 2 and 3. It is interesting that during the simulations we did not encounter samples with non-existing MLE's. This was so because for the cases considered in the tables the probabilities of non-existence are very small as Table 4 shows.
In the simulation study we observe that even for small n (n = 15) the bias ofθ ML is small-the simulated expectations ofθ ML are close to the true values of θ . As n and k increases in such a way that k/n is kept fixed both the bias and standard deviation of θ ML decreases. Moreover, from Tables 1 and 3 we see that for the same values of n the  Table 4 Selected approximate values of probabilities of non-existence of MLE (formulas for these probabilities are taken from Sect. 4) 3 × 10 −7 9 × 10 −19 5 × 10 −6 n = 30 9 × 10 −14 8 × 10 −37 3 × 10 −11 n = 60 9 × 10 −27 6 × 10 −73 7 × 10 −22 bias and standard deviation are smaller when n − k + 1 is larger. This is so because the single experiment terminates at the moment of the (n − k + 1)th component failure and larger n − k + 1 allows to collect more information and thus to obtain a better precision of estimation.
In Table 2, we see a surprising situation. For the same values of n − k + 1 we obtain larger bias when n is larger. This may not agree with our intuition-larger n means that more elements are involved in a single experiment and thus we may expect a better estimation. Yet this is not the case. The reason is that if k = k(n) = [(1 − q)n], where q ∈ (0, 1) is such that the qth quantile of F(θ, ·) is not unique, then due to (38) the behavior of T n−k+1: n is unstable causing worse behavior ofθ ML . Note that the qth quantile of the binomial b(w = 4, θ = 0.5) distribution is not unique when 1 − q = 11/16 and is so when 1 − q = 2/3. Therefore, biases presented in the lefthand side of Table 2 are greater than the corresponding ones given in the right-hand side of this table.
1. For an uncensored sample it is well known that the MLE of θ is equal to the mean. Therefore, based on the whole sample we obtainθ (1) ML = 15.86. 2. Now suppose that we terminate the experiment at the moment of the r = 5th failure, that is after 20 days. Then we have exactly s = 5 air monitors broken. Using the inference based on S(ω) = 5, T 1: 7 (ω) = 8, T 2: 7 (ω) = 8, T 3: 7 (ω) = 10, T 4: 7 (ω) = 10, T 5: 7 (ω) = 20, we getθ (2) ML = 12.80. 3. Finally, let us consider censoring by terminating the experiment at the moment of the r = 3th failure, that is after 10 days. Then we observe s = 4 air monitors breakdowns. Thus, we collect the following data: S(ω) = 4, T 1: 7 (ω) = 8, T 2: 7 (ω) = 8, T 3: 7 (ω) = 10, T 4: 7 (ω) = 10. The MLE based on this data is equal toθ (3) ML = 10.87. We see that the value of the MLE changes significantly when we change the censoring scenario. This unpleasurable feature is due to the fact that n = 7 is very small. Apparently, to obtain more reliable estimates we need to conduct an experiment with a larger number of air monitors.

Conclusions
In this paper, we have focused on maximum likelihood inference of the discrete lifetime distribution of components of a k-out-of-n system in the case when failure times of the components observed up to and including the moment of the breakdown of the system are available. Another problem of interest is the inference in the case when a sample of lifetimes of k-out-of-n systems and numbers of broken components at the moment of the system failure is given. We are currently working on the latter problem and planning to report our findings in a forthcoming paper.
It is also worth pointing out that the new results we obtained for the discrete case are analogous to that known in the literature for the continuous case in the sense that in both the cases under some regularity conditions the MLE's of interest exist almost surely for sufficiently large n and are strongly consistent. Yet, the regularity conditions for the two cases are different and in the proofs different techniques are needed.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. γ q = inf{x ∈ R: F(θ, x) ≥ q} and γ q = sup{x ∈ R: F(θ, x) ≤ q}.
Consequently, for c ∈ N (θ ), where, in the case of M < ∞, we can define h θ (M) = 0. We will show that A n (β) can be taken to be equal to the right-hand side of (31). For this purpose, note that the right-hand side of (31) is a mean of iid rv's with P θ -expectation given by β j=0 g θ ( j) f (θ, j) + h θ (β)F(θ, β).
Set w θ (β) equal to (32). It is easily seen that w θ (β) ∈ [0, ∞). By the strong law of large numbers, we obtain (29), and the proof of part (3) is complete. Now we are ready to prove Theorem 1. We will divide the proof into two parts. In the first part we will consider the simpler case when the qth quantile of F(θ, ·) is unique. In the second part we will deal with the more difficult case when this quantile is not unique.