A new class of fatigue life distributions

# A new class of fatigue life distributions

Marcelo Bourguignon1 Universidade Federal de Pernambuco
Departamento de Estatística, Cidade Universitária, 50740-540 Recife, PE, Brazil
Rodrigo B. Silva2 Gauss M. Cordeiro3
11Corresponding author. Email: m.p.bourguignon@gmail.com
22Email: rodrigobs29@gmail.com
33Email: gauss@de.ufpe.br
###### Abstract

In this paper, we introduce the Birnbaum-Saunders power series class of distributions which is obtained by compounding Birnbaum-Saunders and power series distributions. The new class of distributions has as a particular case the two-parameter Birnbaum-Saunders distribution. The hazard rate function of the proposed class can be increasing and upside-down bathtub shaped. We provide important mathematical properties such as moments, order statistics, estimation of the parameters and inference for large sample. Three special cases of the new class are investigated with some details. We illustrate the usefulness of the new distributions by means of two applications to real data sets.

Keywords: Birnbaum-Saunders distribution; Maximum likelihood estimation; Order statistic; Power series distribution.

journal:

## 1 Introduction

Birnbaum and Saunders (1969) introduced a two-parameter distribution in order to model the failure time distribution for fatigue failure caused under cyclic loading, called the Birnbaum-Saunders distribution (). Since then, this distribution has been extensively used for modelling failure times of fatiguing materials in several fields such as environmental and medical sciences, engineering, demography, biological studies, actuarial, economics, finance and insurance. A random variable following the distribution with shape parameter and scale parameter is defined by

 T=β⎧⎨⎩αZ2+[(αZ2)2+1]1/2⎫⎬⎭2,

where is a standard normal random variable. The cumulative distribution function (cdf) of has the form

 FBS(t;α,β)=Φ(υ),t>0, (1)

where and is the standard normal cumulative function. The corresponding probability density function (pdf) is

 fBS(t;α,β)=κ(α,β)t−3/2(t+β)exp[−τ(t/β)2α2],t>0, (2)

where and . The parameter is the median of the distribution, i.e. . For any . Additionally, the reciprocal property holds, i.e., , see Saunders (1974). The expected value and variance are and , respectively.

The distribution has a close relation to the normal distribution and has been widely studied in the statistical literature. For example, Ng et al. (2006) proposed a simple bias-reduction method to reduce the biases of the maximum likelihood estimators (MLEs) for the model. They also discussed a Monte Carlo EM-algorithm for computing theses estimators. Kundu et al. (2008) investigated the shape of the hazard function.

In the last few years, several classes of distributions were proposed by compounding some useful lifetime and power series distributions. Chahkandi and Ganjali (2009) studied the exponential power series () family of distributions, which contains as special cases the exponential geometric (), exponential Poisson () and exponential logarithmic () distributions, introduced by Adamidis and Loukas (1998), Kus (2007) and Tahmasbi and Rezaei (2008), respectively. Morais and Barreto-Souza (2011) defined the Weibull power series () class of distributions which includes as sub-models the distributions. The distributions can have an increasing, decreasing and upside down bathtub failure rate function. The generalized exponential power series () distributions were proposed by Mahmoudi and Jafari (2012) in a similar way to that found in Morais and Barreto-Souza (2011). In a very recent paper, Silva  emphet al. (2012) introduced the extended Weibull power series () class of distributions which includes as special models the and distributions.

In a similar manner, we introduce the Birnbaum-Saunders power series () class of distributions, which is defined by compounding the Birnbaum-Saunders and power series distributions. The main aim of this paper is to introduce a new class of distributions, which extends the Birnbaum-Saunders distribution, with the hope that the new distribution may have a ’better fit’ in certain practical situations. Additionally, we will provide a comprehensive account of the mathematical properties of the proposed new class of distributions. The hazard rate function of the proposed class can be increasing and upside-down bathtub shaped. We are motivated to introduce the distributions because of the wide usage of (1) and the fact that the current generalization provides means of its continuous extension to still more complex situations.

This paper is unfolds as follows. In Section 2, a new class of distributions is obtained by mixing the and zero truncated power series distributions. The mixing procedure was previously carried out by Morais and Barreto-Souza (2011) and Silva et al. (2012). In Section 3, some properties of the new class are discussed. In Section 4, the estimation of the model parameters is performed by the method of maximum likelihood. In Section 5, we introduce and study three special cases of the class of distributions. In Section 6, two illustrative applications based on real data are provided. Finally, concluding remarks are presented in Section 7.

## 2 The new class

Let be independent and identically distributed (iid) random variables with cdf (1) and pdf (2). We assume that the random variable has a zero truncated power series distribution with probability mass function

 pn=P(N=n)=anθnC(θ),n=1,2,…, (3)

where depends only on , and is such that is finite. Table 1, reported in Morais and Barreto-Souza (2011), summarizes some power series distributions (truncated at zero) defined according to (3) such as the Poisson, logarithmic, geometric and binomial distributions.

By assuming that the random variables and are independent, we define . Then,

 P(X≤x,N=n)=anθnC(θ)[1−(1−Φ(υ))n],x>0,n≥1.

The class of distributions is defined by the marginal cdf of :

 FBSPS(x;θ,α,β)=1−C[θ(1−Φ(υ))]C(θ),x>0. (4)

Hereafter, the random variable following (4) with parameters and is denoted by . The pdf to (4) is

 fBSPS(x;θ,α,β)=θfBS(x)C′[θ(1−Φ(υ))]C(θ),x>0, (5)

where is given in (2).

###### Proposition 1.

The classical distribution with parameters and is a limiting special case of the proposed model when .

###### Proof.

For , we have

 limθ→0+FBSPS(x) =1−limθ→0+∞∑n=1an[θ(1−Φ(υ))]n∞∑n=1anθn =1−limθ→0+1−Φ(υ)+a−11∞∑n=2anθn−1(1−Φ(υ))n1+a−11∞∑n=2anθn−1=Φ(υ).

We now derive two useful expansions for the density function (5). We have and then

 fBSPS(t;θ,α,β) = fBS(t;α,β)∞∑n=1nanθnC(θ)[1−Φ(υ)]n−1. (6)

By using the binomial theorem in equation (6), we can write

 fBSPS(x;θ,α,β) = fBS(x;α,β)∞∑n=1n−1∑k=0ωn,kΦ(υ)k, (7)

where .

###### Proposition 2.

The density function of can be expressed as an infinite linear combination of densities of minimum order statistics of .

###### Proof.

We know that . Therefore,

 fBSPS(x;θ,α,β)=∞∑n=1pnfT(1)(x;α,β,n),x>0,

where is the density function of , for fixed , given by

 fT(1)(t;α,β,n)=nfBS(t)[1−Φ(v)]n−1,t>0.

###### Remark 1.

Let , then the cdf and pdf of are

 FY(y;θ,α,β)=C(θΦ(v))C(θ),y>0 (8)

and

 fY(y;θ,α,β)=θfBS(y)C′(θΦ(v))C(θ),y>0.

Note that

 fY(y;θ,α,β)=∞∑n=1pnfT(n)(y;α,β,n),y>0,

where is the density function of , for fixed , given by

 fT(n)(t;α,β,n)=nfBS(t)Φ(v)n−1,t>0.

The distribution with cdf (8) is called the complementary Birnbaum-Saunders power series () distribution. This class of distributions is a suitable model in a complementary risk problem based in the presence of latent risks which arise in several areas such as public health, actuarial science, biomedical studies, demography and industrial reliability (Basu and Klein, 1982). However, in this work, we do not focus on this alternative class of distributions.

The survival function becomes

 SBSPS(x;θ,α,β)=C[θ(1−Φ(υ))]C(θ),x>0

and the corresponding hazard rate function reduces to

 hBSPS(x;θ,α,β)=θfBS(x)C′[θ(1−Φ(υ))]C[θ(1−Φ(υ))],x>0.

## 3 Quantiles, order statistics and moments

The quantile function of , say , is given by

 t=Q(u)=β⎧⎨⎩αδ2+[(αδ2)2+1]1/2⎫⎬⎭2,

where , is the inverse cumulative function of the standard normal distribution, is the inverse function of and is a uniform random number on the unit interval .

Let be a random sample with density function (5) and define the th order statistic by . The density function of , say , is given by

 fi:m(x;θ,α,β) = 1B(i,m−i+1)fBSPS(x;θ,α,β)FBSPS(x;θ,α,β)i−1[1−FBSPS(x;θ,α,β)]m−i = fBS(x;α,β)B(i,m−i+1)i−1∑k=0(−1)k(i−1k)SBSPS(x;θ,α,β)m+k−i,

where is the survival function.

###### Proposition 3.

The order statistics can be expressed as a linear combination of exponentiated Birnbaum-Saunders () distributions with parameters and .

In fact,

 fi:m(x;θ,α,β)=i−1∑k=0∞∑j=0m+k−i+j∑s=0ωk,j,shs+1(x,α,β),

where

 ωk,j,s=(−1)k+s(i−1k)(m+k+j−is)dj(s+1)B(i,m−i+1)C(θ)m+k−i,
 d0=am+k−i1,dr=1rd0r∑j=1[j(m+k−i)−r+j]aj+1dr−j

and denotes the density function with parameters and .

Many of the important characteristics and features of a distribution are obtained through the moments. The ordinary moments of can be derived from the probability weighted moments (PWMs) (Greenwood et al., 1979) of the distribution formally defined for and non-negative integers by

 τp,r=κ(α,β)∫∞0tp−3/2(t+β)exp[−τ(t/β)2α2]Φr(υ)dt. (9)

The integral (9) can be easily computed numerically in software such as MAPLE, MATHEMATICA, Ox and R, see Cordeiro and Lemonte (2011). An alternative expression for can be derived from Cordeiro and Lemonte (2011) as

 τp,r=βp2rr∑j=0(rj)∞∑k1,…,kj=0A(k1,...,kj)2sj+j∑m=0(−1)m(2sj+jm)×β−(2sj+j−2m)/2I(p+(2sj+j−2m)/2,α).

Here, , , and

 I(p,α)=Kp+1/2(α−2)+Kp−1/2(α−2)2K1/2(α−2),

where the function denotes the modified Bessel function of the third kind with representing its order and the argument. Its integral representation is .

Thus, the th moment of can be obtained from equations (7) and (9) as

 E(Xs) = ∞∑n=1n−1∑k=0ωn,kτs,k.

Expressions for the th moments of the order statistics with cumulative function (4) can be obtained using a result due to Barakat and Abdelkader (2004). We have

 E(Xsi:m) =sm∑j=m−i+1(−1)j−m+i−1(j−1m−i)(mj)∫∞0xs−1SBSPS(x,θ,α,β)jdx =sm∑j=m−i+1(−1)j−m+i−1C(θ)j(j−1m−i)(mj)∫∞0xs−1C[θ(1−Φ(υ))]jdx,

for .

## 4 Maximum likelihood estimation

Here, we discuss maximum likelihood estimation and inference for the distribution. Let be a random sample from and let be the vector of the model parameters. The log-likelihood function for reduces to

 ℓ(θ,α,β) = n{log(θ)−log[C(θ)]+log[κ(α,β)]}−32n∑i=1log(xi)+n∑i=1log(xi+β) − 12α2n∑i=1τ(xiβ)+n∑i=1log{C′[θ(1−Φ(vi))]}.

The components of the score vector are given by

 Uθ(Θ) =nθ−nC′(θ)C(θ)+n∑i=11−ϕ(υi)C′[θ(1−Φ(υi))]∂C′[θ(1−Φ(υi))]∂θ, Uα(Θ) =nα(1+2α2)+1α3n∑i=1τ(xiβ)+θα3n∑i=1ϕ(υi)ρ(xi/β)C′[θ(1−Φ(υi))]∂C′[θ(1−Φ(υi))]∂α and Uβ(Θ) =−n2β+12α2βn∑i=1(xiβ+βxi)−θ2αβn∑i=1ϕ(υi)τ(√xi/β)C′[θ(1−Φ(υi))]∂C′[θ(1−Φ(υi))]∂β,

where is the standard normal density and .

Setting these equations to zero, , and solving them simultaneously yields the MLE of . These equations cannot be solved analytically and statistical software can be used to compute them numerically using iterative techniques such as the Newton-Raphson algorithm.

Often with lifetime data and reliability studies, one encounters censoring. A very simple random censoring mechanism that is often realistic is one in which each individual is assumed to have a lifetime and a censoring time , where and are independent random variables. Suppose that the data consist of independent observations and is such that if is a time to event and if it is right censored for . The censored likelihood for the model parameters is

 L(θ,α,β)∝n∏i=1[fBSPS(xi;θ,α,β)]δi[SBSPS(xi;θ,α,β)]1−δi.

In order to perform interval estimation and hypothesis tests on the model parameters , and , the normal approximation for the MLE can be applied. Note that under certain conditions for the parameters, the asymptotic distribution of is multivariate normal , where is the observed information matrix and . Note also that an () asymptotic confidence interval for the th parameter in is specified by

 ACIi=(ˆθi−z1−γ/2√ˆJθi,θi,^θi+z1−γ/2√ˆJθi,θi),

where stands for the th diagonal element of the inverse of the observed information matrix estimated at , i.e., , for , and is the standard normal quantile.

For testing the goodness-of-fit of the distribution and for comparing it with some sub-models, we can use the likelihood ratio (LR) statistics. The LR statistic for testing the null hypothesis : versus the alternative hypothesis : is given by , where and are the MLEs under the alternative and null hypotheses, respectively. Under the null hypothesis, , , where is the dimension of the subset of interest.

## 5 Special cases

In this section, we investigate some special models of the class of distributions. We offer some expressions for the moments and moments of the order statistics. We illustrate the flexibility of these distributions and provide plots of the density and hazard rate functions for selected parameter values.

### 5.1 Birnbaum-Saunders geometric distribution

The Birnbaum-Saunders geometric () distribution is defined by the cdf (4) with leading to

 FBSG(x;θ,α,β)=1−(1−θ)[1−Φ(v)]1−θ[1−Φ(v)],x>0, (10)

where .

###### Proposition 4.

The distribution of the form (10) is geometric minimum stable.

###### Proof.

The proof follows easily of the arguments given by Marshall and Olkin (1997, p. 647). We omit the details. ∎

The associated density and hazard rate functions are

 fBSG(x;θ,α,β)=(1−θ)fBS(x)[1−θ(1−Φ(v))]−2

and

 hBSG(x;θ,α,β)=fBS(x)[1−Φ(v)]{1−θ[1−Φ(v)]}=hBS(x)1−θ[1−Φ(v)], (11)

for , respectively. From equation (11), we note that is decreasing in . It can also be shown that

 limx→0hBSG(x;θ,α,β)=0andlimx→∞hBSG(x;θ,α,β)=12α2β.
###### Corollary 1.

Let and . Then for all .

###### Proof.

It is straightforward. ∎

Cancho et al. (2012) proposed and studied the geometric Birnbaum-Saunders regression with cure rate. However, they considered instead of . In Figure 1, we plot the density and hazard rate functions of the distribution for selected parameter values. We can verify that this distribution has an upside-down bathtub or an increasing failure rate function depending on the values of its parameters.

Expressions for the density function and the moments of the order statistics from a random sample of the distribution are given by

 fi:m(x)=fBS(x;α,β)B(i,m−i+1)i−1∑k=0(−1)k(i−1k)[(1−θ)Φ(−υ)1−θΦ(−υ)]m+k−i

and

 E(Xsi:m) =sm∑j=m−i+1(−1)j−m+i−1θj(1−θ)−j(j−1m−i)(mj)∫∞0xs−1[1−θΦ(−υ)]jdx,

which are easily determined numerically.

Let be a random sample of size from . The log-likelihood function for the vector of parameters can be expressed as

 ℓ(Θ) = n{log(1−θ)+log[κ(α,β)]}−32n∑i=1log(xi)+n∑i=1log(xi+β) − 12α2n∑i=1τ(xiβ)−2n∑i=1log[1−θ(1−Φ(vi))].

### 5.2 Birnbaum-Saunders Poisson distribution

The Birnbaum-Saunders Poisson () distribution follows by taking in (4), which yields

 FBSP(t;θ,α,β)=1−exp{θ(1−Φ(v))}−1exp(θ)−1.

The density and hazard rate functions of the distribution are

 fBSP(t;θ,α,β)=θfBS(t)exp[θ(1−Φ(v))]exp(θ)−1

and

 hBSP(t;θ,α,β)=θfBS(t)exp{θ(1−Φ(v))}exp{θ(1−Φ(v))}−1.
###### Remark 2.

The limit of the hazard rate function as is .

In Figure 2, we plot the density and hazard rate functions of the distribution for selected parameter values. This distribution can have an upside-down bathtub or an increasing failure rate function depending on the values of its parameters.

The expressions for the density function and the moments of the order statistics from a random sample of the distribution are

 fi:m(x)=fBS(x;α,β)B(i,m−i+1)i−1∑k=0(−1)k(i−1k){exp[θΦ(−υ)]−1eθ−1}m+k−i

and

 E(Xsi:m) =sm∑j=m−i+1(−1)j−m+i−1(exp(θ)−1)j(j−1m−i)(mj)∫∞0xs−1{exp[θΦ(−υ)]−1}jdx.

Let be a random sample of size from . The log-likelihood function for the vector of parameters can be expressed as

 ℓ(Θ) = n{log(θ)−log[exp(θ)−1]+log[κ(α,β)]}−32n∑i=1log(xi)+n∑i=1log(xi+β) − 12α2n∑i=1τ(xiβ)+θn∑i=1[1−Φ(vi)].

### 5.3 Birnbaum-Saunders logarithmic distribution

The cdf of the Birnbaum-Saunders logarithmic () distribution is defined by (4) with , which corresponds to the logarithmic distribution. We obtain

 FBSL(x;θ,α,β)=1−log[1−θ(1−Φ(υ))]log(1−θ)

The associated density and hazard rate functions are

 fBSL(x;θ,α,β)=−θfBS(x)log(1−θ){1−θ[1−Φ(υ)]}

and

 hBSL(x;θ,α,β)=−θfBS(x)log{1−θ[1−Φ(υ)]}{1−θ[1−Φ(υ)]}
###### Remark 3.

The limit of the hazard rate function as is .

In Figure 3, we plot the density and hazard rate functions of the distribution for selected parameter values. We can verify that this distribution can have an upside-down bathtub or an increasing failure rate function depending on the values of its parameters.

The expressions for the density function and the moments of the order statistics from a random sample of the distribution are

 fi:m(x)=fBS(x;α,β)B(i,m−i+1)i−1∑k=0(−1)k(i−1k){log[1−θΦ(−υ)]log(1−θ)}m+k−i

and

 E(Xsi:m) =sm∑j=m−i+1(−1)i−m−1[log(1−θ)]j(j−1m−i)(mj)∫∞0xs−1{log[11−θΦ(−υ)]}jdx.

Let be a random sample of size from . The log-likelihood function for the vector of parameters can be expressed as

 ℓ(Θ) = n{log(θ)−log[−log(1−θ)]+log[κ(α,β)]}−32n∑i=1log(xi)+n∑i=1log(xi+β) − 12α2n∑i=1τ(xiβ)−n∑i=1log{1−θ[1−Φ(υi)]}.

## 6 Applications

In this section, we compare the fits of some special models of the class by means of two real data sets to show the potentiality of the new class. In order to estimate the parameters of these special models, we adopt the maximum likelihood method (as discussed in Section 4) and all the computations were done using the subroutine NLMixed of the SAS software. Obviously, due to the genesis of the distribution, the fatigue processes are by excellence ideally modeled by this distribution. Thus, the use of the and distributions for fitting these two data sets is well justified.

First, we consider a data set from Murthy et al. (2004) consisting of the failure times of 20 mechanical components. The data are: 0.067, 0.068, 0.076, 0.081, 0.084, 0.085, 0.085, 0.086, 0.089, 0.098, 0.098, 0.114, 0.114, 0.115, 0.121, 0.125, 0.131, 0.149, 0.160, 0.485. The MLEs of the parameters (with corresponding standard errors in parentheses), the value of , the Kolmogorov-Smirnov statistic, Akaike information criterion (AIC) and Bayesian information criterion (BIC) for the , , and models are listed in Table 2. Since the values of the AIC and BIC are smaller for the distribution compared with those values of the other models, the new distribution seems to be a very competitive model for these data.