Inference on Reliability of Stress-Strength Models for Poisson Data

Abstract

Researchers in reliability engineering regularly encounter variables that are discrete in nature, such as the number of events (e.g., failures) occurring in a certain spatial or temporal interval. The methods for analyzing and interpreting such data are often based on asymptotic theory, so that when the sample size is not large, their accuracy is suspect. This paper discusses statistical inference for the reliability of stress-strength models when stress and strength are independent Poisson random variables. The maximum likelihood estimator and the uniformly minimum variance unbiased estimator are here presented and empirically compared in terms of their mean square error; recalling the delta method, confidence intervals based on these point estimators are proposed, and their reliance is investigated through a simulation study, which assesses their performance in terms of coverage rate and average length under several scenarios and for various sample sizes. The study indicates that the two estimators possess similar properties, and the accuracy of these estimators is still satisfactory even when the sample size is small. An application to an engineering experiment is also provided to elucidate the use of the proposed methods.

1. Introduction

A stress-strength model, in the simplest terms, considers a unit/system that is subjected to an external stress, modeled by r.v. 𝑋, against which the unit sets its own strength, modeled by r.v. 𝑌, in order to properly operate. The probability that the unit withstands the stress is then given by 𝑅=𝑃(𝑋<𝑌), which is usually called reliability.

A great deal of work has been done about this topic: most of it deals with the computation of reliability, if the distributions of stress and strength are known, or its estimation under various parametric assumptions on 𝑋 and 𝑌, when samples from 𝑋 and 𝑌 are available. A complete review is available in [1]. Many applications of the stress-strength model, for its own nature, are related to engineering or military problems, where it is also referred to as a load-strength model [2]. However, there are also natural applications in medicine or psychology, which involve the comparison of two r.v., representing, for example, the effect of a specific drug or treatment administered to two groups (control and test); here, reliability assumes a wider meaning.

Almost all of these papers consider continuous distributions for 𝑋 and 𝑌, since many practical applications of the stress-strength model in engineering fields presuppose continuous quantitative data. A relatively small amount of work is devoted to discrete or categorical data. Data may be discrete by nature, for example, the number of events occurring in a certain spatial or temporal interval; sometimes discrete data are derived from continuous ones by grouping or discretization or censoring, and then, instead of numerical measurements on 𝑋 and 𝑌, they are presented in a form of ordered categories.

Among the r.v. modeling discrete data, the Poisson can be of interest in several practical applications. The Poisson r.v. is often used to model rare events such as the number of claims in automobile insurance, the number of times a website is accessed, the number of calls to a phone operator, the number of words mistyped per page in a book, and so forth [3, 4]. The distribution of the difference between two independent r.v. each having a Poisson distribution has already attracted some attention [5]. Strackee and van der Gon [6] stated that “in a steady state the number of light quanta, emitted or absorbed in a definite time, is distributed according to a Poisson distribution. In view thereof, the physical limit of perceptible contrast in vision can be studied in terms of the difference between two independent variates each following a Poisson distribution”. Irwin [7] studied the case when the two variables 𝑋 and 𝑌 each have the same expected value; Skellam [8] was the first to discuss the problem when 𝐸(𝑋)=𝜆₁≠𝐸(𝑌)=𝜆₂. Strackee and van der Gon [6] gave tables of the approximate values of the cumulative probability 𝑃(𝑌−𝑋<𝑑) for several combinations of the values of the parameters 𝜆₁ and 𝜆₂. More recently, Karlis and Ntzoufras [9] used the Poisson difference distribution to model the difference in the decayed, missing, and filled teeth index before and after treatment; Karlis and Ntzoufras [10] applied it to model the difference in the number of goals in football games.

In this paper, we examine point and interval estimation for the reliability of the stress-strength model with independent Poisson stress and strength. Although the maximum likelihood (ML) and uniformly minimum variance unbiased (UMVU) estimators of reliability have a known analytical expression, their statistical properties cannot be easily derived and thus need to be assessed through a Monte Carlo simulation study. Confidence intervals for reliability based on approximate expression for variance are also presented, and their performances in terms of coverage rate and average width are empirically investigated.

The paper is laid out as follows: in Section 2 reliability for Poisson stress-strength model and its ML and UMVU estimators are presented and discussed. Section 3 introduces approximate variance estimators and confidence intervals for reliability. Section 4 is devoted to a Monte Carlo (MC) study, which empirically investigates the performance of ML and UMVU estimators, and the corresponding confidence intervals for different combinations of distributional parameters and sample sizes. Section 5 describes an application, and Section 6 gives final remarks.

2. Point Estimators

Let 𝑋 and 𝑌 be independent r.v. modeling stress and strength, respectively, with 𝑋∼Poisson(𝜆₁) and 𝑌∼Poisson(𝜆₂). Then, the reliability 𝑅=𝑃(𝑋<𝑌) of the stress-strength model is given by (see [1]) 𝑅=𝑃(𝑋<𝑌)=^+∞_𝑥=0𝜆^𝑥₁𝑒^−𝜆1𝑥!1−^𝑥_𝑦=0𝜆^𝑦₂𝑒^−𝜆2𝑦!=lim^𝑘_{𝑘→∞𝑥=0}𝜆^𝑥₁𝑒^−𝜆1𝑥!1−^𝑥_𝑦=0𝜆^𝑦₂𝑒^−𝜆2.𝑦! (1)

The terms of the external sum rapidly converge to zero: reliability can be actually computed taking into account only its first terms. As an example, we compute the reliability 𝑅 when 𝜆₁=1 and 𝜆₂=4; the partial sums are reported in Table 1: the value of 𝑅 is already stable at the 7th decimal digit when 𝑘=10.

Table 1 

Partial sums for the computation of 𝑅 for a Poisson stress-strength model (𝜆₁=1,𝜆₂=4).

If two simple random samples 𝐱 of size 𝑛₁ and 𝐲 of size 𝑛₂ from 𝑋 and 𝑌, respectively, are available, reliability can be estimated with the ML estimator, obtained by substituting in (1) the maximum likelihood estimators of the unknown parameters 𝜆₁ and 𝜆₂: 𝑅=^+∞_𝑥=0𝑥^𝑥𝑒^−𝑥𝑥!1−^𝑥_𝑦=0𝑦^𝑦𝑒^−𝑦𝑦!. (2)

Otherwise, one can use the UMVU estimator [1]: 𝑅=^𝑈_𝑥=0𝑛₁𝑥𝑥𝑛₁−1^{𝑛1𝑥−𝑥}𝑛^𝑛1𝑥₁1−^𝑥_𝑦=0𝑛₂𝑦𝑦𝑛₂−1^{𝑛2𝑦−𝑦}𝑛^𝑛2𝑦₂, (3) where 𝑈=min(𝑛₁𝑥,𝑛₂𝑦−1). Note that formula (3) is represented via a finite sum, whereas formula (2) contains a rapidly converging series. The number of calculations that formula (3) performs depends on the sample means and the sample sizes, which jointly define the number of terms of the external sum; in formula (2) the terms of the external sum rapidly converge to zero, so that it may practically need fewer calculations than (3).

In Table 2, the values for the UMVU estimator are reported when 𝑥=1 and 𝑦=2;3;4, for different combinations of sample sizes 𝑛₁ and 𝑛₂. Note that the values of 𝑅 are very close to the value of 𝑅 even for small sample sizes and get closer as the sample sizes increase. These results are pictorially displayed in Figure 1 for 𝑥=1 and 𝑦=4.

Table 2 

Values of the UMVU estimator of 𝑅 for a Poisson stress-strength model (𝑥=1,𝑦=2;3;4) with varying sample sizes.

Figure 1 

Values of the UMVU estimator of 𝑅 for a Poisson stress-strength model (𝑥=1,𝑦=4) with varying sample sizes. The ML estimator is represented by the horizontal line.

Due to the complex expressions involved, the bias of 𝑅 and the variance of either estimators 𝑅 and 𝑅 cannot be analytically derived; a comparison of their performance (in terms of mean square error) can be carried out through MC simulations.

3. Variance Estimators and Confidence Intervals

Whereas the exact value of the variance or the mean square error of either estimator introduced in Section 2 is almost impracticable to derive, an approximate value can be easily supplied recalling the delta method [11]. For the ML estimator 𝑅(𝑥,𝑦), since 𝑋 and 𝑌 are independent estimators of 𝜆₁ and 𝜆₂, the variance of 𝑅 can be approximated as 𝑉𝑅≈𝜕𝑅𝜕𝑥²_{𝑥=𝜆1,𝑦=𝜆2}𝑉(𝑥)+𝜕𝑅𝜕𝑦²_{𝑥=𝜆1,𝑦=𝜆2}𝑉𝑦, (4) with 𝑉(𝑥)=𝜆₁/𝑛₁, 𝑉(𝑦)=𝜆₂/𝑛₂ and 𝜕𝑅𝜕𝑥=^∞_𝑥=0𝑒^−𝑥𝑥^𝑥−1(𝑥!𝑥−𝑥)1−^𝑥_𝑦=0𝑦^𝑦𝑒^−𝑦,𝜕𝑅𝑦!𝜕𝑦=−^∞_𝑥=0𝑥^𝑥𝑒^−𝑥𝑥!^𝑥_𝑦=0𝑒^−𝑦𝑦^𝑦−1𝑦!𝑦−𝑦. (5)

An analogous approximation can be carried out for the variance of the UMVU estimator; remembering that 𝑥!=Γ(𝑥+1), 𝑅 can be rewritten as 𝑅=^𝑈_𝑥=0Γ𝑛₁𝑥+1Γ(𝑥+1)Γ𝑛₁𝑛𝑥−𝑥+1¹−1^{𝑛1𝑥−𝑥}𝑛^𝑛1𝑥₁×1−^𝑥_𝑦=0Γ𝑛₂𝑦+1Γ𝑦+1Γ𝑛₂𝑛𝑦−𝑦+1²−1^{𝑛2𝑦−𝑦}𝑛^𝑛2𝑦₂, (6) and then the two first-order partial derivatives are given by 𝜕𝑅𝜕𝑥=^𝑈_𝑥=01−^𝑥_𝑦=0𝑛₂𝑦𝑦𝑛₂−1^{𝑛2𝑦−𝑦}𝑛^𝑛2𝑦₂×1𝑛₁−1^𝑥𝑛₁𝑛₁𝑥𝑥𝑛₁−1𝑛₁^𝑛1𝑥×Γ^′𝑛₁𝑥+1Γ𝑛₁−Γ𝑥+1^′𝑛₁𝑥−𝑥+1Γ𝑛₁𝑛𝑥−𝑥+1+log¹−1𝑛₁,𝜕𝑅𝜕𝑦=−^𝑈𝑥_𝑥=0𝑦=0𝑛₁𝑥𝑥𝑛₁−1^{𝑛1𝑥−𝑥}𝑛^𝑛1𝑥₁×1𝑛₂−1^𝑦𝑛₂𝑛₂𝑦𝑦𝑛₂−1𝑛₂^𝑛2𝑦×Γ^′𝑛₂𝑦+1Γ𝑛₂−Γ𝑦+1^′𝑛₂𝑦−𝑦+1Γ𝑛₂𝑛𝑦−𝑦+1+log²−1𝑛₂. (7)

The approximate variances of 𝑅 and 𝑅 derived through the delta method can be estimated substituting in (4) the sample means to the unknown parameters and thus getting 𝑣𝑅=𝜕𝑅𝜕𝑥²𝑥𝑛₁+𝜕𝑅𝜕𝑦²𝑦𝑛₂ (8) and an analogous result for 𝑣(𝑅).

The Gamma function Γ(𝑥) and its first derivative, Γ′(𝑥)=∫₀^+∞𝑡^𝑥−1𝑒^−𝑡log𝑡𝑑𝑡, involved in the partial derivatives of 𝑅, have to be numerically computed. In the 𝑅 software environment [12] this task is easily performed through the gamma and digamma functions, the latter providing the ratio Γ′(𝑥)/Γ(𝑥).

Once one has computed 𝑣(𝑅), an approximate (1−𝛼)⋅100% confidence interval for 𝑅 can be built, recalling the asymptotic normality of 𝑅: 𝑅+𝑧_𝛼/2𝑣𝑅,𝑅+𝑧_1−𝛼/2𝑣𝑅, (9) and in an analogous way for 𝑅. Since 𝑅 is bounded in [0,1], special care has to be given when 𝑅 is close to one (close to zero) and/or sample sizes are small: the upper bound may exceed one (the lower bound may fall below zero), and then the CI in (9) will be modified as follows: max0,𝑅+𝑧_𝛼/2𝑣𝑅,min1,𝑅+𝑧_1−𝛼/2𝑣𝑅. (10)

More sophisticated asymptotic confidence intervals for 𝑅 can be built recalling some normalizing transformations, such as logit and arcsine [13].

4. Simulation Study

The simulation study aims at empirically comparing the performance of the ML and UMVU estimators, in terms of bias and mean square error, and the confidence intervals based on them, in terms of the coverage rate and average length. Since the approximation of the variance derived through the delta method (4) holds for large samples, we will investigate to what extent it still holds for small and moderate sample sizes, and how it affects inferential results. In this MC study, the value of the parameter 𝜆₁ of the Poisson distribution for stress 𝑋 is first set equal to a “reference” value, 1, and the parameter 𝜆₂ of the Poisson distribution modeling strength is allowed to vary in order to obtain four different levels of reliability 𝑅, namely, 0.5, 0.6, 0.7, 0.8, and 0.9. Note that a value of 𝜆₂=1.547 is needed in order to get 𝑅=𝑃(𝑋<𝑌)=0.5 while 𝜆₂=𝜆₁=1 lead only to 𝑃(𝑋<𝑌)+0.5⋅𝑃(𝑋=𝑌)=0.5. Then, 𝜆₁ is set equal to greater values (namely 2, 5, and 10), and 𝜆₂ is allowed to vary in order to ensure the five values of reliability 𝑅 above. The corresponding values of 𝜆₂ for each combination of 𝑅 and 𝜆₁ values are reported in Table 3.

Table 3 

Parameter values for the 𝑋 and 𝑌 distributions and corresponding reliability explored in the simulation study of a Poisson stress-strength model.

For each couple (𝜆₁,𝜆₂), a huge number (𝑆=2,000) of samples 𝐱 of size 𝑛₁ and 𝐲 of size 𝑛₂ are drawn from 𝑋∼Poisson(𝜆₁) and 𝑌∼Poisson(𝜆₂) independently. Different and unequal sample sizes are here considered (all the nine possible combinations between the values 𝑛₁=10,20,50, and 𝑛₂=10,20,50). The ML and UMVU estimators are computed on each sample, their approximate variances are calculated, and the corresponding 95% confidence intervals for 𝑅 are built. Some measures of performance for these estimators are supplied. In more detail, the MC root mean square error and the percentage relative bias of the ML estimator are provided: RMSE_MC𝑅=1𝑆^𝑆_𝑠=1𝑅(𝑠)−𝑅²RB_MC𝑅=(1/𝑆)∑^𝑆_𝑠=1𝑅(𝑠)−𝑅𝑅⋅100%, (11) where 𝑅(𝑠) denotes the value of 𝑅 for the 𝑠th sample. Analogous indexes are derived for 𝑅, whose bias is null, and for which we then expect the MC relative bias to be close to zero.

Regarding estimating the variance, the true variance 𝑉(𝑅) is approximated by its MC mean: 𝑉𝑅≈𝑉_MC𝑅=𝐸_MC𝑅−𝑅² (12) with 𝑅=∑^𝑆_𝑠=1𝑅(𝑠)/𝑆, and then the MC relative bias and RMSE of 𝑣(𝑅) are calculated the same way as for 𝑅.

The MC coverage rate of the CIs is simply defined as follows: 1𝑆^𝑆_𝑠=1𝐼𝑅(𝑠)+𝑧_𝛼/2𝑣𝑅(𝑠)≤𝑅≤𝑅(𝑠)+𝑧_1−𝛼/2𝑣𝑅(𝑠), (13) where 𝐼(𝐸) is the indicator function, taking value 1 if 𝐸 is true, 0 otherwise. The length of the confidence interval is then equal to 2𝑧_1−𝛼/2𝑣(𝑅(𝑠)). The same performance indexes are derived for 𝑅.

The simulation results for 𝜆₁=1 are reported in Table 4 (RB and RMSE for ML and UMVU point estimators), Table 5 (RB and RMSE for variance estimators), and Table 6 (coverage rate and average length of confidence intervals).

Table 4 

Simulation results: bias and root mean square error of ML and UMVU estimators.

Table 5 

Simulation results: bias and root mean square error of variance estimators.

Table 6 

Simulation results: coverage rate and average length of confidence intervals based on ML and UMVU estimators.

The simulation results show that the ML estimator 𝑅 always presents a very small bias even for small samples: in absolute value, the MC percentage relative bias is always smaller than 1.841% for all the scenarios considered (whereas the maximum absolute MC percentage relative bias for 𝑅, which is theoretically unbiased, is 0.552%). In 42 scenarios out of 45, 𝑅 underestimates 𝑅. Regarding the RMSE, the ML estimator performs better than UMVU in 27 cases out of 45, worse in 7 cases, and in 11 cases the RMSE is equal at the third decimal digit. However, under each scenario, even for smaller sample sizes, the values of RMSE for the ML and UMVUE estimators are very close. The ML outperforms the UMVU estimator as the value of 𝑅 gets close to 0.5; their performances tend to be alike as 𝑛₁ and 𝑛₂ increase. For both estimators, for fixed sample sizes, the RMSE increases as 𝑅 decreases; for a fixed 𝑅, the RMSE increases, as the sample sizes decrease (as expected). Figure 2 displays the MC distribution of the ML and UMVU estimators in the case 𝑅=0.7, for three values of sample size; it highlights their very similar behaviour.

(a) (𝑛1=𝑛2=10)

(b) (𝑛1=𝑛2=20)

(c) (𝑛1=𝑛2=50).

Figure 2 

MC distribution of 𝑅 (MLE) and 𝑅 (UMVUE) when 𝜆1=1, 𝜆2=2.497 (𝑅=0.7).

Regarding the approximate variance estimators, surprisingly their performance is good even for the moderate sample sizes considered in this study; the percentage relative bias, in absolute value, is never greater than 8%: the worst performance occurs for 𝑛₁=𝑛₂=50 and 𝑅=0.9. Indeed, when both sample sizes equal 50, the RB is greater than for small sample sizes, whereas one would expect that the RB decreases in absolute value when sample sizes increase. The results of further simulations not reported here show that the RB actually decreases to zero for 𝑛₁=𝑛₂=100. For both estimators, the rate of underestimates is almost equal to the rate of overestimates. Under each scenario, and especially when 𝑅=0.5, the value of RB of the variance estimator 𝑣(𝑅) is quite close to the corresponding value of the RB of 𝑣(𝑅), whereas the RMSE of 𝑣(𝑅) is smaller than the RMSE of 𝑣(𝑅) in each of the 45 cases considered. The RB of the two approximate variance estimators does not present a clear trend in terms of 𝑅; while their RMSEs, for each of the couples (𝑛₁,𝑛₂) here explored, seem to present a maximum near 𝑅=0.8 and a minimum for 𝑅=0.5.

The confidence intervals built upon the point estimators and these variance estimators present coverage that is always greater than 87% for UMVU and 90.5% for ML: the lowest value is obtained for 𝑛₁=𝑛₂=10 and 𝑅=0.9. They attain the nominal level (95%) for 𝑛₁=𝑛₂=50; in 31 and 21 cases out of 45, respectively, the coverage rate of the ML and UMVU interval estimators is greater than or equal to 92.5%. Overall, the CIs present better coverage when 𝑅 is close to 0.5. In fact, in this case, the distributions of 𝑅 and 𝑅 tend to be symmetrical and are more finely approximated by the normal distribution; then, the confidence intervals (9), which assume an underlying normal distribution, show a better performance. The CIs based on the ML estimator almost always show a coverage rate greater than those based on the UMVU estimator; moreover, the latter are always a bit wider, unless when 𝑅=0.9. This feature tends to be negligible when the sample sizes are increased. As one would expect, the average length decreases as sample sizes increase, for fixed 𝑅, and as 𝑅 increases, for fixed sample sizes.

The results for 𝜆₁>1, which are not reported here for the sake of brevity, confirm the previous findings. Even if the study is obviously not exhaustive, since only several scenarios have been covered, nevertheless these general features can be outlined.

5. An Example of Application

In this section, we apply the inferential techniques presented in Section 3 to a real dataset. The application is based on the data from an engineering experiment discussed in [3], carried out in an electric company, under several experimental conditions (called “runs”), corresponding to different combinations of 8 factors. The blackening experiment was conducted in a three-layer oven; when each run was completed, 30 masks from each layer in the oven were collected to examine the number of defects in each mask. The total number of defects in the 30 masks from the upper layer for each experimental run is observed (see Table 7). We focus on runs 1 and 2, where 317 and 184 defects, respectively, are counted.

Table 7 

Data for the blackening experiment in [3].

Since the number of defects in a mask is either zero or a positive integer, the appropriate distribution is the Poisson. Denoting with 𝑌 and 𝑋 the variables modeling this number for run 1 and run 2, respectively, we are interested in determining a point estimate and an interval estimator for the probability that the number of defects in run 1 is smaller than in run 2, that is, 𝑃(𝑋<𝑌). Since the sample size is 30 for both variables, then 𝑥=6.13 and 𝑦=10.57; then, supposing that 𝑋 and 𝑌 follow a Poisson distribution, the ML and UMVU estimators and their corresponding approximate variances can be computed according to (2), (3), and (8); the associated confidence intervals can be estimated recalling (10).

The results are presented in Table 8 and show the closeness between the two approaches. All the confidence intervals for 𝑃(𝑋<𝑌) always exclude 0.5, thus meaning that the difference in sample means testifies to the statistical dominance of 𝑌 on 𝑋: the number of defects under run 1 is stochastically larger than the number of defects under run 2.

Table 8 

Results for the application: ML and UMVUE point estimates and asymptotic confidence intervals (𝐿,𝑈) for 𝑅.

6. Conclusions

In this paper, point and interval estimators for the reliability of a Poisson stress-strength model are presented, discussed, and empirically compared through a Monte Carlo simulation study. The results show that the maximum likelihood and uniformly minimum variance unbiased estimators possess similar sampling properties, and the first is slightly preferable to the second in terms of dispersion around the true value of reliability. Moreover, although the variance estimators proposed here are approximate (i.e., biased), being based on the delta method for asymptotically normal r.v., the empirical results emphasize that the estimators' bias is small even for moderate sample sizes, and these estimators can be usefully employed to build approximate confidence intervals, whose coverage is shown to be overall close to the fixed nominal level, especially when the value of reliability is close to 0.5. However, when 𝑅 is close to 1 (or, symmetrically, 0), the intervals can show a poorer performance; then caution is needed when constructing a confidence interval based on a point estimate close to 1 (0). In this case, one can resort, for example, to some variance-stabilizing transformation of the estimate.

Acknowledgments

The author thanks the editor and two anonymous referees for their comments and suggestions on the original paper. Special thanks go to Riccardo Inchingolo for his moral support.