Classical and Bayesian Estimation of Reliability in Multicomponent Stress-Strength Model Based on Weibull Distribution

In this study, we consider a multicomponent system which has k independent and identical strength components X1, . . . , Xk and each component is exposed to a common random stress Y when the underlying distributions are Weibull. The system is regarded as operating only if at least s out of k (1 ≤ s ≤ k) strength variables exceeds the random stress. We estimate the reliability of the system by using frequentist and Bayesian approaches. The Bayes estimate of the reliability has been developed by using Lindley’s approximation and the Markov Chain Monte Carlo methods due to the lack of explicit forms. The asymptotic confidence interval and the highest probability density credible interval are constructed for the reliability. The comparison of the reliability estimators is made in terms of the estimated risks by the Monte Carlo simulations.


Introduction
In the reliability context, the stress-strength model can be described as an assessment of reliability of a system in terms of random variables X representing stress experienced by the system and Y representing the strength of the system available to overcome the stress.If the stress exceeds the strength, then the system will fail.Thus R = P (X < Y ) is the reliability of a system.The main idea was introduced by Birnbaum (1956) and developed by Birnbaum & McCarty (1958).Estimation of R = P (X < Y ) when the random variables X and Y follow a specified distribution has been extensively discussed by many authors in the literature.When the X and Y are independent and follow the generalized exponential, Weibull, three-parameter Weibull and Kumaraswamy distributions, the estimation of R was studied by Kundu & Gupta (2005), Kundu & Gupta (2006), Kundu & Raqab (2009) and Nadar, Kizilaslan & Papadopoulos (2014), respectively.Kotz, Lumelskii & Pensky (2003) provide an excellent review of the development of the stress-strength up to the year 2003.
The reliability in a multicomponent stress-strength model was developed by Bhattacharyya & Johnson (1974).This system consists of k independent and identical strengths component and a common stress, functions when s (1 ≤ s ≤ k) or more of the components simultaneously survive.This model corresponds to the s-out-of-k : G system.Its practical applications range from communication and industrial systems to logistic and military systems.For example, in suspension bridges, the deck is supported by a series of vertical cables hung from the towers.Suppose a suspension bridge consisting of k number of vertical cable pairs.The bridge will only survive if a minimum s number of vertical cable through the deck is not damaged when subjected to stresses due to wind loading, heavy trafic, corrosion, etc.As another example, with a V-8 engine of an automobile it may be possible to drive the car if only four cylinders are firing.However, if less than four cylinders fire, then the automobile cannot be driven.Thus, the functioning of the engine may be represented by a 4-out-of-8 : G system.Other examples include an electrical power station containing eight generating units which produce the right amount of electricity only if at least 6 units are working; the demand of the electricity of a district is fulfilled only if s-out-of-k wind rose are operating at all times; a communication system for a navy can be successful only if 6 transmitters out of 10 are operational to cover a district; a semi-trailer pulled by a truck can be driven safely as long as 6-out-of-8 tires are in good condition.For an extensive reviews of s-out-of-k and related systems see Kuo & Zuo (2003).
Let Y, X 1 , . . ., X k be independent, G(y) be the cumulative distribution function (cdf) of Y and F (x) be the common cdf of X 1 , . . ., X k .The reliability in a multicomponent stress-strength model is given by (1) This reliability can be written as the order statistics of (X 1 , . . ., X k ).This system reliability was considered by Jae & Eun (1981), when the stress and the strength distributions are Weibull with unknown scale parameters and the same known shape parameter.The maximum likelihood estimation (MLE) and the minimum variance unbiased estimation of system reliability were obtained in this study.The system reliability was considered by Hanagal (1999), when (X 1 , . . ., X k ) follow an absolutely continuous multivariate exponential distribution and Y follows an independent exponential distribution.Estimation of system reliability which is given as under the assumption of the strengths of the k components (X 1 , . . ., X k ) are subjected to an independent common stress X k+1 was considered by Hanagal (2003) when (X 1 , . . ., X k+1 ) follow (k + 1) independent of Gamma or Weibull or Pareto distributions.The stress-strength reliability of a system which has n independent components each consisting of m dependent elements was considered by Eryilmaz (2008).The reliability of stress-strength for a general coherent system was studied by Eryilmaz (2010).Recently, estimation of reliability in multicomponent stressstrength for the log-logistic, generalized exponential, generalized inverted exponential, Rayleigh, inverse Rayleigh and Burr Type XII distributions were considered by Rao & Kantam (2010), Rao (2012a), Rao (2012b), Rao (2012c), Rao, Kantam, Rosaiah & Reddy (2013) and Rao, Aslam & Kundu (2014), respectively.In these studies, maximum likelihood (ML), moment estimates and asymptotic confidence interval for the reliability in multicomponent stress-strength were obtained, but Bayesian approach was not taken into consideration.
In this study, we consider the multicomponent stress-strength model which has k independent and identical strength components and a common stress.We assume that the strength variables and stress variable follow Weibull distribution.The system functions if s (1 ≤ s ≤ k) or more of the components simultaneously operate.The estimation of reliability for this system is obtained under the classical and Bayesian frameworks.The Lindley's approximation and Markov Chain Monte Carlo (MCMC) technique are carried out to obtain the Bayes estimates.Moreover, the asymptotic confidence and the highest probability density (HPD) credible intervals are obtained.
A Weibull distribution with the shape parameter σ and scale parameter θ will be denoted by W E(σ, θ).The probability density function (pdf) and the cdf of a random variable X ∼ W E(σ, θ) are given as and F (x; σ, θ) = 1 − e −θx σ , x > 0, θ, σ > 0. (3) The rest of the paper is organized as follows.In Section 2, the MLE and the asymptotic confidence interval of R s,k are obtained when the parameters α, β and σ are unknown.In Section 3, the Bayes estimates of R s,k are obtained by using Lindley's approximation and MCMC method when the parameters α, β and σ have independent Gamma priors.The HPD credible interval for R s,k is also constructed.In Section 4, a simulation study is performed to compare the estimates of R s,k by using Monte Carlo simulations and findings are illustrated by tables and plots.Finally, we conclude the paper in Section 5.

Maximum Likelihood Estimation of R s,k
In our case, we assume that X 1 , . . ., X k be a random sample from Weibull distribution with parameters (σ, α) and Y be a random variable from Weibull distribution with parameters (σ, β).Therefore, for our case R s,k is given by using equations ( 1 In order to obtain the estimators of R s,k , suppose that n systems are put on life-testing experiment.In this case, we obtain the following observed data X i1 , X i2 , . . ., X ik and Y i , i = 1, . . ., n.Then, the likelihood function of the observed sample is given as Revista Colombiana de Estadística 38 (2015) 467-484 and the log-likelihood function is where The MLEs of α, β and σ are given by The MLE of σ, σ is the solution of the following nonlinear equation Therefore, σ can be obtained as a solution of the nonlinear equation of the form H(σ) = σ where Since, σ is a fixed point solution of the nonlinear equation ( 8), its value can be obtained using an iterative scheme as: where σ (j) is the jth iterate of σ.The iteration procedure should be stopped when σ (j) − σ (j+1) is sufficiently small.When we obtain σ, the MLEs of α and β are obtained from equation ( 6).Hence, the MLE of R s,k is obtained from equation (4) by using the invariance property of MLE's

Asymptotic Distribution and Confidence Interval for R s,k
The Fisher information matrix of θ = (θ 1 , θ 2 , θ 3 ) ≡ (α, β, σ) is given as The elements of the matrix are obtained as 2 ) are evaluated by using the formulas 4. 352(1) and 4.358(2 Hence, the other elements of the information matrix are given as The MLE of R s,k , R s,k is asymptotically normal with mean R s,k and variance where I −1 ij is the (i, j)th element of the inverse of the I(θ), see Rao (1965).Then, where and Therefore, an asymptotic 100(1 − γ)% confidence interval of R s,k is given by where z γ/2 is the upper γ/2th quantile of the standard normal distribution and σ R s,k is the value of σ R s,k at the MLE of the parameters.

Bayes Estimation of R s,k
In this section, we assume that all parameters α, β and σ are unknown and have independent Gamma prior distributions with parameters (c i , d i ), i = 1, 2, 3, respectively.The pdf of a Gamma random variable X with parameters (α, β) is Then, the joint posterior density function of α, β and σ is where Then, the Bayes estimator of R s,k under the SE loss function is given by It is not possible to compute equation ( 14) analytically.Two approaches can be applied to approximate equation ( 14), namely, Lindley's approximation and MCMC method.

Lindley's Approximation
Lindley (1980) introduced an approximate procedure for the computation of two integrals.This procedure employed to the posterior expectation of the function U (λ), for a given x, is where Q(λ) = l(λ) + ρ(λ), l(λ) is the logarithm of the likelihood function and ρ(λ) is the logarithm of the prior density of λ.Using the Lindley's approximation, ∂ρ/∂λ j , and σ ij = (i, j)th element in the inverse of the matrix {−L ij } all are evaluated at the MLE of the parameters.

MCMC Method
The joint posterior density function of α, β and σ is given in equation ( 13).It is easily seen that the posterior density functions of α, β and σ are, respectively, where x * and x σ are defined before the equation ( 6).
Therefore, samples of α and β can be easily generated by using Gamma distribution.However, the posterior distribution of σ cannot be reduced analytically to well known distributions and therefore it is not possible to sample directly by standard methods.It is observed that the plot of the posterior distribution is similar to Gaussian distribution.The hybrid Metropolis-Hastings and Gibbs sampling algorithm, which will be used to solve our problem, is suggested by Tierney (1994).This algorithm combines the Metropolis-Hastings with the Gibbs sampling scheme under the Gaussian proposal distribution.
Step 8. Repeat Steps 2-7, N times, and obtain the posterior sample R This sample is used to compute the Bayes estimate and to construct the HPD credible interval for R s,k .The Bayes estimate of R s,k under SE loss function is given as where M is the burn-in period.
The HPD 100(1 − γ)% credible interval of R s,k is obtained by the method of Chen & Shao (1999).

Simulation Study
In this section, we present some numerical results to compare the performance of the estimates of R s,k which is obtained by using different methods for different sample sizes and different priors.The performances of the point estimators are compared by using estimated risks (ERs).The performances of the confidence and credible intervals are compared by using average confidence lengths and coverage probabilities (cps).The estimated risk (ER) of θ, when θ is estimated by θ, is given by under the SE loss function.All of the computations are performed by using Matlab R2010a.All the results are based on 2000 replications.
In the MCMC case, we ran three MCMC chains with fairly different initial values and generated 10000 iterations for each chain.To diminish the effect of the starting distribution, we discard the first half of each sequence and focus on the second half.To provide relatively independent samples for improvement of prediction accuracy, we calculate the Bayesian MCMC estimates by the means of every 5 th sampled values after discarding the first half of the chains (see Gelman, Carlin, Stern & Rubin (2003)).The scale reduction factor estimate R = V ar(ψ) W is used to monitor convergence of MCMC simulations where ψ is the estimate of interest, V ar(ψ) = n−1 n W + 1 n B with the iteration number n for each chain, the between-and within-sequence variances B and W (see Gelman et al. (2003)).In our case, the scale factor value of the all MCMC estimators are found below 1.1 which is an acceptable value for their convergency.Note: The first row represents the average estimates and the second row represents corresponding ERs.The last two columns, the first row represents a 95% confidence interval and the second row represents their lengths and cp's.It is observed that the average ERs for the estimates of R s,k decrease as the sample size increases in all cases and all tables, as expected.The Bayes estimates of R s,k under the SE loss function have smaller ER than that of MLEs in all tables.Moreover, the ERs of the Bayes estimates which is obtained from Lindley's approximation are generally smaller than that of the MCMC method.Also, these ERs are close to each other for sufficiently large sample sizes.The average lengths of the intervals decrease as the sample size increases.The average lengths of the Bayesian credible intervals are smaller than that of the asymptotic confidence intervals.Moreover, the coverage probabilities of the Bayesian credible intervals are closet to the nominal level 95% than asymptotic confidence intervals.
On the other hand, to compare the performance of ML and Bayes estimates, we obtain the graphs for Biases and MSE's based on MLE and Lindley methods.Different sets of hyperparameters α, β and σ are generated from the corresponding Gamma distributions and their averages are used to obtain R s,k .This procedure is repeated 50 times to obtain different true values of R s,k when (s, k) = (1, 3), (2, 4) for given α, β and σ.After that we use the following algorithm for the comparison of the estimates 1.For given α, β and σ, we compute R 1,3 (or R 2,4 ).
2. For given n, we generate a sample from Weibull distribution for strength and stress variables.
3. The MLE and Bayes estimates are computed by using the equations ( 9) and (18).
The Figures 1-6 illustrate the Biases and MSEs of R s,k and R s,k,B for different sample sizes n = 10, 30, 50.It is observed that, the Biases and MSEs of the estimates decrease when the sample size increases, as expected.Moreover, when R s,k is around 0.5 the corresponding MSEs are large and when R s,k is small or large corresponding MSEs take small values for both estimates.

Conclusions
In this paper, we have studied the multicomponent system which hs k independent and identical strength components and each element exposed to a common random stress.We assume that the underlying distributions for both strength and stress variables are Weibull.The reliability of the system is estimated by using ML and Bayesian approaches.The Bayesian estimates are obtained by using Lindley's approximation and MCMC method.The simulation results indicate that the average ERs for the estimates of R s,k and the average lengths of the intervals decrease as the sample size increases.The ERs of the ML estimates are greater than that of Bayes estimates.The ERs of Bayes estimates which are obtained by using Lindley's approximation and MCMC method are close to each other for sufficiently large sample sizes.Due to the heavy computational burden of the MCMC method, Lindley's approximation maybe considered as a better alternative for large sample sizes.

4.
Steps 2-3 are repeated 2000 times, the Biases and MSE's are calculated and are given by Bias

Table 2 :
Estimates of R s,k by using Prior 2 for (α, β, σ) = (3.0387,3.3055, 2.9969).The first row represents the average estimates and the second row represents corresponding ERs.The last two columns, the first row represents a 95% confidence interval and the second row represents their lengths and cp's.

Table 3 :
Estimates of R s,k by using Prior 3 for (α, β, σ) = (2.0522,3.9554,3.0066).Note:The first row represents the average estimates and the second row represents corresponding ERs.The last two columns, the first row represents a 95% confidence interval and the second row represents their lengths and cp's