Reliability of repairable systems with Non-Central Gamma frailty

Brazilian Journal of Biometrics, 42, 182–201 (2024) doi:10.28951/bjb.v42i2.697 ARTICLE Reliability of repairable systems with Non-Central Gamma frailty ID ID Adriane Caroline Teixeira Portela,⋆ ,1 ID Éder Silva de Brito,⋆ ,3 ID Vera Lucia Damasceno Tomazella,2 Carlos Alberto Ribeiro Diniz,2 and ID Paulo Henrique Ferreira4 1 University of São Paulo, São Carlos, Brazil University of São Carlos, São Carlos, Brazil 3 Federal Institute of of Education, Science and Technology of Goiăs, Anăpolis, Brazil 4 Federal University of Bahia, Salvador, Brazil ⋆ Corresponding author. Email: [email protected]; [email protected] 2 Federal (Received: September 30, 2023; Revised: December 22, 2023; Accepted: December 28, 2023; Published: April 15, 2024) Abstract Maintenance actions on industrial equipment are essential to reduce expenses associated with equipment failures. Based on a well-fitted model, it is possible, through the estimated parameters, to predict several functions of interest, such as the cumulative average and reliability functions. In this paper, a new frailty model is proposed to analyze failure times of repairable systems subject to unobserved heterogeneity actions. The Non-Central Gamma distribution is assumed to the frailty random variable effect. The class of minimal repair models for repairable systems is explored considering an approach that includes the frailty term to estimate the unobserved heterogeneity over the systems’ failure process. Classical inferential methods were used to parameter estimation and define the reliability prediction functions. A simulation study was conducted to confirm the properties expected in the estimators. Two real-world data known in literature were used to illustrate the estimation procedures and validate the proposed model as a viable alternative to those already established in the literature. The results obtained highlight the potential of our proposed approach, particularly for industries dealing with such systems, where unquantifiable factors may impact equipment failure times. Keywords: Frailty models; Non-central Gamma distribution; Repairable systems; Minimal repair; Reliability prediction. 1. Introduction The reduction of expenses with actions for the maintenance of machinery(ies) is one of the main issues investigated in many industries that has led to the development of a huge literature dedicated © Brazilian Journal of Biometrics. This is an open access article distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/) Brazilian Journal of Biometrics 183 to the so-called reliability theory. With technological advances and the increasing complexity of systems, the study of their reliability has been of great interest, especially for repairable systems. When considering a system, be it a machine, electronic equipment or software, it is said to be repairable when after a failure its activity can be satisfactorily resumed through repair, without the need for its replacement as a whole, according to the definition given by Ascher & Feingold (1984). Otherwise, a system is non-repairable when after its single failure it is discarded. Therefore, after a failure occurs, it is possible for the repairable system to be restored to operation condition by a repair action. The study of repairable systems has wide application in contexts where the presence of recurring events is identified, such as in industry, where consecutive system failures are often observed. Investments in proper maintenance and repair to extend the lifetime of systems are essential, because in many cases repair is more economically feasible than replacing one or more components. The main challenge when modeling data from repairable systems is how to explain the effect of a repair action taken immediately after a failure. It is generally assumed that repair actions are instantaneous and the repair time is negligible (de Toledo et al., 2015). Among the usual models of repairable systems, the Minimal Repair (MR) model, known in the literature as As Bad as Old (ABAO), is widely used due to the speed of repair and the low associated cost. This repair focuses on correcting only the component originating the failure, consequently leaving the system in the same working condition as before the failure, according to assumptions discussed in different works such as Barlow & Proschan (1996), Wang (2002), and de Toledo et al. (2015). The statistical modeling usually used for these systems is done by counting process, which is a stochastic process that counts the number of occurrences of an event of interest in the time interval [0, t]. We can classify the counting process as Homogeneous Poisson Process (HPP) and Non-Homogeneous Poisson Process (NHPP), being that for the latter, the intensity function is not constant. A particular parametric NHPP case widely used in the literature and that will be used in this study is the Power Law Process (PLP) proposed by Crow (1975), which is attractive for the ease of interpretation of its parameters and simplicity of implementation. In the literature we found some methods for modeling data from repairable systems that do not incorporate the heterogeneity between the systems, and with this it is not possible to measure the most fragile, i.e., those that are more prone to fail (D’Andrea et al., 2017; Vaupel et al., 1979). It is relevant to verify how the heterogeneity of each system can impact on the main measures being investigated, such as increased failure rates (Proschan, 1963). Frailty models are an alternative to estimate the unobserved heterogeneity of systems and are characterized by incorporating a random effect, either additively or multiplicatively, introduced by Vaupel et al. (1979). In the recent reliability literature, there has been increasing interest in data analysis of repairable systems incorporating frailty, with different parametric and non-parametric distributions assumed for frailty. The choice of distribution for an unobservable random variable (frailty) is a problematic part when we are in an approach that considers this term in the model. Incorrect specification of the frailty distribution can lead to errors in the estimates of the parameters of interest, as is exposed in Almeida et al. (2020). The study considered a non-parametric approach to modeling the frailty density. By using Bayesian methods, an a prior Dirichlet process mixture is assigned to frailty, which in turn is more flexible and provides consistent estimates for the PLP model, as well as insights into heterogeneity among systems (Ascher & Feingold, 1984). In particular, Slimacek & Lindqvist (2016) present a new method that seeks to avoid numerical problems often encountered in the estimation of models that incorporate frailty, in which the proposal is to estimate the parameters of a NHPP process with unobserved heterogeneity that does not require parametric assumptions. In addition, examples are illustrated involving the PLP compared 184 Brazilian Journal of Biometrics to the Gamma frailty model and the model that does not consider the frailty term. On the other hand, according to Wienke (2010), a parametric specification of the frailty distribution is a matter of mathematical convenience. Therefore, assigning a mathematically tractable distribution to the frailty variable with adequate properties can result in suitable models to estimate unobserved heterogeneity. In this sense, due to its mathematical simplicity, the Gamma distribution has traditionally been used as a frailty distribution in works in the areas of survival and reliability analysis, as also discussed by Wienke (2010). In order to propose a distribution for the frailty term that is mathematically convenient as the Gamma distribution, which is used in most studies, Tomazella (2003) concluded that the inverse Gaussian distribution has properties equally simple to the Gamma distribution. Additionally, simulation studies were performed to infer on the parameters of interest and a Bayesian approach was proposed considering multiplicative and additive frailty models. In a recent work, to deal with life insurance data, Onchere et al. (2021) propose a more flexible distribution for the frailty term, the Non-central Gamma distribution. According to these authors, this distribution can represent time-varying frailties and fits well to real life assurance data. In addition, they considered the Weibull, log-logistic and log-normal baseline risk functions and their results indicate that the Non-central Gamma-Weibull frailty model is adequate for the data in question. Although the work of Onchere et al. deals with non-recurring events, the idea of a distribution capable of representing time-varying frailties is particularly attractive for models with recurring events. It is reasonable to think about the existence of unobservable factors that can influence the occurrence of recurring events differently in different periods of time. Furthermore, the Noncentral Gamma distribution is flexible and also presents mathematical facilities such as a closed form for the Laplace Transform. In view of the above, it is possible to realize the relevance of incorporating the frailty term with a different distribution in the modeling of repairable systems, seeking for more flexible assumptions that better fit different scenarios. In this way, this paper will explore the models for repairable systems, submitted to minimal repairs with the effect of frailty, where times will be modeled by the PLP process, while for frailty the Non-central Gamma distribution will be considered. So far the approach with this frailty is unknown in the literature in a perspective that considers recurrent events and this is a relevant contribution of our work. This paper is divided as follows. Section 2 reviews the Non-central Gamma distribution, models for repairable systems under minimal repair, and frailty modeling. In Section 3, the Non-central Gamma frailty model is introduced. In Section 4, some properties are verified through a simulation study. In Section 5, applications to real data are made. Finally, in Section 6, final considerations are raised. 2. Background In this section we present the theoretical background used as a basis for our proposed Noncentral Gamma (NCG) frailty model for repairable systems. Initially, in Section 2.1, we define the NCG distribution and outline some of its key features. In Section 2.2, we present the main ideas of modeling failure times for systems subject to minimal repairs. Concluding our review in Section 2.3, we provide a brief characterization of multiplicative frailty models. 2.1 Non-central Gamma distribution According to Knüsel & Bablok (1996) and de Oliveira & Ferreira (2013), the NCG distribution can be defined as a mixture of a discrete distribution and a continuous distribution, where the first distribution is a Poisson and the second is a central Gamma. Onchere et al. (2021), in turn, define 185 Brazilian Journal of Biometrics the NCG distribution as a convolution generated by the sum of Gamma distributions with Poisson weights. Summarizing the ideas and definitions presented by these authors, the NCG distribution can be defined in a general way as a mixture of a Gamma(a + m, b) with weights Poisson(ϕ). In this way, let Z be a random variable with a NCG distribution, the probability density function (pdf ) of Z is defined by z ∞ –ϕ m X e ϕ za+m–1 e– b f (z) = , z > 0, (1) m! Γ (a + m)ba+m m=0 where ϕ > 0 is called the non-centrality parameter and a ≥ 0 and b > 0 comes from the Gamma distribution. The NCG distribution is generally asymmetrical and its skewness is controlled by the noncentrality parameter. Figure 1 shows different behaviors of the density function of the NCG distribution for different parameter values. In this figure, parameter a has a fixed value of zero, since this will be an assumption adopted later in this work. 0.05 0.08 ϕ=5 b = 0.5 ϕ = 10 0.04 b=1 ϕ = 30 b = 1.5 0.06 ϕ = 50 b=2 ϕ = 100 b=4 f(z) f(z) 0.03 0.04 0.02 0.02 0.01 0.00 0.00 0 100 200 300 400 0 100 200 z (a) z (b) Figure 1. (a) NCG pdf for fixed b = 3 and variable ϕ values, and (b) NCG pdf for fixed ϕ = 50 and variable b values. Both cases with a = 0 fixed. A relevant characteristic of the NCG distribution is that it has a closed form for the Laplace Transform, which is given by   bϕs –a . (2) L(s) = (1 + bs) exp – 1 + bs As a direct consequence, we can easily obtain the expectation and variance of the NCG distribution, which are respectively given by E[Z] = b(a + ϕ) and Var[Z] = b2 (a + 2ϕ). The calculations for obtaining the Laplace Transform, the expectation and the variance of the NCG distribution are presented in Appendix 1. 2.2 Review of minimal repair models for repairable systems Minimal repair models are already widely known and discussed in the reliability literature. The tool of theses studies is the time elapsed until the occurrence of a failure, which can be represented by a random variable T. In the context of repairable systems, failure events are recurrent, which leads to obtaining a non-decreasing sequence of positive random variables 0 < T1 ≤ T2 ≤ . . . which represent the failure times. Let N(t) be a number of observed failures in the interval [0, t]. Under the MR assumption, {N(t), t ≥ 0} is considered a NHPP, that is, a counting process with independent increments, N(0) = 186 Brazilian Journal of Biometrics 0 and completely characterized by a non-constant intensity function λ(t) (Gilardoni & Colosimo, 2007). Furthermore, for any t, N(t) has a Poisson distribution with mean Λ(t), where Λ(t) = Rt E[N(t)] = 0 λ(u)du is the cumulative intensity function. Assuming that the intensity function λ(t) can be described by a parametric model with a vector of parameters µ, the failure time history of one or multiple identical systems under MR models can be used to estimate this parameter vector. In a classical inference approach, the likelihood method can be used to obtain the maximum likelihood estimators for the parameters (Tomazella, 2003). Assume that k identical repairable systems are under observation and let 0 < ti,1 < ti,2 < · · · < ti,ni be their failure times, where ti,j is the j-th failure time of the i-th system, with i = 1, . . . , k and j = 1, . . . , ni . Note that ni is the total number of observed failures of the i-th system and consider that for each system i there is still an observation truncation time ti∗ with ti,ni ≤ ti∗ . So, the general likelihood function for the MR model considering the k identical systems is given by    ni k  Y  Y  λ(ti,j ) e–Λ(ti∗ ) . (3) LMR (µ | ti,j ) =   i=1 j=1 Furthermore, the intensity function can be used to compute predictors of the reliability function. Let Ti,n = ti,n be the last observed failure time of the system i and t = Ti,n+1 –ti,n the time until the next failure. Considering the parametric MR model with vector of parameters µ, the general reliability at time t for the system i is given by R(ti,n+1 | ti,n ) = = P[Ti,n+1 – ti,n > t | ti,n ] = P[N(ti,n+1 ) – N(ti,n ) = 0] ! Z ti,n+1 λ(u)du = exp (–Λ(ti,n+1 ) + Λ(ti,n )) . exp – (4) ti,n Note that the function (4) expresses the probability that the system will run without failure for a period of time greater than t after the last failure time ti,n , given the failure history of the process. All the information about MR models that we present here is known in the reliability literature and, particularly, the estimation and prediction of reliability are general for whatever the parametric form of λ(t) is. Among the possible forms, we will present and use in this work the PLP. Introduced by Crow (1975), the PLP is a parametric model widely used in the literature for NHPP, in which its intensity and cumulative intensity functions are respectively given by λ(t | β, η) = β η  β–1 t η and Λ(t | β, η) =  β t , η (5) where η > 0 is the scale parameter and β > 0 is the shape parameter. The PLP is commonly assumed for MR models since its parameters have a physical interpretation in the following sense: the parameter η indicates the time during which exactly one failure is expected to occur, that is, E[N(η)] = 1, while the parameter β indicates whether the system is deteriorating (β > 1) or improving (β < 1) (Dias De Oliveira et al., 2013). To obtain the likelihood and reliability prediction functions considering the PLP parametric model, just replace λ(t) and Λ(t) in the equations (3) and (4) by the intensity and cumulative intensity functions given in (5). In this case, the model vector parameter will be µ = (β, η). This procedure and respective observations can be seen in Brito et al. (2022), Almeida et al. (2020), and Coque Junior (2016). Brazilian Journal of Biometrics 2.3 187 Frailty models As stated in the introduction of this work, frailty models incorporate in the intensity function a latent random variable Z that captures random effects. This variable is drawn from a distribution f (z; θ) and, in this study, it will be included in the intensity function in a multiplicative manner. In this sense, frailty models are considered as natural extensions of the Cox model. Disregarding the presence of covariates, a multiplicative frailty model for repairable systems is described by λ(t | z) = zλ0 (t), where λ0 (.) is the baseline intensity function and z is the frailty term. The choice of f (z) is an important point to consider and some authors have already discussed the impacts of adopting different distributions for the frailty term, e.g. Hougaard (1984) and Tomazella (2003)). Any continuous or discrete distribution with positive support, finite mean and unknown finite variance can be chosen as a frailty distribution. However, some conditions must be satisfied for the baseline intensity function to be identifiable, such as E[Z] = 1 and Var(Z) = θ (Elbers & Ridder, 1982). These restrictions lead to the following interpretations of the frailty term: if zi > 1, the failure intensity is increasing; if zi < 1, the failure intensity is decreasing; and if zi = 1, the model reduces to the baseline intensity function. 3. NCG frailty model for repairable systems In this section we present the proposed NCG frailty model for multiple repairable systems under MR. In addition, we present the likelihood function used to estimate the model parameters as well as the reliability prediction function. As mentioned earlier, to avoid identifiability problems in the model with frailty following a Z distribution, the model requires that its expectation is equal to 1 and its variance is constant. To guarantee this assumption when defining that the frailty Z has an NCG distribution and following the idea presented in Onchere et al. (2021), it is necessary to assume that the parameter a in the equation (1) is equal to zero. Furthermore, considering Var[Z] = θ the parameters b and ϕ of equation (1) will be reparametrized as b = θ/2 and ϕ = 2/θ, so that the expectation is such that E[Z] = 1. In this case, the pdf of the frailty Z in equation (1) can be rewritten as f (z) = ∞ e X m=0 – θ2  m 2 θ m! 2z zm–1 e– θ  m . Γ (m) θ2 (6) Consequently, the Laplace Transform given in equation (2) can be rewritten as L(s) = exp – s 1 + θ2 s ! . Obtaining a closed form for the Laplace Transform is a particularly interesting result for the frailty proposal with NCG distribution, especially for the calculation of the likelihood function. This is because the likelihood function must be obtained in a non-conditional way to frailty Z and the integrals used in these calculations are not trivial. However, if Z has a distribution that assumes a closed form for the Laplace Transform, this problem is easily solved. According to Wienke (2010), the unconditional reliability function is given by R(t) = L(Λ0 (t)) L′ (Λ (t)) and the intensity function unconditional the frailty term is given by λ(t) = –λ0 (t) L(Λ 0(t)) . 0 188 Brazilian Journal of Biometrics Therefore, the unconditional reliability function for the NCG frailty model is given by ! Λ0 (t) R(t) = exp – θΛ (t) 1 + 20 and the unconditional intensity function and the unconditional cumulative intensity function for the NCG frailty model are respectively given by λ(t) =  λ0 (t) 1+  θΛ0 (t) 2 2 and Λ(t) = Λ0 (t) 1+ θΛ0 (t) 2 , (7) where λ0 (t) and Λ0 (t) are the baseline intensity function and the baseline cumulative intensity function, respectively. The functions in (7) can be interpreted as functions of population intensity and cumulative population intensity that depend directly on the baseline functions and the variance of the assumed NCG frailty. The failure intensity function in (7) can assume unimodal forms and not just an increasing or decreasing form (as the PLP intensity function without frailty, for example), which is especially useful for modeling situations where the failures can be concentrated in one period and occur less frequently in other periods. Having defined the frailty model functions, we now want to use observed data to estimate the model parameters. Consider k repairable systems under observation for a period ti∗ , i = 1, . . . , k so that ti,1 , ti,2 , · · · , ti,k are the failure times of the i-th system and ti∗ is the truncation time of the observation of this system. Let us assume that these systems undergo minimal repairs after each failure, so on the classical approach the population parameters can be estimated through the maximum likelihood method. Let us assume that each failure process follows a PLP whose failure intensity and cumulative intensity functions are given by (5) and that there is an unobserved heterogeneity between the failure times modeled by the NCG distribution given in (6). Therefore, the vector of parameters to be estimated in the model is µ = (β, η, θ), where β and η are the parameters of the PLP model and θ is the variance of the NCG frailty. In this case, the likelihood function in (3) can be rewritten as LMR (β, η, θ | ti,j ) = k Y i=1 ( ni Y j=1 β η  "  ti,j β–1 η θ 1+ 2   #–2 ! ti,j β η exp ti∗ β ηβ + θ2 ti∗ β !) , (8) and consequently, the log-likelihood function is given by ni k X X   lMR (β, η, θ | ti,j ) = N log(β) – β log(η) + (β – 1) log(ti,j ) i=1 j=1 –2 ni k X X i=1 j=1 " θ log 1 + 2   # X k ti,j β ti∗ β + , η ηβ + θ ti∗ β i=1 (9) 2 P where N = ki=1 ni is the total number of observed failures across all systems. Given the complexity of the expression (9) and its partial derivatives, it is not possible to obtain the maximum likelihood estimators analytically. Therefore, numerical methods can be used to circumvent this problem and thus obtain the estimates. The confidence intervals of the estimated parameters can be constructed from the asymptotic theory based on the properties of the maximum likelihood estimators and the Normal distribution. Brazilian Journal of Biometrics 189 Finally, based on the reliability predictor function for the classical MR model presented in (4), we can define the reliability predictor function for our proposed NCG frailty model. Given the maximum likelihood estimates (β̂, η̂, θ̂) of the NCG frailty model parameters and let t = Ti,n+1 – ti,n be the time until the last observed failure time ti,n to the next failure time Ti,n+1 , the reliability function is given by R(ti,n+1 | ti,n ) = exp (–Λ(ti,n+1 ) + Λ(ti,n )) , (10) where the cumulative intensity function Λ(t) is the cumulative intensity function of the NCG frailty model given in (7). 4. Simulation study In this section, we conduct a simulation study to assess the efficiency and consistency performance of the maximum likelihood estimator (MLE) defined by maximizing the log-likelihood function in the model presented in (8). The main objective is to show that the asymptotic properties of the proposed maximum likelihood estimators are observed in samples of a simulated population of failure times with frailty effects induced by the NCG distribution. To evaluate the consistency and efficiency of the MLE of the NCG frailty model with PLP base intensity we used three metrics, based on similar studies carried out in the reliability literature: the mean relative error (MRE), the square root of the mean standard error (RMSE) and the coverage probability (CP) of the confidence interval (CI). These three metrics are respectively defined by M 1 X µ̂m MRE(µ̂) = , M µ m=1 where: v u M u1 X (µ̂m – µ)2 RMSE(µ̂) = t M and m=1 CP(µ̂) = M 1 X I[µ̂m ∈ (am , bm )], M m=1 • M is the number of samples to be simulated; • µ represents one of the parameters to be estimated (β, η or θ) and µ̂m is its estimate in the m-th sample; • am = µ̂m – 1.96 × SE(µ̂m ) is the lower bound of the 95% CI of the m-th estimate; • bm = µ̂m + 1.96 × SE(µ̂m ) is the upper bound of the 95% CI of the m-th estimate; • I is the indicator function and SE(µ̂m ) is the standard error of the m-th estimate. By this approach, it is expected that as we increase the sample sizes of observed failure times, the MRE will asymptotically approach 1, the RMSE will approach 0 and the CP will reach an adequate percentage of coverage of the nominal values of the parameters by the asymptotic CI (here we are adopting 95%). The increase in the sample will occur by increasing the number of observed systems. In this case, the number of failures observed in each system is random, however, it is reasonable to expect that the more systems observed, the more failure times are recorded. To generate the failure times for each system, the population unconditional cumulative intensity function of the NCG frailty model defined in (7) will be used. More specifically, let tj and tj+1 be two consecutive failure times such that tj+1 = tj + x. As described in Brito et al. (2022), the cumulative distribution function F(x) for the elapsed time X between two consecutive failures is F(x) = 1 – P[N(tj + x) – N(tj ) = 0]. Under the MR assumption, the counting process N(t) is an NHPP and as discussed in Section 2.2, the cumulative distribution function can be rewritten as F(x) = 1 – e–[Λ(tj +x)–Λ(tj )] , 190 Brazilian Journal of Biometrics where Λ(t) is the cumulative intensity function of the NCG frailty model, given in (7). Thus, to generate each process failure time, simply generate the time between two consecutive failures as a solution to the equation F(x) = u, where u is an observation of the Uniform(0, 1) distribution. This procedure will be repeated until the desired number of failures are generated or until the predefined truncation time is reached. Furthermore, considering k distinct s, this entire procedure will be repeated independently k times. The algorithm below describes the process of generating failure times for k systems, defining a truncation time t∗ of the observations for each system. 1. 2. 3. 4. 5. Set parameter values for β, η and θ; Set i the system index (i in {1, . . . , k}); Set a truncation time ti∗ ; Set ti,0 = 0 and j = 1; Draw ui,j from a Uniform(0, 1) distribution; 6. Solve the equation 1 – e–[Λ(ti,j–1 +x)–Λ(ti,j–1 )] = ui,j , where Λ(·) is given in (7); 7. Set ti,j = ti,j–1 + x; 8. If ti,j > ti∗ stop simulation and discard ti,j . Otherwise, set j = j + 1 and repeat steps 5-7. Different combination scenarios of model parameters were considered for the simulation study, as well as different sample sizes based on the number of systems. Due to the unimodal characteristic of the NCG frailty model’s failure intensity function, for each parameter scenario the truncation time was defined as the time t∗ such that λ(t∗ ) < 0.001, since the model loses influence in the generation of failure times when the intensity asymptotically approaches zero. Two fixed values were defined for the parameter β = (0.8, 1.2), to consider the two cases where the PLP baseline intensity function is increasing and decreasing. For each fixed β, three values were defined for the parameter θ = (0.05, 0.1, 0.2), that include different variances for the assumed NCG frailty. For each combination of parameters β and θ, a parameter η was defined to complete the scenario. Here, we will present the results of six different scenarios that represent the aforementioned model variations. For each defined scenario, 1,000 samples were generated following the previously defined algorithm and, for each sample, the β̂, η̂ and θ̂ estimates of the parameters β, η and θ, respectively, were obtained. From there, the previously defined MRE, RMSE and CP metrics were obtained and the results are summarized and shown in Figures 2 and 3. Additionally, all results are reported in Tables 5-10 in Appendix 2. From Figures 2 and 3, it is possible to see that, in general, the expected behaviors of the three metrics presented are verified in all determined scenarios. It is noticeable that increasing the sample size by increasing the number of observed systems results in CP increasingly close to the nominal value of 95% and RMSE increasingly close to zero for all parameters. It is clear and expected that specifically for the parameter η the RMSE value converges more slowly to zero, given the magnitude of the nominal values of this parameter; however, observing the MRE graphs, the expected convergence of this parameter to the nominal value is verified. In scenarios where there is a greater variance of frailty (θ = 0.2), it can be seen from the MRE graphs that the parameter θ can be slightly underestimated by the model; however, the interval estimates encompass the true nominal values as observed by the high percentage coverage on the CP graphs. We can conclude through the analysis of the simulation studies that the estimators for the proposed model perform well, especially for large samples. 191 Brazilian Journal of Biometrics θ = 0.05 and η = 20 θ = 0.1 and η = 15 θ = 0.2 and η = 10 1.00 0.95 CP 0.90 0.85 MRE Results 2 1 0 10 RMSE 5 0 10 20 30 40 50 10 20 30 40 50 10 20 30 40 50 40 50 Number of Systems beta eta theta Figure 2. Simulation results for NCG frailty model with fixed β = 0.8. θ = 0.05 and η = 20 θ = 0.1 and η = 30 θ = 0.2 and η = 40 1.00 0.95 CP 0.90 0.85 MRE Results 2 1 0 15 RMSE 10 5 0 10 20 30 40 50 10 20 30 40 50 10 20 Number of Systems beta eta theta Figure 3. Simulation results for NCG frailty model with fixed β = 1.2. 30 192 Brazilian Journal of Biometrics 5. Application to real data In this section we present two applications revisiting two known datasets in the reliability literature. The idea is to verify whether in these data it is possible to identify the presence of unobserved heterogeneity in the failure times of the systems and to verify whether the proposed NCG frailty model is adequate for such data. In both examples the parameters were obtained using the maximum likelihood method through the likelihood and log-likelihood functions defined in (8) and (9) respectively, and the confidence intervals were obtained using the asymptotic theory of the normal distribution, as mentioned in Section 3. In addition, we also used the well-known Gamma frailty model to compare the fits obtained with the proposed model in both datasets and we used the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC) as selection criteria between these two models (for more details on these criteria, see, e.g., Burnham & Anderson (2004)). 5.1 Sugarcane harvester data The first application deals with a set of failure times of 9 Brazilian sugarcane harvesters. This dataset was originally presented by Verssani (2018) and a subset of these harvesters was used in the frailty analysis under the assumption of minimal repair presented in D’Andrea et al. (2021). Specifically, the times of failure of the cutting blade system of these harvesters during the period of high sugarcane harvest are observed between the years 2014 and 2015, for a period of 200 days. Note that the truncation time of observations for all of them is the same equal to 200 days. Table 1 shows the observed failure times of each harvester. Note that the truncation time of observations for all of them is the same equal to 200∗ days. Table 1. Harvester cutting blade system failure times System 1 System 2 System 3 System 4 System 5 System 6 System 7 System 8 System 9 10 42 194 200∗ 51 68 85 110 120 146 157 167 83 95 100 116 125 90 103 117 125 136 99 107 121 129 23 24 39 41 72 136 138 165 166 200∗ 4 14 30 48 72 175 200∗ 157 171 173 13 32 44 66 74 93 131 153 157 163 185 200∗ 9 13 14 22 24 40 49 52 61 79 97 106 116 117 132 143 148 171 196 200∗ 62 73 88 107 118 8 9 25 31 40 124 154 158 178 200∗ 2 8 33 44 72 89 129 143 164 189 190 200∗ 60 82 83 101 102 128 87 91 92 116 117 118 182 200∗ 1 29 31 38 129 153 182 200∗ 35 42 59 70 123 151 158 166 181 Substituting the failure times ti,j observed in Table 1 into the likelihood function shown in (8), numerical and computational methods can be used to maximize this function and obtain the desired MLEs. This maximization procedure was performed using the software R (R Core Team, 2021), specifically using its optim function. Despite the algebraic complexity of the likelihood function, Brazilian Journal of Biometrics 193 making it impossible to obtain the MLE analytically, this is not a major problem from a computational point of view, so these estimates were easily obtained. As stated earlier, we also used the Gamma frailty model for repairable systems under minimal repair by D’Andrea et al. (2017) to estimate their parameters and compare with the results obtained by our NCG frailty model. The results of the parameter estimates for the both models are shown in Table 2. In addition, this table lists the maximum values of the log-likelihood function of each model, as well as the AIC and BIC values. Table 2. Estimation results for the frailty NCG and Gamma models applied to sugarcane harvesters data Model β̂ η̂ θ̂ l̂ AIC BIC Gamma 1.107 14.262 0.0378 -463.263 932.527 941.265 NCG 1.115 14.444 0.0349 -463.220 932.439 941.177 It is immediately possible to see the similarity between the estimates of β̂, η̂ and θ̂ for both models. In terms of comparing the models, the AIC and BIC values indicate that the NCG frailty model is slightly superior to the Gamma model for this data set under minimal repair, since they presented lower values. Although the difference in the information criteria for the two models is small, this result shows that our proposed model for NCG frailty is a good alternative to the Gamma model, and may even present slightly more parsimonious estimates. In addition to the point estimates presented in Table 2, the interval estimates for the parameters of the NCG frailty model were also obtained through asymptotic CIs. The 95% CIs for the MLE β̂, η̂ and θ̂ are, respectively, CIβ̂ = (0.81, 1.53), CIη̂ = (8.83, 23.64) and CIθ̂ = (0.006, 0.207). The point estimate of parameter β indicates that the systems are deteriorating over time once β̂ > 1, however this interpretation may not be true since the lower limit of CI is less than 1. Regarding the frailty, despite small, a value greater than zero was obtained for the variance of the frailty variable (θ̂ > 0), indicating the existence of unobserved heterogeneity between the systems and their failure times. In this sense, this result is in line with the other results obtained for this data set in the literature, indicating the existence of unobserved effects on the failure time of the harvesters, since the existence of frailty effects was captured. Finally, we use the equation (10) to estimate the reliability prediction function at the last observed failure time for each harvester. Figure 4 shows the reliability prediction function for the NCG frailty model related to each harvester of this study. It is possible to notice that all harvesters have almost zero reliability after 50 days of the last failure, given the history of each failure process. Furthermore, it is possible to notice that the median expected time to the next failure (i.e., the time such that the expected reliability is 0.5) is close to the estimated η̂ = 14.4 value, for all harvesters. 5.2 Literature systems data The second application deals with the failure times of the set of repairable systems presented in Mettas & Zhao (2005). In the original work, failure times of six systems with different observation truncation times are presented. However, for our application we will use the failure times of only three of these systems, since the other systems have a very low number of observed failures (one or two failures). Table 3 shows the observed failure times of each of the systems. The times marked with “ ∗ ” are the truncation times of the observation of each system. The same procedures as in the previous example were repeated here. Table 4 shows the MLE obtained for the Gamma and NCG frailty models, as well as their estimates of maximum log-likelihood and AIC and BIC values. Once again, it is noticeable that the point estimates are very close for the 194 Brazilian Journal of Biometrics Harvester 1 Harvester 2 Harvester 3 Harvester 4 Harvester 5 Harvester 6 Harvester 7 Harvester 8 Harvester 9 1.00 0.75 0.50 0.25 0.00 1.00 ^ R(t) 0.75 0.50 0.25 0.00 1.00 0.75 0.50 0.25 0.00 0 25 50 75 100 0 25 50 75 100 0 25 50 75 100 Days of operation after T = tn NCG frailty reliability Figure 4. Estimated reliability functions at last failure times tn , for each harvester in the data set, under the fitted NCG frailty model. Table 3. Systems failure times System 1 System 2 System 3 2227.08 2733.23 3524.21 5568.63 5946.30 6018.22 7202.72 8760∗ 3011.11 3121.46 3624.16 4328.32 772.95 1034.46 3758.30 5000∗ 900.99 1289.95 2689.88 3928.82 4704.24 5052.59 5473.17 6200∗ 5886.17 two models. Again, the AIC and BIC selection criteria indicate a slight superiority of our NCG frailty model over the gamma model. Table 4. Estimation results for the frailty NCG and Gamma models applied to literature systems data Model β̂ η̂ θ̂ l̂ AIC BIC Gamma 1.883 1518.044 0.187 -170.477 346.953 350.610 NCG 1.878 1583.923 0.131 -170.348 346.696 350.353 The 95% CIs for the MLE β̂, η̂ and θ̂ of the NCG frailty model are, respectively CIβ̂ = (1.04, 3.38), CIη̂ = (907.92, 2763.25) and CIθ̂ = (0.04, 0.41). In this case, we have more assertive conclusions about the estimates obtained. First, regarding the parameter β, both the point and interval estimates are greater than 1, indicating that the systems are deteriorating over time. Furthermore, the point and interval estimates of parameter θ are greater than zero, again capturing the existence of unobserved heterogeneity between systems and their failure times. 195 Brazilian Journal of Biometrics Once again, we use the equation (10) to estimate the reliability prediction function at the last observed failure time for each system in the studied dataset. The results of the estimated reliability prediction functions for the NCG frailty model are shown in Figure 5. The reliability probability for 2,000 days after the last failure is very close to zero for each system, given the entire history of its failure process. As in this case the estimate θ̂ did not approach zero, the interpretation for the parameter η̂ cannot be done as in the previous example, since the high variance of the frailty variable impacts the interpretation of this parameter in the PLP, as discussed previously. System 1 1.00 0.75 0.50 0.25 0.00 System 2 1.00 ^ R(t) 0.75 0.50 0.25 0.00 System 3 1.00 0.75 0.50 0.25 0.00 0 1000 2000 3000 Days of operation after T = tn NCG frailty reliability Figure 5. Estimated reliability functions at last failure times tn , for each system in the data set, under the fitted NCG frailty model. 6. Final remarks In this article, the main focus was to incorporate and estimate the frailty term in a model for repairable systems subject to unobserved heterogeneity, under the assumption of minimal repair. For the studied model, the likelihood functions and the reliability prediction functions were derived. The times were modeled by the PLP process, while the NCG distribution was considered as distribution for the frailty term. A simulation study considering different scenarios was carried out to verify the quality of the obtained estimators and their behavior. It was found that the results obtained are in line with expectations, since the estimates generated were close to the actual values of the parameters predefined in the study, returning expected values for the analyzed metrics (MRE, RMSE and CP). Two applications were made to real data. We verified that the model has better performance according to the considered evaluation criteria even if they were close to the alternative model that assumes the Gamma distribution for the frailty term. Furthermore, the parameter estimates are in agreement with the results found in the literature in applications for the same datasets. Furthermore, in each application it was possible to predict the reliability of each system, based on their failure process history and the obtained estimates of the model parameters. These results 196 Brazilian Journal of Biometrics show that the approach proposed in this work is effective and can be applied in different industry situations With this, the application in real data showed that the proposed approach is capable of providing useful information for decision making in complex systems, such as forecasting until the first failure. In summary, the results from the simulation and applications on real data corroborate the assertion that the proposed approach constitutes a valuable tool for the analysis of reliability in repairable systems with NCG distribution frailty. Therefore, as up to the present moment the approach with this frailty was unknown in the literature in a framework of repairable systems. The results obtained from theoretical discussions, applications to real data and the simulation study show that the model proposed in this work is relevant and presents a modeling alternative that clearly identifies the existence of unobservable factors, such as heterogeneity over the systems’ failure times. It is worth remembering that such factors directly interfere with the failure intensity function of the systems, which is closely linked to the useful life of the systems. Finally, this study presents an alternative propose when we are dealing with non-central distributions. We emphasize here the importance of one more option for frailty distribution among the models already existing in the literature. The studies developed here are not exhaustive and we emphasize the importance of investigating other approaches, such as other methods of estimating the parameters of interest. Acknowledgments Research developed using computational resources of the Center for Applied Mathematical Sciences to Industry (CeMEAI) financed by FAPESP. The authors also thank CAPES and Petrobras for financial support during the course of this project. Conflicts of Interest The authors declare no conflict of interest. References 1. Almeida, M. P., Paixão, R. S., Ramos, P. L., Tomazella, V., Louzada, F. & Ehlers, R. S. Bayesian non-parametric frailty model for dependent competing risks in a repairable systems framework. Reliability Engineering & System Safety 204, 107145 (2020). 2. Ascher, H. & Feingold, H. Repairable systems reliability: modeling, inference, misconceptions and their causes (Marcel Dekker Inc, New York, 1984). 3. Barlow, R. E. & Proschan, F. Mathematical theory of reliability (SIAM, 1996). 4. Brito, É. S., Tomazella, V. L. & Ferreira, P. H. Statistical modeling and reliability analysis of multiple repairable systems with dependent failure times under perfect repair. Reliability Engineering & System Safety 222, 108375 (2022). 5. Burnham, K. P. & Anderson, D. R. Multimodel inference: understanding AIC and BIC in model selection. Sociological methods & research 33, 261–304 (2004). 6. Coque Junior, M. A. Modelo de confiabilidade para sistemas reparáveis considerando diferentes condições de manutenção preventiva imperfeita. PhD thesis (Universidade de São Paulo, 2016). 7. Crow, L. H. Reliability analysis for complex, repairable systems tech. rep. (Army Materiel Systems Analysis Activityaberdeen Proving Ground Md, 1975). 8. De Oliveira, I. R. C. & Ferreira, D. F. Computing the Noncentral Gamma distribution, its inverse and the noncentrality parameter. Computational Statistics 28, 1663–1680 (2013). Brazilian Journal of Biometrics 197 9. De Toledo, M. L. G., Freitas, M. A., Colosimo, E. A. & Gilardoni, G. L. ARA and ARI imperfect repair models: Estimation, goodness-of-fit and reliability prediction. Reliability Engineering & System Safety 140, 107–115 (2015). 10. Dias De Oliveira, M., Colosimo, E. A. & Gilardoni, G. L. Power law selection model for repairable systems. Communications in Statistics-Theory and Methods 42, 570–578 (2013). 11. D’Andrea, A., Feitosa, C. C., Tomazella, V. & Vieira, A. M. C. Frailty modeling for repairable systems with Minimum Repair: An application to dump truck data of a Brazilian Mining Company. J Math Stat Sci 6, 179–198 (2017). 12. D’Andrea, A. M., Tomazella, V. L., Aljohani, H. M., Ramos, P. L., Almeida, M. P., Louzada, F., Verssani, B. A., Gazon, A. B. & Afify, A. Z. Objective bayesian analysis for multiple repairable systems. Plos one 16, e0258581 (2021). 13. Elbers, C. & Ridder, G. True and spurious duration dependence: The identifiability of the proportional hazard model. The Review of Economic Studies 49, 403–409 (1982). 14. Gilardoni, G. L. & Colosimo, E. A. Optimal maintenance time for repairable systems. Journal of quality Technology 39, 48–53 (2007). 15. Hougaard, P. Life table methods for heterogeneous populations: distributions describing the heterogeneity. Biometrika 71, 75–83 (1984). 16. Knüsel, L. & Bablok, B. Computation of the noncentral gamma distribution. SIAM Journal on Scientific Computing 17, 1224–1231 (1996). 17. Mettas, A. & Zhao, W. Modeling and analysis of repairable systems with general repair in Annual Reliability and Maintainability Symposium, 2005. Proceedings. (2005), 176–182. 18. Onchere, W, Weke, P, Jam, O & Carolyne, O. Non-Central Gamma Frailty with application to life term assurance data. Advances and Applications in Statistics 67, 237–253 (2021). 19. Proschan, F. Theoretical explanation of observed decreasing failure rate. Technometrics 5, 375– 383 (1963). 20. R Core Team. R: A Language and Environment for Statistical Computing R Foundation for Statistical Computing (Vienna, Austria, 2021). https://www.R-project.org/. 21. Slimacek, V. & Lindqvist, B. H. Nonhomogeneous Poisson process with nonparametric frailty. Reliability Engineering & System Safety 149, 14–23 (2016). 22. Tomazella, V. L. D. Modelagem de dados de eventos recorrentes via processo de Poisson com termo de fragilidade. PhD thesis (Universidade de São Paulo, 2003). 23. Vaupel, J. W., Manton, K. G. & Stallard, E. The impact of heterogeneity in individual frailty on the dynamics of mortality. Demography 16, 439–454 (1979). 24. Verssani, B. A. W. Modelo de regressão para sistemas reparáveis: um estudo da confiabilidade de colhedoras de cana-de-açúcar PhD thesis (Universidade de São Paulo, 2018). 25. Wang, H. A survey of maintenance policies of deteriorating systems. European journal of operational research 139, 469–489 (2002). 26. Wienke, A. Frailty models in survival analysis (CRC press, 2010). 198 Brazilian Journal of Biometrics Appendix 1. Calculations for the NCG distribution In this appendix we present some calculations and results concerning the NCG distribution. As the NCG distribution is not widely known and commonly used by professionals and students of Statistics, recording these procedures and results can be useful for future queries about this distribution. The idea here is to obtain the closed form of the Laplace Transform for this distribution and from it calculate the expectation and variance. Let Z be a random variable with NCG distribution, whose pdf is given by the equation (1). It is known that the Laplace Transform of a function f (z) is defined by Z ∞ L(s) = e–sz f (z)dz. (11) 0 Substituting the pdf fZ (z) defined in (1) in the equation (11) and following with algebraic manipulations, the Laplace Transform is rewritten as L(s) z ∞ –ϕ m X e ϕ za+m–1 e– b dz m! Γ (a + m)ba+m 0 m=0 Z ∞ –1 ∞ –ϕ i X (s + b–1 )a+m za+i–1 e–z(s+b ) 1 e ϕ 1 dz = m! ba+m (s + b–1 )a+m 0 Γ (a + m) m=0  a+m ∞ X ϕm 1 –ϕ = e m! 1 + bs = Z ∞ e–sz m=0 = = = by ∞ e –ϕ X 1  ϕ  m m! 1 + bs (1 + bs)a m=0  – ϕ ϕ  e exp 1 + bs (1 + bs)a   1 bsϕ exp – . 1 + bs (1 + bs)a As a consequence, the first and second derivatives of the Laplace transform are respectively given   ϕ i bϕs h b exp – a + L (s) = – 1 + bs 1 + bs (1 + bs)a+1 ′ and L′′ (s) =   b2 bϕs h ϕ  ϕ  ϕ i exp – a + a + 1 + + . 1 + bs 1 + bs 1 + bs 1 + bs (1 + bs)a+2 Finally, the expectation and variance of the random variable Z with NCG distribution can be easily obtained by derivatives of the Laplace Transform (see, e.g. Wienke, 2010), as follows: E[Z] = –L′ (0) = b(a + ϕ), Var[Z] = L′′ (0) – (L′ (0))2 = b2 (a + 2ϕ). Brazilian Journal of Biometrics 199 Appendix 2. Simulation study results tables In this Appendix we present the detailed results of the simulation study discussed in Section 4. Tables 5-7 below refer to the results summarized in Figure 2, while Tables 8-10 refer to the results summarized in Figure 3. Table 5. Simulation study results for scenario β = 0.8, η = 20 and θ = 0.05 RMSE MRE CP k β η θ β η θ β η θ 5 0.209 11.677 0.052 1.094 1.191 1.201 0.912 0.912 0.830 10 0.149 7.681 0.041 1.043 1.090 1.033 0.898 0.928 0.884 15 0.116 5.824 0.036 1.029 1.056 1.010 0.917 0.944 0.886 20 0.100 4.900 0.032 1.018 1.036 0.987 0.904 0.948 0.902 25 0.088 4.228 0.029 1.009 1.023 0.965 0.920 0.946 0.907 30 0.081 3.820 0.027 1.009 1.028 0.946 0.920 0.960 0.919 35 0.075 3.787 0.027 1.005 1.032 0.925 0.921 0.948 0.933 40 0.070 3.397 0.024 1.006 1.022 0.961 0.932 0.943 0.925 45 0.068 3.221 0.024 1.009 1.022 0.974 0.933 0.948 0.918 50 0.066 3.112 0.023 1.006 1.017 0.961 0.932 0.945 0.926 Table 6. Simulation study results for scenario β = 0.8, η = 15 and θ = 0.1 RMSE MRE CP k β η θ β η θ β η θ 5 0.212 9.577 0.067 1.064 1.210 0.919 0.922 0.922 0.913 10 0.148 5.814 0.054 1.028 1.095 0.883 0.908 0.944 0.920 15 0.127 4.738 0.048 1.014 1.078 0.865 0.913 0.945 0.940 20 0.100 3.861 0.041 1.001 1.058 0.856 0.938 0.951 0.955 25 0.098 3.598 0.040 0.999 1.046 0.860 0.922 0.935 0.960 30 0.084 3.105 0.037 0.995 1.043 0.853 0.950 0.946 0.965 35 0.080 2.896 0.035 0.991 1.027 0.851 0.948 0.945 0.976 40 0.078 2.710 0.033 0.998 1.033 0.872 0.937 0.941 0.963 45 0.073 2.557 0.031 0.994 1.029 0.866 0.946 0.945 0.976 50 0.068 2.488 0.030 0.997 1.038 0.872 0.955 0.951 0.973 200 Brazilian Journal of Biometrics Table 7. Simulation study results for scenario β = 0.8, η = 10 and θ = 0.2 RMSE MRE CP k β η θ β η θ β η θ 5 0.234 7.485 0.122 1.048 1.283 0.663 0.901 0.917 0.979 10 0.166 4.846 0.110 1.000 1.176 0.623 0.919 0.919 0.991 15 0.133 3.956 0.105 0.977 1.137 0.609 0.917 0.926 0.997 20 0.118 3.417 0.103 0.965 1.126 0.595 0.936 0.920 1.000 25 0.106 3.075 0.102 0.962 1.136 0.581 0.938 0.917 0.998 30 0.101 2.716 0.096 0.960 1.108 0.601 0.951 0.915 0.999 35 0.090 2.371 0.094 0.953 1.084 0.597 0.949 0.925 1.000 40 0.092 2.319 0.098 0.950 1.099 0.577 0.928 0.930 1.000 45 0.088 2.228 0.094 0.952 1.099 0.585 0.944 0.922 1.000 50 0.084 2.037 0.094 0.949 1.091 0.585 0.936 0.927 0.999 Table 8. Simulation study results for scenario β = 1.2, η = 20 and θ = 0.05 RMSE MRE CP k β η θ β η θ β η θ 5 0.190 6.352 0.015 1.015 1.042 0.971 0.939 0.946 0.933 10 0.126 4.266 0.011 1.007 1.024 0.987 0.955 0.952 0.950 15 0.108 3.756 0.008 1.007 1.021 0.994 0.945 0.937 0.952 20 0.093 3.080 0.007 1.005 1.013 0.997 0.950 0.955 0.948 25 0.081 2.833 0.006 1.003 1.007 0.999 0.958 0.945 0.954 30 0.076 2.499 0.006 1.004 1.011 0.999 0.951 0.944 0.944 35 0.070 2.331 0.006 1.003 1.009 0.997 0.946 0.937 0.938 40 0.066 2.170 0.005 1.006 1.015 1.000 0.944 0.948 0.962 45 0.060 2.016 0.005 1.002 1.006 0.997 0.953 0.949 0.957 50 0.057 1.906 0.004 1.002 1.004 1.002 0.948 0.950 0.958 Table 9. Simulation study results for scenario β = 1.2, η = 30 and θ = 0.1 RMSE MRE CP k β η θ β η θ β η θ 5 0.251 10.924 0.035 1.040 1.075 0.967 0.930 0.926 0.943 10 0.168 7.385 0.024 1.011 1.029 0.973 0.950 0.946 0.954 15 0.142 5.903 0.020 1.007 1.024 0.976 0.938 0.941 0.954 20 0.125 5.391 0.017 1.007 1.024 0.977 0.936 0.919 0.960 25 0.105 4.360 0.015 1.005 1.014 0.982 0.958 0.954 0.950 30 0.099 4.127 0.014 1.005 1.017 0.982 0.954 0.949 0.958 35 0.091 3.794 0.013 1.001 1.013 0.973 0.948 0.940 0.964 40 0.085 3.577 0.012 0.996 1.003 0.976 0.950 0.949 0.973 45 0.078 3.268 0.011 1.001 1.012 0.982 0.952 0.951 0.970 50 0.075 3.171 0.011 1.002 1.010 0.983 0.952 0.939 0.964 Brazilian Journal of Biometrics Table 10. Simulation study results for scenario β = 1.2, η = 40 and θ = 0.2 RMSE MRE CP k β η θ β η θ β η θ 5 0.391 16.962 0.114 1.081 1.130 0.833 0.911 0.908 0.945 10 0.263 11.313 0.095 1.023 1.075 0.794 0.915 0.929 0.975 15 0.210 8.918 0.080 1.011 1.062 0.802 0.936 0.941 0.979 20 0.185 7.822 0.073 1.004 1.038 0.810 0.940 0.931 0.983 25 0.162 7.070 0.069 0.991 1.043 0.791 0.948 0.940 0.987 30 0.145 6.255 0.068 0.977 1.036 0.771 0.963 0.944 0.994 35 0.140 5.703 0.067 0.979 1.032 0.772 0.955 0.941 0.993 40 0.125 5.269 0.059 0.983 1.027 0.797 0.960 0.947 0.995 45 0.129 5.194 0.061 0.983 1.033 0.787 0.943 0.934 0.996 50 0.118 4.980 0.058 0.984 1.030 0.788 0.943 0.933 0.993 201