STATISTICAL ANALYSIS OF PROGRESSIVE STRESS ACCELERATED LIFE TEST FOR THE PRODUCT OF TWO-PARAMETER LAPLACE BS FATIGUE LIFE DISTRIBUTION UNDER INVERSE POWER LAW MODEL∗

Based on the product of two-parameter Laplace Birnbaum-Saunders fatigue life distribution, its failure distribution mode is theoretically derived under the progressive stress accelerated life test with inverse power law model, and then three-parameter generalized Laplace Birnbaum-Saunders fatigue life distribution is introduced. The basic properties of three-parameter generalized Laplace Birnbaum-Saunders fatigue life distribution are analyzed, and the image characteristics of its density function, failure rate function and average failure rate function are investigated. Meanwhile, the point estimate method is given for three parameters, and then the point estimates of parameters are obtained for the product of two-parameter Laplace Birnbaum-Saunders fatigue life distribution under the progressive stress accelerated life test with inverse power law model. In addition, the practical example and simulation examples are illustrated to show the feasibility of the proposed method.


Introduction
The random variable T is assumed to follow two-parameter Birnbaum-Saunders fatigue life distribution BS(α, β) , then its distribution function and density function are respectively where α > 0 is called a shape parameter, β > 0 is called a scale parameter, and ϕ(x), Φ(x) are respectively the distribution function and density function of normal distribution, that is, Birnbaum and Saunders [2] derived the famous two-parameter BS fatigue life distribution by studying the process of material failure mainly caused by crack propagation. It is more suitable than other common life distributions such as Weibull distribution and Lognormal distribution to describe life rules of some fatigue failure products, and become one of common distributions in the reliability engineering. There have been many researches on two-parameter BS fatigue life distribution, as shown in references [3-8, 10, 13-19, 21, 23, 24, 27].
With the in-depth research of BS fatigue life distribution, many scholars further generalize two-parameter BS fatigue life distribution, which is usually called generalized BS fatigue life distribution. The studies all involve generalized BS fatigue life distribution in references [1,9,20,28]. It is worth pointing out that standard normal distribution involved in original two-parameter BS fatigue life distribution is replaced by standard Laplace distribution in these references, namely two-parameter Laplace BS fatigue life distribution. Zhu and Balakrishnan [28] studied this kind of two-parameter Laplace BS fatigue life distribution in detail. They proposed the numerical characteristics of this distribution and researched the image characteristics of its density function and failure rate function. Besides, they studied maximum likelihood estimates of parameters and proved that maximum likelihood estimates exist uniquely.

Failure Mode of Progressive Stress Accelerated
Life Test V (t) = Kt for Two-parameter LBS Fatigue Life Distribution under Inverse Power Law Model The lifetime T of a product is assumed to follow two-parameter Laplace Birnbaum-Saunders fatigue life distribution, denoted as T ∼ LBS(α, β), then its distribution function F T (t) and density function f T (t) are respectively where α > 0 is called a shape parameter, β > 0 is called a scale parameter and The inverse power law model refers to the relationship β = 1 dV c between the scale parameter β ( in hour ) and the voltage ( in volt ) of some products (such as some insulation materials, capacitors, micro motors and electronic devices, etc.) when the voltage is taken as the accelerated stress, according to the physical principle and experimental experience summary, where d > 0 and c > 0 are both constants. If the product is an electronic component, the physical tests show that c is only connected with its type and not related to its specification.
After taking the logarithm of both sides of above equation, β satisfies the logarithmic linear relationship ln β = a + bφ(V ), where a = − ln d , b = −c and φ(V ) = ln V is the function of the stress V.
Most statistical analyses about step stress or progressive stress accelerated life tests are based on the Nelson Assumption ( CE model for short).
Nelson Assumption [11] : The residual life of a product depends only on the accumulated failure part and the stress level at that time, and is not related to the accumulative mode.
Nelson Assumption is essentially a time conversion, that is, if a constant stress is kept on, the products that have not been failed will fail according to the distribution function under that stress, but starting with the accumulated failure part before.
The life T i of a product is assumed to follow two-parameter Laplace BS fatigue life distribution LBS(α, β i ) under a constant stress V i , i = 1, 2 , then its distribution function is Based on Nelson Assumption, we know that is, t2 β2 = t1 β1 , t 1 = β1 β2 t 2 , which means that the working time t 2 of a product under the stress V 2 is equivalent to the working time t 1 = β1 β2 t 2 under the stress V 1 . The general progressive stress accelerated life test V (t) = Kt + V 1 , V 1 > 0 is considered. The life distribution of a product is assumed to follow two-parameter Laplace BS fatigue life distribution LBS(α, β 1 ) under the stress V 1 , and the scale parameter β 1 satisfies the inverse power law model According to W. Nelson [12], the working time t under the stress V (t) = Kt + V 1 , V 1 > 0 is equivalent to working under the constant stress V 1 for the following time Then the life distribution of the product under progressive stress In particular, for V 1 = 0, the life distribution of the product under progressive stress V (t) = Kt is Since the parameter c > 0 is required, the parameter m > 0.5 is also required.
Therefore, if the product whose life follows two-parameter Laplace BS fatigue life distribution is subjected to accelerated life test under progressive stress V (t) = Kt, then its failure distribution is

Image Characteristics and Basic Properties of Three-parameter Generalized Laplace BS Fatigue Life Distribution
Based on the conclusions of previous section, three-parameter generalized Laplace Birnbaum-Saunders fatigue life distribution is introduced. If non-negative random variable X follows three-parameter generalized Laplace BS fatigue life distribution, denoted as X ∼ GLBS(α, β, m), then its distribution function F X (x) and density function f X (x) are respectively where α > 0 is called the first shape parameter, m > 0 is called the second shape parameter, β > 0 is called a scale parameter, and Φ(x) =  It is obvious that the expressions of failure rate function λ(x) and average failure is the reliability function and ϕ(x) = 1 2 e −|x| . The true values of parameters are taken as (α = 0.5, 1, 1.5, β = 1, m = 0.5 ) and ( α = 2, β = 1, m = 0.5, 1, 1.5, 2). The images of failure rate functions λ(x) for GLBS(α, β, m) are shown in Figure 5 and Figure 6. The images of average failure rate functionsλ(x) for GLBS(α, β, m) are shown in Figure 7 and Figure 8. It can be seen that the graphs of λ(x),λ(x) are diversified.  X. Zhu and N. Balakrishnan [28] studied basic properties of two-parameter Laplace BS fatigue life distribution. For example, if X ∼ LBS(α, β) , then cX ∼ LBS(α, cβ) and X −1 ∼ LBS(α, β −1 ). Then the properties of three-parameter GLB-S fatigue life distribution are given as follows. Proof. It is denoted as Y = cX , then for y > 0 , we have That is, cX ∼ GLBS(α, cβ, m).
Then we know In particular, we have Property 3.5. If X ∼ GLBS(α, β, m) and it is denoted as Y = lnX , then we have Proof. The parameters are denoted as µ = ln β and σ = 1 m , then for y < µ, we have and for y ≥ µ, we have It is denoted as Z = Y −µ σ , then the distribution function F Z (z) and density function f Z (z) of Z are respectively If k is odd, we have E(Z k ) = 0. If k is even, we have

Point estimate of the scale parameter β
Due to F X (β) = 0.5, the scale parameter β can be estimated by the sample median. When n is odd, it isβ 1 = X (n+1)/2 . When n is even, it isβ 1 = X (n/2) +X (n/2)+1 2 . Xu et al. [26] suggested using geometric mean of samples as the point estimate of scale parameter β for two-parameter BS fatigue life distribution BS(α, β) . Here for three-parameter GLBS distribution GLBS(α, β, m) , it is denoted as Y = ln X, Y i = ln X i , i = 1, 2, · · · , n andȲ = 1 i . According to Property 3.5, the moment estimate of parameter µ isμ =Ȳ , and then the estimate

Lemma 4.1 ( [25]
). Let X 1 , X 2 , · · · , X n be a simple random sample from the population X of the sample size n , and it is denoted as E(X) = µ, D(X) = σ 2 < +∞. The fourth order central moment v 4 = E(X − EX) 4 of the population X is finite. If the fourth derivative of the function is existent and bounded, then we have Proof.
Let the function be h(x) = e x , and any order derivative of h(x) is e x . By Lemma 4.1, it is known that It is easy to see lim n→+∞ E(β 2 ) = β and lim n→+∞ D(β 2 ) = 0. Thus,β 2 is the approximate unbiased estimate and consistent estimate of β. The true values of parameters are taken as α = 0.5, 1, 1.5, β = 1, m = 0.5, 0.8, 1, 1.2, and the sample size is taken as n = 20(10)50. The means and mean square errors of two point estimatesβ 1 ,β 2 of the scale parameter are calculated by Monte Carlo simulations for 1000 times. The results are shown in Table 1, and it can be seen that (1) the means and mean square errors ofβ 1 ,β 2 are little different, which means that the point estimates of β are very close; (2)β 1 is better for the smaller true value of α , whileβ 2 is better for the larger true value of α ; (3) Considering that the sample quantile is usually highly fluctuant for the smaller sample size n , we recommend usingβ 2 .

Point estimates of the shape parameters α, m
By using Property 3.5, three moment equations can be established as follow.
They can be simplified to (4.1) Then the point estimateα of parameter α can be obtained by solving the equation (4.1). Furthermore, the point estimate of parameter m iŝ m = 1 If the sample data meets the condition 1 < s < 6, then the equation (4.1) has one positive real root at least.
Proof. Let the function be g(α) = +∞ 0 It is denoted as λ = α 2 , then we have Let the function be G(λ) = +∞ 0 Therefore, if the sample data meets the condition 1 < s < 6, then the equation (4.1) has one positive real root at least. The graph of the function g(α) is shown in Figure 9, and it can be seen that g(α) is a strictly monotone decreasing function.
The true values of parameters are taken as α = 0.5, 1, 1.5, β = 0.5, 1, 1.5, m = 0.5, 1, 1.5, and the sample size is taken as n = 50. The number of times that the sample data meets the condition 1 < s < 6 (called Condition A) is calculated by Monte Carlo simulations for 1000 times. The results are shown in Table 2, and it can be seen that (1) For the fixed parameters α, β , the number of times is gradually decreasing with the increase of m ; (2) For the fixed parameters α, m , the number of times is increasing at first and then decreasing with the increase of β ; (3) For the fixed parameters β, m , the number of times is gradually increasing with the increase of α ; (4) When β is around one, the number of times that Condition A is satisfied will be relatively higher.  From above analysis, it can be seen that Condition A is not always satisfied. At the same time, considering that the number of times that Condition A is satisfy is relatively higher for β = 1, the right side of equation (4.1) can be approximately modified as follow.
It is denoted as Y i = ln Xî β2 , i = 1, 2, · · · , n and Y = 1 . The value of s is more likely to be between one and six, and 1 < s < 6 is called Condition B. Thus, equation (4.1) becomes to If Condition A is satisfied, then the estimate of α can be solved by equation (4.1). If Condition A is not satisfied and Condition B is satisfied, then the estimate of α can be solved by equation (4.2).
Example 4.1. Air pollution index means that several routinely monitored air pollution concentrations are simplified into a single conceptual index form, and indicate the degree of air pollution and the state of air quality by classification, which is suitable to represent the status and variation trend of short-term air quality in the city. The monitored pollutants include PM10, ozone, nitrogen dioxide, sulfur dioxide and so on. Filidor Vilca [22] 20, 23, 21, 24, 44, 21, 28, 9, 13, 46, 18, 13, 24, 16, 23, 36, 7, 14, 30, 14, 18, 20, 11, 135, 80, 28, 73, 13. First, the K-S test method is used to conduct the fitting test of distribution, and the null hypothesis H 0 is that the frequency distribution conforms to the theoretical distribution. The test statistic is D n = max |F n (x) − F 0 (x)| . When the actual observed value satisfies D > D n , H 0 is rejected, otherwise H 0 is accepted. That is, when √ nD n < C is satisfied, H 0 is accepted, otherwise H 0 is rejected, where C is the tail quantile of the K-S test.
At the significance level 0.1, we obtain C = 1.224 , D n = 0.0911 and √ nD n = 0.9812 < C . Therefore, according to the K-S test principle, it can be considered that the data in this example follows three-parameter GLBS distribution GLBS(α, β, m).
Then, the method in this paper is used to obtain that two point estimates of parameter β are respectivelyβ 1 = 31.5 andβ 2 = 30.5241 . Due to s = −4.96544 and s = 3.84151 , Equation (4.2) is used to solve the estimatesα = 0.7493 , m = 0.5335.

Point Estimates of parameters in Progressive Stress
Accelerated Life Test V (t) = Kt for Two-parameter LBS Fatigue Life Distribution under Inverse Power Law Model According to Section 2, if the life of a product follows two-parameter Laplace BS fatigue life distribution and the accelerated life test is carried out under progressive stress V (t) = Kt , then the failure distribution of this product is There are n products subjected to progressive stress accelerated life test, and the life of each product follows two-parameter Laplace BS fatigue life distribution. The test will continue until all products fail. Then the failure time of these products are denoted as t 1 , t 2 , · · · , t n , and the order failure time are denoted as t (1) , t (2) , · · · , t (n) . The scale parameter β satisfies the inverse power law model. Therefore, it is easy to obtain the point estimatesα,β 2 ,m of parameters α, β, m by using the method proposed in Section 4, and then the point estimates of parameters c, d areĉ = 2m − 1,d =ĉ + 1 Kĉβĉ