Bayesian and non-Bayesian reliability estimation of multicomponent stress–strength model for unit Weibull distribution

ABSTRACT Mazucheli et al. introduced a new transformed model referred as the unit-Weibull (UW) distribution with uni- and anti-unimodal, decreasing and increasing shaped density along with bathtub shaped, constant and non-decreasing hazard rates. Here we consider inference upon stress and strength reliability using different statistical procedures. Under the framework that stress–strength components follow UW distributions with a known shape, we first estimate multicomponent system reliability by using four different frequentist methods. Besides, asymptotic confidence intervals (CIs) and two bootstrap CIs are obtained. Inference results for the reliability are further derived from Bayesian context when shape parameter is known or unknown. Unbiased estimation has also been considered when common shape is known. Numerical comparisons are presented using Monte Carlo simulations and in consequence, an illustrative discussion is provided. Two numerical examples are also presented in support of this study. Significant conclusion: We have investigated inference upon a stress–strength parameter for UW distribution. The four frequentist methods of estimation along with Bayesian procedures have been used to estimate the system reliability with common parameter being known or unknown.


Introduction
Many physical phenomena from different disciplines are often modelled using some known probability models which include, among others, Burr family of distributions, lognormal, gamma, exponential, Weibull distributions. These probability distributions can model a variety of data exhibiting significant variability. One particular probability model of specific significance is the wellknown Weibull distribution. Over the last five decades or so, several new modifications of this particular model have been proposed and studied by many researchers, such as exponentiated Weibull [1][2][3], extended Weibull [4,5], modified Weibull [6][7][8], odd Weibull [9], Weibull-X class [10], Weibull-G model [11], extended Weibull-G family [12] and so on. These distributions are generally derived by adding some additional parameters to the original probability distribution under consideration. Besides, most of these generalizations share one interesting characteristic that they are based on the support over positive part of the real line. Such inferential trend often leads to insufficient number of distributions with finite support. But at the same time, probability distributions with support on finite range are also of importance in many studies. Physical characteristics of many life test experiments quite often lead to data which may lie in some finite range. Data on fractions, percentages, per capita income growth, fuel efficiency of vehicles, height and weight of individuals, survival times from a deadly disease etc. are likely to lie in some bounded positive intervals. Kumaraswamy [13] studied a two-parameter distribution with support on finite interval and investigated many useful applications of this distribution in meteorological inference. Mazucheli et al. [14,15] derived some useful classical estimates for unknown parameters of a unit-gamma distribution and also introduced unit Birnbaum-Saunders distribution, besides, different estimation methods are used to estimate the model parameters. Further, Mazucheli et al. [16] developed and studied properties of a unit Gompertz distribution and derived inferences for its unknown parameters. Mazucheli et al. [17] initially studied various characteristics of a unit Weibull distribution. Through several applications they showed that this particular model can be treated as an useful alternative to the Kumaraswamy and beta distributions in life test studies. We, here, consider multicomponent reliability estimation under unit Weibull probability distribution. This distribution can assume different shapesmonotonic, unimodal, anti-unimodal based on model parameters and such flexibility makes it quite suitable for fitting various data arising in reliability analysis and industrial experiments. Traditionally, a system consisting more than one component is referred to as multicomponent system, see [18,19]. Under this framework, a system with k independent strength components X 1 , X 2 , . . . , X k properly functions if at least s−out−of −k (s ≤ k) strength variables exceed the stress Y. Such a system is often called as ("s−out−of −k:G") system. Physical aspects of many investigations may lead to multi-component systems and many such examples abound in practice. Such structures are extensively utilized in industrial experiments, military radio networks, bridge construction, building communication networks, etc. In recent past, many scholars have studied multicomponent stress-strength models largely owing to the growing interest on this topic, see for instance, works of [20][21][22][23], Dey and Moala [24] and many others. Seadawy et al. [25], Seadawy et al. [26], Ahmad et al. [27], Riaz et al. [28], Abbasi et al. [29] have also studied some complex structures.
The probability distribution of a unit Weibull distribution, with range in the interval (0, 1), is given by and where α, β govern shape of UW distribution. We write this distribution as UW(α, β). Various frequentist procedures are discussed to estimate the unknown multicomponent reliability function. We further estimate this unknown parametric function using Bayes procedure against proper and improper prior distributions under a well-known loss function. Next, we have constructed asymptotic, two bootstrap and HPD intervals. Further, UMVUE of the multicomponent reliability is discussed as well under the case that common shape is known. To the our knowledge, this estimation problem based on unit Weibull distribution is not discussed before using the aforementioned methods of estimation.
Contribution of this work is many fold as UW distribution has not been discussed much in the literature. We derive inference for R s,k under various classical methods like maximum likelihood estimator (MLE), least square estimator (LSE) and weighted LSE (WLSE) and maximum product of spacing (MPS) estimators. These estimators are evaluated against estimated risks and average bias values computed using various sample sizes. We further estimate R s,k by Bayesian method based on proper and improper priors. Interval estimation is taken care as well. We have constructed asymptotic interval, bootstrap intervals and highest posterior density (HPD) interval. We compare these estimates using their length and coverage probabilities. The uniformly minimum variance unbiased estimator (UMVUE) of R s,k , with known common scale, is derived as well. Our aim is to provide a guideline for selecting the better inference procedure among different classical and Bayesian methods. This would be of deep interest to reliability engineers. The reliability of system is briefly illustrated in Section 2. In the same section, we describe four classical estimation methods, namely, MLE, LSE, WLSE and MPS to estimate multicomponent reliability and construct asymptotic and two bootstrap (boot−t and boot−p) CIs of multicomponent reliability when all parameters of the system are unknown. Further Bayes estimator is discussed when the prior distributions of β, α 1 and α 2 are gamma variables and also independent. In sequel, the HPD interval of multicomponent reliability is also constructed. Subsequently in Section 3, frequentist and Bayesian inference of multicomponent reliability are derived under the assumption that β is known. Simulations results are evaluated in Section 4 to examine numerical performance of proposed procedures. Two real-life examples are studied in Section 5. Finally, paper is concluded in Section 6.

Inference for R s,k with unknown β
Suppose X 1 , X 2 , . . . , X k denote strength components with cdf as given in (2). Likewise variable Y denotes associated stress acting on the system with distribution function G Y (y; α 2 , β). Then we have (see Rao et al. [19]) We next derive some useful inference of this parametric function using different procedures when common shape is unknown or known.

MLE of R s,k
Suppose that X 1 , X 2 , . . . , X k are taken from UW(α 1 , β) distribution and Y follows a UW(α 2 , β) distribution where parameters α 1 , α 2 and β are assumed as unknown. Under this framework, we get that After some simplification, we have We proceed by assuming that n structures are subjected to a test. The likelihood, based on observed data x i1 , x i2 , . . . , x ik and y i , of (α 1 , α 2 , β) is The log-likelihood turns out to be Partially differentiating the loglikelihood, respectively, w.r.t α 1 , α 2 and β and equating to zero yield following results. The MLEα 1 of α 1 iŝ Similarly, from we get the MLEα 2 of α 2 aŝ Finally, ∂l/∂β = 0 is solved for β to obtain its MLEβ satisfying The above likelihood equation is nonlinear and can be solved numerically using some iterative technique. Plugging MLE of β in Equations (6) and (7), we can obtainα 1 andα 2 , respectively. Now MLER s,k of R s,k turns out to be

Asymptotic interval
Here a confidence interval of the multicomponent reliability is given using MLE. The Fisher information of θ = (α 1 , α 2 , β) is The elements of I(θ ) are After some simplification, we get Similarly, we have From the large sample theory, the MLE of reliability is normal with average given by R s,k . The associated variance is computed as where we have and similarly ∂R s,k /∂α 2 is evaluated. Thus 100(1 − η)% CI of the parametric function is given by (R s,k ± q η/2σR s,k ), where q η/2 is associated percentile of normal N(0, 1) variable andσ R s,k is computed at respective MLEs.

Bootstrapping
We develop boot-p and boot-t intervals for R S,k (see Efron [30] and Hall [31] for further details).

Boot-p
(1) Obtain samples (y * 1 , . . . , y * n ) of size n and also simulate (x * i1 , . . . , x * ik ) of size nk where i is a positive integer ranging from 1 to n. Based on it, a sample estimate of reliability is computed asR * s,k . (2) Iterate above stage, say nboot times.

Boot-t
(1) Simulate bootstrap samples (y * 1 , . . . , y * n ) of size n and also generate (x * i1 , . . . , x * ik ) of size nk where i is a positive integer ranging from 1 to n. Based on it, a sample estimate of reliability is computed asR * s,k .

Least square and weighted least square estimates
The LSE method is initially discussed by Swain et al. [32]. Here we obtain this estimator for the multicomponent reliability R s,k . Note that LSEα 1 L of α 1 can be So LSEα l 1 is the solution of the following equation: where X 1 , X 2 , . . . , X nk follow UW(α 1 , β) distribution. The LSEα 2L of α 2 is computed similarly by solving the following equation: Next, weighted least square estimates (WLSE)α w 1 of α 1 is computed by minimizing the following function: and thus the required estimate is determined from the following equation: The WLSE estimateα w 2 of α 2 is similarly computed from the equation presented below: So WLSE of R s,k is given bŷ

Maximum product of spacing (MPS) estimates
Here we obtain another estimator called MPS estimates of the multicomponent reliability as discussed in [33]. Here spacing of a random sample of size nk is defined as also F X (x 0 ) = 0 and F X (x nk+1 ) = 1. The desired estimateα m 1 of α 1 is computed by maximizing the geometric mean ( nk+1 i=1 D i (α 1 )) 1/(nk+1) of spacing. Alternatively we maximize (1/(nk + 1)) nk+1 i=1 ln D i (α 1 ) and compute the desired estimate of α 1 from the following equation: Next using a sample of size n from the Y variable, . . , n, also F Y (y 0 ) = 0 and F Y (y n+1 ) = 1. The MPS estimateα m 2 of α 2 is determined from equation presented as Thus the MPS estimate of R s,k can be obtained aŝ

Bayesian inference
We discuss Bayesian inference of the parametric reliability function. Parameters (α 1 , α 2 , β) are a priori assumed to be independently distributed as gamma variables with hyperparameters (p i , q i ), where p i and q i , respectively, are shape and reverse of scale, index i takes values as 1, 2 and 3. The corresponding pdf is After some simplification, the joint posterior distribution is derived to be The normalizing constant can easily be computed. The required estimator is (10) The above integral is relatively difficult to be solved analytically. However, the corresponding posterior expectation can be approximated numerically. We now use [34] procedure and MH algorithm [35,36] to compute R s,k .

MH algorithm
This iterative procedure is widely used for evaluating Bayes estimates of different parametric functions such as R s,k , among others. The posterior of (α 1 , α 2 , β) given observations is derived in previous section. Thus the marginal posteriors of α 1 and α 2 are seen to be gamma distributed, respectively. While marginal posterior of other parameter β is not evaluated in known form. We simulate samples for β using normal proposal distribution. Accordingly, Bayes estimateR MH s,k of R s,k can be obtained from the following procedure.

2:
In ith iteration obtain sample β utilizing a normal distribution with mean β i−1 and variance σ 2 .
3: Obtain a sample α 1 from the gamma G(a, b) variable where a = (nk + p 1 ) and 6: Simulate a number u from a random variable U which follows a uniform model on (0, 1).

7:
We take the sample as 9: Now iterate stages 2 to 7 to obtain posterior estimates R i s,k , i = 1, 2, . . . , M. Now the required Bayes estimateR MH s,k is given bỹ where M 0 is burn-in samples. The HPD interval of R s,k is constructed from the posterior samples (see Chen and Shao [37]).

Inference for R s,k with known β
We now proceed to derive some more frequentist and Bayesian inference of unknown parametric function under assumption that β = β 0 where β 0 is a known constant. We mention that the expression of the reliability is identically same as considered in previous sections since it does not depend on β 0 .

MLE
Suppose that (x 1 , x 2 , . . . , x k ) denotes an observed value generated from the model as listed in Equation (2) where α 1 being unknown shape and β 0 being known shape. Similarly sample y is drawn when α 2 is unknown. The likelihood is then given by The log-likelihood function can be obtained as We partially differentiate the above function with respect to unknown parameters α 1 and α 2 , respectively. Then solving these likelihood equations, respective MLEs of α 1 and α 2 are obtained aŝ Thus the estimateR s,k of R s,k is evaluated using the invariance property. From large sample theory, we know that this MLE is normally distributed. The associated mean is given by R s,k . Similarly variance is

UMVUE of R s,k
Here UMVUER U s,k of the reliability R s,k is derived. It is enough to find the UMVUE of ξ(α 1 , α 2 ) = α 2 /(α 1 (j + k − i) + α 2 ) as UMVUEs satisfy linearity property. We see that (Y * , Z * ) is a complete sufficient statistic for (α 1 , α 2 ). Further with Y * has a G(nk, λ 1 ) pdf and Z * also has a G(n, is derived using the Lehmann-Scheffe theorem as follows: The conditional distributions involved in the previous equation is derived using [38]. We evaluate following cases: Case (i): The desired UMVUE is now obtained aŝ

Bayesian inference
Bayesian inference for the reliability is discussed here when β is known. It is assumed that α 1 and α 2 are independently distributed as gamma G(p 1 , q 1 ) and G(p 2 , q 2 ), respectively. The corresponding joint posterior distribution is given by The Bayes estimator has the following form: Next we try to evaluate the above expression. We take . Also assume that The above double integral is given by and p 0 = nk + n + p 1 + p 2 . Using hypergeometric series representation the above integral is expressed as when |η| < 1, Re(c) > 0 and Re()¯> 0, also |η| < 1 is the convergence region. We now havê In this particular case, Bayes estimate appears in closedform expression. Such estimates may also be computed from Lindley and MH methods as well. So for comparison purposes, Bayesian inference of the multicomponent reliability is derived from the latter two methods as well.

M-H algorithm
Note that marginal posteriors of α 1 and α 2 are gamma G(nk + p 1 , q 1 + y * ) and G(n + p 2 , q 2 + z * ), respectively. We generate samples using the following algorithm and then apply it to compute the Bayes procedure.
3: Obtain a sample α 2 using G(a * 5: Simulate a number u from a random variable U which follows a uniform model on (0, 1).

6:
We take the sample as α 1 i ← α 1 , α 2 i ← α 2 provided the inequality u ≤ m holds. 7: Compute R i s,k at (λ 1 i , λ 2 i , β 1 ). 8: Now iterate stages 2 to 7 to obtain posterior estimates R i s,k , i = 1, 2, . . . , M. Considering M 0 as discarded samples, the estimated value of the parameterR MH s,k turns out to beR MH The credible intervals are evaluated similarly as computed for unknown β case. (1) Here simulation experiments are performed to compare the performance of various methods proposed for estimating the reliability function. This study is performed against various sample sizes by assuming different values of β, α 1 , α 2 and prior distribution parameters. (2) Evaluation of estimators is done under mean square error (MSE) and average biases (ABs). (3) The performance of asymptotic, boot-t, bootp and HPD intervals is evaluated using interval length and coverage probabilities criteria.      given. All these estimation results are evaluated by considering 3000 replications. (9) We mention that in Tables 3, 4, 7 and 8, for a given n, two values are listed columnwise for each estimate. Among these, the first value denotes average length (AL) of interval of R s,k , immediate lower value denotes the coverage probability (CP) of that estimator. (10) In Tables 3 and 4, we have tabulated AL and coverage probabilities of various intervals. From these tables, it is observed that the average lengths of the intervals decrease with the increase in sample size, as expected. The average interval length of HPD interval is smaller than their counter parts.

Numerical experiments
(11) Next β known case is presented. Tables 5 and 6 contain estimation results of reliability under frequentist and Bayes techniques when β 0 is 2 and 3, respectively.   It is observed that in general, average length of HPD intervals are smaller than asymptotic intervals. We also observe that the length of both intervals tends to become smaller with large sample sizes. Also coverage probabilities of these intervals exhibit satisfactory behaviour.

Data analysis
Data Set I: We present a real-life numerical example in support of estimation procedures considered for evaluating the system reliability based on UW distribution.  The data represent breaking strength of Jute fibre [39]. We have considered 15mm gauge length data for analysis. The diameters of jute fibres are measured with an XSP-8CA digital biological microscope. Consider s = 1 and k = 5, then we observe that Y 1 becomes identical with the sixth observation of considered data. Thus we see that observations X 1k are from first to fifth observations. Similarly Y 2 is the 12th observation and also X 2k are from 7 to 11th data points. Here k is a positive integer ranging from 1 to 5. Data processing for total 30 observations gives n as 5. Observation (X, Y) is given below:   We also update observations X and Y by X ij /(max(X) + 1), Y i /(max(Y) + 1), respectively and then checked using goodness test whether the considered data sets can be fitted to the UW model. The MLEs of respective unknown parameters of all the competing distributions are evaluated numerically. The KS ( Kolmogorov-Smirnov) test statistic is taken into consideration for model evaluation. Table 9 contains associated probability values for all the competing models. It is seen that UW distribution can be considered a good representative model given data. We have also plotted P-P and Q-Q plot (Figure 1) of the given data sets and observe that data fitted satisfactorily. In Figure 2, we plotted curve of the profile likelihood function. This is concave in nature, i.e. unimodal function and maximum at β = 1.0001 for strength data as well as β = 0.7001 for stress data set. In Table 10 (1) Further equivalence testing between UW stress and strength parameters β 1 and β 2 is performed. The likelihood ratio test is applied to derive the inference. (2) H 0 : β 1 = β 2 versus H 1 : β 1 = β 2 . The test statistic is = −2(log 0 − log 1 ) which is χ 2 (1). (3) The estimates of corresponding log likelihoods under data I with H 0 : β 1 = β 2 and alternative hypotheses H 1 : β 1 = β 2 are log 0 = 1.607646 and log 1 = 1.979536, respectively. Here, we obtain = 0.7437795 with probability value 0.388478. (4) We see that null hypothesis β 1 = β 2 indicates favourable result under 5% significance. Thus we conclude that β 1 and β 2 may be equal for the considered numerical example.

Data Set II:
The second data set defines 12 core samples from petroleum reservoirs that were sampled by 4 crosssections. There are a total of 48 observations in this study. Core samples are measured for permeability and each cross-section defined using following components: total area of pores, total perimeter of pores and shape. The data set is available in R language [40], especially on data.frame named rock. Table 11 contains associated estimates for all competing models. We see that UW distribution is a good representative model     for the given data. We have also plotted P-P and Q-Q plots ( Figure 3) of the data. We observe that data fitted satisfactorily by the UW model. In Figure 4, we plotted curve of the profile likelihood function. This is concave in nature, i.e. unimodal function with maximum at β = 5 for strength data and β = 9 for stress data. Following [41], we further explore histogram plot of bootstrap samples for the logit of R s,k based on 3000 replications, see Figure 5. The histogram resembles of a normal distribution with approximate mean 1.727348 and standard deviation 0.5171634. Further associated K-S estimate (0.029392) and P-value (0.06314) also indicate normality of bootstrap samples for the logit of R s,k . In Table 12, estimates for the unknown reliability are evaluated for the complete data case. Associated intervals are also computed. We mention that Bayesian computations are performed against improper prior distribution. We consider s = 1 and k = 7 which is like 1-out-of-7: G system. Let Y 1 denote the 8th breaking stress and X 1k , k = 1, 2, . . . , 7, be the failure times of observations numbered 1 to 7. Further let Y 2 denote the 16th observation and X 2k be the failure times from 9 to 15. We carry on this data processing up to 48th failure, then we get n = 6 data for Y. The obtained data (X, Y) are as follows: (1) Further equivalence testing between UW stress and strength parameters β 1 and β 2 is performed for the data set II. The likelihood ratio test is applied to derive the inference. (2) The estimates of corresponding log likelihoods under data II with H 0 : β 1 = β 2 and alternative hypotheses H 1 : β 1 = β 2 are, log 0 = 58.38494 and, log 1 = 59.55056, respectively. Here, we obtain = 2.331241 with a probability value 0.126804. (3) We see that null hypothesis β 1 = β 2 indicates favourable result under 5% significance. Thus we conclude that β 1 and β 2 may be equal for the considered numerical example.

Conclusion
We have investigated inference upon a stress-strength parameter for UW distribution. The four frequentist methods of estimation have been used to obtain the system reliability with known and unknown β. Besides, asymptotic and two bootstrap CIs are obtained. Further Bayes inferences of R s,k are studied against different approximation techniques. The Bayes credible interval is discussed utilizing posterior samples generated from an MCMC method. Besides, unbiased estimation of reliability is also evaluated for known β case. Our numerical results suggest that Bayes procedure yields better inference compared to MLE and UMVUE. It is also seen that asymptotic, bootstrap and Bayesian intervals exhibit relatively good CP behaviour. However, width of Bayes intervals remains narrower than the corresponding asymptotic and bootstrap CIs. We present two numerical examples to demonstrate and observe the reliability for one system configuration. We found that WLSE provides marginally better results than the corresponding LSE estimates. Further inference obtained from the MPS method shows good behaviour than the WLSE method. We further note that MLE also improves LSE and WLSE. However, MLE and MPS methods are incomparable as in some case MLE performs better than MPS and in other case opposite behaviour is observed. We observe that if proper prior is available then Bayesian procedures have an advantage over the other methods. Overall, better inference of this particular parametric function can be derived under known β case. Although inference is studied when stress-strength components have the UW distributions, however, such studies can be conducted under assumption stress-strength components have different distributions. Moreover, estimation for the multicomponent stress-strength reliability for system with dependent components seems also of interest and importance in practice, which will be studied in the future. In near future, we also try to study such problems under some censoring schemes.