SOME USEFUL PATHWAY MODELS FOR RELIABILITY ANALYSIS Текст научной статьи по специальности «Математика»

SOME USEFUL PATHWAY MODELS FOR RELIABILITY

ANALYSIS

T. PRINCY

Department of Statistics Cochin University of Science and Technology Cochin-682022, Kerala, India princyt@cusat.ac.in

Abstract

In this paper, we first discuss pathway model in general. Then a special case for the real scalar variable is considered. This special case is relevant in reliability problems. In the pathway model, an arbitrary function is introduced so that the hazard function resulting from this model is of a given shape such as a bathtub type hazard function. The model is also derived by using an entropy optimization procedure by introducing a new entropy measure. It is shown that a large number of densities in current use are connected to the pathway model. Certain combinations of pathway densities resulting in hazard functions of desired shapes, multi-component failure situation etc are examined from a reliability point of view. For further use of the proposed model, the unknown parameters are estimated using the method of maximum likelihood estimation. The behaviour of the reliability measure has been observed graphically for arbitrary values of the parameters related to the number of components and operating time.

Keywords: Pathway model, reliability analysis, hazard function, entropy optimization

1. Introduction

The common knowledge used in the literature is that an approximate model for the data at hand can be found, regardless of whether the data comes from the biological, physical, engineering, social, or other fields. The information at hand can be described in the area around the stable condition or along a path that leads there. In (2005) Mathai [1] introduced a pathway model to cover the stable as well as the transitional stages, which describes transitions of rectangular matrix-variate distributions in the real case. It is a mathematical or stochastic model that switches from one functional form to another through pathway parameter so that intermediate steps can be represented. Also the parameters of pathway model connect many families of functions, and as a result, it is possible to identify a suitable member either a given family or between stages of two families. The idea was extended to cover the complex rectangular matrix-variate case in Mathai and Provost [2]. The real scalar version of the pathway model can be stated as the following:

r n

fi(x) = ci|x|Y[1 - a(1 - q)|x|5 ] 1-q (1)

for a > 0, q < 1, y > -1,5 > 0, n > 0,1 - a(1 - q) |x|5 > 0, and zero elsewhere, where c1 can act as the normalizing constant if we wish to create a statistical density out of f1 (x). The support of (1) is on -[a(1 - q)]-5 < x < [a(1 - q)]-It is a finite range model. The functional part in a basic type-1 beta model is of the form xa-1 (1 - x)fi-1, a > 0, fi > 0,0 < x < 1 and zero elsewhere, and hence the model in (1) can be looked upon as a generalized, extended and power-transformed type-1 beta model for q < 1. When q ^ 1 the range of x will go from - ro to ro. Note that (1) is a symmetric model for x < 0 and x > 0. If an asymmetric model is required then different weights can be used for x < 0 and x > 0 situations. These differing weights can be

introduced either by multiplying with constants or through one of the parameters. For q > 1, write 1 — q = — (q — 1), q > 1, then (1) switches into the following model:

r - n

f2(x) = C2|x|7[1 + a(q — 1)1x1*] q—1 (2)

for < x < to, a > 0, q > 1,5 > 0,7 > —1, q > 0. The model in (2) can be looked upon as a generalized, power-transformed and extended type-2 beta family of functions. The functional part of the basic type-2 beta model is xa—1 (1 + x)—(a+ft),0 < x < to, a > 0, ft > 0. Thus, (2) is the generalized and extended version of this basic type-2 beta function. When q ^ 1— in (1) and q ^ 1+ in (2) the models in (1) and (2) go to the generalized and extended gamma family of functions as follows:

f3 (x) = C3 |x|7 e—aq|x|5 (3)

for a > 0, q > 0,5 > 0,7 > —1, — to < x < to. This is the generalized gamma family of functions. Thus, the basic pathway model is (1), and (2) and (3) are available from (1). For q < 1 the family of functions is the generalized and extended type-1 beta family of functions. When q goes to 1, then one goes into the generalized and extended gamma family of functions. When q moves to q > 1 then we go into the generalized and extended type-2 beta family of functions. The parameter q enables us to go to three different families of functions and hence q is called the pathway parameter. The pathway idea in model building situation is to capture the stable situation as well as the unstable neighborhoods by the same model. If the gamma family is the stable situation in a physical problem then the paths leading to this stable situation through the generalized type-1 beta model and generalized type-2 beta model and the transitional stages are captured by the pathway parameter q. If (1) to (3) are to be treated as statistical densities then the following are the normalizing constants:

5[a(1 — q)] 2-1 r( + 1 + ) C1 =-—-qr-, for 7 > —1, a > 0, q < 1, q > 0,5 > 0 (4)

1 2r( )r( +1) r - ' "

C2

5[a(q — 1)] r( "ny) q 7 + 1 ---q-, for q > 0,5 > 0, a > 0, q > 1,7 > — 1, ^r — > 0 (5)

2r( ^ )r( -L — 7+i) 1 q r q—15 ()

C3 = ^rq-T, a > 0, q > 0,5 > 0,7 > —1. (6)

2r( 5 )

Our interest is to consider a special case of the pathway model for x > 0 and when 7 = 5 — 1. In this case the normalizing constants reduce to very simple forms. These forms of the pathway

model are the most relevant in reliability analysis. Our main focus will be on this special case. Before we concentrate on the special case, let us see some other special cases of the general models. Note that (1) for 7 = 0, a = 1,5 = 1, q = 1, q < 1, q > 1, q ^ 1 is Tsallis statistics in non-extensive statistical mechanics. It is claimed that between 1990 and 2010, over 3000 articles are published on this topic of Tsallis statistics. Details of the development may be seen from his book Tsallis [3]. Equation (2) for q > 1, q ^ 1, a = 1,5 = 1, q = 1 is superstatisitcs in statistical mechanics. Dozens of articles are also written in this area since 2003. The basic article in this area is Beck and Cohen [4].

If location and scale parameters are to be incorporated into the models in (1) to (3) then replace |x| by || for some constant ^ and some constant a > 0 in (1) to (3). Note that (3) for 7 = 0,5 = 2 is the normal or Gaussian density. For x > 0, (3) produces the generalized gamma density, Weibull density, gamma density, chisquare density, exponential density, Rayleigh density, Maxwell-Boltzmann density etc. Exponentiation in (2), that is put x = e—cy, c > 0, produces the generalized logistic density of Mathai and Provost [5], logistic density and other special cases of the generalized logistic density, which are relevant in reliability analysis. (2) can produce Cauchy density, Student-t density, F-density etc. For 7 = 5 — 1, (1) and (2) produce many models in reliability analysis. As a limiting form of (2) one can obtain Fermi-Dirac density from (2) and Bose-Einstein density from (1), after exponentiation.

2. Construction of the Pathway Model through Entropy Optimization

Model building in physical situations is often done by optimizing an appropriate entropy measure and then deriving the model from therein. Consider Mathai's entropy, for an earlier version see Mathai and Haubold [6],

1-q+q

/ N fX[ f (X)] n dX - 1

Mq(f ) = Jx[f( )q .,-, q = 1, q < 2 (7)

q - 1

where f (X) is a statistical density where X could be real or complex scalar or matrix variable, and fX denotes the integral over the support of f. A corresponding version for the discrete situation

i-i

can be constructed. One can look upon (7) as an expected value of f q which then corresponds to Kerridge's measure of inaccuracy for q = 1, see Mathai and Rathie [7]. For the real scalar variable case, (7) for q = 1 can be looked upon as a modified Havrda-Charvat entropy, see Mathai and Rathie [7]. In the following discussion we consider the real scalar variable case first. Let us optimize (7) subject to the following moment-type restrictions for real scalar variables:

/ x7 V1 f (x)dx = fixed & / x7^+3f (x)dx = fixed (8)

Jx J x

for 7 > -1, q < 1, $ > 0, and it is assumed that the integrals in (8) exist. Note that for 7 = 0, (8) states that the total integral is unity and that the first moment is given. This is equivalent to the physical law of conservation of energy when we consider energy distribution. It is convenient to use calculus of variation for optimizing (7) subject to the conditions in (8). Then the Euler equation is the following:

d i-q+q 7 (i-q) 7 (i-q) , $ n

df [f - A1 x n f + A2x n +$f ] = 0 (9)

where A1 and A2 are Lagrangian multipliers. Then (9) gives

(1-q) 7 (1-) o

f q = y1 xr q [1 -y2x$] for some constants ^ and , which then gives

t JL

fx = v x7 [1 - p2x$] 1-q (10)

for some v and . Take ^ = a(1 - q) and v as the normalizing constant to obtain the model (1). In (7) if X is a p x 1 vector random variable and if (8) is replaced by the following conditions

r -, (q-1)

[(X - p)'V-1(X - v)]7~ f (X)dX = fixed (11)

X

and

( q - 1 )

[(X - ^)'V-1(X - v)]7~r-+$f (X)dX = fixed (12)

X

where V is p x p real symmetric and positive definite constant matrix, ц is a p x 1 constant vector, a prime denotes transpose, 7 > 0, q > 1, q > 0, then from (2.3) and (10) we have the following density:

f4(X) = c4[(X - y)'V-1 (X - p)]7[1 + a(q - 1){(X - p)'V-1 (X - p)}$]-q-r,q > 1. (13) Then, when q ^ 1, f4(X) goes to f5(X) given by

f5(X) = c5[(X - y)'V-1 (X - y)]7e-aq[(X-y)'V-1 (X-y)]$. (14)

Note that (13) and (14) are also associated with type-2 and gamma distributed random points in Euclidean n-space, p < n, see Mathai [8]. Also, (14) for 7 = 0,5 = 1 is the p-variate Gaussian density with mean value vector ц and covariance matrix V. The quantity

(X — p)'V—1 (X — p) = c > 0 (15)

is known as the ellipsoid of concentration. Hence the constraints can be explained in terms of ellipsoid of concentration. In (7) if X is a p x q, q > p rectangular matrix-variate random variable and if the conditions in (8) are replaced by the following:

i [tr(AXBX')]7V1 f (X)dX = fixed (16)

J X

and

i [tr(AXBX')]7"V+5f (X)dX = fixed (17)

X

where A is p x p and B is q x q constant real positive definite matrices, q > 1, q > 0,5 > 0 then the steps in (9) and (10) give the density

f6(X) = C6[tr(AXBX')]7[1 + a(q — 1){tr(AXBX')}5] -=T,q > 1 (18)

which when q ^ 1 gives

f7(X) = c7 [tr( AXBX' )]7 e—aq[tr(AXBX')]5 (19)

for a > 0, q > 0,5 > 0,7 > —1. Note that (19) for 7 = 0,5 = 1 is the real matrix-variate Gaussian density. By replacing X by X — M, where M is a p x q constant matrix, one can also incorporate a location parameter matrix in the models in (18) and (19). Mathai and Princy [9] illustrate the significance of models (18) and (19) to the real multivariate reliability analysis.

3. A Special Case of Pathway Model and Reliability Analysis

A special case of the pathway model (1),(2),(3) is the case where 7 = 5 — 1, for x > 0. This will

then correspond to a power-transformed basic model. Take the basic model in (1) as the one with

q

7 = 0 and 5 = 1 for x > 0 or c[1 — a(1 — q)x]1—q where a > 0, q < 1 and c is a constant. Make the power transformation here or put x = y5 for some 5 > 0. Then we have

gi(y) = C1 y5—1 [1 — a(1 — q)y5]^,q < 1,q > 0,a > 0,5 > 0 (20)

and Ci is a constant. The cases for q > 1, q ^ 1 are available from (20) as

g2(y) = C2y5—1 [1 + a(q — 1)y5] — , q > 1 (21)

and

g3 (y) = C3y5—1e—aqy5, a > 0, q > 0,5 > 0. (22)

If g1,g2,g3 are to be treated as statistical densities then the normalizing constants are the following:

C1 = 5a(q + 1 — q), q < 1; C2 = 5a(q + 1 — q), q > 1; C3 = 5qa. (23)

The following are the graphs showing the relative positions of g1,g2,g3.

Figure 1: Pathway models for y = 6 — 1 Note that in §2 (y) and §3 (y) the densities of y also belong to the same families. Put x = y then

g2 (x) and g3 (x) go to the following for 0 < x < œ:

g2*(x) = a% + 1 - q)x-5-1 [1 + a(q - 1)x-5] 1-q

-5i A

and

g3„ (x)= aS(n + 1 - q)x-5-1e-anx~

(24)

(25)

Thus, in g2 (y) and g3 (y) both 5 and - 5 with 5 > 0 are admissible with the corresponding change in y5-1. Let us compute the survival functions. For q < 1

_ 1

r [a(1 -q)] 5

S1 (t) = Pr{x > t} = a5(n + 1 - q)x5-1 [1 - a(1 - q)x5] 1-qdx

Jx=t

= [1 - a(1 - q)t5] ^+1,q < 1,n > 0,a > 0,5 > 0.

For q > 1 the survival function S2 (t) is given by the following:

S2(x) = [1 + a(q — 1)t6] q—1 +1,q > 1,n > 0,a > 0,6 > 0

and for q ^ 1

S3 (t) = Pr{x > t} = e—ant6, a > 0, n > 0,6 > 0. The hazard functions for q < 1, q > 1, q ^ 1, denoted by h1 (t), h2 (t), h3 (t) are the following:

(t) = §44 =6 a(n +1—q)t6;1, q < 1

1W S1 (t) 1 — a(1 — q)t6 ^

6a(n + 1 — q)t6—1 ^ .

hi(t) = —t2—;-x , q > 1

iW 1 + a(q — 1)t6 '

h3 (t) = a6n t6—1.

(26)

(27)

(28)

(29)

(30)

(31)

Here (29) to (31) do not show any interesting shapes. A useful shape is the bathtub shape. All different shapes are available from (29) to (31).

4. Exponentiation of the Pathway ModeL in the SpeciaL Case

In this paper our main objective is to construct some useful bathtub shaped models for reliability analysis. For this, let us make the transformation y = ecx, c > 0 in (20). Then the condition

1

1 - a(1 - q)y5 > 0 ^ 1 - a(1 - q)ec5x > 0 or x < — ln[a(1 - q)]

for q < 1. But for q > 1 and q ^ 1, —to < x < to under exponentiation. Then the models reduce to the following:

g4(x) = c«5(n + 1 — q)ec5x [1 — «(1 — q)ec5x] ï-?, q < 1 (32)

for c > 0, « > 0, 5 > 0, n > 0, n + 1 — q > 0, —to < x < -j-f ln[«(1 — q)] and zero elsewhere.

g5(x) = c«5(n + 1 — q)ec5x [1 + «(q — 1)ec5x] — A, q > 1 (33)

for « > 0, n > 0,5 > 0, c > 0, n + 1 — q > 0, —to < x < to.

g6 (x) = c«5nec5x e—«nec5x (34)

for « > 0, n > 0,5 > 0, c > 0, — to < x < to. The hazard functions for (32) to (34) can be seen to the the following:

, .. «c5(n + 1 — q)ec5t

h5(t) = ' „-\~~nr,q < 1. (35)

5W [1 — «(1 — q)ec5t] ^ v '

, «c5(n + 1 — q)ec5t

h6(t)= [i IViqe*],q > (36)

h7(t) = ac5nec5t, « > 0, c > 0,5 > 0, n > 0. (37)

If we make the transformation y = e—cx, c > 0 in (20), we can connect the transformed models to Bose-Einstein density, Logistic and Fermi-Dirac densities. The corresponding transformed models are

g7(x) = c«5(n + 1 — q)e—c5x [1 — «(1 — q)e—c5x] ï-, q < 1 (38)

for c > 0, « > 0,5 > 0, n > 0, n + 1 — q > 0, 1 ln[«(1 — q)] < x < to and zero elsewhere.

g8(x) = c«5(n + 1 — q)e—c5x[1 + «(q — 1)e—c5x] — -F1,q > 1 (39)

for « > 0, n > 0,5 > 0, c > 0, n + 1 — q > 0, —to < x < to.

g9(x) = c«5ne—c5x e—«ne-c5x (40)

for « > 0, n > 0,5 > 0, c > 0, —to < x < to.

4.1. Bose-Einstein density

Consider the limiting case n + 1 = This is not admissible in the density (38). For n + 1 = q can

we re-normalize the function for 0 < x < to and create a density out of it? That is, can

1

gi0(x) = ci0[«(T1^)ec5x — 1]—1 = C7[ea+c5x — 1]—11,0 < x < to (41)

be a density where «(1—q) = e* ?. Let us consider the integral

f TO 1 f to

[ea+c5x — 1]—Tdx = M [eu — 1]—TdM,C = c5, put v = eu

■¡0 C *

1 fTO 1 1 , 1 /"TOr 1 1n ,

= */--=~ / [-t — - ]dv

C Je* vv — 1 C Jea v — 1 v

= fln f-A ) , e* = 1.

C \ea — 1) =

Hence for c10 = C[ln (e*—t)]—1 is a density for 0 < x < to. This density in (41) is the Bose-Einstein density in Physics.

4.2. Logistic and Fermi-Dirac densities

Consider (33) and (39) for q = | and q = 1 and let a = 2. Then (33) and (39) become

e—cSx ecSx

gl1(x) = c5 (1 + e—c5x)2 = c5 (1 + ec5x )2, < x < (42)

This is the logistic density. Hence (33) is a generalized form of the logistic density. Still more general forms can be obtained by deleting the condition 7 = 5 — 1 in (2) and then exponentiating. Consider again the model (39). This can be written as follows:

g8 (y) = ac5(q + 1 — q)e—c5x [1 + a(q — 1)e—c5x ] — ^

_ q „ (n+1—q) = ac5(q + 1 — q)[a(q — 1)] q—1 [e—c5x] q—1 (43)

1 x n

x [1 + , \,ec5x] — q-r. (44)

a(q - 1)

Note that q——j = 1or q + 1 — q = 0is not an admissible value in (44). This is the limiting situation. Can we re-normalize the function in (44) for q + 1 = q and 0 < y < to and create a density from (44)? Consider

1

g12(x) = c12[1 + ea+c5x]—\ = ea. (45)

Consider

/•TO /»TO ^ /»TO

g12(x)dx = c12 / [1 + ea+c5x]—11dx = (1 + e") —1dM, put v = eu

Jo Jo c5 Ja

c12 [ TOr 1 1 IJ c12, z e"

= isL [vT+v]dv = 75"ln(rre«).

c5 ( e

(1+.

Fermi-Dirac density in Physics.

Hence for c12 = c5eX , g12(x) is a density for 0 < x < to and this density is known as ln( )

4.3. Arbitrary function P(x)

Instead of power transformations and exponentiation let us take an arbitrary function P(x) in our basic pathway models in (20) to (22). Then the models become the following:

g13(x) = a5(q + 1 — q)P' (x)[P(x)]5—1 [1 — a(1 — q)(P(x))5 ] ^ (46)

for a > 0,5 > 0, q > 0, q < 1, P'(x) = dxP(x) > 0, P(x) > 0.

g14 (x) = a5(q + 1 — q)P' (x)[P(x)]5—1 [1 + a(q — 1)(P(x))5 ] — ^ (47)

for a > 0,5 > 0, q > 0, q > 1, P'(x) > 0, P(x) > 0.

g15 (x) = a5q P' (x)[P(x)]5—1e—aqtP(x)]5 (48)

for a > 0,q > 0,5 > 0,P'(x) > 0,P(x) > 0. Thus, P(x) must be a positive and increasing function so that P(x) and P'(x) are positive. One such function is the distribution function for an arbitrary density. Let y be a real continuous random variable with density g(y) and distribution function Fy(x) = Pr{y < x}. Then take P(x) = Fy(x) so that P'(x) = g(x) is the density which is positive on the support of g(x). Observe that when Fy(x) is a distribution function then [Fy(x)]5,5 > 0 is again a distribution function for some other random variable. Our problem is the following: Can we select a function P(x) so that the hazard function coming out of the pathway model is of the desired shape? From (29) and (30) note that the hazard function is of the following structure:

h1 (x) = — (q + — — q) dx ln[1 — a(1 — q)(P(x))5],q < 1 (49)

for 1 - a(1 - q) [P(x)]5 > 0,q < 1,a > 0,n > 0,5 > 0, and

h2(x) = (n + -1 q) dx ln[1 + a(q - 1)[P(x)]5],q > 1,

(50)

for a > 0,6 > 0, n > 0. Then for a pre selected hazard function of the desired shape one can solve for (49) and (50) and compute the corresponding P(x). This is the aim. Suppose that the hazard function h1 (x) is of the form

h1(x) =

1

+

1

x - a + e1 b - x + e2

, e1 > 0, e2 > 0, a < x < b

and zero elsewhere. When x = a one has -1 + jr—K-.—

e1 (b-a)+e2

and when x = b it is +

1

bathtub shaped, take b - a large. Then

J h1 (x)dx = ln(x - a + e1 ) - ln(b - x + e2) = ln ^

e2 (b-a)+e1

. It is

x - a + e1 b - x + e2

That is, for example, (49) yields

/ h1 (x)dx = -

x

(n ^ m[1 - a(1 - q)[ P( x )]'] = x-a+S).

That is,

ln[1 - a(1 - q) [P(x)]5] = ln

b - x + e2 x - a + e1

1-q n+1-q

One solution is

[ P( x )]5 =

a(1 - q)

1

b - x + e2 x - a + e1

1-q '

n+1-q

But when we impose the conditions P(x) > 0, P' (x) > 0 the bathtub shape cannot be maintained. The conditions P(x) > 0, P'(x) > 0 can be met only by taking e1 = e2 + (b - a) or by shifting the minimum point to the left-end. Then we get a semi-bathtub shaped hazard function of the following form:

Figure 2: Semi-bathtub shaped hazard function

Consider the case of the arbitrary function P( x) being the distribution function for some random variable. Let P(x) = Fy(x) for a distribution function Fy(t) = Pr{y < t} for some continuous random variable y. Then P' (t) = §y (t) where §y (t) is the density of some random variable y, evaluated at t. Let P(t) = [Fy (t)]. Then the hazard function of (46)

h13(t) =

P'(t) _ a5(n + 1 - q)gy(t)[Fy(t)]

51

Pr{x > t}

1 - a(1 - q) [Fy(t)]5

Take 5 = 1. Then

h13(t) = ,0 < F-") < 1.

1

T. PRINCY RT&A, No 1 (72) SOME USEFUL PATHWAY MODELS FOR RELIABILITY ANALYSIS_Volume 18, March 2023

Take a fast decreasing gy (t) with gy (0) — 0 then such a gy (t) should produce a bathtub shaped hazard function curve. We shall examine a few such cases here. Case (1): Let

gy (t) = 0e-0t, t > 0,0 > 0

and zero elsewhere. In this case the hazard function

a0(n + 1 - q)e-01

hi4(t) =

(1 - a(1 - q))[1 - e-01]

for a(1 - q) < 1, > 0, q < 1, n > 0, n + 1 - q > 0. Some of the plots are given below for the various values of the parameters of h14 (t)

Figure 5: a = 2,0 — 2, n — q = 10

Case (2): Consider a type-1 beta model with the density gy (t) = f (b - t)f-1,0 < t < b, fi > 0. Then the distribution function is Fy (t) = bf - (b - t)f. Again, taking the pathway model for q < 1 and y = 5 - 1 in (1.1), with 5 = 1, produces the hazard function

a(n + 1 - q)f(b - t)f-1 b

h15(t) " 1 - «(1 - q)[bf - (b - t)f],0 < f < b.

When t = 0, h(0) — a (n + 1 - q)fbf 1. When t is nearing b the numerator nears zero and the denominator nears 1 - a(1 - q). Observe that a(1 - q) < 1. Select a and q such that a(1 - q) nears 1 then h15 (t) will be a large quantity. This is plotted for various values of the parameters a, q, f, b.

Figure 6: a = 1, £ = 10, b = 1, n = 20, q = 10 Figure 7: a = 0.0002, £ = ^ b = 50, n = 5, q = 10 Case (3): Consider a Pareto type density for §y (t). Let

§y(t) = Ca(1 + t)—(^+1),0 < t < b,a > 0, C =

(1 + b)a (1 + b)a — 1.

When b = ^ then C = 1 and in this case §y (t) is a type-2 beta density. The hazard function in this case of pathway model of (1.1) for q < 1 with y = 6 — 1 and 6 = 1 is as follows:

a(n + 1 — q)C a(1 + t)—(a+1) h16 (t) 1 — a(1 — q)C [1 — (1 + t)—a ]

for a > 0, q < 1, a > 0,0 < t < b. Note that when t = 0 the denominator is 1 and the numerator is a(n + 1 — q)Ca and when t is large then the numerator is very small and the denominator nears 1 — a(1 — q). Selecting a and q such that a(1 — q) < 1 but close to 1 we can make a bathtub shaped curve. The following are the curves for some selected values of the parameters.

Figure 8: a = 3, a = 1, n = 12, b = 50, q = 10 Figure 9: a = 2, a = n = 1, b = 20, q = 150 Hence this approach does not produce a satisfactory hazard function.

5. Combinations of Pathway ModeLs with Other ModeLs Consider an exponentiated pathway model of (38) for c = 1, of the form

/10(x) = a6(n + 1 — q)e—6x[1 + a(q — 1)e—6x] q > 1, for n > 0,6 > 0, a > 0. Take another function of the power function type.

f11 =

Y

eYx,0 < x < b, y > 0

(51)

(52)

and zero elsewhere. Let f (x) be a convex combination of f10 (x) and f11 (x). Let

f (x) — —x—f10 (x) + —x—f11 (x), a1 > 0, a2 > 0. a1 + a2 a1 + a2

Then the survival function S(t) is the following:

S(t) — Pr{x > t} — — i f10(x)dx + a2 /bf11 (x)dx a1 + a2 t a1 + a2 t

a1 [1 + a(q - 1)e-5t]-*+1 + * [ebY - e*].

a1 + a2

Therefore the hazard function is the following:

a1 + a2 e*b - 1

f (x) _ a1 a5(n + 1 - q)e-5x[1 + a(q - 1)e-5x] q-1 + a2eYTeYx

h( x) —

S( x)

a1 [1 + a(q - 1)e-5x] q-1 +1 + a2iYb^I[ebY - e*x]

,0 < x < b.

(53)

Note that a1 + a2 will be canceled. Hence we may take any positive linear combinations of f10 (x) and f11 (x) to get the numerator of h( x) and the same linear combination to get the denominator. This (53) is plotted for various parameter values. The graphs are for the following combinations of the parameters:

Figure 10

Figure 11

Figure 10: a1 — 1, a2 — 1, b — 40, a — 1,5 — 2, q — 1.5, n — 1, Y — 0.1 Figure 11: a1 — 9, a2 — 1, b — 40, a — 1,5 — 2, q — 1.5, n — 1, Y — 0.1

Figure 12

Figure 13

Figure 12: a1 — 7, a2 — 3, b — 50, a — 1.5 — 2, q — 1.5, n — 1, Y — 0.1 Figure 13: a1 — 7, a2 — 3, b — 70, a — 1,5 — 2, q — 1.5, n — 1, Y — 0.1

Figure 14

Figure 15

Figure 14: a1 = 7, a2 = 3, b = 80, a = 1,6 = 2, q = 1.5, n = 1, Y = 0.1 Figure 15: a1 = 1, a2 = 1, b = 70,6 = 2.2, q = 1.9, n = 1, Y = 0.1

This approach can give hazard functions of the desired shapes by selecting the parameters appropriately.

We shall try combination of /10 (x) with a power function of the type

/12 (x) = ln a ax,0 < x < b, a > 0 y12V 7 ab — 1

(54)

and zero elsewhere. We would like to have a slow rising /12 (x). This can be achieved by taking a near 1 such as a = 1.01. The hazard function in this case is the following, denoting

/*(x) = a1 /10(x) + a2/12(x) and S*(t) = Pr{x > t} in /*(x):

h*(x) =

__n_

/*(x) _ a1 a6(n + 1 — q)e—6x[1 + a(q — 1)e"6x] q—1 + a2 J— S*(x) = a1 [1 — [1 + a(q — 1)e—6x]—^ +1] + a2 [<*—£]

ln a-ax

(55)

For convenience, consider the pathway model /10 (x) = a+xx) ,0 < x < to, a > 0, p > 1 and zero elsewhere. Then

a1 [a(p — 1)( ah — 1)] + a2 [(ln a) ax (1 + a x )p ]

h*(x) =

(56)

a1 (1 + ax)(ab — 1) + a2 (ab — ax)(1 + ax)P '

From here we can get bathtub shaped curves for the hazard function. This is plotted for the following sets of parameters.

Figure 16

Figure 16: a1 = 1, a2 = 1, a = 1, p = 2, b = 100, a = 1.01 Figure 17: a1 = 2, a2 = 1, a = 1, p = 2, b = 100, a = 1.01

Figure 17

Figure 18: a1 = 2, a2 = 1, a = 2, p = 2, b = 100, a = 1.01

6. Moments and Laplace Transforms

Let us consider the special case where y = 5 — 1 in (1),(2),(3) with x > 0. This situation is more relevant to reliability analysis. In this case the pathway models are the following:

K1 (x) = a5(n + 1 — q) x5—1 [1 — a(1 — q) x5 ]T—q (57)

for q < 1, a > 0,5 > 0, n > 0.n + 1 — q > 0,1 — a (1 — q)x5 > 0.

K2 (x) = a5 (n + 1 — q) x5—1 [1 + a(q — 1) x5 ] — ^ (58)

for q > 1, a > 0,5 > 0, n > 0, n + 1 — q > 0, x > 0.

K3 (x) = a5nx5—1e—anx5, a > 0, n > 0,5 > 0. (59)

As illustrated before, in (58) and (59) both x and 1 belong to the same family of distributions, namely type-2 pathway model and generalized gamma model respectively. Consider arbitrary moments E(xh) for a complex number h. This is available from the type-1 beta integral, type-2 beta integral and gamma integral respectively, and they are the following, denoted by E(j) (xh), j = 1,2,3:

_l iMT(JL-

E(1)(xh) = a(n + 1 — q) r(I + 1)r(T— + 1)

xh) = -^--—,--, (60)

) [a (1 — q)]1 +1 r( | + j—q + 2) , ( )

for a < 1, №(| + 1) > 0.

5

E(2,(xh ) = r( 1+^ q— -1 -1) (6D

[ a(q — 1)]1+1 r( q—1)

for q > 1, &(I + 1) > 0, q—Y — h — 1) > 0 or —5 < K(h) < q—T + 5. Thus, only a few moments will exist here. But the above strip of analyticity is sufficient to compute the density via the inverse Mellin transform.

r(h + 1) h

E(3)( xh ) = , a > 0,5 > 0,n > 0,&(- + 1) > 0. (62)

(an)1 5

Let the Laplace transforms with Laplace parameter t be denoted by Lg1 (t), Lg2 (t), Lg3 (t) respectively.

^a(1—q)] — 1 tx 5 ir , N 5lJL-

e—tx„5—1\^„('\ „w5i 1—q

zi[a(1—q)] 5 n

Lg1 (t) = a5(n + 1 — q) / e—txx5—1 [1 — a(1 — q) x5 ] >—q dx

J 0

= n +1 - q) f 1e— [A]1 (1 — u) 1—q du (1 — q) -/0 v ;

Expanding the exponential part and integrating term by term we have the following:

l§1 (t)=(n +1 — q) E (—t)k k !1 ul (1 — u) 1—qdu

g1 () (1 — q) k=0 k! [a(1 — q)]1 J0 ( )

(n + 1 — q)r( 1— + 1) ^ r(k + 1)

E

t

_[a(1 — q)] K

(63)

(1 — q) |=0 k!r(l + 1— + 2)

If 1 = m, m = 2,3,... then we can expand the gammas by using the multiplication formula for gamma functions, namely,

/— 1 1 m — 1

r(mz) = mmz—2r(z)r(z + -)...r(z + m—) (64)

and then (63) can be written as a hypergeometric series of the type m Fm.

Now, consider L§2 (t). Since only a few moments exist the Laplace transform or moment generating function does not exist here. But for truncated case the Laplace transform will exist. Let the right

1 r

tail be truncated out from x = b onward where b < [a(q — 1)] — 6 so that 0 < a(q — 1)x6 < 1. In this case the truncated density is the following:

§* (x) = (nP+{x < t6 x6—1[1+a(q — 1)x6 ]—*,0 < x < b

and zero elsewhere. Then the Laplace transform in this truncated case is the following:

L§*(t) = ('P+,x<t}6 l!e—'xx6—1 [1 + a(q — 1)x6] — *dx

r = uq—T)k[^][a(q — 1)]k¡y^-'e—-dx

= (1^ U q—1 )k^ [^ ]'Y(k6+^ bt) (65)

for | a(qt—1) | < 1 where y(£, c) = /0c xe—xdx is the incomplete gamma function. In the truncated case of §*(x) we can also derive arbitrary moments. Let c = Pr{x < b}. Then

E(xh) = (n + 1c— q)a6 Jbxhx6—1 [1 + a(q — 1)x6] — Adx

(n + 1 — q) 1 f a(q—1)b hM ,-1, = —(-^-U us [1 + u] q—1 du.

c(q — 1 [a(q — 1)] hs J0

Since 0 < u < 1 we can expand the binomial part. (1 + u) q—1 = X^=0( )k— uk. Then

E(xh) = (n+1 — q) 1 h £( jl)k(—t [a(q—1)b\h+kdu ( ) c(q — 1) [a(q — 1)]it* q — 1 )k k! ^

= (n +1 — q) bh6[a(q— 1)b6] _F ( n h +1; h + 2. a(q 1)b6)

= c(q — 1) b (h + 6) 2F1 (q—1, 6 + 1'5 + 2' —a(q — 1)b ).

We already have 0 < a(q - 1)b < 1 and hence the 2F1 series is convergent. Note that for h = 0 the 2F1 reduces to a 1 F0 multiplied by a constant which is equal to n+r—1 which when multiplied by the remaining factor gives 1.

The Laplace transform for §3 is given by the following:

c to

e—txx6—1e—

/•TO 6

,(i) = a6n Jo e—txx6—1e—anx dx

k=0 k!

0

k

(an) ^

kt r(- + 1), for |-r| < 1,6 > 1. (66)

6 (an)s

k

t

7. Component Failures

In engineering fields, a system or network is described as a collection of parts or components. Usually a system is represented by as a network in which the system components are connected together either in series, parallel or a combination of these. Consider a system consisting of k components, these components acting independently. Let the life times of these k components be denoted by x1,..,Xk. Suppose that the system fails if any component fails. Let the system failure time be denoted by x. Then x = min{x1,...., Xk}. Then Pr{x > t} is given by the product of the probabilities Pr{xj > t}, j = 1,..., k because if the smallest is > t then all are > t. The corresponding distribution function of x is

k

P{x < t} = 1 - P{x > t} = 1 - nP{xj > t} = 1 - [P{xj > t}k]

i=1

= 1 - [1 - P{x1 < t}]k, (67)

when xi,..., xn are independently and identically distributed. Then the survival function for x,

denoted by Sk(x), is the following:

kk

Sk(t) = Pr{x > t} = nPr{xj > t} = nSxj(t). (68)

j=1 j=1

Several generalized statistical models are developed by using the formula given in (67) (i.e, the concept of series system) for details see Pascoa et al. [10]. In (1980), Kumaraswamy [11] proposed a two-parameter distribution on (0,1), so called Kumaraswamy distribution. This is contained in the pathway model for q < 1. This type of generalizations contains distributions with unimodal and bathtub shaped hazard functions, see Cordeiro and de Castro [12] and Jones

[13]. These generalized models include Kumaraswamy-Weibull distribution by Cordeiro et al.

[14], Kumaraswamy-Gumbel distribution by Cordeiro et al. [15], Kumaraswamy-generalized gamma distribution by Pascoa et al. [10], Kumaraswamy-log-logistic distribution by Tiago et al. [16], Kumaraswamy-modified Weibull distribution by Cordeiro et al. [17], Kumaraswamy-half Cauchy distribution by Gosh [18], Kw-generalized Rayleigh distribution by Antonio et al. [19] and Kumaraswamy-Gompertz distribution by Rocha et al. [20].

Let the life times be pathway distributed. For convenience let us take the case where the pathway parameter q > 1or1 < q < q + 1. Then from (27)

Sxj(t) = Pr{xj > t} = [1 + a(qj - 1)tbj]bj+1 (69)

where bj = - , bj + 1 = - nj+-1qj. The density of x, denoted by fx (t) is available from Sx (t) by differentiation. That is,

d d k fx (t) = - dtSx (t) = - dt n[1+a- (qj - 1)t'j ]bj+1

k d k : E[-(bj + 1)][1 + «j(qj - 1)t'j]bj [-d-(1 + «j(q; -1)^)] n [1 + «i(qi - 1)t'-'] j = 1 dt i=j=1 kk

= E[ n (1 + «i(qi-1)^)bi+1][1 + «j(q; - 1)t'j][«j(q; -1 j--1 ]. (70)

j=1 i=j=1

Hence the hazard function for x is given by the following:

kk

E (n +1 - qj)[«^t'j-1 ][ n (1 + «i(qi -1)^)bi+1 ][1 + «j(qj -1)^]

hx (t) = j=-^^^- (71)

n(1 + «j(qj - 1)t^j)fcj+1 j=1

This is a very interesting form. Note that in (7.4) all the three families of functions, namely the generalized type-1 beta, type-2 beta and gamma families are involved. For different j, qj can be qj < 1 or qj > 1,1 < qj < nj + 1 or qj ^ 1. Hence we can consider many special cases of various types. For example, let k = 2 and q2 ^ 1. Then

hx (t) = (ni + 1 - qi )

S1 a1t51

-i

[1 + a1 (q1 - 1)t51 ]b1

+ 52 a2 n2152-1.

(72)

The graph for the following parameter values is given below.

Figure 19: Plots of hx (t) for b1 = 1 and different values of the other parameters

8. Estimation of ReLiabiLity Function

Let x1,..., xn denote a random sample of size n from the pathway model (21) with parameters a, q,6, n. From (21), the logarithmic likelihood function is

n n n 6

L(a,6,n,q;y) = nlna + nln6 + nln(n + 1 — q) + (6 — 1) ^lnyi---- ^ln(1 + a(q — 1)y6).

q1

i=1

i=1

(73)

First of all we differentiate (73) with respect to all unknown parameters and equate these differential equations to zero. The MLEs of the unknown parameters are obtained on solving these differential equations simultaneously. Let a, q, nj and 6 be the MLEs of the a, q, n, and 6 respectively.

Theorem 1. The MLE of S2 (t) is given by

S2(t) = [1 + a (q - 1)t5 ]-^+1

(74)

Using the invariance property of the MLEs, we can easily establish the above result. The first derivative of the logarithmic likelihood function relative to the parameters is non-linear, and analytical solutions are difficult to obtain. A constrained optimization method can be used to solve such kinds of equations. This optimization problem can be carried out using constrOptim () or optim () function in R software.

8.1. Simulation Study

Random samples were generated from the pathway model using its distribution function. We considered a random sample of sizes n = 200,400,600 and the procedure was repeated 1000 times. The maximum likelihood estimate was computed using the optim () function in R. The results are given in the following tables, bias and MSE can be observed to decrease as the sample size increases. Similarly, the parameters of the distribution corresponding to q < 1 can also estimated by using this method. We are able to estimate the parameters of all the models proposed in this paper using the same procedure.

Table 1: ML estimate, Bias, MSE's for the parameter value (q,a,n,5) = (1.2,2,3,1)

n q Abs.Bias MSE a Abs.Bias MSE

200 1.2256 0.0256 0.00065 2.1298 0.1298 0.0168

400 1.2134 0.0237 0.00018 2.0237 0.0237 0.0005

600 1.2109 0.0108 0.00011 2.0079 0.0079 0.00006

n n Abs.Bias MSE S Abs.Bias MSE

200 3.1186 0.1186 0.0140 1.015 0.0158 0.0002

400 3.1148 0.1148 0.0131 1.0063 0.0066 0.00004

600 3.0894 0.0894 0.0079 1.0049 0.0049 0.00002

Table 2: ML estimate, Bias, MSE's for the parameter value (q,a, n, 5) = (1.1,2,2,5)

n q Abs.Bias MSE a Abs.Bias MSE

200 1.1360 0.0360 0.0013 2.1013 0.1013 0.0102

400 1.1190 0.0190 0.0003 2.0496 0.0496 0.0024

600 1.1113 0.0113 0.00012 2.0433 0.0433 0.0018

n n Abs.Bias MSE S Abs.Bias MSE

200 2.0831 0.0831 0.0069 5.0953 0.0953 0.00908

400 2.0397 0.0397 0.0015 5.0466 0.0466 0.00217

600 2.0148 0.0148 0.0002 5.0344 0.0344 0.00118

Table 3: ML estimate, Bias, MSE's for the parameter value (q,a, n, 5) = (1.2,4,3,5)

n q Abs.Bias MSE a Abs.Bias MSE

200 1.2423 0.0423 0.0017 4.3068 0.3068 0.0941

400 1.2243 0.0243 0.0005 4.0791 0.0791 0.0062

600 1.2155 0.0155 0.0002 4.0746 0.0746 0.0055

n n Abs.Bias MSE S Abs.Bias MSE

200 3.1272 0.1272 0.0161 5.0793 0.0793 0.0062

400 3.1246 0.1246 0.0155 5.0331 0.0331 0.0010

600 3.0505 0.0505 0.0025 5.0247 0.0247 0.00061

9. Reliability for Arbitrary Values of the Parameters

In this section we discussed some simple numerical results for illustrating the theory developed in this paper. The reliability measurement of the system was obtained for arbitrary values of parameters related to the number of components, running time, etc. The reliability behavior of the system was graphically observed to identify the best possible configuration of the components with enhanced reliability of the system.

Table 4: Reliability measure for a = 0.5,5 = 0.1, n = 2 and various values of q at time t = 10

No of components q = 1.2 q = 1.3 q = 1.4 q = 1.5 q = 1.6

1 0.3439755 0.3752348 0.4072686 0.4400367 0.4735025

2 0.1183192 0.1408012 0.1658677 0.1936323 0.2242046

3 0.0406989 0.0528335 0.0675527 0.0852053 0.1061614

4 0.0139994 0.019825 0.0275121 0.0374935 0.0502677

5 0.0048155 0.007439 0.0112048 0.0164985 0.0238019

6 0.0016564 0.0027914 0.0045634 0.00726 0.0112702

7 0.0005698 0.0010474 0.0018585 0.0031946 0.0053365

8 0.000196 0.000393 0.0007569 0.0014058 0.0025268

9 0.0000674 0.0001475 0.0003083 0.0006186 0.0011965

10 0.0000232 0.0000553 0.0001255 0.0002722 0.0005665

Table 5: Reliability measure for q = 1.6,5 = 0.1, q = 2 and various values of a at time t = 10

No. Of Components (n) a = 0.02 a = 0.03 a = 0.04 a = 0.05 a = 0.06

1 0.9656186 0.9490587 0.9329014 0.9171343 0.9017451

2 0.9324194 0.9007123 0.870305 0.8411353 0.8131442

3 0.9003615 0.8548288 0.8119088 0.771434 0.7332488

4 0.8694059 0.8112827 0.7574309 0.7075086 0.6612035

5 0.8395145 0.7699549 0.7066083 0.6488803 0.596237

6 0.8106509 0.7307323 0.6591959 0.5951104 0.5376538

7 0.7827796 0.6935078 0.6149648 0.5457961 0.4848267

8 0.7558666 0.6581796 0.5737015 0.5005684 0.4371901

9 0.7298789 0.6246511 0.535207 0.4590884 0.394234

10 0.7047846 0.5928305 0.4992953 0.4210457 0.3554986

Table 6: Reliability measure for q = 1.6, a = 0.06, q = 2 and various values of 5 at time t = 10

No of Components 5 = 0.2 5 = 0.3 5 = 0.4 5 = 0.5 5 = 0.6

1 0.8785595 0.8505635 0.8170973 0.7775798 0.7316045

2 0.7718669 0.7234583 0.667648 0.6046303 0.5352452

3 0.678131 0.6153472 0.5455333 0.4701483 0.3915878

4 0.5957784 0.5233919 0.4457538 0.3655778 0.2864874

5 0.5234268 0.445178 0.3642242 0.2842659 0.2095955

6 0.4598616 0.3786522 0.2976066 0.2210394 0.153341

7 0.4040158 0.3220677 0.2431735 0.1718758 0.112185

8 0.3549519 0.273939 0.1986964 0.1336471 0.082075

9 0.3118464 0.2330025 0.1623543 0.1039213 0.0600465

10 0.2739756 0.1981835 0.1326593 0.0808071 0.0439303

Table 7: Reliability measure for q = 1.6, a = 0.06,5 = 0.6 and various values of q at time t = 10

No of Components q = 3 q = 4 q = 5 q = 6 q = 7

1 0.5852359 0.4681505 0.3744899 0.2995675 0.2396344

2 0.342501 0.2191649 0.1402427 0.0897407 0.0574247

3 0.2004439 0.1026022 0.0525195 0.0268834 0.0137609

4 0.117307 0.0480333 0.019668 0.0080534 0.0032976

5 0.0686522 0.0224868 0.0073655 0.0024125 0.0007902

6 0.0401778 0.0105272 0.0027583 0.0007227 0.0001894

7 0.0235135 0.0049283 0.001033 0.0002165 0.0000454

8 0.0137609 0.0023072 0.0003868 0.000064 0.0000109

9 0.0080534 0.0010801 0.0001449 0.0000194 0.0000026

10 0.0047131 0.0005057 0.0000543 0.0000058 0.0000006

The following are the graphical representations of reliability relative to the number of components n.

Figure 20

Figure 21

Figure 22

Figure 23

The results obtained for arbitrary values of the parameters indicate that reliability of a series system of 10 identical components keeps on decreasing with the increasing number of components. Several distributions are obtained from our paper for particular parameter values and various arbitrary functions P(x). As q ^ 1 the models (20) and (21) become the Weibull distribution, so it can be considered as an extended version of the Weibull distribution.

10. Conclusions

An arbitrary function is introduced in the pathway model, for constructing the hazard functions of desired shapes. In the present study, we conclude that reliability continues to decline as the number of components increases. It is recommended to use the smallest number of components in a series system for better performance. However, the performance of these systems can be enhanced by using components that follow the extended form of Weibull failure laws.

References

[1] Mathai, A. M. (2005): A pathway to matrix-variate gamma and normal densities, Linear Algebra and its Applications, 396: 317-328.

[2] Mathai, A. M. and Provost, S. B. (2006). Some complex matrix variate statistical distributions in rectangular matrices. Linear Algebra and its Applications, 410: 198-216.

[3] Tsallis, C. Introduction to Nonextensive Statistical Mechanics: Approaching a Complex World, Springer, New York, 2009.

[4] Beck, C. and Cohen, E. G. D. (2003). Superstatistics. Physica A, 322: 267-275.

[5] Mathai, A. M. and Provost, S. B. (2006a). On q-logistic and related distributions. IEEE Transactions on Reliability, 55: 237-244.

[6] Mathai, A. M. and Haubold, H. J. (2007). Pathway model, superstatistics, Tsallis statistics and generalized measure of entropy. Physica A, 375: 110-122.

[7] Mathai, A. M. and Rathie, P. N. Basic Concepts in Information Theory and Statistics: Axiomatic Foundations and Applications, Wiley Halsted, New York and Wiley Eastern, New Delhi, 1975.

[8] Mathai, A. M. An Introduction to Geometrical Probability: Distributional Aspects with Applications, Gordon and Breach, Newark, USA, 1999.

[9] Mathai, A. M. and Princy, T. (2017). Analogues of reliability analysis for matrix-variate cases. Linear Algebra and its Applications, 532: 287-311.

[10] Pascoa, A. R. M., Ortega, E. M. M. and Cordeiro, G. M. (2011). The Kumaraswamy generalized gamma distribution with application in survival analysis. Statistical Methodology, 8: 411-"433.

[11] Kumaraswamy, P. (1980). Generalized probability density function for double-bounded random processes. Journal of Hydrology, 462: 79-88.

[12] Cordeiro, G. M., Castro, M. de. (2011). A new family of generalized distributions, Journal of Statistical Computation and Simulation. 81: 883-898.

[13] Jones, M. C. (2009). Kumaraswamy's distribution: a beta-type distribution with some tractability advantages, Statistical Methodology. 6: 70-81.

[14] Cordeiro, G. M., Ortega, E. M. M. and Nadarajah, S. (2010). The Kumaraswamy Weibull distribution with application to failure data, Journal of the Franklin Institute, 347: 1399-1429.

[15] Cordeiro, G. M., Nadarajah, S. and Ortega, E. M. M. (2011). The Kumaraswamy Gumbel distribution. Statistical Methods and Applications. 21: 139-168.

[16] Tiago Viana Flor de Santana, Edwin M. M. Ortega, Gauss M. Cordeiro and Giovana O. Silva, (2012). The Kumaraswamy-log-logistic distribution. Journal of Statistical Theory and Applications. 11: 265-291.

[17] Cordeiro, G. M., Ortega, E. M. M. and Silva, G. O. (2012). The Kumaraswamy modified Weibull distribution: theory and applications. Journal of Statistical Computation and Simulation. 84: 1387-1411.

[18] Ghosh, I. (2014). The Kumaraswamy-half-Cauchy distribution: properties and applications. Journal of Statistical Theory and Applications. 13: 122-134.

[19] Antonio E. Gomes, Cibele Q. da-Silva, Gauss M. Cordeiro and Edwin M.M. Ortega (2014). A new lifetime model: the Kumaraswamy generalized Rayleigh distribution. Journal of Statistical Computation and Simulation, 84: 290-309.

[20] Rocha,R., Nadarajah, S., Tomazella, V., Louzada, F. and Eudes, A. (2015). New defective models based on the Kumaraswamy family of distributions with application to cancer data sets. Statistical methods in medical research, 0962280215587976.