A Family of Lifetime Distributions

Abstract

A four-parameter family of Weibull distributions is introduced, as an example of a more general class created along the lines of Marshall and Olkin, 1997. Various properties of the distribution are explored and its usefulness in modelling real data is demonstrated using maximum likelihood estimates.

1. Introduction

Probability distributions are often used in survival analysis for modeling data, because they offer insight into the nature of various parameters and functions, particularly the failure rate (or hazard) function. Throughout the last decades, a considerable amount of research was devoted to the creation of lifetime models with more than the classical increasing and decreasing hazard rates; apparently, the motivation for this trend was to provide with more freedom of choice in the description of complex practical situations (see e.g., [1–9], and the references therein). In this paper a general class of models is introduced, by adding an extra parameter to a distribution in the sense of Marshall and Olkin [10], and subsequently used in developing a four-parameter modified Weibull extension distribution, with various failure rate curves that compete well with other alternatives in fitting real data. Specifically, Xie et al. [11] generalized the Chen [12] distribution by adding the lacking scale parameter, thus creating a three-parameter Weibull distribution; although the variety of shapes of the reliability curves was not enriched, the resulting model provided better fit to real data. The proposed distribution extends the Xie et al. [11] distribution by adding a shape parameter; it will be seen that compared to the previous and other models, the cost of the addition is balanced by the improvement in fitting real data.

The paper is organized as follows. Section 2 includes the general class of models and some properties. The proposed four-parameter Weibull model is introduced in Section 3 and some properties and reliability aspects are studied. The parameters are estimated by the method of maximum likelihood and the observed information matrix is obtained; the fit of the proposed distribution to two sets of real data is examined against three and two parameter competitors.

2. The Class of Distributions

It is possible to generalize a distribution by adding a shape parameter, in the sense of Marshall and Olkin [10]. Thus, starting with a distribution with survival function 𝑠₀, the survival function of the proposed family with the additional parameter 𝑝 is given by𝑠(𝑥)=ln1−1−𝑝𝑠₀(𝑥)ln𝑝,𝑥∈ℝ,𝑝∈ℝ₊−{0}, (1) and when 𝑝→1, then 𝑠→𝑠₀. The probability density and hazard functions are readily found to be𝑓(𝑥)=𝑝−1𝑓₀(𝑥)1−1−𝑝𝑠₀(𝑥,ℎ(𝑥)=𝑠)ln𝑝𝑝−1₀(𝑥)ℎ₀(𝑥)1−1−𝑝𝑠₀(𝑥)ln1−1−𝑝𝑠₀(𝑥,) (2) where 𝑓₀ and ℎ₀ are the probability density and hazard functions corresponding to the distribution with survival function 𝑠₀, and 𝑥∈ℝ, 𝑝∈ℝ₊−{0}. Since,lim_{𝑥→−∞}ℎ(𝑥)=𝑝−1𝑝ln𝑝lim_{𝑥→−∞}ℎ₀(𝑥),lim_𝑥→+∞ℎ(𝑥)=lim_𝑥→+∞ℎ₀(𝑥) (3) it follows from (2) that𝑝−1ℎ𝑝ln𝑝⁰(𝑥)≤ℎ(𝑥)≤ℎ₀(𝑥)ℎ,𝑥∈ℝ,𝑝≥1,₀(𝑥)≤ℎ(𝑥)≤𝑝−1ℎ𝑝ln𝑝⁰(𝑥),𝑥∈ℝ,𝑝∈(0,1]. (4) Therefore, ℎ(𝑥)/ℎ₀(𝑥) with 𝑥∈ℝ is increasing for 𝑝≥1 and decreasing for 𝑝∈(0,1]. When 𝑠₀(0)=1, the hazard function at the origin, ℎ(0), behaves quite differently than the corresponding functions for the Weibull and gamma distributions; for both these families, the distribution can be exponential, or ℎ(0)=0, or ℎ(0)=∞, so that ℎ(0) is discontinuous in the shape parameter. This is not the case for the hazard functions in (2), and therefore the proposed family may be useful in fine-tuning the distribution with survival function 𝑠₀.

3. A Weibull Extension Model

The survival function of the modified Weibull extension distribution, introduced by Xie et al. [11] and studied further by Tang et al. [13], is given by𝑠₀(𝑥)=𝑒^{−𝛼𝜆(𝑒𝛽(𝑥/𝛼)−1)}, (5) for 𝑥,𝛼,𝛽,𝜆∈ℝ₊−{0}; hereinafter we shall be referring to the distribution with survival function given by (5) as the XTG distribution for brevity. By substituting (5) in (1), the survival function of the proposed distribution is obtained in the form𝑠()=𝑒𝑥;𝜃ln1−1−𝑝^{−𝛼𝜆(𝑒𝛽(𝑥/𝛼)−1)}ln𝑝, (6) where 𝜃=(𝛼,𝛽,𝜆,𝑝) and 𝑥,𝛼,𝛽,𝜆,𝑝∈ℝ₊−{0}. From (2) or (6), the pdf is readily obtained as𝑓()=()𝑥;𝜃𝛽𝜆𝑝−1𝑥/𝛼^𝛽−1𝑒^{(𝑥/𝛼)𝛽−𝛼𝜆(𝑒𝛽(𝑥/𝛼)−1)}1−1−𝑝𝑒^{−𝛼𝜆(𝑒𝛽(𝑥/𝛼)−1)}ln𝑝, (7) where 𝛼 is ascale parameter and 𝛽, 𝜆, 𝑝 are shape parameters. It can be shown that the pdf is monotone decreasing, unimodal or even roller-coaster type; the different shapes of the pdf are illustrated in Figure 1 for selected values of the parameters.

Figure 1 

Probability density functions of the EXTG distribution for 𝛼=2,𝛽=0.5,𝜆=0.03,𝑝=0.2 (solid line), 𝛼=12,𝛽=4,𝜆=2,𝑝=0.2 (dashed line), and 𝛼=2,𝛽=0.7,𝜆=0.3,𝑝=0.2 (dotted line).

Clearly, for 𝑝→1 the proposed distribution reduces to the XTG distribution therefore the proposed model can be viewed as an extension of the XTG model (which is asymptotically related to the usual two-parameter Weibull distribution) and if, in addition, 𝛼=1, then (7) defines the Chen [12] distribution; hereinafter we shall be referring to this distribution as the extended XTG distribution (EXTG distribution for brevity). Furthermore, it can be shown that for 𝑝∈(0,1) (7) is a compound of the logarithmic and the XTG distributions. Indeed, by incorporating the results of Barlow and Proschan [14] and Arnold et al. [15], consider the lifetime 𝑋=min(𝑋₁,𝑋₂,…,𝑋_𝑍) of a “series-system” of 𝑍 identical components, where failure occurs if at least one component ceases to function. If the lifetimes of the components are iid random variables with survivals given by (5) and the distribution of their number 𝑍 is logarithmic, independently of the 𝑋’s, with pmf𝑝𝑧;𝑝=−1−𝑝^𝑧𝑧ln𝑝, (8) for 𝑧∈ℕ−{0}, 𝑝∈(0,1), then the distribution of 𝑋∣𝑍 has pdf𝑓𝑥𝑥∣𝑧;𝛽,𝜆=𝛽𝜆𝑧𝛼^𝛽−1𝑒^{(𝑥/𝛼)𝛽−𝛼𝜆(𝑒𝛽(𝑥/𝛼)−1)}, (9) for 𝑥,𝛼,𝛽,𝜆∈ℝ₊−{0}, and the distribution of 𝑋 is the EXTG with pdf given by (7).

The calculations of the 𝑟th raw moments of the EXTG distribution involve the use of standard numerical integration procedures available in every mathematical package; for 𝑝∈(0,1) they can be expressed in the form,𝐸𝑋^𝑟1=−ln𝑝^∞_𝑘=11−𝑝^𝑘𝑘^∞₀𝛼^𝑟𝑥ln𝛼𝜅𝜆+1^𝑟/𝛽𝑒^−𝑥𝑑𝑥. (10)

By straightforward reversal of the cdf, obtained from (6) using that 𝐹(𝑥;𝜃)=1−𝑆(𝑥), the quantile function is calculated to be𝐹⁻¹𝑞1=𝛼ln1−𝛼𝜆ln1−𝑝^1−𝑞1−𝑝^1/𝛽, (11) for 𝑝∈(0,1); hence the median is 𝑀=𝛼[ln{1+(1/𝛼𝜆)ln(1+√𝑝)}]^1/𝛽.

3.1. Failure Rate and Mean Residual Life Functions

From (6) and (7) the failure rate (also known as hazard rate) function of the EXTG distribution isℎ()=()𝑥;𝜃𝛽𝜆𝑝−1𝑥/𝛼^𝛽−1𝑒^{(𝑥/𝛼)𝛽−𝛼𝜆(𝑒𝛽(𝑥/𝛼)−1)}1−1−𝑝𝑒^{−𝛼𝜆(𝑒𝛽(𝑥/𝛼)−1)}ln1−1−𝑝𝑒^{−𝛼𝜆(𝑒𝛽(𝑥/𝛼)−1)}. (12) It can be shown that for 𝛽≥1, the EXTG is an IFR distribution [16]. However, for 𝛽<1 it can be IFR, DFR, and BTFR distribution, although it is not easy to determine analytically the ranges of the parameter values; the IFR, DFR, and BTFR characteristics are depicted in Figure 2 for selected values of the parameters. Given that there is no failure prior to 𝑥₀, the residual life is the period from time 𝑥₀ until the time of failure. The mean residual lifetime, for 𝑝∈(0,1), is𝑚𝑥₀;𝜃=𝐸𝑋−𝑥₀∣𝑋≥𝑥₀=^∞_𝑥0𝑥−𝑥₀𝑓(𝑥;𝜃)1−𝐹𝑥₀1;𝜃𝑑𝑥=−ln1−1−𝑝𝑒^{−𝛼𝜆(𝑒0𝛽(𝑥/𝛼)−1)∞}_𝑘=11−𝑝^𝑘𝑘×^∞_{𝛼𝑘𝜆(𝑒0𝛽(𝑥/𝛼)−1)}𝛼𝑥ln𝛼𝑘𝜆+1^1/𝛽−𝑥₀𝑒^−𝑥𝑑𝑥. (13) Other reliability aspects of the distribution can be obtained numerically. For example, the renewal function, which is important for maintenance, can be calculated approximately either by the well-known method of considering its limit at infinity and the first and second raw moments given by (10), or by applying the method of the linear combination of the cdf and the hazard function; see Cui and Xie [17], Jiang [18] and the references therein.

Figure 2 

Hazard rate functions of the EXTG distribution for 𝛼=0.05,𝛽=0.45,𝜆=0.002,𝑝=20 (solid line), 𝛼=0.2,𝛽=0.2,𝜆=3,𝑝=2 (dashed line), and 𝛼=10,𝛽=0.9,𝜆=0.5,𝑝=200 (dotted line).

3.2. Inference

Assuming a random sample of n observations, 𝑦_obs=(𝑥_𝑖;𝑖=1,…,𝑛), from (7), the log-likelihood is given byℓ𝜃;𝑦_obs+=𝑛ln𝑝−1𝛽𝜆−lnln𝑝+𝛼𝜆+1−𝛽ln𝛼𝛽−1^𝑛_𝑖=1ln𝑥_𝑖+^𝑛_𝑖=1𝑥_𝑖𝛼^𝛽−𝛼𝜆^𝑛_𝑖=1𝑒^{(𝑥𝑖/𝛼)𝛽}−^𝑛_𝑖=1ln1−1−𝑝𝑒^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)}, (14)

and by differentiating with respect to 𝜃 the gradients are𝜕ℓ𝜕𝛼=𝑛𝜆−𝑛𝛼⁻¹𝛽−1−𝛼⁻¹𝛽^𝑛_𝑖=1𝑥_𝑖𝛼^𝛽−𝜆^𝑛_𝑖=1𝑒^{(𝑥𝑖/𝛼)𝛽}𝑥1+𝛽^𝑖𝛼^𝛽×−𝜆1−𝑝_𝑛𝑖=1𝑒^{(𝑥𝑖/𝛼)𝛽}1+𝛽𝑥_𝑖/𝛼^𝛽−1𝑒^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)}1−1−𝑝𝑒^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)},𝜕ℓ𝜕𝛽=𝑛𝛽⁻¹+^𝑛_𝑖=1𝑥1+^𝑖𝛼^𝛽𝑥−𝛼𝜆^𝑖𝛼^𝛽𝑒^{(𝑥𝑖/𝛼)𝛽}𝑥ln^𝑖𝛼−𝛼𝜆1−𝑝^𝑛_𝑖=1𝑥_𝑖/𝛼^𝛽ln𝑥_𝑖/𝛼𝑒^{(𝑥𝑖/𝛼)𝛽−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)}1−1−𝑝𝑒^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)},𝜕ℓ𝜕𝜆=𝑛𝛼+𝜆⁻¹−𝛼^𝑛_𝑖=1𝑒^{(𝑥𝑖/𝛼)𝛽}×−𝛼1−𝑝_𝑛𝑖=1𝑒^{(𝑥𝑖/𝛼)𝛽}−1𝑒^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)}1−1−𝑝𝑒^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)},𝜕ℓ𝜕𝑝=𝑛𝑝−1⁻¹+𝑝ln𝑝⁻¹−^𝑛_𝑖=1𝑒^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)}1−1−𝑝𝑒^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)}. (15) The mle of 𝜃, 𝜃, is obtained by solving simultaneously the four nonlinear normal equations, 𝜕ℓ/𝜕𝛼=0, 𝜕ℓ/𝜕𝛽=0, 𝜕ℓ/𝜕𝜆=0 and 𝜕ℓ/𝜕𝑝=0 by any iterative numerical method such as the Newton-Raphson, quasi-Newton, or Nelder-Mead procedures.

If the true parameter vector is an interior point of the parameter space, then 𝜃 can be treated as being approximately trivariate normal with mean 𝜃 and covariance matrix the inverse of the Fisher’s expected information matrix 𝐽(𝜃)=𝐸(𝐼;𝜃), where 𝐼(𝜃) is the observed information matrix with elements 𝐼_𝑖𝑗=−𝜕²ℓ/𝜕𝜃_𝑖𝜕𝜃_𝑗, 𝑖,𝑗=1,…,4 and the expectation is to be taken with respect to the distribution of 𝑋. By differentiating the normal equations in (15) the elements in the upper triangular part of 𝐼(𝜃) are found to be𝐼₁₁=−𝑛𝛼⁻²𝛽−1−𝛼⁻¹𝛽^𝑛_𝑖=1𝑥_𝑖𝛼^𝛽𝛼⁻¹𝛽+1−𝜆𝑒^{(𝑥𝑖/𝛼)𝛽}×𝛽𝑥_𝑖𝛼^𝛽+𝛽+1+𝜆1−𝑝^𝑛_𝑖=1𝑒^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)}1−1−𝑝𝑒^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)}×𝑒^{(𝑥𝑖/𝛼)𝛽}𝜆𝑒^{(𝑥𝑖/𝛼)𝛽}𝛽²𝑥_𝑖𝛼^2𝛽−1+2−𝛼⁻¹𝛽𝑥_𝑖𝛼^𝛽𝛽𝑥_𝑖𝛼^𝛽−𝜆𝑒+𝛽+11−𝑝^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)}𝑒^{(𝑥𝑖/𝛼)𝛽}𝑥1+𝛽^𝑖𝛼^𝛽−1²1−1−𝑝𝑒^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)},𝐼₁₂=𝑛𝛼⁻¹+𝛼^𝑛−1_𝑖=1𝑥_𝑖𝛼^𝛽𝑥1+𝛽ln^𝑖𝛼+𝜆𝑒^{(𝑥𝑖/𝛼)𝛽}×𝑥1+1+𝛽ln^𝑖𝛼𝑥+𝛽^𝑖𝛼^𝛽+𝜆1−𝑝^𝑛_𝑖=1𝑥_𝑖/𝛼^𝛽𝑒^{(𝑥𝑖/𝛼)𝛽}𝑒^{𝑖𝛽−𝛼𝜆(𝑒(𝑥/𝛼)−1)}1−1−𝑝𝑒^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)}×𝑥1+1+𝛽ln^𝑖𝛼𝑥+𝛽^𝑖𝛼^𝛽𝑥−𝛼𝜆ln^𝑖𝛼𝑒^{(𝑥𝑖/𝛼)𝛽}𝑥1+𝛽^𝑖𝛼^𝛽×𝑒−11+1−𝑝^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)}1−1−𝑝𝑒^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)},𝐼₁₃=−𝑛+^𝑛_𝑖=1𝑒^{(𝑥𝑖/𝛼)𝛽}𝑥1+𝛽^𝑖𝛼^𝛽+1−𝑝^𝑛_𝑖=1𝑒^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)}𝑒^{(𝑥𝑖/𝛼)𝛽}1+𝛽𝑥_𝑖/𝛼^𝛽−11−1−𝑝𝑒^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)2}×1−1−𝑝𝑒^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)}−𝛼𝜆𝑒^{(𝑥𝑖/𝛼)𝛽},𝐼−1₁₄=−^𝑛_𝑖=1𝑒^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)}1−1−𝑝𝑒^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)2},𝐼₂₂=𝑛𝛽⁻²−^𝑛_𝑖=1𝑥_𝑖𝛼^𝛽𝑥1−𝛼𝜆1+^𝑖𝛼^𝛽𝑒^{(𝑥𝑖/𝛼)𝛽}ln²𝑥_𝑖𝛼+𝛼𝜆1−𝑝^𝑛_𝑖=1𝑥_𝑖/𝛼^𝛽ln²𝑥_𝑖/𝛼𝑒^{(𝑥𝑖/𝛼)𝛽−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)}1−1−𝑝𝑒^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)2}×𝑥1+^𝑖𝛼^𝛽𝑥−𝛼𝜆^𝑖𝛼^𝛽𝑒^{(𝑥𝑖/𝛼)𝛽}−𝑥1−𝑝1+^𝑖𝛼^𝛽𝑒^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)},𝐼₂₃=𝛼^𝑛_𝑖=1𝑥_𝑖𝛼^𝛽𝑒^{(𝑥𝑖/𝛼)𝛽}𝑥ln^𝑖𝛼+𝛼1−𝑝^𝑛_𝑖=1𝑥_𝑖/𝛼^𝛽ln𝑥_𝑖/𝛼𝑒^{(𝑥𝑖/𝛼)𝛽−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)}1−1−𝑝𝑒^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)2}×1−𝛼𝜆𝑒^{(𝑥𝑖/𝛼)𝛽}−1−1−𝑝𝑒^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)},𝐼₂₄=−𝛼𝜆^𝑛_𝑖=1𝑥_𝑖/𝛼^𝛽ln𝑥_𝑖/𝛼^𝛽𝑒^{(𝑥𝑖/𝛼)𝛽−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)}1−1−𝑝𝑒^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)2},𝐼₃₃=𝑛𝜆⁻²−𝛼²1−𝑝^𝑛_𝑖=1𝑒^{(𝑥𝑖/𝛼)𝛽}−1²𝑒^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)}1−1−𝑝𝑒^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)2},𝐼₃₄=−𝛼^𝑛_𝑖=1𝑒^{(𝑥𝑖/𝛼)𝛽}−1𝑒^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)}1−1−𝑝𝑒^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)2},𝐼₄₄=𝑛𝑝−1⁻²−1+ln𝑝𝑝ln𝑝⁻²−^𝑛_𝑖=1𝑒^{−2𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)}1−1−𝑝𝑒^{−𝛼𝜆(𝑒𝑖𝛽(𝑥/𝛼)−1)2}. (16)

The latter is a consistent estimator of 𝐽(𝜃) and can be used for constructing asymptotic confidence intervals for the parameters. However, if any of the true parameter values is zero then the asymptotic distribution of the maximum likelihood estimators is a mixture distribution [19]; in this case obtaining the asymptotic confidence intervals becomes quite difficult and shall not be pursued here.

3.3. Examples

In this section two sets of real data are considered in order to test the goodness of fit of the proposed model. The first set of data consists of times to first failure of fifty devices [20]. The second set of data involves forty-four observations obtained from a life test concerning failure times (in hours) of all subsystems of a machine, that is, engine, hydraulic and air-conditioning subsystems, brakes, transmissions, tyres and wheels, body and chassis [21, 22]; in both cases, the data were grouped and the empirical hazard rate was estimated by many methods to indicate a BT shape. In addition to the EXTG, the XTG distribution, the two-parameter Chen [12] distribution, and the three-parameter model introduced by Dimitrakopoulou et al. [3], were fitted to the datasets; for brevity, hereinafter we shall be referring to the latter two models as the Chen and DAL distributions respectively. The fit of each distribution was examined by the Akaike information criterion (AIC) and the Kolmogorov-Smirnov (K-S) goodness-of-fit test using maximum likelihood estimates; the estimates, the maximized log-likelihoods, the values of the AIC, and the values of the K-S statistic with the associated 𝑃-values are presented in Table 1. Furthermore, the values of the likelihood ratio test statistic for testing 𝐻₀∶𝑝=1, calculated from the first and the second set of data, were 8.7939 (𝑃=0.003) and 6.8366 (𝑃=0.0089), respectively; the analogous computations for testing 𝐻₀∶𝑝=𝛼=1 were 11.8371 (𝑃=0.0027) and 6.479 (𝑃=0.0392). All the results indicate that the EXTG distribution describes these data better than the other models; these findings are also supported by the empirical and fitted survivor functions, plotted in Figure 3.

Table 1 

Parameter estimates, values of the log-likelihood (LL) and Akaike information criterion (AIC), and Kolmogorov-Smirnov (K-S) statistic obtained from the fit of each of the four distributions, to the times of first failure of devices (data set 1) and the machine’s subsystems lifetimes (data set 2).

(a)

(b)

Figure 3 

Reliability curves of the empirical distribution (starred line), the EXTG distribution (solid line), the XTG distribution (dashed line), the Chen distribution (dotted line), and the DAL distribution (dot-dashed line) for the times of first failure of fifty devices (a) and the forty-four failure times of a machine's subsystem (b).

Acknowledgment

The authors would like to thank a referee for useful comments and suggestions.