Parameter Estimation for Type III Discrete Weibull Distribution: A Comparative Study

Abstract

The type III discrete Weibull distribution can be used in reliability analysis for modeling failure data such as the number of shocks, cycles, or runs a component or a structure can overcome before failing. This paper describes three methods for estimating its parameters: two customary techniques and a technique particularly suitable for discrete distributions, which, in contrast to the two other techniques, provides analytical estimates, whose derivation is detailed here. The techniques’ peculiarities and practical limits are outlined. A Monte Carlo simulation study has been performed to assess the statistical performance of these methods for different parameter combinations and sample sizes and then give some indication for their mindful use. Two applications of real data are provided with the aim of showing how the type III discrete Weibull distribution can fit real data, even better than other popular discrete models, and how the inferential procedures work. A software implementation of the model is also provided.

1. Introduction

Most reliability studies assume that time is continuous, and continuous probability distributions such as exponential, gamma, Weibull, normal, and lognormal are commonly used to model the lifetime of a component or a structure. These distributions and the methods for estimating their parameters are well known. In many practical situations, however, lifetime is not measured with calendar time: for example, when a machine works in cycles or on demands and the number of cycles or demands before failure is observed; or when the regular operation of a system is monitored once per period, and the number of time periods successfully completed is observed. Moreover, reliability data are often grouped into classes or truncated according to some censoring criterion. In all these cases, lifetime is modeled as a discrete random variable (r.v.). Indeed, not much work has been done in discrete reliability. Generally, most reliability concepts for continuous lifetimes have been adapted to the discrete case; in particular, discrete analogues of continuous distributions have been introduced [1]. In this context, geometric and negative binomial distributions are the corresponding discrete alternatives for the exponential and gamma distributions, respectively. Yet, discrete lifetime distributions can be defined without any reference to a continuous counterpart. Bracquemond and Gaudoin [2] provided an exhaustive survey on discrete lifetime distributions.

The (two-parameter) continuous Weibull distribution is one of the most widely used stochastic distributions for modeling the life of a component; it is very flexible since, for different chosen shape parameters, it models either increasing or decreasing failure rates [3]. As a discrete alternative to the continuous Weibull distribution, three main forms have been proposed. The first one was introduced in Nakagawa and Osaki [4] and is referred to as type I discrete Weibull; it mimics the cumulative distribution function of the continuous Weibull distribution. The second one (type II) was proposed and studied in Stein and Dattero [5]; it mimics the hazard rate of its continuous counterpart. The third, which is referred to as type III and this paper considers, was introduced in Padgett and Spurrier [6]: their approach does not start from the continuous Weibull distribution but tries to generalize the notions of hazard rate and mean residual life to the discrete case [7]. From a different perspective and with a different objective, Roy and Dasgupta [8] proposed a discretization method for continuous random variables for computing reliability in complex stress-strength models, with a specific application to the Weibull r.v.; however, the support of the discrete r.v. thus provided is not the set of nonnegative integers.

Type III discrete Weibull (henceforth simply discrete Weibull) r.v. is not similar in functional form to any of the functions describing a continuous Weibull distribution. It is defined by the following cumulative distribution: 𝐹𝑥;𝑐,𝛽=𝑃𝑋≤𝑥;𝑐,𝛽=1−𝑒^{−𝑐∑𝑥+1𝑗=1𝑗𝛽}𝑥=0,1,… (1) with 𝑐>0 and 𝛽≥−1 (and not 𝛽∈ℝ, as stated in Padgett and Spurrier [6]) or equivalently by the following probability mass function: 𝑃𝑋=0;𝑐,𝛽=𝑝0;𝑐,𝛽=1−𝑒^−𝑐𝑃𝑋=𝑥;𝑐,𝛽=𝑝𝑥;𝑐,𝛽=𝑒^{−𝑐∑𝑥𝑗=1𝑗𝛽}1−𝑒^{−𝑐(𝑥+1)𝛽}𝑥=1,2,… (2) which can be more compactly rewritten as follows: 𝑃𝑋=𝑥;𝑐,𝛽=𝑝𝑥;𝑐,𝛽=𝑒^{−𝑐∑𝑥𝑗=1𝑗𝛽}1−𝑒^{−𝑐(𝑥+1)𝛽}𝑥=0,1,2,…, (3) letting ∑^𝑥_𝑗=1𝑗^𝛽=0 if 𝑥=0. The corresponding survival function is 𝑆(𝑥)=𝑃𝑋≥𝑥;𝑐,𝛽=𝑒^{−𝑐∑𝑥𝑗=1𝑗𝛽}𝑥=0,1,…, (4) and the hazard rate (or failure rate) function is ℎ(𝑥)=𝑝(𝑥)𝑆(𝑥)=1−𝑒^{−𝑐(𝑥+1)𝛽}𝑥=0,1,…. (5)

The complexity of the expression of the probability mass function (3) has somehow hindered the use and diffusion of this discrete model. A plot of the probability mass function of the discrete Weibull for three combinations of its parameters is given in Figure 1. Note that for 𝛽=0, the discrete Weibull reduces to a geometric r.v. with parameter 1−exp(−𝑐), characterized then by a constant failure rate; for 𝛽>0, the distribution has an increasing failure rate, whereas for 𝛽<0, it has a decreasing failure rate. The first and second moments cannot in general be expressed in a closed form but can only be expressed as an infinite series 𝔼(𝑋)=^+∞_𝑥=0𝑥𝑃𝑋=𝑥;𝑐,𝛽=^+∞_𝑥=1𝑒^{−𝑐∑𝑥𝑗=1𝑗𝛽},𝔼𝑋²=^+∞_𝑥=0𝑥²𝑃=𝑋=𝑥;𝑐,𝛽_+∞𝑥=12𝑥𝑒^{−𝑐∑𝑥𝑗=1𝑗𝛽}−𝔼(𝑋). (6)(7)

Figure 1 

Probability mass function of type III discrete Weibull distribution for three combinations of parameters.

The first moment (6) is finite for 𝛽>−1 or for 𝛽=−1 and 𝑐>1; the second moment (7) is finite for 𝛽>−1 or for 𝛽=−1 and 𝑐>2 (see the appendix).

This study first describes and discusses procedures for estimating parameters 𝑐 and 𝛽 of the discrete Weibull r.v., one of which is an original adaptation of a technique used for some discrete distributions, emphasizing their practical limits (Section 2). An extensive Monte Carlo study, implemented through an adhoc package developed under the R environment, assesses and compares the performance of these estimators, in terms of bias and variability, for different combinations of the parameters and sample sizes (Section 3). The estimation procedures are applied to two datasets taken from the literature (Section 4), and, finally, the summarizing remarks conclude the paper (Section 5).

2. Parameters Point Estimation

Focusing on the point estimation of the parameters of the discrete Weibull r.v., based on an observed simple random sample 𝐱=(𝑥₁,𝑥₂,…,𝑥_𝑛) of size 𝑛, three techniques are now described: the method of proportion, which is strictly related to the specific features of the distribution function of the discrete Weibull r.v.; the maximum likelihood method; and the method of moments. Below, we assume that both parameters are unknown. Special attention will be devoted to detecting samples leading to implausible estimates or to no estimate at all with either technique.

Method of Proportion. This method was originally introduced in Khan et al. [9] for type I discrete Weibull and was extended by Jazi et al. [10] to discrete inverse Weibull. It relies upon the estimation of a probability (or two or more probabilities, according to the number of parameters involved) using the corresponding sample proportion(s). For the present model, we have 𝑃(𝑋=0)=1−exp(−𝑐), which can be estimated through the proportion of 0’s in the sample, 𝑦/𝑛, where 𝑦=∑^𝑛_𝑖=1𝐼_𝑥𝑖=0 denotes the number of 0’s in the sample, 𝐼_𝐴 the indicator function, which equals 1 if 𝐴 is true and 0 if 𝐴 is false. Thus, an estimate of 𝑐 is 𝑐_𝑃̂𝑦=−log1−𝑛. (8) Similarly, the probability 𝑃(𝑋=1)=exp(−𝑐)[1−exp(−𝑐2^𝛽)] can be estimated by the proportion of 1’s in the sample, 𝑧/𝑛, with 𝑧=∑^𝑛_𝑖=1𝐼_𝑥𝑖=1. Consequently, also substituting for 𝑐 its estimate (8) (see the appendix), 𝛽 is estimated by 𝛽_𝑃=loglog1−𝑦/𝑛−𝑧/𝑛/log1−𝑦/𝑛−1.log2 (9) Though the method of proportion easily provides an analytical expression for both estimates, an apparent weakness is that this method does not exploit all the information contained in the sample, since the estimates involve only the rates of 0’s and 1’s. Moreover, it fails to provide feasible estimates for 𝛽 if there are no 0’s in the sample. In this case, 𝑐_𝑃̂=0, which is a boundary value for 𝑐, and 𝛽_𝑃 in (9) cannot be computed. This is particularly common for small samples, and for small values of 𝑐. If in the sample there are no 1’s (𝑧=0), then by formula (9), the estimate 𝛽_𝑃 is not available again. The method fails in computing an estimate for 𝛽 even if the sample contains only 0’s and 1’s (𝑦+𝑧=𝑛).

The method of proportion can also lead to “implausible” estimates of 𝛽, that is, estimates that do not belong to its parameter space: 𝛽_𝑃<−1. This happens when log1−𝑦/𝑛−𝑧/𝑛1log1−𝑦/𝑛−1<2, (10) which is equivalent to say 𝑦1−𝑛−𝑧𝑛>𝑦1−𝑛^3/2, (11) or in terms of 𝑧𝑧𝑛<𝑦1−𝑛−𝑦1−𝑛^3/2. (12) The probability of an “implausible” estimate could be theoretically derived from the last inequality, remembering that 𝑦 and 𝑧 are the first two marginal distributions of a trinomial r.v. with parameters (𝑛, 𝑝₀=𝑃(𝑋=0), 𝑝₁=𝑃(𝑋=1), 𝑝₂=1−𝑝₀−𝑝₁). Figure 2 shows the combinations of 𝑦/𝑛 and 𝑧/𝑛 leading to such implausible estimates as a subset of the triangular region (0≤𝑦/𝑛≤1; 0≤𝑧/𝑛≤1−𝑦/𝑛). For example, in the sample 𝐱=(0,0,0,1,2,2,3,5,6,8), with 𝑦/𝑛=0.3 and 𝑧/𝑛=0.1, the method of proportion provides, through (8) and (9), the estimates 𝑐_𝑃̂=0.357 and 𝛽_𝑃=−1.210.

Figure 2 

“Plausible” and “implausible” regions for the estimate of 𝛽 through the method of proportion.

Maximum Likelihood Method. Having defined the log-likelihood function as 𝑙(𝑐,𝛽;𝐱)=log∏^𝑛_𝑖=1𝑃(𝑋=𝑥_𝑖;𝑐,𝛽), we obtain 𝑙𝑐,𝛽;𝐱=−𝑐^𝑛𝑥_𝑖=1^𝑖_𝑗=1𝑗^𝛽+^𝑛_𝑖=1log1−𝑒^{−𝑐(𝑥𝑖+1)𝛽}. (13) The maximum likelihood estimates of 𝑐 and 𝛽(𝑐̂_ML,𝛽_ML) are defined as the values that maximize the log-likelihood function 𝑐̂_ML,𝛽_ML=argmax𝑙𝑐,𝛽;𝐱. (14) The two first derivatives of 𝑙(𝑐,𝛽;𝐱) are quite easily computed, 𝜕𝑙𝑐,𝛽;𝐱=𝜕𝑐^𝑛_𝑖=1𝑥_𝑖+1^𝛽𝑒^{−𝑐(𝑥𝑖+1)𝛽}1−𝑒^{−𝑐(𝑥𝑖+1)𝛽}−^𝑥𝑖_𝑗=1𝑗^𝛽=^𝑛_𝑖=1𝑒^{−𝑐(𝑥𝑖+1)𝛽}∑^𝑥𝑖+1_𝑗=1𝑗^𝛽−∑^𝑥𝑖_𝑗=1𝑗^𝛽1−𝑒^{−𝑐(𝑥𝑖+1)𝛽},𝜕𝑙𝑐,𝛽;𝐱𝜕𝛽=𝑐^𝑛_𝑖=1𝑒^{−𝑐(𝑥𝑖+1)𝛽}∑^𝑥𝑖+1_𝑗=1𝑗^𝛽log𝑗−∑^𝑥𝑖_𝑗=1𝑗^𝛽log𝑗1−𝑒^{−𝑐(𝑥𝑖+1)𝛽}, (15) but, as already noted in Padgett and Spurrier [6], the solution to the maximization of 𝑙(𝑐,𝛽;𝐱) (with the constraints that 𝑐 and 𝛽 belong to in their natural parameter spaces) cannot be derived in a closed form, by equating the expressions in (15) to zero, but can be obtained only numerically, for example, using one of the functions nlm, optim, and Rsolnp in the R programming environment, which allows the user to solve nonlinear constrained or unconstrained minimization/maximization problems. Even this method cannot be successfully applied to every possible sample; in particular, the method fails in providing a solution if 𝑦+𝑧=𝑛 (i.e., if the sample contains only 0’s and 1’s). In this case, in fact, it can be easily proved that the first-order 𝛽 derivative of the log-likelihood is never null; the log-likelihood does not have an absolute maximum in the parameter space. Let us consider, for example, the following sample: 𝐱=(0,1). The corresponding likelihood function is given by 𝐿𝐱;𝑐,𝛽=1−𝑒^−𝑐𝑒^−𝑐1−𝑒^−𝑐2𝛽 (16) and the log-likelihood 𝑙𝑐;𝛽=−𝑐+log1−𝑒^−𝑐+log1−𝑒^−𝑐2𝛽. (17) The first-order derivatives are𝜕𝑙𝑐,𝛽;𝑥_𝑖𝑒𝜕𝑐=−1+^−𝑐1−𝑒^−𝑐+2^𝛽𝑒^−𝑐2𝛽1−𝑒^−𝑐2𝛽,𝜕𝑙𝑐,𝛽;𝑥_𝑖=𝑒𝜕𝛽^−𝑐2𝛽𝑐2^𝛽log21−𝑒^−𝑐2𝛽. (18a)(18b)The derivative (18b) is never null but tends to zero for 𝛽→+∞; when 𝛽→+∞ the third addend in the derivative (18a) goes to zero, and solving the equation −1+exp(−𝑐)/[1−exp(−𝑐)]=0 one can obtain 𝑐̂^∗_ML=log2. Note that this is the same estimate the method of proportion supplies.

Method of Moments. Through this method, the parameter estimates are obtained solving, in terms of 𝑐 and 𝛽, the equations 𝔼(𝑋)=𝜇₁=𝑚₁ and 𝔼(𝑋²)=𝜇₂=𝑚₂, where 𝑚₁ and 𝑚₂ are the first and second sample moments: 𝑚₁=∑^𝑛_𝑖=1𝑥_𝑖/𝑛, 𝑚₂=∑^𝑛_𝑖=1𝑥²_𝑖/𝑛. Since they cannot be solved analytically, as suggested in Khan et al. [9], one can numerically minimize, with respect to 𝑐 and 𝛽, the following quadratic loss function: ℒ𝑐,𝛽=𝑚₁−𝜇₁²+𝑚₂−𝜇₂². (19) The task can be carried out, for example, again using the functions nlm, optim, or solnp in the R environment [11]. The solution is the couple (𝑐_𝑀̂,𝛽_𝑀), which should correspond to the value ℒ(𝑐_𝑀̂,𝛽_𝑀)=0. Note that the minimization should be subject to the natural constraints 𝑐>0 and 𝛽>−1. Without these constraints, the minimization algorithm may get stuck in implausible intermediate-step solutions. Indeed, the optim function under the R environment seems to provide excellent results in terms of convergence to the optimal solution, even without setting the constraints on 𝑐 and 𝛽. As to the launch values for 𝑐 and 𝛽, required by any minimization function, one can set 𝛽=𝛽₀=0 and 𝑐=𝑐₀=log[(1+𝑥)/𝑥]. This pair of values is always computable and feasible (unless the sample contains all 0’s) and ensures at the first iteration that 𝑚₁=𝔼(𝑋): in fact, recalling (6), 𝔼(𝑋;𝑐₀,𝛽₀)=∑^+∞_𝑖=1exp(−𝑐₀∑^𝑖_𝑗=11)=∑^+∞_𝑖=1exp(−𝑖𝑐₀)=exp(−𝑐₀)/[1−exp(−𝑐₀)]=𝑥.

When the sample contains only 0’s and 1’s, the method of moments is not applicable. To see this, first consider the probability mass function of the discrete Weibull r.v., and note that by letting 𝛽 tend to +∞ in (3), it degenerates into a r.v. that takes only two values: 0 with probability 1−exp(−𝑐) and 1 with probability exp(−𝑐). It is then clear that if the sample is made up of a fraction 𝑟 of 0’s and a fraction (1−𝑟) of 1’s, the equality of both first and second moments computed on the sample and on the original r.v. holds for 1−exp(−𝑐)=𝑟; that is, 𝑐=−log(1−𝑟) and 𝛽→+∞. In this case, the loss function ℒ(𝑐,𝛽) does not admit an absolute minimum, it only admits an inferior limit.

For some samples, the numerical minimization procedure can require a huge computation time, much larger than that required by the maximum likelihood method. This is due to the iterated calculation of the first and second moments, which is itself numerical and particularly time consuming for the negative values of 𝛽 (in this case, in fact, the convergence of the series in (6) and (7) is slower).

Given the complexity of the estimators derived through the methods listed in this section (and only the method of proportion provides an analytical expression for them), not as much can be analytically derived about their statistical properties for finite sample size, that is, bias and variability. When the method of proportion can be applied, it provides a consistent estimator for 𝑐; 𝛽_𝑃 is consistent as well, but nothing can be said about their unbiasedness (see the analogous discussion in Khan et al. [9] for type I discrete Weibull). For large samples, the general properties of the estimators derived from the maximum likelihood method and the method of moments can be recalled. In the next section, an extensive simulation study is presented, which was performed to investigate the performance of the estimation methods and outline some practical advice for their use.

3. Simulation Study

The estimators presented in the previous section were investigated through an extensive Monte Carlo study; they were compared in terms of bias (𝐵), defined as 𝐵(𝜃_𝐸)=𝔼(𝜃_𝐸)−𝜃, and root mean square error (RMSE), defined as RMSE(𝜃_𝐸)=𝔼[(𝜃_𝐸−𝜃)²], where 𝜃 denotes one of the two parameters (𝑐 or 𝛽) and 𝜃_𝐸 denotes one of the three corresponding estimators, according to the method indicated by the subscript 𝐸 (method of proportion, 𝑃; maximum likelihood method, ML and method of moments, 𝑀).

3.1. Simulation Design

In this study, several parameter combinations and sample sizes (𝑛=20,50,100) were considered. The values (𝑐,𝛽) were chosen in order to explore a large spectrum of the discrete Weibull distribution, in particular to comprise increasing, constant, and decreasing failure rates. At the same time, the parameters were set in order to keep the discrete nature of the distribution reasonable: values entailing a nonnegligible probability for a large number of integer values were deliberately excluded (in this case, a continuous r.v. should be preferred to model failure data); parameter values ensuring a nonnegligible probability for just the first integers were avoided as well (one assumes the component that is monitored is likely to last for more than 1 or 2 cycles). Table 1 shows the combinations of values of 𝑐 and 𝛽 examined in the simulation study along with the corresponding expected value, standard deviation, and 99% quantile. A note about the computation of the expected value and standard deviation is due: they are calculated numerically—see formulas (6) and (7)—considering the first 𝑛_max integers, with 𝑛_max=2𝐹⁻¹(1−𝜖) and 𝜖 as small as possible (here, 𝜖=0.0001, which actually proved to be a satisfactory practical choice for the considered scenarios).

Table 1 

Parameter combinations and corresponding expected value, standard deviation, and 0.99-quantile for the simulation study.

3.2. Software Implementation

The simulation study was based on 5,000 Monte Carlo replications and carried out under the R environment [11]. In particular, the necessary code was developed to implement the probability mass function, the cumulative distribution function and its inverse, and the random generation for the discrete Weibull model; to compute the first and second moments; and to to realize the algorithms corresponding to each estimation method. This code, structured as an R package, DiscreteWeibull, is freely available in the CRAN repository [12].

3.3. Simulation Results

Table 2 shows the bias and root mean square error for both estimators derived from each of the three methods, under each combination of the two parameters and for each sample size. The results took into account the applicability limits of the three methods, computing the Monte Carlo averages over the feasible samples only. For the smallest sample size (𝑛=20) and, more rarely, for 𝑛=50, the method of proportion met a certain percentage of infeasible samples under several scenarios, the highest (about 24%) with 𝑛=20, 𝑐=0.1, and 𝛽=0.25. The maximum likelihood method and the method of moments encountered only a few infeasible samples under some scenarios, whose rate was in any case always smaller than 1%.

Table 2 

Simulation results: bias (B) and root mean square error (RMSE) for the estimators of 𝑐 and 𝛽.

Note that strictly speaking, to compare the behavior of the estimators under different scenarios, one should use their relative bias and relative root mean square error as performance indexes, that is, the bias and the root mean square error divided by the MC expected value of the estimator. This way of proceeding would provide a more correct “double” reading of the results, assessing the performance of a single estimator moving through different scenarios and comparing the performance of the three estimators under a fixed scenario. Our choice of using “absolute” instead of “relative” indexes was partially dictated by the necessity of including negative and null values for 𝛽, which would represent difficulty when computing the relative indexes (the denominator could take very small values, dramatically amplifying the corresponding absolute value of the index and making the reading arduous).

Trying to synthesize the results, let us start with the largest sample size (𝑛=100) and the estimators of parameter 𝑐. It is quite evident that all three estimators show a modest bias, regardless of the values of 𝑐 and 𝛽. The maximum likelihood method looks the best in this sense, while the method of moments shows the largest bias in absolute value for high values of 𝑐 and negative values of 𝛽, and it is often negatively biased. As to the root mean square error, the method of moments and the maximum likelihood method present very similar values under each scenario; the method of proportion shows a larger value, especially for small values of 𝑐. The magnitude of RMSE seems to be much more affected by the value of 𝑐 rather than by 𝛽 for all the estimators.

As to the estimators of 𝛽, on average, the method of proportion provides the least biased estimator, unless 𝑐 is too small, while the maximum likelihood method and, to a larger extent, the method of moments are significantly biased. While 𝑀 and ML estimators always overestimate the true value of the parameter, the 𝑃 estimator seems to underestimate it for lower values. In absolute value, the bias of the ML estimator increases as 𝑐 increases for a fixed 𝛽; it increases as 𝛽 increases for a fixed 𝑐. The bias of the 𝑀 estimator increases as 𝑐 increases for a fixed 𝛽; it increases as 𝛽 decreases for a fixed 𝑐. In terms of variability, the 𝑀 and ML estimators perform in a similar way for each parameter combination, while the 𝑃 estimator is much more variable and presents very different values moving through the scenarios: its RMSE is from 2 to 6 times larger than the corresponding value of its competitors. The worst scenarios for the estimator, in this sense, are those connected with 𝑐=0.1, and this is quite reasonable. In fact, in this case, the probability 𝑃(𝑋=0) equals 0.095, which corresponds to the expected fraction of zeros in the sample, and this small value leads to a modest performance of the 𝑃 method, which uses only proportions of 0’s and 1’s of the samples and discards the other information.

Decreasing the sample size to 50 and 20, the bias in absolute value and the root mean square error of the estimators for each scenario obviously increase. However, the behaviors and trends exposed for 𝑛=100 hold still. The bias of the estimators of 𝛽 derived by the method of moments and the maximum likelihood method becomes substantial especially for large values of 𝑐 and negative values of 𝛽; in these cases, the method of proportion is less biased, but its RMSE is still much larger than those of its competitors. For high values of 𝑐 (viz., equal to or larger than 0.8) and 𝑛=20, the bias in absolute value of the estimators of 𝑐 for the method of proportion and the method of moments tends to become much more substantial than that for the method of maximum likelihood.

Looking at the overall values of RMSE for the estimators of 𝑐 and 𝛽 obtained through the three methods, it is also evident that much more uncertainty is associated with the point estimation of the second parameter.

Figure 3 shows the MC distributions of the estimators of 𝑐 and 𝛽 for 𝑛=50 and three combinations of parameters. From the analysis, it emerges that the positive bias of 𝛽_ML and 𝛽_𝑀 is due to the presence of a certain number of samples providing estimates much larger than the true value of 𝛽, while the MC medians of both estimators are very close to it. These boxplots also emphasize the presence of implausible values for the estimates of 𝛽 yielded by the method of proportion (𝛽_𝑃<−1).

(a) 𝑐(̂𝑐=0.1,𝛽=0.75)

(b) 𝛽(𝑐=0.1,𝛽=0.75)

(c) 𝑐(̂𝑐=0.4,𝛽=0.25)

(d) 𝛽(𝑐=0.4,𝛽=0.25)

(e) 𝑐(̂𝑐=0.8,𝛽=−0.25)

(f) 𝛽(𝑐=0.8,𝛽=−0.25)

Figure 3 

Boxplot of Monte Carlo distributions of 𝑐̂ and 𝛽 for three parameter combinations (𝑛=50 and 5,000 MC runs). The dashed horizontal lines indicate the true value of the parameters 𝑐 and 𝛽.

Trying to synthesize all the results presented here—which are not, however, exhaustive—the method of proportion, despite its straightforward analytical derivation of the estimators, performs worse overall (especially in terms of variability) than the method of moments and the maximum likelihood method, and it is competitive only for some specific scenarios, namely, for medium/high values of 𝑐, where it can better exploit the information contained in the sample.

4. Applications to Real Data

The methods that have been illustrated and empirically investigated in the previous sections are applied to two datasets. The first one, called J1 [13], contains the number of failures of software observed over 62 weeks, and its frequency distribution is reported in Table 3. Assuming that the statistical distribution underlying the data is a type III discrete Weibull, it is possible to compute the estimates yielded by the three estimators described in Section 3, which are reported in Table 4. Note that the method of moments and the maximum likelihood method provide very similar estimates for 𝑐 and 𝛽; the method of proportion supplies an estimate of 𝛽 that is quite different from the other two methods (moreover, it is negative), while the estimate of 𝑐, 𝑐_𝑃̂ is quite close to 𝑐_𝑀̂ and 𝑐̂_ML. If we want to test the goodness of fit of the discrete Weibull model for the data, we can use the chi-square statistic 𝑇=∑^𝑘_𝑖=1(𝑓_𝑖−𝑛𝑝_𝑖)²/(𝑛𝑝_𝑖), where 𝑓_𝑖 denotes the category frequencies and 𝑝_𝑖 denotes the probability of an observation falling into the 𝑖th category under the study model, such that ∑^𝑘_𝑖=1𝑝_𝑖=1 and ∑^𝑘_𝑖=1𝑓_𝑖=𝑛. If the model holds, 𝑇 is asymptotically distributed as a chi-square with 𝑘−1−𝑝 degrees of freedom, where 𝑝 is the number of parameters to be estimated (in this case, 𝑝=2). Before computing the empirical value of 𝑇, under each model we have to group the categories in such a way that all the expected frequencies 𝑛𝑝_𝑖 are not smaller than 5. Table 5 reports these groupings. The goodness-of-fit chi-square test accepts the type III Weibull model with the ML and 𝑀 parameter estimates (statistic’s value of 3.8843 and 𝑃 value of 0.2742 for ML, 3.8514 and 0.2780 for 𝑀), while it refuses the model with the 𝑃 parameter estimates (statistic’s value of 10.2323 and 𝑃-value of 0.0060) at a 5% nominal level of significance. Again, there is a notable difference between the results for the method of proportion and the maximum likelihood method and the method of moments.

Table 3 

Empirical frequency distributions for the first example dataset.

Table 4 

Estimates of the parameters of the discrete Weibull for the first example dataset.

Table 5 

Empirical and theoretical frequency distributions for the first example dataset.

The second example considers the data regarding accidents to 647 women working on H. E. Shells during five weeks [14], whose frequency distribution is reported in Table 6. Even if these are not strictly failure data, in the meaning explained in Section 1, nevertheless they are count data, and the discrete Weibull r.v. can be used to model them. The parameter estimates yielded by the three methods are reported in Table 7. The estimates of 𝑐 are again very close to each other; all the estimates of 𝛽 are negative and quite small in absolute value. Note that all the pairs (𝑐,̂𝛽) fall outside the scenarios explored in the simulation study; nevertheless, the empirical distribution of the data (many 0’s and 1’s) is favorable to the method of proportion, and the large sample size (𝑛=647) should ensure that all three methods are reliable. Testing the goodness-of-fit of the type III discrete Weibull model with the ML parameters, which requires grouping the last two categories (then 𝑘=5), the value of the chi-square statistic with 𝑘−1−𝑝=2 degrees of freedom is 3.9595, and the corresponding 𝑃-value is 0.1381. Note that this model fits the data better than the natural negative binomial model considered by Greenwood and Yule [14] and the generalized Poisson model proposed by Consul and Jain [15].

Table 6 

Empirical frequency distribution for the second example dataset.

Table 7 

Estimates of the parameters of the discrete Weibull for the second example dataset.

5. Conclusion

This paper examined three estimators for the parameters of the type III discrete Weibull random variable, which represents an alternative distribution to the geometric and negative binomial for modeling discrete reliability data, and can ensure increasing and decreasing failure rates. One of the methods presented (method of proportion) has recently been introduced in discrete models, and here it is newly phrased; it provides a closed form for the estimates of both parameters. The other two methods (maximum likelihood method and method of moments) are standard approaches for estimating parameters, but due to the complex expression of the probability mass function, they provide the estimates as a numerical solution to a minimization/maximization problem. It follows that not as much can be said about the statistical properties of the estimators for a finite sample size; then, an extensive Monte Carlo simulation study was carried out to assess their behavior. Far from giving a definitive solution to the problem, the study highlighted that the method of proportion, when applicable, can provide reliable estimates even for small sample sizes only under specific parameter configurations, whereas under other configurations it may provide poor results, especially in terms of the accuracy of the estimator of the second parameter. The other two standard methods can be usefully adopted under most scenarios; caution is necessary since some parameter combinations may lead to nonnegligible bias of the corresponding estimators. The paper stressed the potential practical applicability limits of each method, also in terms of computational burden, and illustrated their use through two examples with real data.

Appendix

A. Existence of First- and Second-Order Moments

Let us denote with 𝑎_𝑛, 𝑛≥1, the general element of the series in (6), 𝑎_𝑛=exp[−𝑐∑^𝑛_𝑗=1𝑗^𝛽]. 𝔼(𝑋) exists finite if and only if the series is convergent. If 𝛽≥0, the series is convergent according to the ratio criterion, since the limit 𝑙₁=lim_𝑛→∞𝑎_𝑛+1/𝑎_𝑛=lim_𝑛→∞exp[−𝑐(𝑛+1)^𝛽] is 0 if 𝛽>0 and is exp(−𝑐)<1 if 𝛽=0.

If −1<𝛽<0, note that ∑^𝑛_𝑗=1𝑗^𝛽>∫^𝑛₁𝑗^𝛽𝑑𝑗=𝑛^𝛽+1/(𝛽+1)−1/(𝛽+1), and then 𝑒^{−𝑐∑𝑛𝑗=1𝑗𝛽}<𝑒^{−𝑐[𝑛1+𝛽/(1+𝛽)−1/(1+𝛽)]}=𝑒^{𝑐/(1+𝛽)}𝑒^{−𝑐𝑛1+𝛽/(1+𝛽)}. (A.1)

The series with the minoring term is convergent (e.g., using the comparison criterion with the harmonic converging series) and thus the original series in (6), too.

If 𝛽=−1, note that ∑^𝑛_𝑗=11/𝑗>∫^𝑛₁1/𝑗𝑑𝑗=log𝑛 and then 𝑒^{−𝑐∑𝑛𝑗=11/𝑗}<𝑒^{−𝑐log𝑛}=𝑛^−𝑐. (A.2)

The series with the minoring term is convergent for 𝑐>1 and thus the original series in (6) is too. Moreover, since ∫^𝑛₁1/𝑗𝑑𝑗=log𝑛>∑^𝑛_𝑗=21/𝑗 for 𝑛≥2𝑒^{−𝑐∑𝑛𝑗=11/𝑗}=𝑒^{−𝑐(1+∑𝑛𝑗=21/𝑗)}=𝑒^−𝑐𝑒^{−𝑐∑𝑛𝑗=21/𝑗}>𝑒^−𝑐𝑒^{−𝑐log𝑛}=𝑒^−𝑐𝑛^−𝑐. (A.3)

The series with the minoring term is divergent for 𝑐≤1 and thus the original series in (6) is too.

If 𝛽<−1, then lim_𝑛→∞𝑎_𝑛 is greater than zero, and the series is divergent.

Let us now denote with 𝑏_𝑛, 𝑛≥1, the general element of the series in (7), 𝑏_𝑛=2𝑛exp[−𝑐∑^𝑛_𝑗=1𝑗^𝛽]; then 𝔼(𝑋²)=∑^∞_𝑛=1𝑏_𝑛−𝔼(𝑋). If 𝛽≥0, the series is convergent according to the ratio criterion, since the limit 𝑙₂=lim_𝑛→∞𝑏_𝑛+1/𝑏_𝑛=lim_𝑛→∞(𝑛+1)/𝑛⋅exp[−𝑐(𝑛+1)^𝛽] is 0 if 𝛽>0 and is exp(−𝑐)<1 if 𝛽=0; thus 𝔼(𝑋²) is finite.

If −1<𝛽<0, 2𝑛𝑒^{−𝑐∑𝑛𝑗=1𝑗𝛽}<2𝑛𝑒^{−𝑐[𝑛1+𝛽/(1+𝛽)−1/(1+𝛽)]}=2𝑒^{𝑐/(1+𝛽)}𝑛𝑒^{−𝑐𝑛1+𝛽/(1+𝛽)}. (A.4) The series with the minoring term is convergent (e.g., using the comparison criterion with the harmonic converging series) and thus the original series in (7) is too. Then 𝔼(𝑋²) is finite.

If 𝛽=−1, 2𝑛𝑒^{−𝑐∑𝑛𝑗=11/𝑗}<2𝑛𝑒^{−𝑐log𝑛}=2𝑛^−𝑐+1. (A.5) The series with the minoring term is convergent for 𝑐>2 and thus the original series in (7) is too. Then 𝔼(𝑋²) is finite. Moreover, for 𝑛≥2, 2𝑛𝑒^{−𝑐∑𝑛𝑗=11/𝑗}=2𝑛𝑒^{−𝑐(1+∑𝑛𝑗=21/𝑗)}=2𝑒^−𝑐𝑛𝑒^{−𝑐∑𝑛𝑗=21/𝑗}>2𝑒^−𝑐𝑛𝑒^{−𝑐log𝑛}=2𝑒^−𝑐𝑛^−𝑐+1. (A.6)

The series with the minoring term is divergent for 𝑐≤2, and thus the original series in (7) is too. 𝔼(𝑋²) is not finite.

For 𝛽<−1, 𝔼(𝑋²) is not finite, since 𝔼(𝑋) is not finite and 𝔼(𝑋²)≥𝔼(𝑋).

B. Derivation of 𝛽_𝑃

Equating the probability 𝑃(𝑋=1)=exp(−𝑐)[1−exp(−𝑐2^𝛽)] to the sample proportion 𝑧/𝑛, we get exp−𝑐2^𝛽𝑧=1−𝑛exp(𝑐), (B.1) and after taking the natural logarithm of both sides of the equation twice, the second time after changing signs, we obtain ()(𝑐𝛽=log−log1−𝑧/𝑛exp)/𝑐.log2 (B.2) Substituting 𝑐 with its estimate (8) in (B.2) we get ()/𝛽=loglog1−𝑧/𝑛1−𝑦/𝑛/log1−𝑦/𝑛,log2 (B.3) and further simplifying it 𝛽=loglog1−𝑦/𝑛−𝑧/𝑛/log1−𝑦/𝑛−1,log2 (B.4) which represents the estimate 𝛽_𝑃.

Acknowledgments

The author would like to thank the editor and an anonymous reviewer for their suggestions on the early version of the paper.