Generalization of the Weibull distribution: the odd Weibull family

Use load(“LungCancerData.RData”) to load the posted Lung Cancer Data. Use “time” for survival time. A failure density function (model) is assigned to each group (find in the table below).( Odd Weibull (OW))

1. Study the parameters of the assigned time to failure density function. Each parameter must be investigated separately while you fix the other parameters. Visualize the model while you change one parameter at a time.

2. Estimate the parameters of the model using MLE (Maximum Likelihood Estimation) for the Lung Cancer dataset as follows.

a) For the entire data.

b) For Female Patients.

c) For Male Patients.

3. Estimate a set of two parameters at a time by fixing the other parameters and repeat this until you study all possible combinations of two parameters.

4. Using the estimated parameters in part 2, find and visualize the reliability and hazard for groups a, b, and c in question 2.

Presentation must include brief theory, simulations, observations, and conclusions. Please submit a report (pdf), your presentation (ppt), and the code.

Mathematical Sciences (2019) 13:105–114
https://doi.org/10.1007/s40096-019-0283-7
ORIGINAL RESEARCH
A new uniform distribution with bathtub‑shaped failure rate
with simulation and application
Jamal N. Al abbasi1 · Mundher A. Khaleel2 · Moudher Kh. Abdal‑hammed3 · Yue Fang Loh4 · Gamze Ozel5
Received: 11 October 2018 / Accepted: 8 April 2019 / Published online: 20 April 2019
© The Author(s) 2019
Abstract
Bathtub failure rate shape is widely used in industrial and medical applications. In this paper, a three-parameter lifetime
distribution, so-called the generalized Weibull uniform distribution that extends the Weibull distribution, is proposed and
studied. This distribution has bathtub-shaped or decreasing failure rate function which enables it to fit real lifetime data sets.
Various structural properties of the new distribution are derived, including explicit expressions for the quantile function,
moments, moment-generating function and order statistics. Parameter estimations are provided by a maximum likelihood
estimation, and the performance of the maximum likelihood estimation is evaluated using a simulation study. An application
to real-life data demonstrates that the proposed distribution can be very useful in fitting real data.
Keywords Weibull distribution · Reliability · Failure rate · Bathtub shape · Quantile function
Mathematics Subject Classification 62E15 · 62F10 · 62P30
Introduction
In recent years, many classical distributions have been generalized by adding more shape parameters since numerous
application in the field of engineering, financial, biomedical
and environmental sciences indicated that classical distributions are not suitable to explain the data sets. Hence, there
is continual need for the extension of these distributions to
make effective progress in the application of these areas.
The Weibull distribution with two parameters is a flexible
* Gamze Ozel
gamzeozl@hacettepe.edu.tr
1
Department of Statistics, Al Nahrain University, Baghdad,
Iraq
2
Department of Mathematics, Faculty of Computer Science
and Mathematics, University of Tikrit, Tikrit, Iraq
3
Department of Business Administration, Faculty
of Administration and Economics, University of Tikrit,
Tikrit, Iraq
4
Department of Actuarial Science and Applied Statistics,
Faculty of Business and Information Science, UCSI
University, Kuala Lumpur, Malaysia
5
Department of Statistics, Hacettepe University, Ankara,
Turkey
distribution to model different types for lifetime data that
exhibit monotone-shaped failure rates. The Weibull distribution can be specified through its cumulative distribution
function (cdf) given by
𝜔
G(Z) = 1 − e−𝜐Z , Z, 𝜐, 𝜔 > 0.
(1)
where 𝜐, 𝜔 refer to the scale and shape parameters, respectively. According to Eq. (1), the probability density function
(pdf) of the Weibull distribution is given as follows:
𝜔
(2)
From (1) and (2), the failure rate function (FRF) is obtained
as
g(Z) = 𝜐𝜔Z 𝜔−1 e−𝜐Z .
h(Z) = 𝜐𝜔Z 𝜔−1 .
(3)
Indeed, when 𝜔 > 1, the FRF is increasing, whereas 𝜔 < 1, it is decreasing, and when 𝜔 = 1, it is constant. The Weibull distribution, on the other hand, is incapable to model data that exhibit bathtub-shaped or unimodal-shaped FRFs. The later can be considered as a common situation in industrial and medical applications, for example, electric machine life cycles, human mortality, life cycle of cancer patients, etc. Furthermore, the Weibull distribution may not be a suitable distribution for many of the engineering applications requiring load-bearing capability in the area of material science like strength of brittle materials which include most 13 Vol.:(0123456789) 106 Mathematical Sciences (2019) 13:105–114 ceramics, glass, optical fiber and some polymers, and in other numerous applications. However, the strength values are limited where the Weibull distribution requires an unlimited value. Therefore, it is necessary to find generalizations of the Weibull distribution. Several trials have been conducted to define new generalizations of the Weibull distribution by adding several additional shape parameter(s), for instance, see [7, 13, 15, 19]. Broadly speaking, for bringing more flexibility to baseline distribution can be added an extra parameter, so it can be more benefit in exploring tail properties and exhibiting bathtub-shaped failure rate function. Some models for bathtub and unimodal-shaped FRFs are summarized in Table 1. Note that the generalization is also complex and causes estimation problems, especially when the number of parameters is four or more. In this study, the main purpose is to obtain the more flexible distribution with a few parameters. Recently, a family of distributions has been proposed by Bourguignon [3] called as the Weibull-G “W-G” family by using the Weibull G(t;𝜼) distribution as in (1) and replacing the argument t with G(t;𝜼) , ̄ ̄ where G(t;𝜼)) = 1 − G(t;𝜼)). Its cdf is defined as −𝜐 F(t;𝜐, 𝜔, 𝜼) = 1 − e [ G(t;𝜼) ̄ G(t;𝜼) ]𝜔 , t ∈  ⊆ ; 𝜐, 𝜔 > 0,
(4)
where G(t;𝜼) refers to a parent cdf, which relies on a parameter vector 𝜼 . The corresponding pdf to (4) is
] [ ]𝜔
G(t;𝜼)𝜔−1 −𝜐 G(t;𝜼)
̄
f (t;𝜐, 𝜔, 𝜼) = 𝜐𝜔g(t;𝜼)
.
e G(t;𝜼)
𝜔+1
̄
G(t;𝜼)
[
Assuming that the parent distribution is uniform in the interval (0, 𝜑), 𝜑 > 0, then, the cdf is given by
t
,
𝜑
G(t;𝜑) =
(6)
0 < t < 𝜑 < ∞. Then, the corresponding pdf becomes 1 , 𝜑 g(t;𝜑) = (7) 0 < t < 𝜑 < ∞. Inserting (6) in (4) yields a three-parameter Weibull uniform (WU) distribution with the cdf given by [ ( F(t) = 1 − exp − 𝜐 t 𝜑−t )𝜔 ] , 0 < t < 𝜑 < ∞, 𝜐, 𝜔 > 0.
(8)
This distribution is also known as the Kies distribution [8].
The pdf of the WU distribution is given by
f (t) =
𝜐𝜔𝜑
(𝜑 − t)2
(
t
𝜑−t
)𝜔−1
[ (
exp − 𝜐
t
𝜑−t
)𝜔 ]
.
(9)
More recently, Cordeiro et al. [6] proposed a generalized
Weibull-G “GW-G” family of distributions by substituting the
argument t with − loge {1 − G(t;𝜼)} in Eq. (1) and defined the
cdf of their family of distributions by
𝜔
F(t;𝜐, 𝜔, 𝜼) = 1 − e−𝜐[− loge {1−G(t;𝜼)}] ,
t ∈  ⊆ ;
𝜐, 𝜔 > 0.
(10)
Then, the pdf of the GW-G family is given by
(5)
Table 1  Some models for bathtub-shaped or unimodal-shaped FRFs
Distribution
Distribution function, F(t)
Failure rate function, h(t)
Weibull uniform (WU), Kies [8]
)𝜔 ]
t
,
1 − exp − 𝜐 𝜑−t
𝜐𝜔𝜑
(𝜑−t)2
Exponentiated Weibull (EW), Mudholkar and Srivastava [11]
(
𝜔 )𝜆
1 − e−𝜐t ,
Additive Weibull (AW), Xie and Lai [18]
Modified Weibull (MW), Lai et al. [9]
[
(
𝜔
00
−𝜐t𝜔 e𝜆t
Odd Weibull (OW), Cooray [4]
1−e
,
t>0
[
( 𝜐t𝜔
)𝜆 ]−1
1− 1+ e −1
,t > 0
Flexible Weibull (FW), Bebbington et al. [2]
1 − e−e
SZModified Weibull (SZMW), Sarhan and Zaindin [14]
Kumaraswamy Weibull (KwW), Cordeiro et al. [5]
𝜐t− 𝜔t
𝜔
1 − e−𝜐t −𝛾t ,
t>0
[
]𝛾
(
𝜔 )𝜆
−𝜐t
1− 1− 1−e
, t>0
Weibull Pareto (WP), Alzaatreh et al. [1]
1−e
New Weibull Pareto (NWP), Tahir et al. [17]
1−e
Weibull power (WPo), Tahir et al. [16]
1−e
13
t>0
,
[
( )]𝜔
−𝜐 loge 𝜑t
−
[( )𝜐
−𝜐
t
𝜑
[
]𝜔
−1
t𝜆
𝜑𝜆 −t𝜆
]𝜔
,
,
,
t>𝜑>0
{
t
𝜑−t
}𝜔−1
𝜔
𝜔
𝜐𝜔𝜆t𝜔−1 e−𝜐t (1−e−𝜐t )𝜆−1
1−(1−e−𝜐t𝜔 )𝜆
𝜔−1
𝜆−1
𝜐𝜔t
+ 𝛾𝜆t
𝜐(𝜔 + 𝜆t)t𝜔−1 e𝜆t
𝜔
𝜔
𝜔−1 𝜐t
𝜐𝜔𝜆t
e (e𝜐t − ]1)𝜆−1
[
−1
𝜔
× 1 + (e𝜐t − 1)𝜆
(
)
𝜔
𝜐 + t𝜔2 e𝜐t− t
𝜐𝜔t𝜔−1 + 𝛾
𝜔
𝜔
𝜐𝜔𝛾𝜆t𝜔−1 e−𝜐t (1−e−𝜐t )𝜆−1
1−(1−e𝜐t𝜔 )𝜆
𝜐𝜔
t
[
𝜐 loge
( )]𝜔−1
t
𝜑
[( )𝜐
]𝜔−1
−1
[
]𝜔
t>𝜑>0
𝜐𝜔t𝜐−1
𝜑𝜐
0 0.
(14)
Here, the shapes of the pdf and FRF of the GWU distribution are discussed. The derivative loge f (t) with respect to t
is given by
(15)
So, the roots of Eq. (15) represent the critical points of the
pdf as follows:
{
)}𝜔
(
)
(
t
t
+ loge 1 −
+ (𝜔 − 1) = 0.
𝜐𝜔 − loge 1 −
𝜑
𝜑
(16)
Herein, there may be more than one root of Eq. (16). If t = t0
that is a point of inflexion, a local maximum or a local mini2
mum relying on whet her d 2 loge f (t;𝜐, 𝜔, 𝜑) = 0 ,
dt
d2 log f (t;𝜐, 𝜔, 𝜑) < 0 , d2 log f (t;𝜐, 𝜔, 𝜑) > 0 , where
e
e
dt 2
dt 2
d2 log f (t;𝜐, 𝜔, 𝜑) is given by
e
dt 2
{
(
)}{
(
)
}
{
(
)}𝜔
t
t
t
loge 1 −
loge 1 −
− (𝜔 − 1) − 𝜐𝜔 − loge 1 −
𝜑
𝜑
𝜑
)
}
(
{
t
+ (𝜔 − 1) − (𝜔 − 1)
× loge 1 −
𝜑
)2 {
( )}2
(
𝜑2 1 − 𝜑t
loge 𝜑−t
𝜑
It is difficult to determine the behavior of the density function from Eq. (16); despite this, in Fig. 1, some of the possible
shapes are illustrated of the pdf when 𝜑 = 2.5 and for different values of 𝜐 and 𝜔. For fixed 𝜔 < 1, the density is always U-shaped, and as 𝜐 tends to 0 (∞), the density function f(t) tends to J (reversed-J) shape, for fixed 𝜔 > 1, the right (left)skewed increases as 𝜐 decreases (increases), and when 𝜔 = 1,
the density function is monotonically decreasing if 𝜐 > 1,
constant if 𝜐 = 1 and monotonically increasing if 𝜐 < 1. Some calculation indicates the limiting behaviors of the pdf, and the limit of the f (t;𝜐, 𝜔, 𝜑) as T → 𝜑− is 0. Further, the limit of f(t) as T → 0+ is given by ⎧∞ ⎪ lim f (t;𝜐, 𝜔, 𝜑) = ⎨ 0 t→0+ ⎪ 𝜑𝜐 ⎩ if 0 < 𝜔 < 1 if 𝜔 > 1
if 𝜔 = 1
Shape of failure rate function
Based on Eq. (14), the cumulative hazard function H(t) is simply given by
13
108
Mathematical Sciences (2019) 13:105–114
)}𝜔
{
(
t
.
H(t;𝜐, 𝜔, 𝜑) = 𝜐 − loge 1 −
𝜑
The FRF has other important system’s quantity characterizing life phenomena and characteristic of interest of a random
variable defined by
h(t;𝜐, 𝜔, 𝜑) =
f (t;𝜐, 𝜔, 𝜑)
,
R(t;𝜐, 𝜔, 𝜑)
h(t) of the GWU distribution takes the following form:
h(t;𝜐, 𝜔, 𝜑) =
𝜐𝜔
d
H(t;𝜐, 𝜔, 𝜑) = (
)
dt
𝜑 1− t
{
)}𝜔−1
(
t
.
− loge 1 −
𝜑
𝜑
{ O b v i(o u s l y,)} b e c a u s e
− loge 1 − 𝜑t
𝜔−1
through
the
(17)
term
, the shape of FRF relies only on the
parameter 𝜔, also there is no direct effect on the shapes by the
remaining two parameters. In order to discuss the shapes of
h(t), with respect to t, the derivative of loge h(t) is given by
⎡
⎤
⎢
⎥
𝜔−1
d
1
loge h(t;𝜐, 𝜔, 𝜑) = �
� ⎢1 −
� ⎥. (18)
�
dt
𝜑 1 − 𝜑t ⎢
loge 1 − 𝜑t ⎥
⎣
⎦
So, the roots of Eq. (18) represent the critical points of the
FRF as follows:
)
(
t
− (𝜔 − 1) = 0
loge 1 −
(19)
𝜑
One can obtain the turning point of the FRF t0 by solving
Eq. (19) as follows:
]
[
t0 = 𝜑 1 − e𝜔−1
(20)
It can be shown from Eq. (17) that when 𝜔 ≥ 1, h(t) tends
to increase in t, and when 0 < 𝜔 < 1, h(t) tends to decrease for t < t0 and increase for t > t0. This implies a bathtub shape for
the FRF. However, at a minimum value t = t0 for the h(t), the
value of t0 decreases as 𝜔 increases, and h(t0 ) can be obtained
after simplification using Eqs. (17) and (20) as
h(t0 ;𝜐, 𝜔, 𝜑) =
Fig. 1  Plots of f(t) for selected parameter values of GWU distribution
13
𝜐𝜔(1 − 𝜔)𝜔−1
.
𝜑e𝜔−1
Although bathtub-shaped FRF are common in reliability
and survival analysis, these properties are related to several
lifetime distributions, as the GWU distribution.
The limiting behaviors of the FRF can be easily shown from
Eq. (17), while the limit of the as h(t) as T → 𝜃 − is ∞. Further,
the limit of h(t) as T → 0+ is given by
Mathematical Sciences (2019) 13:105–114
⎧∞
⎪
lim h(t;𝜐, 𝜔, 𝜑) = ⎨ 0
t→0+
⎪ 𝜑𝜐
⎩
109
if 0 < 𝜔 < 1 if 𝜔 > 1
if 𝜔 = 1
Figure 2 displays several of the possible shapes of
h(t;𝜐, 𝜔, 𝜑) for selected values of (𝜐, 𝜔, 𝜑).
Figure 2 reveals that the proposed distribution can produce
flexible FRF shapes such as increasing, decreasing, bathtub,
J, snf reversed-J. This fact implies that the GWU distribution
can be very useful to fit different data sets with various shapes.
Now, suppose Z is a random variable having the Weibull
distribution in Eq. (1) with parameter 𝜐 and 𝜔, the mgf of Z
is given by
(
)
v
∞ sv Γ 1 +
∑
𝜔
(22)
MZ (s) =
v
v!𝜐 𝜔
v=0
Transformation
Since the relationship between the random variable Z has
the Weibull distribution in Eq. (1) and the random[ variable
]
T that follows GWU distribution given by T = 𝜑 1 − e−Z ,
the first non-central moment of T is given by
[
]
[
]
𝜇1� = 𝜑 1 − E(e−Z ) = 𝜑E 1 − MZ (− 1)
[
]
The random variable T = 𝜑 1 − e−Z , when a random variable Z is Weibull-distributed with scale parameter 𝜐 > 0 and
shape parameter 𝜔 > 0, follows the GWU(v, w, 𝜑) distribution
that helps us to simulate the random variable T in easy way
by first[ simulating
] the random variable Z and then calculate
T = 𝜑 1 − e−Z which has the GWU distribution.
�
�
�
�
⎡
⎡∞
v ⎤
v+1 Γ 1 + v ⎤
∞ (− 1)v Γ 1 +
(−
1)
�
�
⎢
⎢
𝜔 ⎥
𝜔 ⎥
𝜇1� = 𝜑⎢1 −
n
n
⎥ = 𝜑⎢
⎥
𝜔
𝜔
v!𝜐
v!𝜐
v=0
⎢
⎥
⎥
⎢ v=1
⎣
⎦
⎦
⎣
Moment
Assume Z is a random variable possessing the moment-generating function (mgf) MZ (s) = E(esZ ). One can express MZ (s)
in a Maclaurin series expansion as
MZ (s) =
∞ v �
∑
s 𝜇v
v=0
(21)
v!
v
so that 𝜇v′ is the coefficient of sv! in expansion (21).
Setting s = − 1 in (22), we have MZ (− 1), and hence, the first
non-central moment of T becomes
(23)
So, the rth non-central moment of the GWU distribution
can be obtained directly from Eq. (22) for s = − r as
[
]
[
]
𝜇r� = 𝜑r 1 − E(e−rZ ) = 𝜑r 1 − MZ (−r)
)
(
v
∞ (− 1)v+1 r v Γ 1 +
∑
𝜔
.
= 𝜑r
n
v!𝜐 𝜔
v=1
(24)
Further, the central moments 𝜇r of T are obtained from
Eq. (24) as follows:
r
∑
( )
r
𝜇r =
(− 1)r 𝜇i� 𝜇r−i ,
i
i=0
(25)
For r = 2, 3, 4 , Eq. (25) gives 𝜇2 , 𝜇3 and 𝜇4 , which are of
the most interest because of the relations between variance,
kurtosis and skewness of the random variables Z. Note that
𝜇0� = 1 and 𝜇1� = 𝜇. Table 2 displays the numerical values for
mean, variance, skewness and kurtosis for several values of
parameters for the GWU distribution.
Moment‑generating function
The result in (25) can be used in (21) to rewrite the mgf as
(
)
w
∞
∞ (− 1)w+1 vw Γ 1 +
v
v
∑
∑
𝜔
𝜑s
(26)
MT (s) =
w
v!
w!𝜐 𝜔
v=0
w=1
Fig. 2  Plots of h(t;𝜐, 𝜔, 𝜑) for selected parameter values of GWU distribution
For this section, Eq. (26) represents
[ ] the main result. The
characteristic function Φ(s) = E eisT of T can be √
obtained
from this equation with is in place of s, where i = − 1.
13
110
Mathematical Sciences (2019) 13:105–114
Table 2  Mean, variance, skewness and kurtosis for some values of 𝜐
and 𝜔 with 𝜑 = 1 (*:undefined)
𝜐
𝜔
Mean
Variance
Skewness
Kurtosis
1.0
1.2
1.5
2.0
3.0
4.0
5.0
8.0
10.0
1.0
1.2
1.5
2.0
3.0
5.0
10.0
1.0
1.2
1.5
2.0
5.0
10.0
1.0
1.2
1.5
2.0
5.0
10.0
1.0
1.5
2.0
5.0
10.0
1.0
1.5
2.0
5.0
10.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.2
1.2
1.2
1.2
1.2
1.2
1.2
1.5
1.5
1.5
1.5
1.5
1.5
2.0
2.0
2.0
2.0
2.0
2.0
5.0
5.0
5.0
5.0
5.0
10.0
10.0
10.0
10.0
10.0
0.500
0.454
0.400
0.333
0.250
0.200
0.167
0.111
0.091
0.512
0.471
0.422
0.361
0.284
0.203
0.123
0.527
0.492
0.449
0.395
0.250
0.170
0.546
0.517
0.482
0.438
0.313
0.237
0.592
0.563
0.543
0.480
0.435
0.611
0.596
0.586
0.553
0.528
∗
∗
∗
0.389
0.338
0.293
0.258
0.188
0.158
0.430
0.434
0.430
0.413
0.372
0.303
0.208
0.446
0.450
0.449
0.439
0.355
0.271
0.460
0.465
0.467
0.464
0.412
0.350
0.475
0.484
0.489
0.493
0.486
0.473
0.479
0.483
0.492
0.496
∗
∗
∗
∗
0.231
0.245
0.241
0.209
0.187
∗
∗
0.881
0.156
0.223
0.254
0.226
− 0.035
0.009
0.062
0.127
0.252
0.259
− 0.061
− 0.024
0.022
0.080
0.221
0.268
− 0.129
− 0.090
− 0.062
0.029
0.094
− 0.157
− 0.138
− 0.124
− 0.078
− 0.042
∗
∗
∗
∗
∗
0.232
0.240
0.226
0.210
∗
∗
∗
∗
0.199
0.242
0.242
0.121
0.120
0.126
0.143
0.229
0.260
0.127
0.123
0.123
0.131
0.197
0.247
0.148
0.136
0.130
0.126
0.137
0.159
0.151
0.146
0.133
0.127
Setting 𝜆 = 12 in (27) gives the median of the GWU distribution as follows:
} ]]
[
[ {
− loge (2) 𝜔
t 1 = 𝜑 1 − exp −
2
𝜐
(28)
Further, setting 𝜆 = 14 and 𝜆 = 34 gives the lower quartile
and upper quartile, respectively. Simulating the GWU random variable is straightforward. If U is a random number
following a uniform distribution  (0, 1), then the random
variable
} ]]
[
[ {
− loge (1 − u) 𝜔
T = Q(U) = 𝜑 1 − exp −
𝜐
(29)
will follow a GWU distribution with parameters 𝜐 , 𝜔 and 𝜑.
Order statistics
The pdf fj;n (t) of the jth orders T(j) is given by
(
)
n−1
fj;n (t) = n
f (t)F(t)j−1 {1 − F(t)}n−j , j = 1, 2, … , n.
j−1
(30)
Using
j−1
F(t)
)}𝜔 ]j−1
{
(
−𝜐 − loge 1− 𝜑t
[
= 1−e
=
j−1
∑
(
)
(
{
)}𝜔
t
j−1
𝓁 −𝓁𝜐 − loge 1− 𝜑
(− 1) e
𝓁
(31)
𝓁=0
and
n−j
{1 − F(t)}
)}𝜔
{
(
−(n−j)𝜐 − loge 1− 𝜑t
=e
(32)
and substituting (31) and (32) into (30), we get
)}𝜔−1
(
{
(
) j−1
n−1 ∑
𝜐𝜔
t
fj;n (t) = n
(−1)𝓁 (
) 𝜐 − loge 1 −
j − 1 𝓁=0
𝜑
𝜑 1 − 𝜑t
{
(
)}𝜔
−𝜐(n−j+𝓁+1)𝜐 − loge 1− 𝜑t
× e
Quantile function and simulation
The quantile function of the GWU distribution is obtained
by inverting (12) as
} ]]
[ {
− loge (1 − 𝜆) 𝜔
Q(𝜆) = 𝜑 1 − exp −
𝜐
[
13
(27)
(33)
)
(
) j−1 (
(− 1)𝓁
n−1 ∑ j−1
f (t;𝜐n,𝓁 , 𝜔, 𝜑)
fj;n (t) = n
𝓁
j − 1 𝓁=0
(n − j + 𝓁 + 1)
(34)
where f (t;𝜐n,𝓁 , 𝜔, 𝜑) is the GWU PDF with parameters
𝜐n,𝓁 = (n − j + 𝓁 + 1)𝜐 , 𝜔 and 𝜑.
Therefore, the rth non-central moment of the jth OS of
GWU distribution is given by
Mathematical Sciences (2019) 13:105–114
𝜇r� = n
j;n
(
×
n−1
j−1
)
) j−1 (
∑ j−1
𝓁
(
v
∞ (−1)v+1 r v Γ 1 +
∑
𝜔
𝓁=0
(−1)𝓁
𝜑r
(n − j + 𝓁 + 1)
)
111
Setting (40) and (41) to zero and simplifying, we obtain
n
�
̂
𝜐=
(35)
∑n �
i=1
t
− loge 1 − t i
(n)
v
v=1 v!{(n − j + 𝓁 + 1)𝜐} 𝜔
So, the pdf of the largest order statistics T(n) is given by
fn;n (t) = n
n−1
∑
(
𝓁=0
)}w−1
(
{
)
n−1
𝜐𝜔
t
𝓁
(−1) (
) − loge 1 −
𝓁
𝜑
𝜑 1 − 𝜑t
{
(
)}𝜔
−𝜐(𝓁+1) − loge 1− 𝜑t
× e
(36)
fn;n (t) = n
n−1
∑
(
n−1
𝓁
𝓁=0
)
(− 1)𝓁
f (t;𝜐𝓁 , 𝜔, 𝜑)
(𝓁 + 1)
(37)
n
=
∑n �
𝜔
̂
i=1
n
�
(38)
where f (t;𝜐n , 𝜔, 𝜑) is the pdf of the GWU distribution with
parameters 𝜐n = n𝜐 , 𝜔 and 𝜑.
t
− loge 1 − t i
(n)
��𝜔̂
Consider a random sample t1 , t2 , … , tn from GWU(𝜐, 𝜔, 𝜑)
distribution, and assume Θ = (𝜐, 𝜔, 𝜑)T is a vector of the distribution parameter with known 𝜑 (since t < 𝜑), then the loglikelihood function denoted by 𝓁(Θ) can be written as ( ) n ( ) ∑ t 𝓁(Θ) = n loge (𝜐) + loge (𝜔) − loge (t(n) ) − loge 1 − i − (𝜔 − 1) t(n) i=1 ( { ( )} )}𝜔 n n { ∑ ∑ t t × loge − loge 1 − i −𝜐 − loge 1 − i t(n) t(n) i=1 i=1 (39) Since 𝜑 is(supposed known,)the score vector can be defined 𝜕 𝜕 𝓁(Θ), 𝜕𝜔 𝓁(Θ) 𝜕𝜐 T U𝜔 (Θ) = n n ∑ − 𝜐 i=1 { � ��𝜔̂ t − loge 1 − i t(n) (43) U𝜐𝜔 (Θ) = − n 𝜐2 n ∑ i=1 (44) { )}𝜔 )} { ( ( t t loge − loge 1 − i − loge 1 − i t(n) t(n) (45) U𝜔𝜔 (Θ) = − n ∑ n −𝜐 2 𝜔 i=1 )}𝜔 [ { ( )}]2 ( t t loge − loge 1 − i − loge 1 − i t(n) t(n) { (46) For large n, the bivariate normal distribution with zero mean vector 0, and variance-covariance matrix√I2 (𝚯)−1 can ̂ − 𝚯) , be used to approximate the distribution of n(𝚯 where I(𝚯) = limn→∞ n−1 In (𝚯) is the expected. The estimated asymptotic bivar iate nor mal distr ibution ̂ can be used to construct approximate N2 (0, I(𝚯)−1 ) of 𝚯 confidence intervals for the distribution parameters. A large sample intervals for 𝜐 and 𝜔 √ 100(1 − 𝜸)% confidence √ 𝜐 ± z 𝛾 var(̂ ̂ ± z 𝛾 var(𝜔) ̂ , respectively, where are ̂ 𝜐) and 𝜔 2 2 the var(⋅)’s are the diagonal elements of In (𝚯)−1 , and z 𝜸 is 2 𝜸 the upper 2 th percentile of the standard normal distribution. where the components cor- responding to the GWU(𝜐, 𝜔, 𝜑) distribution by differentiating Eq. (39) determined the parameters of new model. The elements of the score model are given by U𝜐 (Θ) = i=1 �� ̂ of 𝜔 . The solution of Eq. (43) represents the MLE 𝜔 𝜐 of 𝜐 can be found by substituting the estimate The MLE ̂ 𝜔 ̂ in Eq. (42). The elements of the observed information matrix In (Θ) = {Irs } for r;s = 𝜐, 𝜔 are found by f1;n (t) = f (t;𝜐n , 𝜔, 𝜑) Maximum likelihood estimation n � � ��� � �� � � n � t t − loge − loge 1 − i × loge − loge 1 − i t(n) t(n) i=1 U𝜐𝜐 (Θ) = − where f (t;𝜐𝓁 , 𝜔, 𝜑) is the pdf of the GWU distribution with parameters 𝜐𝓁 = (𝓁 + 1)𝜐, 𝜔 and 𝜑, and the pdf of the smallest order statistics T(1) is given by: by U(Θ) = (42) ��𝜔̂ ( )}𝜔 t − loge 1 − i t(n) (40) { ( )} ( )}𝜔 n n { ∑ t t n ∑ + −𝜐 − loge 1 − i loge − loge 1 − i 𝜔 i=1 t(n) t (n) i=1 )} ( { ti × loge − loge 1 − t(n) (41) Simulation study A simulation study is conducted to examine the performance of the MLE method to estimate the GWU distribution parameters. Monte Carlo technique is employed to simulate the GWU random number generation depending on the inversion method, i.e., varieties of the GWU distribution are generated using (29). The simulation study is conducted for nine-parameter incorporation, 𝜐 = 0.5, 1.0, 1.5 and 𝜔 = 0.5, 1.0, 1.5 and consider sample sizes with n = 30, 50 and 100. The process is repeated 5000 times for each parameter incorporation and each sample size. 13 112 1 The average bias AB = 5000 Mathematical Sciences (2019) 13:105–114 and i=1 1 ∑5000 mean square error MSE = 5000 i=1 (estimate − actual)2 of ∑5000 ̂ + 2q, AIC = − 2(𝚯) (estimate − actual) ̂ + q ln(n), BIC = − 2(𝚯) the MLEs are used to measure the performance. It can be observed in Table 3 that for all the different sets when we increase n, the bias for 𝜐 and 𝜔 decreases. The MSEs tend to toward zero when we increase the sample size n. The simulation results validate the performance of the maximum likelihood estimators. ̂ refers to the log-likelihood function that evaluwhere (𝚯) ates at the MLE, while q is the number of parameters and n is the sample size. Voltage data The data have been obtained from Meeker and Escobar [10]. They represent the lifetimes of 30 electronic components taken from power-line voltage spikes during electric storms. The data have been previously used by Nadarajah et al. [12] and Tahir et al. [16] (Table 4). The total time on test (TTT) transform plot illustrates that data set has a bathtub-shaped FRF. So, the GWU distribution may be suitable for modeling this data set. Table 5 provides the MLEs of the unknown parameters with corresponding standard errors for different distributions. The values of these statistics in data sets are smaller for the GWU Application Voltage data were employed to illustrate the flexibility of the GWU distribution, in addition to compare the behavior of the new model with other generalization of Weibull distribution. The model selection is achieved by using the maximized log-likelihood, Akaike information criterion AIC and Bayesian information criterion (BIC) given by: Table 3  MLE, bias and MSE for various parameter values when 𝜑 = 1 Sample size Actual values Bias n 𝜐 𝜔 B(̂𝜐) B(𝜔) ̂ MSE(̂𝜐) MSE(𝜔) ̂ 30 0.5 1.0 1.5 0.5 1.0 1.5 0.5 1.0 1.5 0.5 1.0 1.5 0.5 1.0 1.5 0.5 1.0 1.5 0.5 1.0 1.5 0.5 1.0 1.5 0.5 1.0 1.5 0.5 0.5 0.5 1.0 1.0 1.0 1.5 1.5 1.5 0.5 0.5 0.5 1.0 1.0 1.0 1.5 1.5 1.5 0.5 0.5 0.5 1.0 1.0 1.0 1.5 1.5 1.5 0.01219 0.04255 0.07695 0.00979 0.03798 0.08354 0.01502 0.03740 0.06931 0.00828 0.03238 0.04621 0.00687 0.02107 0.05013 0.00641 0.02259 0.04742 0.01051 0.02847 0.02957 0.00561 0.01278 0.02433 0.00513 0.01473 0.02716 0.11761 0.03608 0.02466 0.04648 0.04534 0.04657 0.06757 0.07139 0.06745 0.10939 0.02617 0.01424 0.02832 0.02710 0.02772 0.03924 0.04102 0.04072 0.12317 0.02028 0.00824 0.01289 0.01104 0.01188 0.01681 0.01706 0.01876 0.02083 0.04569 0.10246 0.01468 0.04540 0.10968 0.01519 0.04537 0.10517 0.01239 0.02506 0.05540 0.00904 0.02529 0.05355 0.00943 0.02542 0.05776 0.00644 0.01220 0.02356 0.00472 0.01197 0.02496 0.00458 0.01194 0.02552 0.01009 0.00588 0.00624 0.02446 0.02448 0.02503 0.05393 0.05480 0.05508 0.00586 0.00318 0.00320 0.01350 0.01335 0.01328 0.03009 0.02978 0.03020 0.00326 0.00147 0.00152 0.00632 0.00638 0.00616 0.01469 0.01430 0.01429 50 100 13 MSE Mathematical Sciences (2019) 13:105–114 Table 4  Lifetimes of 30 electronic components from Meeker and Escobar [10] Table 5  MLEs of the parameters (the standard errors are given in parentheses) to the voltage data and the values of the − log likelihood, AIC and BIC statistics 275 106 88 Dist. 113 13 300 247 147 300 28 23 212 143 181 300 300 30 300 23 likelihood Estimates 65 300 300 − Log 10 2 80 300 261 245 173 293 266 AIC BIC 135.29 274.59 278.79 300.03 139.22 282.43 286.63 0.005 (0.001) – 169.92 347.84 353.44 7.701 (0.219) 0.004 (0.0000) – 171.33 352.71 358.32 0.337 (0.057) 2.874 – 300.00 152.58 311.15 315.34 0.006 (0.055) – – – 184.31 368.63 371.43 𝜐 𝜔 𝜆 𝛾 𝜑 GWU​ 0.603 (0.037) 0.624 (0.025) – – 300.03 WU 0.377 (0.102) 0.214 (0.029) – – BW 0.079 (0.030) 0.066 (0.036) 7.936 (0.765) KuW 0.052 (0.024) 0.229 (0.091) WP 0.772 (0.252) W 1.265 (0.252) Bold values show the best model which is suggested in this paper Fig. 3  Voltage data set. a TTT- transform plot represents to the voltage data set. b The empirical and estimated cdf of the distribution. c Probability plots. d The expected and observed number of observations 13 114 distribution compared to those values of other fitted distributions. The approximate 95% two-sided asymptotic confidence interval of the parameters 𝜐 , 𝜔 and 𝜑, respectively, is: (0.391,0.815) , (0.417, 0.831) and (2030.72, 369.34). It is concluded that the proposed model provides a goodness of fit for the given data set, which leads to ability to choose the GWU distribution as the best model (Fig. 3). Conclusion The GWU distribution was introduced with the new three parameters. The new model can have bathtub and increasing shapes for the FRF. The new model displays the lifetime data on a bathtub shape as a versatility of FRF. Many structural properties of the GWU distribution have been studied and derived, encompassing the moment, moment-generating function and order statistics. The MLE was used to estimate the asymptotic interval and point estimates of the parameters. One set of real-life data was employed to demonstrate the flexibility of GWU distribution. The encourage results of new model were compared with the results of the others distributions in the literature. Discover revealed that the GWU distribution gives us a better fit than the other lifetime models. It is expected that the new model can be employed in other several sciences. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat​iveco​ mmons​.org/licen​ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. References 1. Alzaatreh, A., Famoye, F., Lee, C.: Weibull-pareto distribution and its applications. Commun. Stat. Theory Methods 42, 1673– 1691 (2013) 2. Bebbington, M., Lai, C.D., Zitikis, R.: A flexible Weibull extension. Reliab. Eng. Syst. Saf. 92, 719–726 (2007) 3. Bourguignon, M., Silva, R.B., Cordeiro, G.M.: The Weibull-G family of probability distributions. J. Data Sci. 12, 53–68 (2014) 13 Mathematical Sciences (2019) 13:105–114 4. Cooray, K.: Generalization of the Weibull distribution: the odd Weibull family. Stat. Model. 6, 265–277 (2006) 5. Cordeiro, G.M., Ortega, E.M.M., Nadarajah, S.: The Kumaraswamy Weibull distribution with application to failure data. J. Frankl. Inst. 347, 1399–1429 (2010) 6. Cordeiro, G.M., Ortega, E.M.M., Ramires, T.G.: A new generalized Weibull family of distributions: mathematical properties and applications. J. Stat. Distrib. Appl. 2, 13 (2015) 7. Famoye, F., Lee, C., Olumolade, O.: The beta-Weibull distribution. J. Stat. Theory Appl. 4, 121–136 (2005) 8. Kies, J.A.: The strength of glass. Technical Report 5093 Naval Research Laboratory Washington, DC (1958) 9. Lai, C.D., Xie, M., Murthy, D.N.P.: A modified Weibull distribution. IEEE Trans. Reliab. 52, 33–37 (2003) 10. Meeker, W.Q., Escobar, L.A.: Statistical Methods for Reliability Data. Wiley, New York (1998) 11. Mudholkar, G.S., Srivastava, D.K.: Exponentiated Weibull family for analyzing bathtub failure-rate data. IEEE Trans. Reliab. 42, 299–302 (1993) 12. Nadarajah, S., Cordeiro, G.M., Ortega, E.M.M.: General results for the beta-modified Weibull distribution. J. Stat. Comput. Simul. 81, 1211–1232 (2011) 13. Phani, K.: A new modified Weibull distribution function. J. Am. Ceram. Soc. 70, 182–184 (1987) 14. Sarhan, A.M., Zaindin, M.: Modified Weibull distribution. Appl. Sci. 11, 123–136 (2009) 15. Silva, G.O., Ortega, E.M.M., Cordeiro, G.M.: The beta modified Weibull distribution. Lifetime Data Anal. 16, 409–430 (2010) 16. Tahir, M., Alizadeh, M., Mansoor, M., Cordeiro, G.M., Zubair, M.: The Weibull-power function distribution with applications. Hacet. J. Math. Stat. 45, 245–265 (2016a) 17. Tahir, M.H., Cordeiro, G.M., Alzaatreh, A., Mansoor, M., Zubair, M.: A new Weibull-Pareto distribution: properties and applications. Commun. Stat Simul. Comput. 45, 3548–3567 (2016b) 18. Xie, M., Lai, C.D.: Reliability analysis using an additive Weibull model with bathtub-shaped failure rate function. Reliab. Eng. Syst. Saf. 52, 87–93 (1996) 19. Xie, M., Tang, Y., Goh, T.N.: A modified Weibull extension with bathtub-shaped failure rate function. Reliab. Eng. Syst. Saf. 76, 279–285 (2002) Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Likelihood theory, numerical methods, simulation methods April 8, 2015 1 Numerical maximization of the likelihood For many statistical models, closed form expressions for the maximum likelihood estimators can not be derived. For example, suppose that t1 , t2 , . . . , tn is a sample from a Weibull distribution with parameters a and b. The density function for the Weibull distribution is given by a f (t) = b  a−1 t a t e−( b ) . b (1) The Weibull distribution can be derived from the assumption that the rate of mortality λ(t) = (1/b)(t/b)a−1 . When the parameter a = 1, the rate of mortality is constant and the models becomes equivalent to the exponential model. For other values of a, the rate of mortality either increases or deceases with age. Note that the parameterization used here (and by R) is somewhat differerent from the parameterization used in notat 2 in ST0103. Given this model, the likelihood function for the observed data is thus    a n Y a ti a−1 − tbi L(a, b) = e b b i=1 Y  ti a−1 P a n −an e− (ti /b) , =a b b (2) (3) and the log likelihood l(a, b) = n ln a − na ln b − (a − 1) X ln ti − X (ti /b)a . (4) Noting that ∂ ∂ a ln(ti /b) (ti /b)a = e = ea ln(ti /b) ln(ti /b) = (ti /b)a ln(ti /b), (5) ∂a ∂a and setting the partial derivatives of the log likelihood function l(a, b) equal to zero yields X X n − n ln b − ln ti − (ti /b)a ln(ti /b) = 0 (6) a and na a X + (ti /b)a . (7) b b These equations are non-linear in the unknowns a and b and no closed form solution can be found. The maximum likelihood estimates must therefore instead be computed using numerical methods. This can by first defining a function in R which computes the likelihood (or the log likelihood) for given values of the unknown parameters a and b and for a given set of data. We shall then use a special function in R which finds the parameter values which maximises the likelihood. This function (the R function optim) works by making repeated calls to a given − 1 function for different parameter values. Based on the evaluated value of likelihood at these parameter values, the optimisation algorithm usually finds it’s way at least to a local maximum in the parameter space. Let us start by defining a function in R which computes the likelihood. This should be a function of the unknown parameters a and b and since the likelihood also depends on the observed sample t1 , t2 , . . . , tn we need an additional argument containing the observations. Our R function may then compute either the likelihood of the log likelihood. In practice, optimisation on the log likelihood works best. By default, optim minimises the given function. Thus if we let our function compute the negative log likelihood and minimise this we maximise the likelihood. Based on equation (4) we can define the likelihood function in R as follows l