# Discrimination between accelerated life models via Approximate Bayesian Computation

Mohamed Rabhi, Anis Ben Abdessalem, Laurent Saintis, Bruno Castanier, Rodrigue Sohoin

► To cite this version:

Mohamed Rabhi, Anis Ben Abdessalem, Laurent Saintis, Bruno Castanier, Rodrigue Sohoin. Discrimination between accelerated life models via Approximate Bayesian Computation. Quality and Reliability Engineering International, 2023, 10.1002/qre.3283. hal-03993798

HAL Id: hal-03993798

https://univ-angers.hal.science/hal-03993798

Submitted on 17 Feb 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Discrimination between accelerated life models via Approximate Bayesian Computation

Mohamed Rabhi<fnref n="1" />,<fnref n="2" /> |

<figure><img src="image_2.png" /></figure>

Anis Ben Abdessalem<fnref n="1" /> |

Laurent Saintis<fnref n="1" /> |

Bruno Castanier<fnref n="1" /> |

Rodrigue Sohoin<fnref n="2" />

## Correspondence

Mohamed Rabhi, Univ Angers, LARIS, SFR MATHSTIC, Angers, France

Email:

mohamed.rebhi@etud.univ-angers.fr

## Abstract

Accelerated life testing (ALT) is widely used in high-reliability product estimation to get relevant information about an item's performance and its failure mechanisms. To analyse the observed ALT data, reliability practitioners need to select a suitable accelerated life model based on the nature of the stress and the physics involved. A statistical model consists of (i) a lifetime distribution that represents the scatter in product life and (ii) a relationship between life and stress. In practice, several accelerated life models could be used for the same failure mode and the choice of the best model is far from trivial. For this reason, an efficient selection procedure to discriminate between a set of competing accelerated life models is of great importance for practitioners. In this paper, accelerated life model selection is approached by using the Approximate Bayesian Computation (ABC) method and a likelihood-based approach for comparison purposes. To demonstrate the efficiency of the ABC method in calibrating and selecting accelerated life model, an extensive Monte Carlo simulation study is carried out using different distances to measure the discrepancy between the empirical and simulated times of failure data. Then, the ABC algorithm is applied to real accelerated fatigue life data in order to select the most likely model among five plausible models. It has been demonstrated that the ABC method outperforms the likelihood-based approach in terms of reliability predictions mainly at lower percentiles particularly useful in reliability engineering and risk assessment applications. Moreover, it has shown that ABC could mitigate the effects of model misspecification through an appropriate choice of the distance function.

## KEYWORDS

accelerated life models, accelerated life tests, Approximate Bayesian Computation, likelihood-based approach, model selection, reliability, percentiles

<footnote n="1">¹Univ Angers, LARIS, SFR MATHSTIC, Angers, France</footnote>

<footnote n="2">²Liebherr-Aerospace Toulouse SAS, Toulouse, France</footnote>
# 1 INTRODUCTION

Designing components with high reliability and longer lifetime presents a challenging task for industry because the necessary testing time under normal use conditions is often excessive. Such circumstances lead to consider accelerated life testing (ALT) as the most suitable technique to overcome this issue through the extrapolation of the failure data collected from severe stress levels to the nominal stress level.¹ It is commonly assumed that the failure modes and failure mechanisms remain the same for all the applied stress levels to be able to extrapolate. Roughly speaking, the components are subjected to higher stress conditions than those for which they were designed to operate forcing them to fail more rapidly. It is commonly assumed that for all stress levels, the failure times are governed by the same parametric lifetime distribution. The acceleration life models are formulated by expressing the parameters of the reliability model with relationships, called acceleration models, which express their dependencies on the stress levels. The statistical physics-based models explain the relationship between the applied stress levels and the system failure time by using the parameters of the failure physics in conjunction with the statistical parameters to obtain realistic models.¹⁻³ A detailed description of ALT can be found in Meeker and Escobar.⁴

In ALT, researchers are widely divided between two assumptions: (1) the stress independence of the distribution shape parameter is commonly accepted because the resulting models are typically easier to use and the data set suggests that from a practical standpoint, such a constraint is sometimes valid since extrapolation is made assuming that the failure mechanism remains the same for all the accelerated stress levels. Therefore, the tested units fail in the same manner across different loads. Carlsson et al.⁵ made this assumption for the service life assessment of solar thermal components, and Gaertner et al.⁶ for the life prediction at normal use conditions of cathode ray tubes. However, the second assumption (2) is that the shape parameter does depend on stress, that is, in the case of more than one parameter of the lifetime distribution being dependent on the stress, each one can be related to the stress by an acceleration model. In the literature, some authors suggested checking the validity of the stress independence assumption by a graphical procedure. Whitman et al.⁷ showed that the slopes of the probabilities lines for the different stress levels are nearly parallel. Also, Barlow et al.⁸ showed that ALT data of Kevlar/Epoxy strands follow Weibull distribution with a stress-dependent shape parameter. In the same perspective, Giuseppe et al.⁹ studied the estimation of fatigue reliability of structural components via a two-parameter Birnbaum-Saunders (BS) model with stress-dependent parameters (both scale and shape parameters).

Selecting a suitable acceleration life model is still quite challenging for practitioners. In practice, several competing models could potentially fit the data distribution reasonably well, regarding the physics involved. Birnbaum and Saunders¹⁰ showed that most two-parameter lifetime distributions fit well the fatigue data and produce similar fitting qualities in the central tendency region of the reliability curve. Nevertheless, large discrepancies could be seen in lower and upper percentiles when comparing these distributions. This is a challenging issue for industry interested in low percentiles predictions. Box and Draper¹¹ stated that ‘All models are wrong but some are useful’, this means that there is no true model for all the situations but there is only a model that performs better than the other competing models for a specific goal. Hence, the challenging task consists in selecting the ‘most useful’ model and estimating its parameters from ALT observations. This is of a paramount importance to make good reliability predictions at normal use conditions, especially for lower lifetime percentiles that are important for early failure predictions, costly product failure claims avoiding and warranty analysis.¹² Thus, models should be compared in order to select the one which presents the highest performance. A comprehensive review of reliability model selection was presented by Settani et al.¹³

Several classical methods have been proposed in the literature to deal with model selection issues. The largely used ones are presented hereafter. The classical techniques commonly used are the weighted least squares method,¹⁴⁻¹⁶ the graphical analysis method,¹⁷,¹⁸ the likelihood ratio test (see Kalbfleisch et al.¹⁹) and the Information Criteria such as the Akaike Information Criterion (AIC)²⁰ and the Bayesian Information Criterion (BIC).²¹ Ghaly et al.²² used likelihood ratio test to compare between exponentiated inverted Weibull and inverted Weibull distributions using real life data set. These approaches require the estimation of the model parameters that maximize the likelihood function and an additional term to penalise overly complex models. The major limitation of these approaches is that they undertake separately the competing models regardless the uncertainties related to the parameters estimates which could lead to an inaccurate selection of the ‘best’ model as stated by Gupta et al.²³

Contrarily to the classical approaches, Santana et al.²⁴ proposed a novel Bayesian formulation based on the residuals evaluation to examine the benefits of the uncertainty evaluation in the model selection procedure. Other Bayesian approaches based on goodness of fit indices evaluation are used to discriminate between the competing models. The idea is to select the model that presents the best point estimate of the indices to be evaluated.²⁵ Ling et al.²⁶ proposed
an optimization method for parameter estimation and model selection based on the minimization of the Kolmogorov-Smirnov distance between the observed probability of failure and the expected probability of failure associated to the competing distributions.

Further methods based on Bayes theory have been proposed in the literature for model selection. The Reversible-Jump Markov-Chain Monte Carlo (RJ-MCMC) algorithm is one of the largely used Bayesian methods in model selection issues.²⁷ Several studies have addressed the strengths and weakness of RJ-MCMC and the other trans-dimensional setups of the algorithm.²⁸,²⁹ Indeed, their implementations require the construction of across-model proposals between the state spaces. Hence, the change of models dimensions introduces additional difficulty since the definition of the between-model transition matrices may no longer be intuitive. This could result in slowing the state space exploration. In some circumstances, the proposed between-model jumps tend to be rejected and the chain becomes trapped in one model, even though that model has low posterior probability.³⁰

Bayes factors (BFs)³¹ have been considered for a long time as the standard tools for performing Bayesian model comparison. However, they provide only a relative comparison of competing models, not the absolute values of their posterior probabilities. In fact, from a decision-making point of view, it can be perceived as the posterior probability of the null hypothesis. The main drawback of BFs is their sensitivity to the choice of priors. In fact, the approach proposed by Jeffreys³¹ is not available for comparing models with improper priors since these are defined only up to a multiplicative constant, leaving the BF undetermined.³² This issue was addressed by the intrinsic Bayes factors (IBF) proposed by Berger et al.³³ The IBF uses a part of the data as a training sample to convert the non-informative priors to proper posterior distributions. Then, the IBF are computed with the remainder of the data and using the posterior distribution estimated from the training data as a prior distribution. Nevertheless, the main challenges of IBF are that they manage only pairwise comparisons to discriminate numerous models.

So far, the common limitation of the classical Bayesian approaches is that they require the definition of a likelihood function to measure the level of agreement between the observed and simulated data particularly for high-dimensional models and low-sized samples. In some circumstances, the likelihood function could not be available in a closed form and may encounters numerical problems such the non-convergence or the convergence to the wrong roots raised in Refs. 34–37 Another aspect concerns the effects of model misspecification on the desired estimates of some life characteristics. Balakrishnan et al.³⁸ showed that the effect of model misspecification between Gamma and Weibull distributions should not be ignored for one-shot device testing data. The authors evaluated bias, root mean square errors (MSEs), coverage probabilities and average widths of confidence intervals for the mean lifetime as well as the reliability at a given mission time and suggested implementing a specification test to improve the estimation accuracy. Furthermore, Khakifirooz et al.³⁹ considered generalised Gamma, Lognormal and Weibull distributions for studying model misspecification effects under accelerated lifetime testing with censored data. Castilla et al.⁴⁰ studied the effect of model misspecification between Gamma, Weibull, Lognormal and BS lifetime distributions on the design of optimal constant-stress accelerated life-test for one-shot device.

In this work, and in order to overcome the above-mentioned issues, a likelihood-free method called also Approximate Bayesian Computation (ABC) which has been widely used in the literature (see, Beaumont et al.⁴¹) to deal with model selection and parameter estimation is adopted. Contrary to the classical Bayesian inference, which requires to write down a likelihood function, the ABC algorithms relax the need for an explicit likelihood function. In a nutshell, it consists in approaching the posterior distributions by generating samples from prior distributions of the model parameters, evaluating the level of agreement between the ALT data and a simulated data through an appropriate distance metric and propagate good particles until the convergence condition is satisfied. Nevertheless, a well-known problem of the ABC algorithm in its basic form is the low acceptance rate of samples. It decreases dramatically along the iterations which increases exponentially the computational time. This has been reported by Bonassi et al.⁴² To overcome this shortcoming and alleviate the computational burden, an ABC algorithm based on an elegant ellipsoidal Nested Sampling technique (ABC-NS) and a re-weighting scheme introduced in Ref. 43 is used. It promises drastic speedups compared to other variants⁴³ and provides a good alternative for parameter estimation as shown in Ref. 44 Additionally, it has been demonstrated that the ABC-NS can mitigate the effects of model misspecification through an appropriate choice of the distance function which can be tailored according to the user's interest.

The paper is organized as follows. In Section 2, the principle of the Bayesian inference and ABC is set out with a brief description of the algorithm implementation. Section 3 introduces the selected accelerated life models used in this study. In Section 4, a simulation study is carried out to demonstrate the performance of the proposed likelihood-free Bayesian procedure compared to a likelihood-based approach. In Section 5, the procedure is applied on the basis of real fatigue life data. Finally, Section 6 provides some concluding remarks with an emphasis on future research challenges.
## 2 | BAYESIAN INFERENCE

### 2.1 | Bayes formalism

Bayesian inference has proven very efficient to deal with model selection and parameter estimation in many disciplines. It consists of assessing the relative plausibility of a set of competing models $M = \{M_k, k = 1, ..., \kappa\}$ as well as the parameters $\theta$ within each model (for simpler presentation a subscript $k$ on $\theta$ is omitted) using a combination of one's prior knowledge and a set of available data, $D$. Both levels of inference (i.e. parameter estimation and model selection) can be addressed via the sequential application of Bayes' theorem:

$$ \mathbb{P}(\theta|D, M_k) = \frac{\mathbb{P}(D|\theta, M_k)\mathbb{P}(\theta|M_k)}{\mathbb{P}(D|M_k)} \quad (1) $$

$$ \mathbb{P}(M_k|D) = \frac{\mathbb{P}(D|M_k)\mathbb{P}(M_k)}{\mathbb{P}(D)} \quad (2) $$

Evaluation of Equation (1) requires the definition of the prior, $\mathbb{P}(\theta|M_k)$, and the likelihood, $\mathbb{P}(D|\theta, M_k)$. The prior is a subjective probability distribution which describes our prior belief for the parameter values. The likelihood encodes the information contained in the data according to model, $M_k$, with parameters, $\theta$. The denominator of Equation (1) is a normalising constant which ensures that $\mathbb{P}(\theta|D, M_k)$ integrates to unity. Successful evaluation of Equation (1) gives the posterior parameter distribution, which describes the probability of parameter vector, $\theta$, given the data, $D$, and the chosen model structure, $M_k$.

With regard to Equation (2), $\mathbb{P}(M_k)$ is the prior probability of $M_k$ that expresses the modeller's judgement on the initial relative plausibility of $M_k \in M$, before observing data $D$. Prior model probability $\mathbb{P}(M_k)$ can be specified depending on the existing prior knowledge about the credibility of model $M_k$, or it can be given a uniform probability, $\mathbb{P}(M_k) = 1/\kappa$, if not referring to any information. $\mathbb{P}(D)$ is a normalising constant and $\mathbb{P}(D|M_k)$ is equal to the evidence term on the denominator of Equation (2). $\mathbb{P}(M_k|D)$ is a distribution describing the relative probability of different competing model structures conditional on the data, $D$.

### 2.2 | Likelihood-free method

To deploy ABC methods, we need to be able to simulate data from the model, and require a distance function, $\rho(D, D^*)$, where $D^*$ is the simulated data generated through the model. This distance function provides a measure of discrepancy between the observed and the simulated data. After calculating $\rho(D, D^*)$, we accept the particle $\theta$ where $\rho(D, D^*) \le \varepsilon$, where $\varepsilon$ is the tolerance threshold. This leads to a modification of original Bayes theorem (see, Equation 1, the model is dropped for simplicity):

$$ \mathbb{P}(\theta|\rho(D, D^*) \le \varepsilon) = \frac{\mathbb{P}(\rho(D, D^*) \le \varepsilon|\theta)\mathbb{P}(\theta)}{\mathbb{P}(\rho(D, D^*) \le \varepsilon)} \quad (3) $$

For the prior probability distribution, $\mathbb{P}(\theta)$, we will assume a uniform prior over the possible parameter set, $\theta$. By being able to produce simulations from the model, we are able to perform inference for the parameters of interest, subject to data $D$. The tolerance threshold $\varepsilon$ determines the level of approximation as $\varepsilon \to 0$, $\mathbb{P}(\theta|\rho(D, D^*) \le \varepsilon) \to \mathbb{P}(\theta|D)$.

As already mentioned, the ABC method can address both levels of inference. Conversely to the classical Bayesian inference, the ABC algorithm brings a straightforward solution to approximate $\mathbb{P}(M_k|D)$ (see Refs. 45, 46). Briefly, a model $M_k$ is sampled from the prior and the model is stored only if the difference between simulated data and the observed data is less than a prefixed threshold. The posterior probability $\mathbb{P}(M_k|D)$ is then approximated by the acceptance frequency of $M_k$. When no prior information is available for the preference of the models, one assumes that the prior on the models follows a uniform distribution. As a result, comparing $\mathbb{P}(M_k|D)$ is the same as comparing $\mathbb{P}(D|M_k)$, and the model with the largest value of $\mathbb{P}(D|M_k)$ should be selected. In the literature, $\mathbb{P}(D|M_k)$ is known as the Bayesian evidence or the marginal likelihood under model $M_k$, and it has also been widely used for model selection (e.g. see Refs. 47, 48). The most elementary implementation of ABC is ABC rejection sampling (ABC-RS),$^{49,50}$ (see Algorithm 1). This algorithm requires a huge computational time to converge because it samples blindly over the model and parameter spaces. To overcome this issue, many variants have been proposed in the literature.$^{51,52}$ In the
<table><caption>ALGORITHM 1 ABC-RS for model choice.</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>1:</td><td>for i = 1 to N do</td></tr><tr><td>2:</td><td>repeat</td></tr><tr><td>3:</td><td>Draw M* from the prior P(M) (where M is a vector that indexes the models)</td></tr><tr><td>4:</td><td>Sample θ* from the prior p(θ|M*)</td></tr><tr><td>5:</td><td>Simulate a candidate data set D* ~ F(·|θ*, M*)</td></tr><tr><td>6:</td><td>until ρ(D*, D) ≤ ε<sub>T</sub></td></tr><tr><td>7:</td><td>Store the pair (M<sup>(i)</sup> = M*, θ<sup>(i)</sup> = θ*)</td></tr><tr><td>8:</td><td>end for</td></tr><tr><td colspan="2">Ensure: Model posterior probabilities, marginal distributions</td></tr></tbody></table>

present study, the ABC algorithm coupled with an elegant nested sampling technique (ABC-NS) introduced in Ref. 43 is employed.

## 2.3 | Steps of the ABC-NS algorithm

The iterative process is given in Algorithm 2 for completeness (see Ref. 43 for more details). The algorithm starts by selecting a model from the candidate models supposed to be equally probable. A proposal particle is then sampled from the prior over the model parameters, simulating the data using the probability model and keep the pair (model, particle) if the constraint $ρ(D, D*) ≤ ε_1$ is satisfied (here, D for observed data and D* for simulated data, ρ(·) is the function measuring the discrepancy). In this study, the initial tolerance threshold is chosen to target an acceptance rate roughly equal to ~ 50%. The simulations are repeated until one gets N particles distributed over the candidate models (N is equal to 1000 for all the simulations). Then, weights are calculated for each particle and normalised according to each model using the kernel shown in step 10. The first iteration of the ABC-NS algorithm is similar to the ABC-RS algorithm but with a large tolerance threshold value to speed-up the algorithm.

The subsequent tolerance threshold is defined based on the discrepancy values ranked in descending order (highest on top, see, step 12) as the $(α_0N)$th value where $α_0$ is equal to 0.3. The dropped particles represent 30% from the total number of particles. After that, we normalise the weights of the remaining particles (see, step 14) and a weight of zero is assigned to the dropped particles. From the remaining particles, a number of $(β_0N)$ are randomly selected and propagated to the next population ($β_0$ is equal to 0.6). In this way, the algorithm can visit any part of the parameter space and cannot be trapped in a local minima to ensure a good exploration of the model and parameter spaces. The remaining particles associated to each model are then enclosed in an ellipsoid in which the mass centre $μ$ and the covariance matrix $C$ are estimated based only on the remaining particles; one denotes this ellipsoid by $E_t^{(t)} = (μ_t^{(t)}, C_t^{(t)})$. The aim of the elliptical bound is to restrict the parameter space around the most interesting part of the parameter space which improves the acceptance ratio and efficiency through the iterations. The volume of the generated ellipsoid could be enlarged by a factor $f_0$ as mentioned in step 17 to ensure that the particles on the borders will be inside. It is worth noting that the ellipsoidal sampling was firstly proposed in Ref. 53 to improve the efficiency of the nested sampling algorithm which has been widely used for Bayesian inference. Finally, the population is replenished by re-sampling particles inside the enlarged ellipsoid for each model (see step 18) following the scheme in Ref. 54 In the subsequent steps, the threshold is updated adaptively in the same way and samples selection are subjected to more stringent threshold. The priors on the models are also updated when one of the competing model is eliminated. Through the populations and as $ε → 0$, a larger number of particles are selected for the most likely model(s) and the samples for the parameters better reflect the real posterior distribution. Several stopping criteria could be used to stop the algorithm. In this work, the algorithm is stopped when the difference between two successive tolerance threshold values falls below a pre-specified value.

The hyper-parameters applied in this study are those suggested in Ref. 43 through a sensitivity analysis. To sum up, the ABC-NS hyper-parameters are selected in such a way as to guarantee the best trade-off between the computational time and the statistical efficiencies by ensuring relatively high acceptance rates over iterations (≥50%) and precise estimates.
ALGORITHM 2 ABC-NS for model selection and parameter estimation.

**Require:** data: *D*, initial tolerance: *ε*₁, distance function: *ρ*(·), prior distributions: *P*(*Θ*), tuning parameters: (*α*₀, *β*₀, *f*₀), *N*: number of particles, *T*: number of populations

1: **At iteration t = 1**

2: **for** *i* = 1 to *N* **do**

3: **repeat**

4: Draw *M*\* from the prior *P*(*M*)

5: Sample *Θ*\* from the prior *P*(*Θ*|*M*\*)

6: Simulate a candidate data set *D*\* ~ *F*(·|*Θ*\*, *M*\*)

7: **until** *ρ*(*D*, *D*\*) ≤ *ε*₁

8: Store *M*<sub>*t*</sub><sup>(*i*)</sup> = *M*\*, *Θ*<sub>*t*</sub><sup>(*i*)</sup> = *Θ*\*, *e*<sub>*t*</sub><sup>(*i*)</sup> = *ρ*(*D*, *D*\*)

9: **end for**

10: Evaluate weights for all particles: *ω*<sub>*i*</sub> = (1 - (e<sub>*i*</sub>/ε<sub>1</sub>)<sup>2</sup>), *i* = 1, ..., *N*

11: Sort *e*<sub>*t*</sub><sup>(*i*)</sup> in descending order and store them.

12: Compute the next tolerance threshold *ε*₂ = *e*<sub>*t*</sub><sup>(*α*₀*N)</sup>

13: Drop particles with *ρ*(*D*, *D*\*) ≥ *ε*₂, {ω<sub>*k*</sub><sup>*α*₀*N*</sup>}, *k* = 0 ▷ A weight of zero is assigned to the dropped particles

14: Normalise the weights such that Σ<sup>(1-α₀)N</sup>*l*<sub>1</sub> ω<sub>l</sub> = 1

15: Select *A*<sub>*t*</sub> = *β*₀*N* particles from the remaining ones based on the weights

16: For each model, define the ellipsoid by its mass centre and covariance matrix *E*<sub>*t*</sub><sup>(*k*)</sup> = {*μ*<sub>*t*</sub><sup>(*k*)</sup>, *C*<sub>*t*</sub><sup>(*k*)</sup>}, *k* = 1, ..., *κ*

17: Enlarge the ellipsoids by *f*₀ ▷ The same notation for the updated ellipsoids is kept

18: **for** *t* = 2 to *T* **do**

19: **for** *j* = 1 to (1 - *β*₀)*N* **do**

20: **repeat**

21: Draw *M*\* from the prior *p*(*M*)▷ updated through the iterations

22: Sample *Θ*\* inside the ellipsoid *E*<sub>*t*-1</sub><sup>*</sup>

23: Simulate a candidate data set *D*\* ~ *F*(·|*Θ*\*, *M*\*),

24: **until** *ρ*(*D*, *D*\*) ≤ *ε*<sub>*t*</sub>

25: Store *M*<sub>*t*</sub><sup>(*j*)</sup> = *M*\* and add *Θ*<sub>*t*</sub><sup>(*j*)</sup> = *Θ*\* to the population of particles of {*Θ*(*M*\*)}, *e*<sub>*t*</sub><sup>(*j*)</sup> = *ρ*(*D*, *D*\*)

26: **end for**

27: Store the new particles in *S*<sub>*t*</sub>

28: Obtain the new particle set for each model, *N*<sub>*new*</sub> = [*A*<sub>*t*-1</sub>; *S*<sub>*t*</sub>] with their correspondent distance value *e*<sub>*t*</sub><sup>*</sup>

29: Sort *e*<sub>*t*</sub><sup>*</sup> and define *ε*<sub>*t*+1</sub> = *e*<sub>*t*</sub><sup>(*α*₀*N)</sup>

30: Compute new weights for all particles as in step (10) and normalise them

31: Define the new set of selected particles *A*<sub>*t*</sub> as in step (15)

32: Update the hyper-parameters of the ellipsoid using *A*<sub>*t*</sub>, *E*<sub>*t*</sub> = {*μ*<sub>*t*</sub><sup>(*k*)</sup>, *C*<sub>*t*</sub><sup>(*k*)</sup>}, *k* = 1, ..., *κ* ▷ The same factor is used for all the simulations

33: **end for**

**Ensure:** Model posterior probabilities, Marginal distributions

In fact, the stopping criterion threshold value, ξ, is selected by the user according to the required precision on the posterior estimates. Ben Abdessalem et al.⁴³ showed that, for a fixed ξ, the ellipsoid enlargement factor, *f*₀, impacts slightly the acceptance rate of the ABC-NS approach over the populations. Furthermore, the effect of *f*₀ on the posterior estimates is negligible. As for the initial tolerance threshold value *ε*₁, it is calibrated, particularly for each case study, by running some few preliminary numerical simulations. Indeed, it should be set large in order to explore all the parameters' space. Finally, the best trade-off between *α*₀ and *β*₀ is, respectively, 0.3 and 0.6 as these values ensure a stable convergence and a high acceptance rate.

At the last population, the algorithm produces a Markov chain on (*M*, *Θ*) for which the marginal distribution is *p*(*M*|*D*). Let us suppose that (*M*<sub>*k*</sub><sup>(*i*)</sup>, *Θ*<sub>*k*</sub><sup>(*i*)</sup>) with (*i* = 1, ..., *N*, *k* = 1, ..., *κ*) are obtained at the last population, the posterior model
probability can then be estimated by

$$ \hat{P}(M_k|D) = \frac{1}{N} \sum_{i=1}^{N} \mathbb{I}(M^{(i)} = M_k), \quad k = 1, \dots, \kappa \qquad (4) $$

## 2.4 | Selection of distance functions

To implement the ABC-NS algorithm, one needs to select an effective distance function. Various types of distances have been proposed in the literature to measure the discrepancy between the empirical data (times to failure) and the data simulated by the reliability model. In this section, the computational forms of the distances used in this study are given.⁵⁵

* **The Cramer-von Mises distances** The Cramer-von Mises (CM) distance measures the distance between the cumulative distribution function (CDF) of the derived distribution against the dataset's cumulative histogram. It is defined by

$$ W_n^2 = n \int_{-\infty}^{+\infty} [F(t) - S_n(t)]^2 \Psi(t) dF(t) \qquad (5) $$

where $F(t)$ is the predicted cumulative probability from the theoretical CDF and $S_n(t)$ is the empirical cumulative distribution function (ECDF) of the observed data. $\Psi(t)$ is a weight function such that $\Psi(t) = 1$. For given ordered sample, $t_1, t_2, \dots, t_n$, from a statistical distribution, the distance can be written as follows:

$$ \rho_{\text{CM}} = \frac{1}{12n} + \sum_{i=1}^{n} \left[ z_i - \frac{2i-1}{2n} \right]^2 \qquad (6) $$

where $z_i = F(t_i)$ is the predicted value for the ith observation. In the following, the same alphabets denote the same meaning.

* **The Anderson-Darling distance** The Anderson-Darling (AD) distance first introduced by Anderson and Darling⁵⁶ to place more weight at the tails of the distribution.⁵⁷ It is defined by the following equation:

$$ A_n^2 = n \int_{-\infty}^{+\infty} [F(t) - S_n(t)]^2 \Psi(t) dF(t) \qquad (7) $$

where $\Psi(t) = [F(t)(1-F(t))]^{-1}$ is a non-negative weight function. Equation (7) can be written for a finite data sample as

$$ \rho_{\text{AD}} = -n - \frac{1}{n} \sum_{i=1}^{n} (2i-1)[\ln(z_i) + \ln(1-z_{n+1-i})] \qquad (8) $$

Two modified distances of the AD distance called ADL and ADR which give more weight to the left and right tail, respectively, have been tested. They are computed by the following formulas:

$$ \rho_{\text{ADL}} = -\frac{3n}{2} + 2 \sum_{i=1}^{n} z_i - \frac{1}{n} \sum_{i=1}^{n} (2i-1) \ln z_i \qquad (9) $$

$$ \rho_{\text{ADR}} = \frac{n}{2} - 2 \sum_{i=1}^{n} z_i - \frac{1}{n} \sum_{i=1}^{n} (2i-1) \ln(1-z_{n+1-i}) \qquad (10) $$

Similarly, the $\rho_{AD2R}$ and $\rho_{AD2L}$ metrics, which assign larger weights to the tails of the distribution, are considered. They are computed by the following formulas:

$$ \rho_{\text{AD2L}} = 2 \sum_{i=1}^{n} \log z_i + \frac{1}{n} \sum_{i=1}^{n} \frac{2i-1}{z_i} \qquad (11) $$
$$ \rho_{\text{AD2R}} = 2 \sum_{i=1}^{n} \log(1 - z_i) + \frac{1}{n} \sum_{i=1}^{n} \frac{2i - 1}{1 - z_{n+1-i}} \quad (12) $$

# 3 | FORMULATION OF THE ACCELERATION LIFE MODELS

## 3.1 | Choice of the parametric lifetime distributions

The accelerated life models consists in choosing a parametric distribution and an acceleration model. In practice, the ALT data (times-to-failure) may be represented by several lifetime distributions. In this study, three parametric distributions named Weibull, BS and Lognormal distributions will be used. It should be noted that a candidate distribution is chosen for one of the following reasons: (i) there is a physical statistical argument that theoretically matches a failure mechanism to a life distribution model, (ii) a particular model has previously been used successfully and confirmed by experience for the same or a similar failure mechanism, (iii) a convenient model provides a good empirical fit to all the failure data. In this study, our choices are motivated by the failure mechanism which is induced by fatigue. The specifications of the parametric distributions used in this study with some relevant life characteristics are defined below.

* **The two-parameter Weibull distribution** It is a popular distribution used for modelling the lifetime data because of its flexibility within multiple types of failure mechanisms and its easiness of computing probabilities without the need for numerical integration (Cao et al.⁵⁸). The probability density function and the reliability function of the two-parameter Weibull distribution are given respectively by

$$ f(t, \eta, \beta) = \left(\frac{\beta}{\eta}\right) \left(\frac{t}{\eta}\right)^{\beta-1} \exp \left\{ -\left(\frac{t}{\eta}\right)^{\beta} \right\}, \quad (\eta, \beta > 0); \quad t \ge 0 \quad (13) $$

$$ R(t, \eta, \beta) = \exp \left\{ - \left( \frac{t}{\eta} \right)^{\beta} \right\} \quad (14) $$

where $\eta$ and $\beta$ are the scale and shape parameters, respectively. $t$ is the lifetime of the component. The percentiles, denoted $t_p$, is the time by which a specified fraction $p$ of the population fails is given by

$$ t_p = F^{-1}(p) = \eta(-\ln(1-p))^{1/\beta}, \quad 0 < p < 1 \quad (15) $$

where $F^{-1}(\cdot)$ is the inverse CDF. The mean time to failure (MTTF) is the expected life $E(T)$, it is defined by

$$ \text{MTTF} = E(T) = \int_{0}^{\infty} R(t)dt = \eta\Gamma(1 + \frac{1}{\beta}) \quad (16) $$

where $\Gamma(\cdot)$ denotes the gamma function.

* **The two-parameter Lognormal distribution** The Lognormal distribution is widely used for modelling failure times of materials and structures induced by fatigue. When $Y$ has a Lognormal distribution with parameters $\mu$ and $\sigma^2$, then $X = \ln(Y)$ has a normal distribution with mean equal to $\mu$ and a variance $\sigma^2$.⁴ The probability density function of Lognormal distribution is

$$ f(t) = \frac{1}{\sqrt{2\pi}\sigma t} \exp\left(-\frac{[\ln t - \mu]^2}{2\sigma^2}\right) = \frac{1}{\sigma t} \phi\left[\frac{\ln t - \mu}{\sigma}\right], \quad (17) $$

where $\phi$ is the standard normal pdf. The percentiles is given by the following equation:

$$ t_p = e^{\sigma\Phi^{-1}(p)+\mu}, \quad 0 < p < 1 \quad (18) $$
where $\Phi^{-1}(\cdot)$ is the inverse standard normal CDF. The MTTF is defined by

$$MTTF = e^{\left(\frac{\sigma^2}{2} + \mu\right)}. \qquad (19)$$

* **The two-parameter Birnbaum-Saunders distribution** The BS model is based on a physical argument of cumulative damage that produces fatigue in materials, which is exerted by cyclical stress. The failure follows from the development and growth of a dominant crack in the material. The probability density function and the reliability function of the BS model are given, respectively, by

$$f(t, \alpha, \nu) = \frac{t^{-3/2}(t + \nu)}{2\sqrt{2\pi}\alpha\sqrt{\nu}} \exp\left\{-\frac{1}{2\alpha^2}\left(\frac{t}{\nu} + \frac{\nu}{t} - 2\right)\right\}, \quad (\nu, \alpha > 0); \quad t \ge 0 \qquad (20)$$

$$R(t, \alpha, \nu) = 1 - \Phi\left\{\frac{1}{\alpha}\left(\sqrt{\frac{t}{\nu}} - \sqrt{\frac{\nu}{t}}\right)\right\} \qquad (21)$$

where $\alpha$ and $\nu$ are the shape and scale parameters, respectively. $t$ is the lifetime of the component and $\Phi(\cdot)$ is the CDF of the normal distribution. By definition, the BS random variable $T \sim \text{BS}(\alpha, \beta)$, is linked to the standard normal random variable $Z$ (where $Z \sim N(0,1)$) by the following relationship:

$$Z = \frac{1}{\alpha} \left( \sqrt{\frac{T}{\nu}} - \sqrt{\frac{\nu}{T}} \right) \qquad (22)$$

The percentiles and the MTTF are given, respectively, by the following equations:

$$t_p = \frac{\nu}{4} \left\{ \alpha \Phi^{-1}(p) + \sqrt{4 + (\alpha \Phi^{-1}(p))^2} \right\}^2, \quad 0 < p < 1 \qquad (23)$$

$$MTTF = \nu \left( 1 + \frac{\alpha^2}{2} \right) \qquad (24)$$

## 3.2 | Choice of the acceleration model

This is a key element of ALT, as it expresses the relationship between the life characteristics and the applied stress. Combined with the life distribution, it provides a single model, known as the accelerated life model, which the ALT data analysis can be based on. It has been assumed in this study, that the parameters of the selected distributions are expressed as stress function through the inverse power law (IPL) relationship which is commonly used for non-thermal accelerated stresses. The choice is also motivated by the failure mechanism and the nature of the applied stress. The IPL is given by Equation (25):

$$L(S) = e^{\gamma_0} x_i^{-\gamma_1} \qquad (25)$$

where $L(S)$ is the life characteristic, $x_i = \frac{S}{S_0}$ is the ratio of severe stress by operating stress, $\gamma_0$ and $\gamma_1$ are the IPL parameters of the acceleration model estimated from the ALT data.

From Equation (25), a linear extrapolation can be obtained using a logarithmic transformation on the accelerating variables as follows:

$$\ln(L(S)) = \gamma_0 - \gamma_1 \ln(x_i) \qquad (26)$$

Following this, the dependence of the scale parameter $\eta$ of the Weibull distribution at the ith stress level can be expressed as follows:

$$\eta(S_i) = \exp\{\gamma_0 - \gamma_1 \ln(x_i)\} \qquad (27)$$
<table><caption>TABLE 1 Accelerated life models.</caption><thead><tr><th>Parametric distribution</th><th>Acceleration model: stress dependency</th><th>Model size</th><th>Model</th></tr></thead><tbody><tr><td rowspan="2">Weibull</td><td>- Only scale parameter</td><td>3</td><td>M<sub>1</sub></td></tr><tr><td>- Both scale and shape parameters</td><td>4</td><td>M<sub>2</sub></td></tr><tr><td rowspan="2">Birnbaum-Saunders</td><td>- Only scale parameter</td><td>3</td><td>M<sub>3</sub></td></tr><tr><td>-Both scale and shape parameters</td><td>4</td><td>M<sub>4</sub></td></tr><tr><td>Lognormal</td><td>- Only the mean parameter</td><td>3</td><td>M<sub>5</sub></td></tr></tbody></table>

In the same way, the shape parameter could be linked to the applied stress by the same formula:

$$ \beta(S_i) = \exp\{\gamma'_0 - \gamma'_1 \ln(x_i)\} \quad (28) $$

The same IPL has been used for the BS and Lognormal distributions. In this study, for illustrative purpose, only five competing accelerated life models are considered. The models are obtained by coupling the IPL to the two-parameter Weibull, the BS and the Lognormal distributions as listed in Table 1.

To summarize, let denote by $\theta_1 = (\beta, \gamma_0, \gamma_1)$ the parameter vector of $M_1$, $\theta_2 = (\gamma'_0, \gamma'_1, \gamma_0, \gamma_1)$ the parameter vector of $M_2$, $\theta_3 = (\alpha, \gamma_0, \gamma_1)$ the parameter vector of $M_3$, $\theta_4 = (\gamma'_0, \gamma'_1, \gamma_0, \gamma_1)$ the parameter vector of $M_4$ and $\theta_5 = (\sigma, \gamma_0, \gamma_1)$ the parameter vector of $M_5$.

## 3.3 | Formulation of the log-likelihood function of the accelerated life models

The likelihood-based approach uses the log-likelihood functions for the different competing models. One assumes that the components are tested only under three accelerated stress levels. Hence, $t_{i,j}$ ($i = 1, 2, 3$ and $j = 1... n_i$) stands for the set of the ALT data in the $i$th group and $n_i$ is the number of failures observed at $x_i$. Let $\rho_{LH}$ denotes the log-likelihood function of $M_1$. It is given by the following equation⁵⁹:

$$ \rho_{LH}(\theta_1) = N[\ln \beta - \beta \ln \gamma_0] + -\gamma_1 \beta \sum_{i=1}^{3} \ln x_i + (\beta + 1) \sum_{i=1}^{3} \sum_{j=1}^{n_i} \ln t_{ij} + \sum_{i=1}^{3} \sum_{j=1}^{n_i} \left( \frac{t_{ij}}{\gamma_0 x_i^{-\gamma_1}} \right) \quad (29) $$

Similarly, the log-likelihood functions of models $M_2$, $M_3$, $M_4$ are given, respectively, by the following equations⁹:

$$ \rho_{LH}(\theta_2) = N \left[ \ln \gamma'_0 + \sum_{i=1}^{3} \left( \gamma'_0 \gamma_1 x_i^{-\gamma'_0} - \gamma'_1 \right) \right] \ln x_i + \sum_{i=1}^{3} \sum_{j=1}^{n_i} \left( \gamma'_0 x_i^{-\gamma'_1} - 1 \right) \ln t_{ij} - \sum_{i=1}^{3} \sum_{j=1}^{n_i} \left[ \frac{t_{ij}}{\gamma_0 x_i^{-\gamma_1}} \right]^{\gamma'_0 x_i^{-\gamma'_1}} \quad (30) $$

$$ \rho_{LH}(\theta_3) = N \left[ \ln\left(\frac{1}{2\sqrt{2\pi}}\right) - \ln \alpha \right] - \frac{n}{\alpha^2} - \frac{1}{2\alpha^2} \left[ \sum_{i=1}^{3} \sum_{j=1}^{n_i} \left( \frac{t_{ij}}{\gamma_0 x_i^{-\gamma_1}} + \frac{\gamma_0 x_i^{-\gamma_1}}{t_{ij}} \right) \right] + \sum_{i=1}^{3} \sum_{j=1}^{n_i} \ln(t_{ij} + \gamma_0 x_i^{-\gamma_1}) + \frac{\gamma_1}{2} \sum_{i=1}^{3} n_i \ln x_i \quad (31) $$

$$ \rho_{LH}(\theta_4) = N \left[ \ln \left( \frac{1}{2\sqrt{2\pi}} \right) - \left( \log \gamma'_0 + \frac{1}{2} \ln \gamma_0 \right) \right] - \frac{1}{2\gamma'_0^2} \left[ \sum_{i=1}^{3} \sum_{j=1}^{n_i} \left( \frac{t_{ij}}{\gamma_0 x_i^{-(2\gamma'_0+\gamma_1)}} + \frac{\gamma_0 x_i^{2\gamma'_0-\gamma_1}}{t_{ij}} \right) \right] + \sum_{i=1}^{3} \sum_{j=1}^{n_i} \ln(t_{ij} + \gamma_0 x_i^{-\gamma_1}) + (\gamma'_1 + \frac{\gamma_1}{2}) \sum_{i=1}^{3} n_i \ln x_i + \frac{1}{\gamma'_0^2} \sum_{i=1}^{3} n_i \ln x_i^2 - \frac{3}{2} \sum_{i=1}^{3} \ln t_{ij} \quad (32) $$

$$ \rho_{LH}(\theta_5) = -\frac{N}{2} \ln(2\pi) - N \ln(\sigma) - \sum_{i=1}^{3} \sum_{j=1}^{n_i} \ln t_{ij} - \sum_{i=1}^{3} \sum_{j=1}^{n_i} \left[ \frac{(\ln t_{ij} - e^{\gamma_0 x_i^{-\gamma_1}})^2}{2\sigma^2} \right] \quad (33) $$
<table><caption>TABLE 2 Parameter bounds for the considered accelerated life model.</caption><thead><tr><th>Parameter</th><th>Type of distribution</th><th>Lower bound</th><th>Upper bound</th></tr></thead><tbody><tr><td>β</td><td>Uniform</td><td>0.1 × β̃</td><td>10 × β̃</td></tr><tr><td>γ<sub>0</sub></td><td>Uniform</td><td>0.1 × γ̃<sub>0</sub></td><td>10 × γ̃<sub>0</sub></td></tr><tr><td>γ<sub>1</sub></td><td>Uniform</td><td>0.1 × γ̃<sub>1</sub></td><td>10 × γ̃<sub>1</sub></td></tr></tbody></table>

## 4 | SIMULATION STUDY

In this section, a simulation study is conducted to investigate the performance of the ABC method to deal with model selection and parameter estimation of accelerated life models. First, a comparison between the different distances in terms of the bias and the MSE under different sample sizes is given. Then, a qualitative comparison in terms of reliability predictions of some relevant life characteristics is presented. Finally, a comparison of the PCS when very close competing models are considered under different sample sizes is given.

Let us consider a three-level constant-stress ALT without censoring. The accelerated data are collected under three stress levels: $S_1 = 320$ MPa, $S_2 = 360$ MPa and $S_3 = 400$ MPa. The operating stress level $S_0$ is supposed to be 200 MPa. It is assumed that the lifetime $T$ of an item at use condition follows a Weibull distribution with parameters: $\beta = 1.5$, $\gamma_0 = 10$ and $\gamma_1 = 5$. The IPL is assumed between the scale parameter and the stress level only for illustrative purpose. In this manner, the lifetime data for each stress level can be obtained by combining the CDF of the parametric Weibull distribution and the IPL following Equations (34) and (35):

$$F(T) = y = 1 - \exp \left\{ - \left( \frac{T}{\eta(S)} \right)^{\beta} \right\} \quad (34)$$

$$T = T(S) = \eta(S)[-\ln(1 - y)]^{1/\beta} \quad (35)$$

where $y$ is a uniform random number in the interval $[0, 1]$.

Sample sizes of $n = 5, 10, 20, 30, 50, 100$ and $200$ have been considered to investigate the effects of the times to failure data. In this study, the model parameters were given uniform priors spanning a range of one order of magnitude above and below crude estimates obtained by least square method¹⁸ as illustrated in Table 2 for $M_1$.

The ABC-NS algorithm is used here to deal with parameter estimation. It is implemented with the following hyper-parameters:

* the initial threshold $\epsilon_1$ is set to 1000.

* the number of particles per population is set to 1000.

* the proportion of dropped particles and surviving particles are set respectively, to 30 and 60%.

* The ellipsoid enlargement factor $f_0$, is set to 1.1.

* The stopping criterion threshold value is equal to $10^{-6}$.

The quality of the estimates from the different distances is studied by considering the bias, the MSE and the deficiency (DEF) given by

$$\text{Bias}(\hat{\theta}) = \bar{\theta} - \theta \quad (36)$$

such that

$$\bar{\theta} = \frac{1}{1000} \sum_{i=1}^{1000} \hat{\theta}_i \quad (37)$$

$$\text{MSE}(\hat{\theta}) = \frac{1}{1000} \sum_{i=1}^{1000} (\hat{\theta}_i - \theta)^2 \quad (38)$$

The DEF is a natural measure of joint efficiencies of the estimators $\hat{\beta}$, $\hat{\gamma}_0$ and $\hat{\gamma}_1$. DEF is formulated as follows:

$$\text{DEF}(\theta) = \text{MSE}(\hat{\beta}) + \text{MSE}(\hat{\gamma}_0) + \text{MSE}(\hat{\gamma}_1) \quad (39)$$
By replicating the experiment 1000 times, we can obtain the estimated MSEs, bias and the deficiencies ($\widehat{DEF}$) for the considered distances and sample sizes. From Table 3, the following observations can be drawn:

* As expected, the MSE of all the estimators become smaller as the sample sizes increase. The numerical results indicate that the ML estimator behave slightly better than the ABC in terms of the MSE.

* The likelihood-free method provides overall good estimates mainly when the $\rho_{AD}$ and $\rho_{ADR}$ are used showing an ABC estimator to be competitive with the ML estimator.

* By analysing the bias which is the average tendency of the simulated data to over-or under-estimate the real parameters (positive values indicate a model under-estimation and negative values indicate an over-estimation), it is clear that the $\rho_{LH}$ metric provides the best estimates, seconded by the $\rho_{ADR}$ and the $\rho_{AD}$ distances.

* From the same table, one can see that the latter distances provide an estimates close to the ML estimates. This mean that these distances could be used as a robust alternative estimator when for instance the likelihood function is difficult to formulate or challenging to be minimized which is a commonly encountered in practice mainly when one deals with small sample sizes or when it is challenging to derive a likelihood function.

## 4.1 | Estimation of the reliability lifetime characteristics

Put the focus now on the estimation of some important life characteristics widely used in reliability engineering applications. Some of them are: the MTTF, the percentile $B_{100p}$, ($0 \le p \le 1$) by which time 100$p\%$ of the product population will fail is used as the warranty life. Generally, the percentile life for airplane parts is $p = 0.1$ or $p = 0.01$ and is $p = 0.05$ or $p = 0.10$ for machine parts.⁶⁰ In Table 4, these life characteristics using the true parameters with the assumed ALT model are given. The objective here is to compare the predictive capability of the ABC method using the considered distance functions over the likelihood-based method. Furthermore, one investigates the effect of the sample size on the predicted lifetimes.

The results summarized in Figures 1 and 2 display bias and $log(MSE)$ of the average values of the considered percentiles obtained from the different distances and the true values for different sample sizes. The following observations can be drawn from both figures:

* It can be seen that the differences on the percentiles estimated by the different approaches and distances can be significant.

* One can see that the $\rho_{AD}$ distance outperforms the other distances to get precise estimates of the lower extreme percentiles ($B_1, B_5, B_{10}$) mainly for small sample sizes.

* Another important point should be highlighted is that, regardless the sample size, $\rho_{CM}$ and $\rho_{LH}$ distances always overestimate the lower percentiles ($B_1, B_5$ and $B_{10}$) by several orders of magnitude.

* For $\rho_{LH}$ distance, the estimates present the highest MSE for extremely small sample sizes ($n = 5$). However, the MSE using $\rho_{LH}$ decreases for large sample sizes and performs better than the other distance functions.

* Only the $\rho_{AD2L}$ and $\rho_{AD2R}$ yield conservative estimates of the percentiles as they ensure negative bias associated to lower percentiles. In the same way, one can see that the other distances provide always non-conservative predictions as the bias are greater than 0.

* It is clear that the MSE of each low percentiles for the considered distances tend to zero for large sample sizes (for instance, for $n = 5$ the $log(MSE)$ of $B_1$ estimated for $\rho_{AD2L}$ distance is 25% greater than that estimated for $n = 200$ using the same distance function).

* For large sample sizes, the MSE deviations between the considered distances start to widen contrary to small sample sizes. Indeed, $\rho_{AD2L}$ yields the highest values (as could be seen for $n = 20, n = 30, n = 50, n = 100, n = 200$).

* One can see that, for a given sample size, the MSE magnitudes increase for larger percentiles regardless the used distance function.

* To conclude, it is recommended to use the ABC method coupled with the $\rho_{AD}$ distance to precisely estimate the lower lifetime percentiles as it presents the lesser bias and the smaller MSE mainly for small sample sizes. Furthermore, $\rho_{LH}$ distance presents the less precise predictions of the low percentiles regardless the sample size. However, this overestimation is less pronounced for large sample sizes.
<table><caption>TABLE 3 Bias, MSE and DEF values using different sample sizes.</caption><thead><tr><th rowspan="2">Sample size</th><th rowspan="2">Parameter/ DEF</th><th colspan="2">ρ<sub>LH</sub></th><th colspan="2">ρ<sub>CM</sub></th><th colspan="2">ρ<sub>AD</sub></th><th colspan="2">ρ<sub>ADL</sub></th><th colspan="2">ρ<sub>ADR</sub></th><th colspan="2">ρ<sub>AD2L</sub></th><th colspan="2">ρ<sub>AD2R</sub></th></tr><tr><th>Bias</th><th>MSE</th><th>Bias</th><th>MSE</th><th>Bias</th><th>MSE</th><th>Bias</th><th>MSE</th><th>Bias</th><th>MSE</th><th>Bias</th><th>MSE</th><th>Bias</th><th>MSE</th><th>Bias</th></tr></thead><tbody><tr><td rowspan="5">5</td><td>β</td><td>0.2570</td><td>0.02507</td><td>0.1249</td><td>0.4901</td><td>-0.1802</td><td>0.1540</td><td>0.1004</td><td>0.2538</td><td>-0.0645</td><td>0.2070</td><td>-0.1219</td><td>0.4694</td><td>-0.0184</td><td>0.3250</td></tr><tr><td>γ<sub>0</sub></td><td>-0.0819</td><td>1.3295</td><td>-0.0382</td><td>1.6784</td><td>-0.0145</td><td>1.5025</td><td>-0.0194</td><td>1.9836</td><td>-0.0331</td><td>1.3976</td><td>0.0355</td><td>2.9310</td><td>-0.0504</td><td>1.6164</td></tr><tr><td>γ<sub>1</sub></td><td>-0.0878</td><td>3.8053</td><td>-0.0441</td><td>4.8572</td><td>-0.0561</td><td>4.3259</td><td>-0.0328</td><td>5.6693</td><td>-0.0914</td><td>4.0394</td><td>0.0377</td><td>8.5061</td><td>-0.1389</td><td>4.6503</td></tr><tr><td>D<sub>EF</sub></td><td>5.3855</td><td></td><td>7.0257</td><td></td><td>5.9824</td><td></td><td>7.9067</td><td></td><td>5.6439</td><td></td><td>11.9064</td><td></td><td>6.5917</td><td></td></tr><tr><td>β</td><td>0.1154</td><td>0.0692</td><td>0.0267</td><td>0.0930</td><td>-0.1045</td><td>0.0589</td><td>-0.0763</td><td>0.0785</td><td>-0.0540</td><td>0.0703</td><td>-0.1832</td><td>0.1658</td><td>-0.0712</td><td>0.1090</td></tr><tr><td rowspan="5">10</td><td>β</td><td>-0.0023</td><td>0.6406</td><td>0.0034</td><td>0.8345</td><td>0.0167</td><td>0.7539</td><td>0.0153</td><td>1.0015</td><td>0.0131</td><td>0.6836</td><td>0.0867</td><td>1.5569</td><td>-0.0180</td><td>0.7051</td></tr><tr><td>γ<sub>1</sub></td><td>0.0185</td><td>1.8343</td><td>0.0105</td><td>2.3838</td><td>0.0108</td><td>2.1552</td><td>0.0203</td><td>2.8739</td><td>0.0024</td><td>1.9467</td><td>0.0702</td><td>4.5236</td><td>-0.0567</td><td>2.0191</td></tr><tr><td>D<sub>EF</sub></td><td>2.5441</td><td></td><td>3.3113</td><td></td><td>2.9679</td><td></td><td>3.9539</td><td></td><td>2.7006</td><td></td><td>6.2463</td><td></td><td>2.8332</td><td></td></tr><tr><td>β</td><td>0.0561</td><td>0.0290</td><td>0.0221</td><td>0.0410</td><td>-0.0501</td><td>0.0282</td><td>-0.0268</td><td>0.0405</td><td>-0.0319</td><td>0.0318</td><td>-0.1598</td><td>0.1035</td><td>-0.0888</td><td>0.0620</td></tr><tr><td>γ<sub>0</sub></td><td>-0.0066</td><td>0.3479</td><td>-0.0088</td><td>0.4385</td><td>0.0030</td><td>0.3952</td><td>-0.0002</td><td>0.5216</td><td>-0.0006</td><td>0.3645</td><td>0.0667</td><td>0.8086</td><td>-0.0157</td><td>0.3811</td></tr><tr><td rowspan="5">20</td><td>γ<sub>1</sub></td><td>0.0010</td><td>0.9961</td><td>-0.0068</td><td>1.2558</td><td>-0.0003</td><td>1.1330</td><td>0.0030</td><td>1.5013</td><td>-0.0068</td><td>1.0417</td><td>0.0249</td><td>2.3014</td><td>-0.0251</td><td>1.0842</td></tr><tr><td>D<sub>EF</sub></td><td>1.3729</td><td></td><td>1.7353</td><td></td><td>1.5564</td><td></td><td>2.0635</td><td></td><td>1.4380</td><td></td><td>3.2134</td><td></td><td>1.5272</td><td></td></tr><tr><td>β</td><td>0.0337</td><td>0.0179</td><td>0.0060</td><td>0.0260</td><td>-0.0390</td><td>0.0194</td><td>-0.0290</td><td>0.0264</td><td>-0.0250</td><td>0.0212</td><td>-0.1532</td><td>0.0770</td><td>-0.0838</td><td>0.0429</td></tr><tr><td>γ<sub>0</sub></td><td>-0.0224</td><td>0.2087</td><td>-0.0268</td><td>0.2679</td><td>-0.0208</td><td>0.2390</td><td>-0.0255</td><td>0.3156</td><td>-0.0198</td><td>0.2229</td><td>0.0209</td><td>0.5222</td><td>-0.0235</td><td>0.2336</td></tr><tr><td>γ<sub>1</sub></td><td>-0.0288</td><td>0.6044</td><td>-0.0418</td><td>0.7704</td><td>-0.0389</td><td>0.6912</td><td>-0.0435</td><td>0.9073</td><td>-0.0375</td><td>0.6436</td><td>-0.0530</td><td>1.5006</td><td>-0.0329</td><td>0.6756</td></tr><tr><td rowspan="5">30</td><td>D<sub>EF</sub></td><td>0.8310</td><td></td><td>1.0643</td><td></td><td>0.9497</td><td></td><td>1.2492</td><td></td><td>0.8877</td><td></td><td>2.0998</td><td></td><td>0.9522</td><td></td></tr><tr><td>β</td><td>0.0184</td><td>0.0107</td><td>0.0058</td><td>0.0165</td><td>-0.0231</td><td>0.0125</td><td>-0.0132</td><td>0.0168</td><td>-0.0176</td><td>0.0136</td><td>-0.1315</td><td>0.0517</td><td>-0.0899</td><td>0.0290</td></tr><tr><td>γ<sub>0</sub></td><td>0.0186</td><td>0.1198</td><td>0.0191</td><td>0.1601</td><td>0.0237</td><td>0.1426</td><td>0.0254</td><td>0.1959</td><td>0.0213</td><td>0.1297</td><td>0.0914</td><td>0.3607</td><td>0.0035</td><td>0.1357</td></tr><tr><td>γ<sub>1</sub></td><td>0.0333</td><td>0.3363</td><td>0.0328</td><td>0.4503</td><td>0.0350</td><td>0.3996</td><td>0.0418</td><td>0.5548</td><td>0.0311</td><td>0.3621</td><td>0.0701</td><td>0.9982</td><td>0.0186</td><td>0.3804</td></tr><tr><td>D<sub>EF</sub></td><td>0.4669</td><td></td><td>0.6269</td><td></td><td>0.5546</td><td></td><td>0.7675</td><td></td><td>0.5055</td><td></td><td>1.4106</td><td></td><td>0.5450</td><td></td></tr><tr><td rowspan="5">50</td><td>β</td><td>0.0122</td><td>0.0051</td><td>0.0065</td><td>0.0072</td><td>-0.0093</td><td>0.0056</td><td>-0.0033</td><td>0.0078</td><td>-0.0076</td><td>0.0061</td><td>-0.0926</td><td>0.0282</td><td>-0.0691</td><td>0.0183</td></tr><tr><td>γ<sub>0</sub></td><td>-0.0106</td><td>0.0635</td><td>-0.0116</td><td>0.0846</td><td>-0.0096</td><td>0.0753</td><td>-0.0114</td><td>0.0992</td><td>-0.0096</td><td>0.0706</td><td>0.0300</td><td>0.1701</td><td>-0.0215</td><td>0.0713</td></tr><tr><td>γ<sub>1</sub></td><td>-0.0158</td><td>0.1820</td><td>-0.0175</td><td>0.2424</td><td>-0.0173</td><td>0.2158</td><td>-0.0183</td><td>0.2846</td><td>-0.0170</td><td>0.2024</td><td>-0.0138</td><td>0.4867</td><td>-0.0190</td><td>0.2022</td></tr><tr><td>D<sub>EF</sub></td><td>0.2505</td><td></td><td>0.3343</td><td></td><td>0.2967</td><td></td><td>0.391</td><td></td><td>0.2792</td><td></td><td>0.6851</td><td></td><td>0.2918</td><td></td></tr><tr><td>β</td><td>0.0069</td><td>0.0024</td><td>0.0007</td><td>0.0033</td><td>-0.0057</td><td>0.0027</td><td>-0.0049</td><td>0.0035</td><td>-0.0027</td><td>0.0031</td><td>-0.0783</td><td>0.0200</td><td>-0.0467</td><td>0.0093</td></tr><tr><td rowspan="5">200</td><td>β</td><td>-0.0101</td><td>0.0293</td><td>-0.0076</td><td>0.0391</td><td>-0.0075</td><td>0.0347</td><td>-0.0085</td><td>0.0484</td><td>-0.0069</td><td>0.0313</td><td>0.0211</td><td>0.0902</td><td>-0.0204</td><td>0.0353</td></tr><tr><td>γ<sub>0</sub></td><td>-0.0178</td><td>0.0824</td><td>-0.0151</td><td>0.1108</td><td>-0.0159</td><td>0.0981</td><td>-0.0174</td><td>0.1374</td><td>-0.0151</td><td>0.0887</td><td>-0.0248</td><td>0.2510</td><td>-0.0227</td><td>0.0985</td></tr><tr><td>γ<sub>1</sub></td><td>0.1141</td><td></td><td>0.1532</td><td></td><td>0.1355</td><td></td><td>0.1893</td><td></td><td>0.1231</td><td></td><td>0.3612</td><td></td><td>0.1431</td><td></td></tr><tr><td>D<sub>EF</sub></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>
<table><caption>TABLE 4 Percentiles and MTTF using true model parameters.</caption><thead><tr><td>Life characteristic</td><td>True value</td></tr></thead><tbody><tr><td>B<sub>1</sub></td><td>1025.8059</td></tr><tr><td>B<sub>5</sub></td><td>3040.7815</td></tr><tr><td>B<sub>10</sub></td><td>4913.5654</td></tr><tr><td>MTTF</td><td>19884.2883</td></tr></tbody></table>

In the same way, we shall compare the MTTF life characteristics. Figure 3 displays the log-transformed bias and MSE of the estimated MTTF using different distances and for different sample sizes. From Figure 3, it is obvious that for small sample sizes (n = 5), the predictions provided by the $ρ_{LH}$ distance are slightly better as the bias is closer to 0 and the MSE value is the lowest. However, as we increase the sample size, the advantage of the ML method becomes less pronounced. For different sample sizes, not much differences were observed between all the distance functions except for $ρ_{AD2L}$ that provides the worst predictions since it presents the highest bias and MSE values of the MTTF. Furthermore, it is important to mention that for n = 200, the $ρ_{AD2R}$ outperforms $ρ_{LH}$ distance by providing a better prediction of the MTTF with roughly the same magnitude of the MSE.

## 4.2 | Proportion of correct selection

In this subsection, we shall compute the proportion of correct selection (PCS) using different sample sizes and for all the distances by considering very close accelerated lifetime models and nested ALT models. The idea here is to compare the performance of each distance to recover the ‘parent model’. Since it is not possible to consider all possible models in a single study, the choice is made to generate the simulated data sets through $M_1$ and to consider only $M_1, M_2, M_3$ and $M_5$ (already presented in Section 3.2) for the discrimination studies.

This choice is made to, first, keep a reasonable computational time and, second, to present the results for the main scenarios: (1) discrimination between similar models ($M_1$ vs. $M_3$ and $M_1$ vs. $M_5$) and (2) discrimination between nested models ($M_1$ vs. $M_2$). The computation of the PCS matrices was performed using 1000 replications.

### 4.2.1 | Discrimination between close ALT models

The matrices exhibiting the PCS for $M_1$ versus $M_3$ and $M_1$ versus $M_5$ discrimination studies are shown, respectively, in Figures 4 and 5.

One can clearly see, from both figures, that the PCS increases when the sample size increases regardless the distance function. It is apparent from both discrimination studies that the ABC approach using $ρ_{ADL}$ and $ρ_{AD2L}$ distances is performing the best for small sample sizes compared to the likelihood-based distance, through exhibiting the highest proportion of ‘parent model’ selection. This is quite interesting, as these distances could be efficiently used to mitigate the risk of model misspecification by inferring and discriminating very close acceleration life models when observed failure times are limited.

It is important to mention that the evolution of $M_1$ versus $M_5$'s PCS for different sample sizes shares the same trend as that presented in $M_1$ versus $M_3$'s PCS matrix.

Furthermore, from both figures, one can see that $ρ_{ADR}$ and $ρ_{AD2R}$ distances exhibit the lowest PCS for both cases. This can be explained by the fact that as these distances give more weight to the right-tail region where the available data are limited (the three competing ALT models are right-skewed). It will be difficult to discriminate between overlapping models. One can conclude that these distances are not suitable to discriminate between very close acceleration life models.

### 4.2.2 | Discrimination between nested ALT models

In the following study, one investigates the performance of the ABC approach and the information criteria (IC) to recover the correct model by considering nested ALT models for different sample sizes. The PCS matrices are estimated for different tolerance thresholds ($Δε$) using both ABC algorithm and IC.
<figure><img src="image_3.png" /><figcaption>FIGURE 1 Bias of the low percentiles predicted using different distances and for different sample sizes.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIGURE 2 Mean squared errors of low percentiles predicted using different distances and for different sample sizes.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIGURE 3 Mean squared error (MSE) and bias of the predicted MTTF using different distances and for different sample sizes.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 4 Proportion of correct selection between $M_1$ and $M_3$ accelerated life models when the data are generated from model $M_1$.</figcaption></figure>

Two popular IC have been used: the AIC and the BIC. They are given by

$$ \text{AIC} = -2 \ln(\mathcal{L}(\hat{\theta})) + 2k, \quad (40) $$

where k stands for the number of model parameters and $\mathcal{L}(\hat{\theta})$ denotes the likelihood function evaluated at the estimated parameters values.

$$ \text{BIC} = -2 \ln(\mathcal{L}(\hat{\theta})) + k \ln(N), \quad (41) $$

where N represents the number of observed data.
<figure><img src="image_7.png" /><figcaption>FIGURE 5 Proportion of correct selection for $M_1$ and $M_5$ accelerated life models discrimination when the data are generated from model $M_1$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 6 Proportion of correct selection for $M_1$ and $M_2$ accelerated life models discrimination for $\xi = 10^{-2}$.</figcaption></figure>

It should be mentioned that for a set of candidate models, the selected model is the one that presents the lowest IC value. Hereafter, $\rho_{LH}$ distance is not considered in PCS matrix since the competing models do not have the same dimensionality. Hence, for the same data sets generated from model $M_1$, the estimated PCS matrix for $\Delta\epsilon \le 10^{-2}$ is illustrated in Figure 6. Similarly, a PCS matrix was estimated for $\Delta\epsilon \le 10^{-6}$. However, it is not presented in this paper since the proportions of parent model selection tend towards zero for all the considered distances regardless the sample sizes. The following observations can be drawn:

* Overall, we notice that, for small tolerance thresholds, all the considered distances largely favour the complex model ($M_2$).

* The proportions of $M_2$ selection obtained for $\Delta\epsilon = 10^{-2}$ are less pronounced than those for $\Delta\epsilon = 10^{-6}$. This proves that tightening the expected tolerance threshold reduces the selection probability of the simple model.

* It can be seen that increasing sample sizes for a given tolerance thresholds favours the selection of $M_2$ regardless the considered distance and the tolerance threshold. This could be explained by the fact that $M_2$ is more flexible and offer a better fitting quality to the data.

* Conversely to the ABC algorithm, AIC and BIC ensure the selection of simpler model (i.e. $M_1$) regardless the tolerance threshold. The PCS for IC increases slightly for large sample sizes. Indeed, for $n = 5$, AIC and BIC present almost the same PCS value, respectively, 77.5 and 78.5%, and for $n = 200$, the PCS increases to, respectively, 86 and 86.1%.

* Independently of the considered tolerance thresholds, the IC penalises complex models and leads to the selection of the simple model. This could be justified by the fact that models with low number of parameters generalise better.

To sum up, the ABC-NS algorithm coupled with appropriate distances ensures, depending on the user's objective, the flexibility of selecting the complex model and infer its parameters to better fit the considered data (low tolerance threshold) or generalising for unseen data (large tolerance threshold).
<table><caption>TABLE 5 Definition of the priors on the parameters of the accelerated life models.</caption><thead><tr><th>Model</th><th>Priors</th><th></th><th></th></tr></thead><tbody><tr><td>M<sub>1</sub></td><td>β ~ U(0, 20)</td><td>γ<sub>0</sub> ~ U(0, 100)</td><td>γ<sub>1</sub> ~ U(0, 50)</td></tr><tr><td>M<sub>2</sub></td><td>γ'<sub>0</sub> ~ U(0, 10)</td><td>γ'<sub>1</sub> ~ U(-10, 10)</td><td>γ<sub>0</sub> ~ U(0, 100)</td><td>γ<sub>1</sub> ~ U(0, 50)</td></tr><tr><td>M<sub>3</sub></td><td>β ~ U(0, 5)</td><td>γ<sub>0</sub> ~ U(0, 100)</td><td>γ<sub>1</sub> ~ U(0, 50)</td></tr><tr><td>M<sub>4</sub></td><td>γ'<sub>0</sub> ~ U(-10, 10)</td><td>γ'<sub>1</sub> ~ U(0, 20)</td><td>γ<sub>0</sub> ~ U(0, 100)</td><td>γ<sub>1</sub> ~ U(0, 50)</td></tr><tr><td>M<sub>5</sub></td><td>σ ~ U(0, 3)</td><td>γ<sub>0</sub> ~ U(0, 10)</td><td>γ<sub>1</sub> ~ U(-10, 10)</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>FIGURE 7 Evolution of the model posterior probabilities over some selected populations using the $\rho_{AD}$ distance.</figcaption></figure>

## 5 | REAL DATA APPLICATION

In this section, the ABC-NS algorithm is employed to discriminate between five competing accelerated life models (see Table 1) using real data. It is constituted of 304 fatigue failure times using 6061-T6 aluminium coupons under three different stress levels: $S_1 = 31,000$ Psi with $n_1 = 101$ observations, $S_2 = 26,000$ Psi with $n_2 = 102$ observations and $S_3 = 21000$ Psi with $n_3 = 101$ observations. All tests have been conducted until failure (see Ref. 61 for more details). As in the first part, a comparison between likelihood-based approach and likelihood-free method applied for model selection is conducted.

For the ABC-NS algorithm, vague uniform priors based on crude estimates from a least square method are employed. The estimates are obtained, for each model, through a preliminary analysis that consists of using the data collected at the stress levels. The priors on the parameters of the considered five accelerated life models are given in Table 5. Furthermore, equal prior probabilities are assumed for the five competing models, that is, $P(M_i) = \frac{1}{5}$. The other hyper-parameters of ABC-NS algorithm are kept the same as for the simulation study.

### 5.1 | Results

In this subsection, we investigate the performance of the ABC-NS algorithm using the $\rho_{AD}$ distance for illustrative purpose. This choice is motivated by the fact that this distance provides the best predictions at low lifetime percentiles and very efficient when very close models are in competition. The posterior probabilities over the populations are shown in Figure 7.

A few important observations are worth mentioning: first, the parsimony principle is well embedded in the ABC-NS algorithm; the algorithm tries, first, to select the simpler models ($M_1, M_3$ and $M_5$). Then, when these latter are no longer able to fit the data very well, the algorithm switches to a more complex models ($M_2$ and $M_4$). Second, one can see that $M_1$ is eliminated first, then $M_5$ pursued by $M_3$. After that, the ABC-NS algorithm oscillates between $M_2$ and $M_4$ and finishes by converging to model $M_2$. This model has been selected around the 100th population with a posterior probability equal to 1. The algorithm continues to run until reaching the stopping condition. In other words, from the 100th
<figure><img src="image_10.png" /><figcaption>FIGURE 9 Acceptance rate over the populations.</figcaption></figure>

<figure><img src="image_11.png" /></figure>

to the 136th population, the algorithm tends to refine the posterior estimates. Figure 8 displays the marginal distributions of the selected model parameters. As one can see, they are very well peaked with low variability, the value of the distance function at the last population is equal to $ρ_{AD} = 3.142$.

Figure 9 shows the acceptance rate over the populations computed by dividing the number of particles needed to replenish each population by the total number of simulated particles. One can clearly see that the ABC-NS algorithm ensures a relatively high acceptance rates. It is worth mentioning that for the first populations, the acceptance rate decreases because the algorithm explores large model and parameter spaces, then it steadily rises as the sampling spaces tighten and the least likely models are eliminated.

For comparison purpose, the algorithm is now employed with the other distance functions in order to select the adequate model(s) and to compare the fitting qualities. Table 6 gives the selected model for the different distances. Again, model $M_2$ is selected as the most likely model to fit the accelerated data when using $ρ_{CM}$, $ρ_{ADL}$ and $ρ_{AD2L}$ distances (rank 1 in bold in the table). However, for $ρ_{ADR}$ and $ρ_{AD2R}$ distances, the most likely model seems to be model $M_4$ (rank 1 in bold in the table). The obtained results show that the selected model may vary depending on the selected distance function which measures the level of agreement between the empirical and the simulated data and also expresses a desired objective of the user.

Table 7 shows that the differences on the percentiles estimated by the different selected models and distances are significant. One important point that should be highlighted is that the $ρ_{AD2L}$ and $ρ_{ADL}$ distances tend to be more ‘pessimistic’ to predict lifetime at the left-tail. This lies to the fact that both distances give more weights to the left observations which rises the risk of over-fitting. In other words, they provide a good fitting quality at the left-tail region when observed data are available, but they are poor in terms of predictions (for unseen data). Both $ρ_{ADR}$ and $ρ_{AD2R}$ provide probably an
<table><caption>TABLE 6 Model ranking for the considered distance functions.</caption><thead><tr><th rowspan="2">Distance function</th><th colspan="5">Competing accelerated life models</th></tr><tr><th>M<sub>1</sub></th><th>M<sub>2</sub></th><th>M<sub>3</sub></th><th>M<sub>4</sub></th><th>M<sub>5</sub></th></tr></thead><tbody><tr><td>ρ<sub>CM</sub></td><td>3</td><td>1</td><td>4</td><td>2</td><td>5</td></tr><tr><td>ρ<sub>AD</sub></td><td>3</td><td>1</td><td>4</td><td>2</td><td>5</td></tr><tr><td>ρ<sub>ADL</sub></td><td>3</td><td>1</td><td>4</td><td>2</td><td>5</td></tr><tr><td>ρ<sub>ADR</sub></td><td>4</td><td>2</td><td>3</td><td>1</td><td>5</td></tr><tr><td>ρ<sub>AD2L</sub></td><td>2</td><td>1</td><td>4</td><td>3</td><td>5</td></tr><tr><td>ρ<sub>AD2R</sub></td><td>4</td><td>3</td><td>2</td><td>1</td><td>5</td></tr></tbody></table>

<table><caption>TABLE 7 Prediction of some life characteristics under normal stress.</caption><thead><tr><th rowspan="2">Distance function</th><th rowspan="2">Selected model</th><th colspan="4">Life characteristics/(cycles number)</th></tr><tr><th>B<sub>1</sub></th><th>B<sub>5</sub></th><th>B<sub>10</sub></th><th>MTTF</th></tr></thead><tbody><tr><td>ρ<sub>LH</sub></td><td>M<sub>2</sub></td><td>3475.5130</td><td>6277.591</td><td>8150.623</td><td>16408.188</td></tr><tr><td>ρ<sub>CM</sub></td><td>M<sub>2</sub></td><td>2051.3905</td><td>4504.1864</td><td>6374.7216</td><td>16725.2032</td></tr><tr><td>ρ<sub>AD</sub></td><td>M<sub>2</sub></td><td>2575.4527</td><td>5201.1226</td><td>7094.1467</td><td>16586.8824</td></tr><tr><td>ρ<sub>ADL</sub></td><td>M<sub>2</sub></td><td>2002.9660</td><td>4452.2793</td><td>6335.5730</td><td>16910.1358</td></tr><tr><td>ρ<sub>ADR</sub></td><td>M<sub>4</sub></td><td>6019.4383</td><td>7939.3868</td><td>9255.9937</td><td>17834.3071</td></tr><tr><td>ρ<sub>AD2L</sub></td><td>M<sub>1</sub></td><td>1734.9376</td><td>4121.2831</td><td>6039.1011</td><td>17699.7480</td></tr><tr><td>ρ<sub>AD2R</sub></td><td>M<sub>4</sub></td><td>8949.3723</td><td>10768.5230</td><td>11905.1746</td><td>17734.2997</td></tr></tbody></table>

<table><caption>TABLE 8 Comparison between the competing models using two information criteria.</caption><thead><tr><th>Model</th><th>Maximum Log-likelihood <i>L</i></th><th>AIC</th><th>BIC</th></tr></thead><tbody><tr><td>M<sub>1</sub></td><td>-1793.0338</td><td>3592.0675</td><td>3603.2186</td></tr><tr><td>M<sub>2</sub></td><td><b>-1782.7183</b></td><td><b>3573.4365</b></td><td><b>3588.3046</b></td></tr><tr><td>M<sub>3</sub></td><td>-1806.5762</td><td>3619.1523</td><td>3630.3034</td></tr><tr><td>M<sub>4</sub></td><td><b>-1784.6265</b></td><td><b>3577.2530</b></td><td><b>3592.1212</b></td></tr><tr><td>M<sub>5</sub></td><td>-1818.1496</td><td>3642.2994</td><td>3653.4504</td></tr></tbody></table>

over-estimated values of lower percentiles as one can see from the same table. It is obvious that these distances are not recommended in predicting the 'safe life' at lower percentiles. Another reason to explain the big difference in terms of lifetime predictions is that both distances select $M_4$ as the most likely model. Despite that the selected models are very close, they yield quite different values. It is obvious from the obtained results that the choice of the model but also the distance used to infer the accelerated life model play a crucial role with regard to prediction of tail percentiles. Based on the obtained results and for the purpose of specifying a safe life, one can say that $M_2$ provides the best compromise between a good fitting quality and a better generalization (i.e. prediction of unseen data). Finally, the disparity between the predicted MTTF given by the used distances is less pronounced as one can see from Table 7.

## 5.2 | Model selection through likelihood-based approach

Table 8 gives an estimate of AIC and BIC criterion for the competing models. One can see that the most likely model using both criteria is $M_2$ where it presents the lowest IC values (in bold in the table). In addition, it can be seen that $M_4$ provides an acceptable fit to the data as the difference between the IC values is not significant. The obtained results from the likelihood-based approach is coherent with the likelihood-free method. It is obvious that the $M_2$ and $M_4$ are the most likely models to fit the data (the values are in bold in the table). The advantage of the likelihood-free method is that it can be tailored through an appropriate choice of the distance function according to a specific user objective.
# 6 | CONCLUDING REMARKS

In this paper, a comparison between likelihood-free and likelihood-based approaches has been conducted to deal with accelerated life model selection and parameter estimation. Several distance functions have been tested for the first approach to get more insights. A focus on the ability of each approach to precisely estimate the model parameters, to predict some lower percentiles and to discriminate between very close competing accelerated life models have been investigated. It has been shown that each approach has its own strength and its own drawbacks. The likelihood-based approach is a good estimator but presents some weaknesses in predicting lower percentiles and also to discriminate between very close models when the available data are limited. However, the likelihood-free approach coupled with an appropriate distance (selected based on the user's goal) may be a good alternative to overcome these issues as it has been demonstrated. In addition, the ABC method is a single step procedure as all the models are considered simultaneously and the least likely models are ruled out progressively over the populations. The presented results indicated that the likelihood-free Bayesian approach applied for the first time to discriminate between accelerated life models is a promising option for practical applications regarding reliability of systems/components, mainly when a good predictions of lower percentiles lifetime is of major concern. The method can be tailored to our needs through an appropriate selection of a distance function/metric.

The overall conclusion here is that the likelihood-free method may be of considerable value in engineering reliability problems for selecting a 'useful' model. It is flexible in the sense that different distance functions can be used and can be considered as a promising alternative to overcome the drawbacks of the likelihood-based approaches. As it has been demonstrated, it has a good prediction capability especially for lower percentiles. Although, we have discussed only the complete data, but similar results can be obtained for censored data also. The work is in progress and it will be reported later. The improvement of the ABC-NS algorithm using reinforcement learning is another avenue under investigation to gain more in efficiency.

## ACKNOWLEDGEMENTS

We would like to thank the University of Angers and Liebherr Aerospace Toulouse for providing facilities throughout the project. This article has benefited from the contributions of Mr Anis Ben Abdessalem, Mr Laurent Saintis, Prof Bruno Castanier and Mr Rodrigue Sohoin. We thank the journal and the reviewers for the thoughtful comments and constructive suggestions and for the attention given to this work. This research received no specific grant from any funding agency in the public, commercial or not-for-profit sectors.

## DATA AVAILABILITY STATEMENT

The data and the codes that support the findings of this study are available upon request.

## ORCID

Mohamed Rabhi https://orcid.org/0000-0002-6815-0553

## REFERENCES

1. Pan R. A bayes approach to reliability prediction utilizing data from accelerated life tests and field failure observations. *Qual Reliab Eng Int.* 2009;25:229-240.

2. Fan TH, Yu CH. Statistical inference on constant stress accelerated life tests under generalized gamma lifetime distributions. *Qual Reliab Eng Int.* 2013;29(5):631-638.

3. Freels JK, Pignatiello JJ, Warr RL, Hill RR. Bridging the gap between quantitative and qualitative accelerated life tests. *Qual Reliab Eng Int.* 2015;31(5):789-800.

4. Meeker WQ, Escobar LA. *Statistical Methods for Reliability Data*. John Wiley and Sons; 2014.

5. Carlsson B, Möller K, Köhl M, et al. The applicability of accelerated life testing for assessment of service life of solar thermal components. *Sol Energy Mater. Sol Cells*. 2004;84:255-274.

6. Gaertner G, Raasch D, Barratt D, Jenkins S. Accelerated life tests of CRT oxide cathodes. *Appl Surf Sci.* 2003;215:72-77.

7. Whitman CS. Accelerated life test calculations using the method of maximum likelihood: an improvement over least squares. *Microelectron Reliab.* 2003;43:859-864.

8. Barlow R, Toland R, Freeman T. A Bayesian analysis of the stress rupture life of Kevlar/epoxy spherical pressure vessels. In: *Accelerated Life Testing and Experts Opinions in Reliability*. Elsevier Science Ltd.; 1988;203-236.

9. D'Anna G, Giorgio M, Riccio A. Estimating fatigue reliability of structural components via a Birnbaum-Saunders model with stress dependent parameters from accelerated life test. *Compos Part B*. 2017;119:206-214.

10. Birnbaum ZW, Saunders SC. A new family of life distributions. *J Appl Probab.* 1969;6(2):319-327.
11. Box GE, Draper NR. *Empirical Model-Building and Response Surfaces*. Wiley Series in Probability and Mathematical Statistics. Wiley; 1987.

12. Shanshan L, Zhanwen N, Guodong W, Qu L, Zhen H. Lower percentile estimation of accelerated life tests with nonconstant scale parameter. *Qual Reliab Eng Int.* 2017;33:1437-1446.

13. Settanni E, Newnes LB, Thenent NE, Bumblauskas D, Parry G, Goh YM. A case study in estimating avionics availability from field reliability data. *Qual Reliab Eng Int.* 2016;32(4):1553-1580.

14. Prabhakar Murthy DN, Bulmer M, Eccleston JA. Weibull model selection for reliability modelling. *Reliab Eng Syst Saf.* 2004;86(3):257-267.

15. Jukic D, Bensic M, Scitovski R. On the existence of the nonlinear weighted least squares estimate for a three-parameter Weibull distribution. *Comput Stat Data Anal.* 2008;52(9):4502-4511.

16. Li H, Zuo H, Liu R, Liu J, Jing C. Evaluation of aileron actuator reliability with censored data. *Chin J Aeronaut.* 2015;28(4):1087-1103.

17. Jiang R, Murthy DN. Reliability modeling involving two Weibull distributions. *Reliab Eng Syst Saf.* 1995;47(3):187-198.

18. Zhang T, Dwight R. Choosing an optimal model for failure data analysis by graphical approach. *Reliab Eng Syst Saf.* 2013;115:111-123.

19. Kalbfleisch JG. *Probability and Statistical Inference*. 2nd ed. Springer-Verlag; 1979:1-60.

20. Akaike H. A new look at statistical model identification. *IEEE Trans Autom Control.* 1974;19(6):716-722.

21. Schwarz G. Estimating the dimension of a model. *Ann Statist.* 1978;6(2):461-464.

22. Ghaly A, Aly HM, Salah RN. Different estimation methods for constant stress accelerated life test under the family of the exponentiated distributions. *Qual Reliab Eng Int.* 2016;32:1095-1108.

23. Gupta A, Mukherjee B, Upadhyay SK. Weibull extension model: a Bayes study using Markov chain Monte Carlo simulation. *Reliab Eng Syst Saf.* 2008;93(10):1434-1443.

24. Santana D, Filho CF, Márcio IS, Martins A. A novel Bayesian approach to reliability modeling: the benefits of uncertainty evaluation in the model selection procedure. *Qual Reliab Eng Int.* 2018;34:1127-1141.

25. Al-Garni AZ, Jamal A. Artificial neural network application of modeling failure rate for boeing 737 tires. *Qual Reliab Eng Int.* 2011;27(2):209-219.

26. Ling J, Pan J. A new method for selection of population distribution and parameter estimation. *Reliab Eng Syst Saf.* 1998;60:247-255.

27. Green PJ. Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. *Biometrika*. 1995;80(4):711-732.

28. Dellaportas P, Jonathan J, Ntzoufras I. On Bayesian model and variable selection using MCMC. *Stat Comput.* 2002;12:27-36.

29. Godsill SJ. On the relationship between Markov chain Monte Carlo methods for model uncertainty. *J Comput Graph Stat.* 2001;10(2):230-248.

30. Green PJ, Hastie DI. Model choice using reversible jump Markov chain Monte Carlo. *Statistica Neerlandica.* 2012;66:309-338.

31. Jeffreys H. *Theory of Probability*. 1st ed. Oxford University Press; 1939:10-56.

32. Sansó B, Pericchi LR, Moreno E, Racugno W. On the robustness of the intrinsic Bayes factor for nested models. *IMS Lect Notes Monogr Ser.* 1996;6:155-173.

33. Berger JO, Pericchi LR. The intrinsic Bayes factor for model selection and prediction. *J Am Stat Assoc.* 1996;91(433):109-122.

34. Puthenpura S, Sinha NK. Modified maximum likelihood method for the robust estimation of system parameters from very noisy data. *Automatica*. 1986;22(2):231-235.

35. Vaughan DC. The generalized secant hyperbolic distribution and its properties. *Commun Stat-Theory Methods*. 2002;31(2):219-238.

36. Acitas S, Kasap P, Senoglu B, Arslan O. Robust estimation with the skew distribution. *Pak J Stat.* 2013;29(8):409-430.

37. Acitas S, Kasap P, Senoglu B, Arslan O. One-step m-estimators: Jones and Faddy's skewed t-distribution. *J Appl Stat.* 2013;40(3):1545-1560.

38. Balakrishnan N, Ling MH. Model misspecification analyses of Weibull and gamma models based on one-shot device test data. *IEEE Trans.* 2017;66(3):641-650.

39. Khakifirooz M, Tseng ST, Fathi M. Model misspecification of generalized gamma distribution for accelerated lifetime-censored data. *Technometrics*. 2020;62(3):357-370.

40. Balakrishnan N, Castilla E, Ling MH. Optimal designs of constant-stress accelerated life-tests for one-shot devices with model misspecification analysis. *Qual Reliab Eng Int.* 2022;38:989-1012.

41. Beaumont MA, Zhang W, Balding DJ. Approximate Bayesian computation in population genetics. *Genetics*. 2002;162:2025-2035.

42. Bonassi V, West M. Sequential Monte Carlo with adaptive weights for approximate Bayesian computation. *Bayesian Anal.* 2015;10:171-187.

43. Ben Abdessalem A, Dervilis N, Wagg D, Worden K. Model selection and parameter estimation of dynamical systems using a novel variant of approximate Bayesian computation. *Mech Syst Signal Process.* 2019;122:364-386.

44. Rabhi M, Ben Abdessalem MA, Saintis L, Castanier B, Sohoin RK. Parameters estimation of accelerated lifetime testing models using an efficient approximate Bayesian computation method. In: ESREL. 2021.

45. Toni T, Welch D, Strelkowa N, Ipsen A, Stumpf MPH. Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems. *J R Soc Interface*. 2009;6:187-202.

46. Barnes CP, Filippi S, Stumpf MP, Thorne T. Considerate approaches to constructing summary statistics for ABC model selection. *Stat Comput.* 2012;22(6):1181-1197.

47. Everitt RG, Johansen AM, Rowing E, Evdemon-Hogan M. Bayesian model comparison with un-normalised likelihoods. *Stat Comput.* 2017;27(2):403-422.

48. Ong VM, Nott DJ, Tran MN, Sisson SA, Drovandi CC. Variational Bayes with synthetic likelihood. *Stat Comput.* 2018;28(4):971-988.

49. Pritchard JK, Seielstad MT, Perez-Lezaun A, Feldman MW. Population growth of human Y chromosomes: a study of Y chromosome microsatellites. *Mol Biol Evol.* 1999;16(12):1791-1798.
50. Sunnaker M, Busetto AG, Numminen E, Corander J, Foll M, Dessimoz C. Approximate Bayesian computation. *PLOS Comput Biol.* 2013;9(1):e1002803.

51. Chiachio M, Beck JL, Chiachio J, Rus G. Approximate Bayesian computation by subset simulation. *SIAM J Sci Comput.* 2014;36(3):A1339-A1358.

52. Picchini U, Anderson R. Approximate maximum likelihood estimation using data-cloning ABC. *Comput Stat Data Anal.* 2017;105:166-183.

53. Mukherjee P, Parkinson D, Liddle AR. A nested sampling algorithm for cosmological model selection. *Astrophys J.* 2006;638:L51-L54.

54. Shaw JR, Bidges M, Hobson MP. Efficient Bayesian inference for multimodal problems in cosmology. *Mon Notice R Astron Soc.* 2006;378:1-7.

55. Luceno A. Fitting the generalized pareto distribution to data using maximum goodness-of-fit estimator. *Comput Stat Data Anal.* 2006;51:904-917.

56. Anderson TW, Darling DA. A test of goodness-of-fit. *J Am Stat Assoc.* 1954;49:765-769.

57. Farrel PJ, Stewart KR. Comprehensive study of tests for normality and symmetry: extending the Spiegelhalter test. *J Stat Comput Simul.* 2006;76:803-816.

58. Cao QV, McCarty SM. New methods for estimating parameters of Weibull functions to characterize future diameter distributions in forest stands. In: 13th Biennial Southern Silvicultural Research Conference. 2006;338-340.

59. Bai DS, Chung SW. An accelerated life test model with the inverse power law. *Reliab Eng Syst Saf.* 1989;24:223-230.

60. Park SH, Kim JH. Lifetime estimation of led lamp using gamma process model. *Microelectron Reliab.* 2016;57:71-78.

61. Birnbaum ZW, Saunders SC. Estimation for a family of life distributions with applications to fatigue. *J Appl Probab.* 1969;6(2):328-347.

## AUTHOR BIOGRAPHIES

<figure><img src="image_12.png" /></figure>

**Mohamed Rabhi** earned a PHD in Reliability engineering from University of Angers – France where he worked on the development of a methodology for quantifying the impact of preponderant stresses on the reliability of complex systems (2023) and a Mechanical Engineer from the National Engineering School of Tunis (2017). He currently holds a position of Reliability Specialist at Liebherr Aerospace Toulouse where he is in charge of the methodological development of accelerated life testing of complex systems.

<figure><img src="image_13.png" /></figure>

**Anis B. Abdessalem** (ABA) is a Lecturer in the Department of Mechanical Engineering at the University of Angers and a member of the dependability and decision support group. His background is structural dynamics, nonlinear system identification, reliability, stochastic modelling and Bayesian inference. ABA obtained a PhD in Physical, Mathematics and Information Sciences for Engineers (2007–2011) from INSA Rouen. From 2015 to 2017, he joined the Dynamics Research Group at the University of Sheffield to work on the EPSRC Programme Grant ‘Engineering Nonlinearity’ (EK/K003836). His current research sits between reliability analysis of complex systems, model selection, verification and validation in situations where data are sparse in multiple engineering disciplines.

<figure><img src="image_14.png" /></figure>

**Laurent Saintis** is an Associate Professor at the University of Angers in the Department of Quality, Innovation and Reliability Engineering. He is actually head of dependability research team of the LARIS lab. He works principally on reliability testing and complex systems modelling for dependability evaluation. Since 2015, he is mainly involved in the design of reliability test plan for electronic medical devices, and is responsible for the RECOME project for the laboratory. His main teaching themes are dependability, reliability assessment, maintainability design and value analysis.

<figure><img src="image_15.png" /></figure>

**Bruno Castanier** is a Professor at Université d'Angers in the Department of Quality, Innovation and Reliability Engineering. He conducts his research in the domain of mathematical models for reliability assessment and maintenance optimization. He has published more than 30 papers in the main journal of the domains and more than 100 contributions in international conferences. His research activity is strongly connected with industry. He has participated as leader or contributor to 15 different research contracts. He was the General Chair of the 31st European Reliability and Safety Conference – ESREL 2021 – in Angers, France.
<figure><img src="image_16.png" /></figure>

**Rodrigue Sohoin** earned a PHD in Reliability Engineering from the University of Angers - France where he worked on the definition of a methodology for reliability estimation and qualification of mechanical systems in the development stage. He currently holds a position of Reliability Specialist at Liebherr Aerospace Toulouse where he is in charge of reliability assessment for complex systems.

**How to cite this article:** Rabhi M, Abdessalem AB, Saintis L, Castanier B, Sohoin R. Discrimination between accelerated life models via Approximate Bayesian Computation. *Qual Reliab Eng Int.* 2023;1-25.

https://doi.org/10.1002/qre.3283