A New Extended Birnbaum–Saunders Model: Properties, Regression and Applications

Abstract

We propose an extended fatigue lifetime model called the odd log-logistic Birnbaum–Saunders–Poisson distribution, which includes as special cases the Birnbaum–Saunders and odd log-logistic Birnbaum–Saunders distributions. We obtain some structural properties of the new distribution. We define a new extended regression model based on the logarithm of the odd log-logistic Birnbaum–Saunders–Poisson random variable. For censored data, we estimate the parameters of the regression model using maximum likelihood. We investigate the accuracy of the maximum likelihood estimates using Monte Carlo simulations. The importance of the proposed models, when compared to existing models, is illustrated by means of two real data sets.

1. Introduction

Motivated by problems of vibration in commercial aircraft that caused fatigue in the materials, Birnbaum and Saunders [1] proposed the two-parameter Birnbaum–Saunders (BS) model, also known as the fatigue life distribution, with shape parameter $a>0$ and scale parameter $b>0$, i.e., $BS(a,b)$. This distribution can be used to model lifetime data and it is widely applicable to model failure times of fatiguing materials. A random variable W having the $BS(a,b)$ distribution is defined by

$$\begin{array}{c}\hfill W=b{\left[\frac{aZ}{2}+{\left\{{\left(\frac{aZ}{2}\right)}^2+1\right\}}^{1/2}\right]}^2,\end{array}$$

(1)

where Z is a standard normal random variable. Its cumulative distribution function (cdf) is given by

$$\begin{array}{c}\hfill G_{a,b}\left(x\right)=\Phi \left(\nu \right),x>0,\end{array}$$

(2)

where $\nu =a^{-1}\rho (x/b)$, $\rho \left(q\right)=q^{1/2}-q^{-1/2}$ and $\Phi (·)$ is the standard normal cumulative function. The parameter b is the median of the distribution, i.e., $G_{a,b}\left(b\right)=\Phi \left(0\right)=1/2$. For any $k>0$, $kW\sim BS(a,kb)$. Some distributions have been proposed to extend the fatigue lifetime BS model. For example, Cordeiro et al. [2] defined the McDonald–Birnbaum–Saunders distribution, Ortega et al. [3] introduced the $\beta$-Birnbaum–Saunders distribution and Hashimoto et al. [4] proposed the Poisson Birnbaum–Saunders model with long-term survivors. Further, Pescim et al. [5] studied the Kummer beta generalized Birnbaum–Saunders distribution, Cordeiro et al. [6] introduced the negative binomial Birnbaum–Saunders model, Cordeiro et al. [7] proposed the gamma-Birnbaum–Saunders distribution, Ortega et al. [8] constructed the odd Birnbaum–Saunders regression model and Ortega et al. [9] studied a four-parameter odd Birnbaum–Saunders geometric distribution.

The probability density function (pdf) corresponding to Equation (2) is given by

$$\begin{array}{c}\hfill g_{a,b}\left(x\right)=\kappa (a,b)x^{-3/2}(x+b)exp\left[-\frac{\tau (x/b)}{2a^2}\right],x>0,\end{array}$$

(3)

where $\kappa (a,b)={\left(2a{e}^{a^2}\sqrt{2\pi b}\right)}^{-1}$ and $\tau \left(q\right)=q+q^{-1}$. The fractional moments of Equation (3) are determined as $E\left(W^p\right)=b^{p}I(p,a)$ ([10]), where

$$\begin{array}{c}\hfill I(p,a)=\frac{{\theta }_{p+1/2}\left(a^{-2}\right)+{\theta }_{p-1/2}\left(a^{-2}\right)}{2{\theta }_{1/2}\left(a^{-2}\right)},\end{array}$$

(4)

and ${\theta }_{\nu }\left(z\right)=0.5{\int }_{-\infty }^{\infty }exp\{-zcosh\left(t\right)-\nu t\}dt$ denotes the modified Bessel function of the third kind with $\nu$ representing its order and z the argument. A discussion of this function can be found in [11].

Recently, some attempts have been made to define new classes of lifetime distributions that provide great flexibility in modeling skewed data in practice. For any baseline cdf G(x) ($x\in \mathbb{R}$), Gleaton and Lynch [12] defined the cdf $F\left(x\right)$ of the odd log-logistic-G (“OLL-G” for short) class by

$$\begin{array}{c}\hfill F(x;\mathbf{\tau },\alpha )={\int }_0^{\frac{G(x;\mathbf{\tau })}{1-G(x;\mathbf{\tau })}}\frac{\alpha t^{\alpha -1}}{{(1+t^{\alpha })}^2}dt=\frac{G{(x;\mathbf{\tau })}^{\alpha }}{G{(x;\mathbf{\tau })}^{\alpha }+\overline{G}{(x;\mathbf{\tau })}^{\alpha }}.\end{array}$$

(5)

The pdf corresponding to Equation (5) is given by

$$\begin{array}{c}\hfill f(x;\mathbf{\tau },\alpha )=\frac{\alpha g(x;\mathbf{\tau })G{(x;\mathbf{\tau })}^{\alpha -1}\overline{G}{(x;\mathbf{\tau })}^{\alpha -1}}{{\left[G{(x;\mathbf{\tau })}^{\alpha }+\overline{G}{(x;\mathbf{\tau })}^{\alpha }\right]}^2}.\end{array}$$

(6)

In addition, some authors have investigated the OLL-G class in different settings. For instance, da Cruz et al. [13] proposed the log-odd log-logistic Weibull regression model in survival analysis, Braga et al. [14] studied the odd log-logistic Student t distribution, da Cruz et al. [15] considered the bivariate odd-log-logistic-Weibull regression model for oral health-related quality of life, Cordeiro et al. [16] proposed the generalized odd log-logistic family and, recently, Prataviera et al. [17] introduced the heteroscedastic odd log-logistic generalized gamma regression model for censored data. We use a similar method to define a new distribution by compounding the Birnbaum–Saunders and Poisson distributions whose procedure is described below.

The method of compounding discrete and continuous distributions started fifteen years ago to model complex situations under series and parallel structures. Many compounded classes of distributions can be constructed by choosing a positive integer discrete random variable N with probability mass function (pmf) $P(N=n)$. Let $Z_1,\cdots ,Z_N$ be independent identically distributed (iid) random variables with common cdf given in Equation (5). The minimum and maximum random variables $min\{Z_1,\cdots ,Z_N\}$ and $max\{Z_1,\cdots ,Z_N\}$ can generate several models that arise in series and parallel systems with identical components and have many industrial and biological applications. For example, the time to the failure of a series system with an unknown number of protected components or the cancer recurrence of an individual after undergoing a certain treatment can be modeled by the generated distribution of the minimum random variable. In a dual mechanism, the time to the failure of a parallel system with an unknown number of protected components or a disease manifestation if it occurs only after an unknown number of factors have been active can be modeled by the generated distribution of the maximum random variable.

The odd log-logistic-G- Poisson (OLLG-P for short) family is defined by the maximum random variable $M_N=max\{Z_1,\cdots ,Z_N\}$, under the assumption that the $Z_i$s are iid OLLG random variables independent of N, by taking N as a truncated Poisson random variable with pmf

$$\begin{array}{c}\hfill P(N=n)=\frac{{\beta }^n}{({\mathrm{e}}^{\beta }-1)n!},n=1,2,\dots ,\beta >0.\end{array}$$

(7)

Then, the cdf of the OLLG-P family can be expressed as

$$\begin{array}{c}\hfill F(x;\mathbf{\tau },\alpha ,\beta )=\sum _{n=1}^{\infty }P(M_N\le x|N=n)P(N=n)=\frac{1}{({\mathrm{e}}^{\beta }-1)}\left[exp\left\{\frac{\beta G{(x;\mathbf{\tau })}^{\alpha }}{G{(x;\mathbf{\tau })}^{\alpha }+\overline{G}{(x;\mathbf{\tau })}^{\alpha }}\right\}-1\right].\end{array}$$

(8)

The pdf of the OLLG-P family reduces to

$$\begin{array}{c}\hfill f(x;\mathbf{\tau },\alpha ,\beta )=\frac{\alpha \beta g(x;\mathbf{\tau })G{(x;\mathbf{\tau })}^{\alpha -1}\overline{G}{(x;\mathbf{\tau })}^{\alpha -1}exp\left\{\frac{\beta G{(x;\mathbf{\tau })}^{\alpha }}{G{(x;\mathbf{\tau })}^{\alpha }+\overline{G}{(x;\mathbf{\tau })}^{\alpha }}\right\}}{({\mathrm{e}}^{\beta }-1){\left[G{(x;\mathbf{\tau })}^{\alpha }+\overline{G}{(x;\mathbf{\tau })}^{\alpha }\right]}^2},\end{array}$$

(9)

where $\alpha >0$; $\beta >0$; $G(x;\mathbf{\tau })$ and $g(x;\mathbf{\tau })$ denote the cdf and pdf of the baseline distribution, respectively; and $\mathbf{\tau }$ is the baseline parameter vector.

In this paper, we propose a new lifetime model named the odd log-logistic Birnbaum–Saunders–Poisson (OLLBSP) distribution from Equations (8) and (9) to extend the BS model. We obtain some of its structural properties and discuss maximum likelihood estimation of the model parameters.

Regression models can be proposed in different forms in survival and reliability analysis. Among them, the location-scale regression model and Cox proportional-hazards model are distinguished and frequently used. We also construct a location-scale regression model based on the OLLBSP distribution, called the log-odd log-logistic Birnbaum–Saunders–Poisson (LOLLBSP) regression model for censored data, as a feasible alternative to the log-BS regression model. We consider a classic analysis for the LOLLBSP regression model. The inferential part is carried out using the asymptotic distribution of the maximum likelihood estimators (MLEs).

The paper is outlined as follows. In Section 2, we define the OLLBSP distribution and provide plots of its density and hazard rate functions. In Section 3, we present the maximum likelihood method to estimate the model parameters and a simulation study. In Section 4, we derive a useful linear representation for the OLLBSP density function and obtain its ordinary moments. In Section 5, we define the LOLLBSP regression model for censored data and estimate the model parameters by maximum likelihood. In Section 6, we prove empirically the potentiality of the new distribution for fatigue life data and the flexibility and relevance of the proposed regression model by means of two applications to real data sets. Finally, Section 7 ends with some concluding remarks.

2. The OLLBSP Distribution

Most generalized BS distributions have been proposed in reliability to provide a better fit to certain data sets than the traditional two- and three-parameter BS models. Providing a new lifetime distribution is always precious for statisticians. The pdf and cdf of the OLLBSP distribution are defined (for $x>0$) by substituting Equations (2) and (3) into Equations (8) and (9)

$$\begin{array}{c}\hfill F\left(x\right)=\frac{1}{({\mathrm{e}}^{\beta }-1)}\left[exp\left\{\frac{\beta \Phi {\left(\nu \right)}^{\alpha }}{\Phi {\left(\nu \right)}^{\alpha }+{\left[1-\Phi \left(\nu \right)\right]}^{\alpha }}\right\}-1\right]\end{array}$$

(10)

and

$$\begin{array}{c}\hfill f\left(x\right)=\frac{\alpha \beta \kappa (a,b)x^{-3/2}(x+b)exp\left[-\frac{\tau (x/b)}{2a^2}\right]\Phi {\left(\nu \right)}^{\alpha -1}{\left[1-\Phi \left(\nu \right)\right]}^{\alpha -1}exp\left\{\frac{\beta \Phi {\left(\nu \right)}^{\alpha }}{\Phi {\left(\nu \right)}^{\alpha }+{\left[1-\Phi \left(\nu \right)\right]}^{\alpha }}\right\}}{({\mathrm{e}}^{\beta }-1){\left[\Phi {\left(\nu \right)}^{\alpha }+{\left[1-\Phi \left(\nu \right)\right]}^{\alpha }\right]}^2},\end{array}$$

(11)

respectively. Evidently, the density function in Equation (11) does not involve any complicated function. Hereafter, we denote by $X\sim$OLLBSP$(\alpha ,\beta ,a,b)$ a random variable having the density in Equation (11).

Figure 1 displays some possible shapes of the density function in Equation (11) for selected parameter values. It is evident that the OLLBSP distribution is much more flexible than the BS distribution presenting one additional shapes, namely bi-modality. Some real situations require bimodal density behavior and a distribution that presents this shape is a good tentative model. Further, there are not many distributions having this shape.

The hazard rate function (hrf) corresponding to Equation (11) is given by

$$\begin{array}{c}\hfill h\left(x\right)=\frac{\alpha \beta \kappa (a,b)x^{-3/2}(x+b)exp\left[-\frac{\tau (x/b)}{2a^2}\right]\Phi {\left(\nu \right)}^{\alpha -1}{\left[1-\Phi \left(\nu \right)\right]}^{\alpha -1}exp\left\{\frac{\beta \Phi {\left(\nu \right)}^{\alpha }}{\Phi {\left(\nu \right)}^{\alpha }+{\left[1-\Phi \left(\nu \right)\right]}^{\alpha }}\right\}}{\left\{exp\left(\beta \right)-exp\left[\frac{\beta \Phi {\left(\nu \right)}^{\alpha }}{\Phi {\left(\nu \right)}^{\alpha }+{\left[1-\Phi \left(\nu \right)\right]}^{\alpha }}\right]\right\}{\left[\Phi {\left(\nu \right)}^{\alpha }+{\left[1-\Phi \left(\nu \right)\right]}^{\alpha }\right]}^2}.\end{array}$$

(12)

Plots of the hrf $h\left(x\right)$ for some parameter values are displayed in Figure 2. We note that these plots illustrate the four classic types of hazard shapes. Further, we note the increasing–decreasing–increasing shape, which is important because there are fewer lifetime distributions that yield complicated shapes similar to this. Further, we require a mixture of two or more lifetime distributions to model this shape. Some data sets present this behavior, usually in regression study, when we have different levels of explanatory variables. Thus, probably, this model will have a better fit to certain data sets than other known models, as we show in the Application Section.

The new distribution is simulated as follows: if U is a uniform random variable $(0,1)$, then

$$\begin{array}{c}\hfill X=Q_{BS}\left(\frac{V^{1/\alpha }}{V^{1/\alpha }+{[1-V]}^{1/\alpha }}\right),\end{array}$$

(13)

has the density in Equation (11), where $Q_{BS}(·)$ is the BS quantile function (qf) and

$$\begin{array}{c}\hfill V=\frac{log\left\{U[exp(\beta )-1]+1\right\}}{\beta }.\end{array}$$

(14)

3. Inference and Estimation

In this section, we discuss the maximum likelihood method for inference and estimation of the model parameters of the OLLBSP distribution. Let $x_1,\cdots ,x_n$ be a random sample of size n from the OLLBSP$(\alpha ,\beta ,a,b)$ distribution. The log-likelihood function for the vector of parameters $\mathbf{\theta }={(\alpha ,\beta ,a,b)}^T$ is given by

$$\begin{array}{ccc}\hfill l\left(\mathbf{\theta }\right)& =& nlog\left[\alpha \beta \kappa (a,b)\right]-nlog({\mathrm{e}}^{\beta }-1)-\frac{3}{2}\sum _{i=1}^{n}log\left(x_i\right)+\sum _{i=1}^{n}log(x_i+b)-\frac{1}{2a^2}\sum _{i=1}^n\tau \left(\frac{x_i}{b}\right)+\hfill \\ & & (\alpha -1)\sum _{i=1}^{n}log[\Phi \left({\nu }_i\right)]+(\alpha -1)\sum _{i=1}^{n}log[1-\Phi \left({\nu }_i\right)]+\beta \sum _{i=1}^n\frac{{\Phi }^{\alpha }\left({\nu }_i\right)}{{\Phi }^{\alpha }\left({\nu }_i\right)+{[1-\Phi \left({\nu }_i\right)]}^{\alpha }}+\hfill \\ & & 2\sum _{i=1}^{n}log\{{\Phi }^{\alpha }\left({\nu }_i\right)+{[1-\Phi \left({\nu }_i\right)]}^{\alpha }\}.\hfill \end{array}$$

(15)

Maximization of Equation (15) can be performed using well-established routines such as nlm or optimize in the R statistical package. Setting these equations to zero, $U_{\alpha }\left(\mathbf{\theta }\right)=0$, $U_{\beta }\left(\mathbf{\theta }\right)=0$, $U_a\left(\mathbf{\theta }\right)=0$, $U_b\left(\mathbf{\theta }\right)=0$ and solving them simultaneously gives the MLE $\widehat{\mathbf{\theta }}$ of $\mathbf{\theta }$. These equations cannot be solved analytically and statistical or mathematical software can be used to evaluate them numerically using iterative techniques such as the Newton–Raphson algorithm.

For interval estimation and tests of hypotheses on the parameters in $\mathbf{\theta }$, we require the $4\times 4$ observed information matrix $\mathbf{J}=\mathbf{J}\left(\mathbf{\theta }\right)=\left\{j_{r,s}\right\}$, whose elements $j_{r,s}$ (for $r,s=\alpha ,\beta ,a,b$) can be evaluated numerically. Under standard regularity conditions, the approximate multivariate normal $N_4(0,\mathbf{J}{(\widehat{\mathbf{\theta }})}^{-1})$ distribution can be used to construct approximate confidence intervals for the model parameters.

The likelihood ratio (LR) statistics are useful for comparing the new distribution with some special models. For example, we may use the LR statistic to check if the fit using the OLLBSP distribution is statistically “superior” to the fit using the BS distribution for a given data set.

We assess the finite sample performance of the MLEs in the OLLBSP distribution by varying the sample size n using the inverse method with Equation (13). The sample sizes are taken as $n\in \{50,100,200,300\}$. We perform a Monte Carlo simulation study with 1000 replications under three setups for parameters of the OLLBSP distribution: $a=2.5$, $b=1.5$, $\alpha =2$ and $\beta =0.75$ (Setup 1); $a=0.75$, $b=3$, $\alpha =1$ and $\beta =2$ (Setup 2); and $a=1.5$, $b=3$, $\alpha =0.75$ and $\beta =1.25$ (Setup 3).

We conclude from the figures in

Table 1. Simulation study.

	Sample Size	Parameter	Mean	Bias	MSE
		a	2.7479	0.2479	0.5496
	$N=50$	b	1.2505	−0.2495	0.2038
		$\alpha$	2.1559	0.1559	0.3803
		$\beta$	1.4741	0.7241	1.4878
		a	2.7228	0.2228	0.5063
	$N=100$	b	1.3640	−0.1360	0.1263
		$\alpha$	2.1356	0.1356	0.3608
Setup 1		$\beta$	1.1648	0.4148	0.7094
		a	2.6960	0.1960	0.4276
	$N=200$	b	1.4476	−0.0524	0.0709
		$\alpha$	2.1111	0.1111	0.3104
		$\beta$	0.9102	0.1602	0.3297
		a	2.6507	0.1507	0.3695
	$N=300$	b	1.4466	−0.0534	0.0551
		$\alpha$	2.1084	0.1084	0.2584
		$\beta$	0.9014	0.1514	0.2706
		a	0.7487	−0.0013	0.2457
	$N=50$	b	3.0076	0.0076	0.4509
		$\alpha$	0.9691	−0.0309	0.5599
		$\beta$	2.1168	0.1168	1.4317
		a	0.7544	0.0044	0.1970
	$N=100$	b	2.9874	−0.0126	0.3606
		$\alpha$	0.9754	−0.0246	0.4350
Setup 2		$\beta$	2.1233	0.1233	1.3254
		a	0.7616	0.0116	0.0920
	$N=200$	b	2.9967	−0.0033	0.3039
		$\alpha$	0.9928	−0.0072	0.1912
		$\beta$	2.1087	0.1087	1.1643
		a	0.7525	0.0025	0.0447
	$N=300$	b	2.9996	−0.0004	0.2216
		$\alpha$	0.9853	−0.0147	0.0876
		$\beta$	2.0678	0.0678	0.8310
		a	1.5922	0.0922	0.7011
	$N=50$	b	3.0035	0.0035	0.5183
		$\alpha$	0.8068	0.0568	0.2529
		$\beta$	1.3356	0.0856	0.4881
		a	1.5724	0.0724	0.4521
	$N=100$	b	3.0179	0.0179	0.2984
		$\alpha$	0.7955	0.0455	0.1701
Setup 3		$\beta$	1.2694	0.0194	0.2892
		a	1.5195	0.0195	0.2394
	$N=200$	b	3.0461	0.0461	0.1839
		$\alpha$	0.7637	0.0137	0.0910
		$\beta$	1.2415	−0.0085	0.1530
		a	1.5229	0.0229	0.1487
	$N=300$	b	3.0199	0.0199	0.1186
		$\alpha$	0.7651	0.0151	0.0573
		$\beta$	1.2527	0.0027	0.1108

that the average estimates of the parameters tend to be closer to the true parameters when n increases. This fact supports that the asymptotic normal distribution provides an adequate approximation to the finite sample distribution of the MLEs. The normal approximation can be oftentimes improved by using bias adjustments to these estimators. Approximations to the their biases in simple models may be determined analytically. Bias correction typically does a very good job for correcting the MLEs. However, it may also increase the mean squared errors (MSEs). Whether bias correction is useful in practice depends basically on the shape of the bias function and on the variance of the MLE.

4. Structural Properties

4.1. Linear Representation

Expansions for Equations (10) and (11) can be derived using the concept of exponentiated distributions to obtain some structural properties of the OLLBSP distribution. Cordeiro et al. [2] defined the exponentiated Birnbaum–Saunders (EBS) distribution with positive parameters $\alpha$, $\beta$ and c, say $Y\sim EBS(a,b,c)$, if its cdf and pdf are given by $H(y;a,b,c)=\Phi {\left(\nu \right)}^c$ and $h(y;a,b,c)=c{g}_{a,b}\left(y\right)\Phi {\left(\nu \right)}^{c-1}$, respectively, where $\Phi \left(\nu \right)$ is defined in Equation (2) and $g_{a,b}\left(y\right)$ is the BS$(a,b)$ pdf. The properties of some exponentiated distributions have been studied by several authors (see, for example, Mudholkar and Srivastava [18] for exponentiated Weibull (EW) and [19] for the exponentiated generalized gamma (EGG) distributions).

By using Taylor expansion, we can write Equation (10) as

$$\begin{array}{c}\hfill F\left(x\right)=\frac{1}{{\mathrm{e}}^{\beta }-1}\sum _{i=1}^{\infty }\frac{{\beta }^i}{i!}\frac{\Phi {\left(\nu \right)}^{i\alpha }}{{\left[\Phi {\left(\nu \right)}^{\alpha }+\overline{\Phi }{\left(\nu \right)}^{\alpha }\right]}^i}.\end{array}$$

(16)

For p real non-integer and $\left|z\right|<1$, the convergent power series

$$\begin{array}{c}\hfill {(1-z)}^p=\sum _{j=0}^{\infty }{(-1)}^j\left(\genfrac{}{}{0pt}{}{p}{j}\right)z^j\end{array}$$

(17)

holds, where the binomial coefficient is defined for any real. By applying this equation twice to $\Phi {\left(\nu \right)}^{i\alpha }={\{1-[1-\Phi \left(\nu \right)]\}}^{i\alpha }$ and interchanging sums, we obtain the convergent power series (for $\alpha >0$)

$$\begin{array}{c}\hfill \Phi {\left(\nu \right)}^{i\alpha }=\sum _{k=0}^{\infty }{a}_k\left(i\alpha \right)\Phi {\left(\nu \right)}^k,\end{array}$$

(18)

where (for $k\ge 0$)

$$a_k\left(i\alpha \right)=\sum _{r=k}^{\infty }{(-1)}^{k+r}\left(\genfrac{}{}{0pt}{}{i\alpha }{r}\right)\left(\genfrac{}{}{0pt}{}{r}{k}\right).$$

By using the binomial expansion and Equation (17), we obtain (for $i\ge 1$)

$$\begin{array}{c}\hfill {\left[\Phi {\left(\nu \right)}^{\alpha }+\overline{\Phi }{\left(\nu \right)}^{\alpha }\right]}^i=\Phi {\left(\nu \right)}^{i\alpha }+\sum _{r=1}^i\sum _{j=0}^{\infty }{d}_{r,j}(i,\alpha )\Phi {\left(\nu \right)}^{(i-r)\alpha +j},\end{array}$$

(19)

where $d_{r,j}(i,\alpha )={(-1)}^j\left(\genfrac{}{}{0pt}{}{i}{r}\right)\left(\genfrac{}{}{0pt}{}{r\alpha }{j}\right)$.

Further, using Equation (18) in the two powers of the right-hand side and interchanging sums, we have

$$\begin{array}{c}\hfill {\left[\Phi {\left(\nu \right)}^{\alpha }+\overline{\Phi }{\left(\nu \right)}^{\alpha }\right]}^i=\sum _{k=0}^{\infty }{b}_k(i,\alpha )\Phi {\left(\nu \right)}^k,\end{array}$$

(20)

where $b_k(i,\alpha )=a_k\left(i\alpha \right)+{\sum }_{r=1}^i{\sum }_{j=0}^{\infty }{d}_{r,j}(i,\alpha )a_k((i-r)\alpha +j)$. Hence, from Equations (18) and (20), we obtain the ratio of two convergent power series

$$\begin{array}{c}\hfill \frac{\Phi {\left(\nu \right)}^{i\alpha }}{{\left[\Phi {\left(\nu \right)}^{\alpha }+\overline{\Phi }{\left(\nu \right)}^{\alpha }\right]}^i}=\frac{{\sum }_{k=0}^{\infty }{a}_k\left(i\alpha \right)\Phi {\left(\nu \right)}^k}{{\sum }_{k=0}^{\infty }{b}_k(i,\alpha )\Phi {\left(\nu \right)}^k}=\sum _{k=0}^{\infty }{c}_k(i,\alpha )\Phi {\left(\nu \right)}^k,\end{array}$$

(21)

where the coefficient $c_k(i,\alpha )$ can be determined recursively as

$$c_k(i,\alpha )=\frac{1}{b_0(i,\alpha )}\left[a_k\left(i\alpha \right)-\frac{1}{b_0(i,\alpha )}\sum _{r=1}^k{b}_r(i,\alpha )c_{k-r}(i,\alpha )\right].$$

Then, we can write

$$\begin{array}{c}\hfill F\left(x\right)=\sum _{k=0}^{\infty }{p}_k{H}_k\left(x\right),\end{array}$$

(22)

where (for $k\ge 0$)

$$\begin{array}{c}\hfill p_k=\frac{1}{{\mathrm{e}}^{\beta }-1}\sum _{i=1}^{\infty }\frac{{\beta }^i{c}_k(i,\alpha )}{i!}\end{array}$$

and $H_k\left(x\right)=\Phi {\left(\nu \right)}^k$ denotes the EBS cdf with parameters $\alpha$, $\beta$ and k.

Consequently, the pdf of X can be expressed as

$$\begin{array}{c}\hfill f\left(x\right)=\sum _{k=0}^{\infty }{p}_{k+1}{h}_{k+1}\left(x\right),\end{array}$$

(23)

where $h_{k+1}\left(x\right)$ denotes the $EBS(k+1,a,b)$ density function. Based on Equation (23), some mathematical properties of the OLLBSP distribution can be obtained by knowing those of the EBS distribution.

4.2. Moments

The ordinary moments of X can be derived from the probability weighted moments (PWMs) ([20]) of the BS distribution formally defined for p and r non-negative integers by

$$\begin{array}{c}\hfill {\tau }_{p,r}={\int }_0^{\infty }{x}^p\Phi {\left(\nu \right)}^r{g}_{a,b}\left(x\right)dx.\end{array}$$

(24)

The integral in Equation (24) can be computed numerically in software such as MAPLE, MATLAB, MATHEMATICA, Ox and R. Cordeiro and Lemonte [21] derived an expression for ${\tau }_{p,r}$ given by

$$\begin{array}{c}\hfill {\tau }_{p,r}=\frac{{\beta }^p}{2^r}\sum _{j=0}^r\left(\genfrac{}{}{0pt}{}{r}{j}\right)\sum _{k_1,\dots ,k_j}^{\infty }A(k_1,\dots ,k_j)\sum _{m=0}^{2s_j+j}{(-1)}^m\left(\genfrac{}{}{0pt}{}{2s_j+j}{m}\right)I\left(p+\frac{(2s_j+j-2m)}{2},\alpha \right),\end{array}$$

(25)

where $s_j=k_1+\dots +k_j$, $A(k_1,\dots ,k_j)={\alpha }^{-2s_j-j}{\overline{a}}_{k_1},\dots ,{\overline{a}}_{k_j}$,${\overline{a}}_k={(-1)}^k{2}^{(1-2k)/2}{\left[\sqrt{\pi }(2k+1)\right]}^{-1}$ and $I(p+(2s_j+j-2m)/2,\alpha )$ is obtained from Equation (4). These algebraic platforms currently have the ability to deal with analytic expressions of formidable size and complexity.

The sth moment of X can be expressed from Equation (23) as

$$\begin{array}{c}\hfill {\mu }_s^{\prime }=\sum _{k=0}^{\infty }(k+1)p_{k+1}{\tau }_{s,k},\end{array}$$

(26)

where ${\tau }_{s,k}$ is given by Equation (25). Equation (26) can be computed numerically in any symbolic mathematical software by taking in the sum a large number of summands instead of infinity.

The central moments (${\mu }_r$) and cumulants (${\kappa }_r$) of X are ${\mu }_r={\sum }_{k=0}^r{(-1)}^k\left(\genfrac{}{}{0pt}{}{r}{k}\right){\mu }_1^{\prime k}{\mu }_{r-k}^{\prime }$ and ${\kappa }_r={\mu }_r^{\prime }-{\sum }_{k=1}^{r-1}\left(\genfrac{}{}{0pt}{}{r-1}{k-1}\right){\kappa }_k{\mu }_{r-k}^{\prime }$, respectively, where ${\kappa }_1={\mu }_1^{\prime }$.

Plots of the skewness and kurtosis of X, for fixed values of a, b and $\beta$, as functions of $\alpha$ are displayed in Figure 3a,b, respectively. Plots of the skewness of X for a, b and $\alpha$ fixed, as functions of $\beta$ are displayed in Figure 3c,d, respectively. We can note that it is possible to obtain negative skewness for some parameter values.

5. The LOLLBSP Regression Model with Censored Data

If X is a random variable following the OLLBSP density function in Equation (11), the random variable $Y=log\left(X\right)$ has the LOLLBSP distribution (also referred to as the odd log-logistic sinh-normal Poisson distribution). The density function of Y obtained by replacing $b={\mathrm{e}}^{\mu }$ and performing two re-parameterizations reduces to (for $y\in \mathbb{R}$)

$$\begin{array}{c}\hfill f_Y\left(y\right)=\frac{2\beta \alpha {\xi }_{1}exp(-\frac{2}{a^2}{\xi }_2^2){\{\Phi \left(\frac{2}{a}{\xi }_2\right)[1-\Phi \left(\frac{2}{a}{\xi }_2\right)]\}}^{\alpha -1}}{a\sigma \sqrt{2\pi }[exp\left(\beta \right)-1]{\{{\Phi }^{\alpha }\left(\frac{2}{a}{\xi }_2\right)+{[1-\Phi \left(\frac{2}{a}{\xi }_2\right)]}^{\alpha }\}}^2}exp\left\{\frac{\beta {\Phi }^{\alpha }\left(\frac{2}{a}{\xi }_2\right)}{{\Phi }^{\alpha }\left(\frac{2}{a}{\xi }_2\right)+{[1-\Phi \left(\frac{2}{a}{\xi }_2\right)]}^{\alpha }}\right\},\end{array}$$

(27)

where

$${\xi }_1=cosh\left(\frac{y-\mu }{\sigma }\right)\mathrm{and}{\xi }_2=sinh\left(\frac{y-\mu }{\sigma }\right).$$

The parameter $\mu \in \mathbb{R}$ is a location parameter, $\sigma >0$ is a scale parameter and $\beta$, $\alpha$ and a are positive shape parameters. We refer to Equation (27) as the non-standard LOLLBSP model. Thus,

$$\mathrm{if}X\sim \mathrm{OLLBSP}(\alpha ,\beta ,a,b),\mathrm{then}Y=log\left(X\right)\sim \mathrm{LOLLBSP}(\alpha ,\beta ,a,\sigma ,\mu ).$$

The survival function of Y is given by

$$\begin{array}{c}\hfill S\left(y\right)=\left\{exp\left(\beta \right)-exp\left[\frac{\beta {\Phi }^{\alpha }\left(\frac{2}{a}{\xi }_2\right)}{{\Phi }^{\alpha }\left(\frac{2}{a}{\xi }_2\right)+{[1-\Phi \left(\frac{2}{a}{\xi }_2\right)]}^{\alpha }}\right]\right\}{[exp\left(\beta \right)-1]}^{-1}.\end{array}$$

(28)

The pdf of the standardized random variable $Z=(Y-\mu )/\sigma$ is given by

$$\begin{array}{ccc}\hfill f_Z\left(z\right)& =& \frac{2\beta \alpha cosh\left(z\right)exp[-\frac{2}{a^2}{sinh}^2\left(z\right)]{\{\Phi [\frac{2}{a}sinh\left(z\right)]\{1-\Phi [\frac{2}{a}sinh\left(z\right)]\}\}}^{\alpha -1}}{a\sqrt{2\pi }[exp\left(\beta \right)-1]{\{{\Phi }^{\alpha }[\frac{2}{a}sinh\left(z\right)]+{\{1-\Phi [\frac{2}{a}sinh\left(z\right)]\}}^{\alpha }\}}^2}\times \hfill \\ & & exp\left\{\frac{\beta {\Phi }^{\alpha }[\frac{2}{a}sinh\left(z\right)]}{{\Phi }^{\alpha }[\frac{2}{a}sinh\left(z\right)]+{\{1-\Phi [\frac{2}{a}sinh\left(z\right)]\}}^{\alpha }}\right\}.\hfill \end{array}$$

(29)

Plots of the density function in Equation (27) for selected parameter values are given in Figure 4. These plots show great flexibility for different values of the shape parameters $\alpha$ and $\beta$. They indicate that the density function in Equation (27) is very flexible, especially for modeling bimodal and asymmetric data and hence can be used in several applications.

In many practical applications, the lifetimes are affected by explanatory variables and parametric regression models are adopted to estimate univariate survival functions for censored data. A parametric model that provides a good fit to lifetime data tends to yield more precise estimates of the quantities of interest. Based on the standardized LOLLBSP density, we propose a linear location-scale regression model linking the response variable $y_i$ and the explanatory variable vector ${\mathbf{v}}_i^T=(v_{i1},\cdots ,v_{ip})$ by

$$y_i={\mathbf{v}}_i^T\mathbf{\tau }+\sigma z_i,i=1,\dots ,n,$$

(30)

where the random error $z_i$ has the density function in Equation (29), $\mathbf{\tau }={({\tau }_1,\cdots ,{\tau }_p)}^T$, $a>0$, $\alpha >0$, $\beta >0$ and $\sigma >0$ are unknown parameters. The parameter ${\mu }_i={\mathbf{v}}_i^T\mathbf{\beta }$ is the location of $y_i$. The location parameter vector $\mathbf{\mu }={({\mu }_1,\cdots ,{\mu }_n)}^T$ is given by a linear model $\mathbf{\mu }=\mathbf{v}\mathbf{\beta }$, where $\mathbf{V}={({\mathbf{v}}_1,\cdots ,{\mathbf{v}}_n)}^T$ is a known model matrix. The LOLLBSP regression model in Equation (30) opens new possibilities for fitting many different types of data.

Consider a sample $(y_1,{\mathbf{v}}_1),\cdots ,(y_n,{\mathbf{v}}_n)$ of n independent observations, where each random response is defined by $y_i=min\{log\left(t_i\right),log\left(c_i\right)\}$. We assume non-informative censoring such that the observed lifetimes and censoring times are independent. Let F and C be the sets of individuals for which $y_i$ is the log-lifetime or log-censoring, respectively. Conventional likelihood estimation techniques can be applied here. The log-likelihood function for the vector of parameters $\mathbf{\theta }={(\alpha ,\beta ,a,{\mathbf{\tau }}^T,\sigma )}^T$ from the model in Equation (30) has the form $l\left(\mathbf{\theta }\right)={\displaystyle \sum _{i\in F}}{l}_i\left(\mathbf{\theta }\right)+{\displaystyle \sum _{i\in C}}{l}_i^{\left(c\right)}\left(\mathbf{\theta }\right)$, where $l_i\left(\mathbf{\theta }\right)=log\left[f\left(y_i\right)\right]$, $l_i^{\left(c\right)}\left(\mathbf{\theta }\right)=log\left[S\left(y_i\right)\right]$, $f\left(y_i\right)$ is the density in Equation (27) and $S\left(y_i\right)$ is the survival function in Equation (28) of $Y_i$. The total log-likelihood function for $\mathbf{\theta }$ reduces to

$$\begin{array}{cc}\hfill l\left(\mathbf{\theta }\right)& =qlog\left(\frac{2\beta \alpha }{a\sigma \sqrt{\pi }}\right)-nlog[exp\left(\beta \right)-1]+\sum _{i\in F}log\left({\xi }_{i1}\right)-\frac{2}{a^2}\sum _{i\in F}{\xi }_{i2}^2\hfill \\ & +(\alpha -1)\sum _{i\in F}log\left[\Phi \left(\frac{2}{a}{\xi }_{i2}\right)\right]+(\alpha -1)\sum _{i\in F}log\left[1-\Phi \left(\frac{2}{a}{\xi }_{i2}\right)\right]\hfill \\ & -2\sum _{i\in F}log\left\{{\Phi }^{\alpha }\left(\frac{2}{a}{\xi }_{i2}\right)+{\left[1-\Phi \left(\frac{2}{a}{\xi }_{i2}\right)\right]}^{\alpha }\right\}+\beta \sum _{i\in F}\frac{{\Phi }^{\alpha }\left(\frac{2}{a}{\xi }_{i2}\right)}{{\Phi }^{\alpha }\left(\frac{2}{a}{\xi }_{i2}\right)+{\left[1-\Phi \left(\frac{2}{a}{\xi }_{i2}\right)\right]}^{\alpha }}\hfill \\ & +\sum _{i\in C}log\left\{exp\left(\beta \right)-exp\left[\frac{\beta {\Phi }^{\alpha }\left(\frac{2}{a}{\xi }_{i2}\right)}{{\Phi }^{\alpha }\left(\frac{2}{a}{\xi }_{i2}\right)+{[1-\Phi \left(\frac{2}{a}{\xi }_{i2}\right)]}^{\alpha }}\right]\right\},\hfill \end{array}$$

(31)

where q is the observed number of failures and

$${\xi }_{i1}={\xi }_{i1}\left(\mathbf{\theta }\right)=cosh\left(\frac{y_i-{\mathbf{v}}_i^T\mathbf{\beta }}{\sigma }\right),{\xi }_{i2}={\xi }_{i2}\left(\mathbf{\theta }\right)=sinh\left(\frac{y_i-{\mathbf{v}}_i^T\mathbf{\beta }}{\sigma }\right),$$

for $i=1,\dots ,n$. The MLE $\widehat{\mathbf{\theta }}$ of the vector of unknown parameters can be calculated by maximizing the log-likelihood of Equation (31). We use the NLMixed procedure in SAS to calculate the estimate $\widehat{\mathbf{\theta }}$. Initial values for a are taken from the fit of the LBS regression model with $\alpha =1$. From the fitted LOLLBSP model, we obtain the estimated survival function for $y_i$ given by

$$\begin{array}{c}\hfill S(y_i;\widehat{a},\widehat{\alpha },\widehat{\beta },\widehat{\sigma },{\widehat{\mathbf{\tau }}}^T)=\left\{exp\left(\widehat{\beta }\right)-exp\left[\frac{\widehat{\beta }{\Phi }^{\widehat{\alpha }}\left(\frac{2}{a}{\widehat{\xi }}_{i2}\right)}{{\Phi }^{\widehat{\alpha }}\left(\frac{2}{a}{\widehat{\xi }}_{i2}\right)+{[1-\Phi \left(\frac{2}{\widehat{a}}{\widehat{\xi }}_{i2}\right)]}^{\widehat{\alpha }}}\right]\right\}{[exp\left(\widehat{\beta }\right)-1]}^{-1},\end{array}$$

(32)

where ${\widehat{\xi }}_{i2}=sinh\left(\frac{y_i-{\mathbf{v}}_i^T\widehat{\mathbf{\beta }}}{\widehat{\sigma }}\right).$ Under standard regularity, the asymptotic distribution of $(\widehat{\mathbf{\theta }}-\mathbf{\theta })$ is multivariate normal $N_{p+4}(0,K{\left(\mathbf{\theta }\right)}^{-1})$, where $K\left(\mathbf{\theta }\right)$ is the expected information matrix. The asymptotic covariance matrix $K{\left(\mathbf{\theta }\right)}^{-1}$ of $\widehat{\mathbf{\theta }}$ can be approximated by the inverse of the $(p+4)\times (p+4)$ observed information matrix $-\ddot{\mathbf{L}}\left(\mathbf{\theta }\right)$. The approximate multivariate normal distribution $N_{p+4}(0,-\ddot{\mathbf{L}}{\left(\mathbf{\theta }\right)}^{-1})$ for $\widehat{\mathbf{\theta }}$ can be used in the classical way to construct approximate confidence regions for some parameters in $\mathbf{\theta }$.

6. Applications

In this section, we present two applications to prove empirically the flexibility of the OLLBSP distribution. We compare the results of the fits of the OLLBSP model and other models under some goodness-of-fit statistics. In fact, we compare the fits of the OLLBSP distribution with six other lifetime models, namely:

• BS distribution [1];
• Marshall–Olkin extended-BS (MOEBS) distribution [22];
• Odd log-logistic BS (OLLBS) distribution [7];
• Beta BS (BBS) distribution [21];
• Kumaraswamy BS (KwBS) distribution [23]; and
• Exponentiated generalized BS (EGBS) distribution [24].

The MLEs of the model parameters are computed and the goodness-of-fit statistics for these models fitted to the data sets are compared. These estimates are evaluated using the BFGS script as well as the goodness-of-fit measures. The MLEs of the parameters and these measures including the Cramér-von Mises ($W^{*}$) and Anderson–Darling ($A^{*}$) statistics defined by Chen and Balakrishnan [25] and Kolmogorov–Smirnov (K-S) statistics are calculated to compare the fitted models. In general, the smaller the values of these statistics, the better is the fit to the data. All computations are performed using the R statistical software [26].

6.1. Data Set 1: Breaking Stress of Carbon Fibres (In Gba)

As a first application, we consider the data set consisting of $n=100$ observations on breaking stress of carbon fibers (in Gba) reported by [27]. A brief summary of these data are: mean $=2.621$, median $=2.700$, standard deviation $=1.013885$, skewness $=0.36264$, kurtosis $=0.0431$, minimum $=0.390$, and maximum $=5.560$. The data have positive skewness and kurtosis. These values indicate that the empirical distribution is skewed to the right and leptokurtic.

Table 2. MLEs of the parameters for the fitted models to the breaking stress of carbon fibres data.

Model	a	b	$\mathit{\alpha }$	$\mathit{\beta }$	${\mathit{W}}^{*}$	${\mathit{A}}^{*}$
BS$(a,b)$	0.4621	2.3661			0.2978	1.6181
	(0.0326)	(0.1064)
MOEBS$(a,b,\beta )$	2.8026	0.0250		9909.4580	0.0798	0.5580
	(0.4486)	(0.0078)		(435.3655)
OLLBS$(a,b,\alpha )$	217.3455	2.4788	530.5732		0.2461	1.2847
	(0.1095)	(0.0557)	(0.1012)
OLLBSP$(a,b,\alpha ,\beta )$	0.3173	1.2653	0.3993	4.7372	0.0711	0.3945
	(0.0977)	(0.3080)	(0.1697)	(1.6564)
BBS$(a,b,\alpha ,\beta )$	97.3898	0.0075	629.2049	466.4159	0.1040	0.5363
	(24.6425)	(0.0011)	(236.0238)	( 203.7381)
KWBS$(a,b,\alpha ,\beta )$	12.9322	0.0054	26.0008	2.1588	0.1486	0.7541
	(1.6490)	(0.0006)	(5.2036)	(0.7928)
EGBS$(a,b,\alpha ,\beta )$	1.7685	69.7519	105.9274	0.3928	0.1603	0.8570
	(0.0736)	(0.0859)	(0.4150)	(0.0726)

lists the MLEs of the parameters and their standard errors (in parentheses) from the fitted models and the measures $W^{*}$, $A^{*}$, K-S and its p-value. The results indicate that the OLLBSP model has the smallest values of these statistics and large p-value for the K-S distance among the fitted models, and therefore it could be chosen as the best model to the current data.

Plots of the estimated pdf and the empirical and estimated cdf of the OLLBSP model are displayed in Figure 5. It is clear from these plots that the OLLBSP model provides the best fit to the histogram of the breaking stress data.

6.2. Data Set 2: Ambient Temperature Data

The data set consists of failure times (X) of eight components at three different temperatures. The data are obtained from [28]. The original sample size is $n=24$. The following variables are associated with each component: $x_i$, observed time (in years); and $v_{i1}$, temperature (temperature category: 100, 120 and 140), for $i=1,2,\dots ,24$. We analyze these data using the LOLLBSP regression model. First, we consider the regression model

$$\begin{array}{c}\hfill y_i={\tau }_0+{\tau }_1{v}_{i1}+\sigma z_i,i=1,\dots ,24,\end{array}$$

(33)

where the random errors $z_i$’s ($i=1,\dots ,n$) are independent random variables having the density function in Equation (29) and $y_i=log\left(x_i\right)$ for $i=1,\dots ,24$. These data were originally analyzed by Cordeiro et al. [2] using the log-McDonald–Birnbaum–Saunders (LMcBS) regression model. The LMcBS distributions includes as special cases two important distributions: the log-beta Birnbaum–Saunders (LBBS) model and the log-Birnbaum–Saunders (LBS) model. Next, we compare the proposed regression model with other known regression models.

The values of the AIC (Akaike Information Criterion), CAIC (Corrected Akaike’s Information Criterion) and BIC (Bayesian Information Criterion) statistics to compare the LOLLBSP, LMcBS, LBBS and LBS regression models are listed in

Table 3. AIC, BIC and CAIC statistics.

	Statistic
Model	AIC	CAIC	BIC
LOLLBSP	61.1	66.0	68.2
LMcBS	76.7	67.7	68.9
LBBS	73.2	78.2	80.3
LBS	66.4	67.6	69.9

. The LOLLBSP regression model outperforms the LMcBS, LBBS and LBS models, irrespective of the criteria and thus the proposed regression model can be used effectively in the analysis of these data.

Table 4. MLEs, standard errors (SEs) and [p-values] for the LOLLBSP regression model fitted to the ambient temperature data.

Model	$\mathit{\beta }$	$\mathit{\alpha }$	a	$\mathit{\sigma }$	${\mathit{\beta }}_0$	${\mathit{\beta }}_1$
LOLLBSP	1.5548	0.1062	0.1031	3.9382	6.6303	−0.02936
	(0.9542)	(0.01876)	(0.005682)	(1.06512)	(0.2890)	(0.001231)
					[$<0.001$]	[$<0.001$]

lists the MLEs of the parameters for these regression models fitted to the current data using the NLMixed procedure in SAS. Iterative maximization of the logarithm of the likelihood function in Equation (31) starts with initial values for $\mathbf{\tau }$ and $\sigma$ taken from the fit of the LBS regression. We note from the fitted LOLLBSP regression model that the explanatory variable $v_1$ is significant at $1\%$ and that there is a significant difference among the levels of the temperature for the failure times.

A graphical comparison among the levels of the explanatory variable $v_1$ is given in Figure 6a. These plots provide the Kaplan–Meier (KM) estimate and the estimated survival function of the LOLLBSP regression model given by Equation (32) for different levels of $v_1$. The plot of the LOLLBSP hrf in Figure 6b reveals that it has an increasing–decreasing–increasing shape. As discussed in Section 2, the OLLBSP distribution is adequate for this type of hrf (see Figure 6c). Based on these plots, it is clear that the LOLLBSP regression model gives a superior fit to the ambient temperature data.

7. Concluding Remarks

We propose the odd log-logistic Birnbaum–Saunders–Poisson (OLLBSP) model to extend the Birnbaum–Saunders (BS) distribution pioneered by Birnbaum and Saunders [1]. We provide a mathematical treatment for the new distribution including expansions for the cumulative and density functions. We derive explicit expressions for the ordinary moments. The estimation of the model parameters is approached by the maximum likelihood method. We verify by Monte Carlo simulations that the normal distribution provides an adequate first-order approximation to the finite sample distribution of the maximum likelihood estimators. A new regression model is defined based on the logarithm of the OLLBSP random variable. We provide two applications to real data sets to illustrate the usefulness of the proposed models in the literature of the fatigue life distributions.