# A Bimodal Weibull Distribution: Properties and Inference

Roberto Vilaa1 Mehmet Niyazi Çankayaa2,a3 a1 Departamento de Estatística,
Universidade de Brasília, Brazil  a2 Faculty of Applied Sciences, Department
of International Trading and Finance, Uşak University,  a3 Faculty of Art and
Sciences, Department of Statistics, Uşak University, Uşak, Turkey

###### Abstract

Abstract

Modeling is a challenging topic and using parametric models is an important
stage to reach flexible function for modeling. Weibull distribution has two
parameters which are shape $\alpha$ and scale $\beta$. In this study,
bimodality parameter is added and so bimodal Weibull distribution is proposed
by using a quadratic transformation technique used to generate bimodal
functions produced due to using the quadratic expression. The analytical
simplicity of Weibull and quadratic form give an advantage to derive a bimodal
Weibull via constructing normalizing constant. The characteristics and
properties of the proposed distribution are examined to show its usability in
modeling. After examination as first stage in modeling issue, it is
appropriate to use bimodal Weibull for modeling data sets. Two estimation
methods which are maximum $\log_{q}$ likelihood and its special form including
objective functions $\log_{q}(f)$ and $\log(f)$ are used to estimate the
parameters of shape, scale and bimodality parameters of the function. The
second stage in modeling is overcome by using heuristic algorithm for
optimization of function according to parameters due to fact that converging
to global point of objective function is performed by heuristic algorithm
based on the stochastic optimization. Real data sets are provided to show the
modeling competence of the proposed distribution.

Keywords. Bimodal distribution; Estimation; Modelling; $q$-inference; Weibull
distribution

##  I Introduction

Weibull distribution is very popular and used extensively in the applied field
of science. There are two main parameters which are shape $\alpha$ and scale
$\beta$ in probability density (p.d.) function of
Weibull$(\boldsymbol{\theta})$, where $\boldsymbol{\theta}=(\alpha,\beta)$.
The special values of shape and scale parameters drop to exponential and
Rayleigh distributions Weibull . Weibull distribution is used extensively and
different forms of Weibull were generated by EPdist ; Caste  and references
therein. It has many applications in survival analysis, bioscience Martinez ;
OrtegaBetaWei , management of maintenance, reliability engineering reliability
. The exponential kernels in Weibull and Gamma distributions are
$\exp(-z^{\alpha})$ and $\exp(-z)$, respectively, which shows that Weibull is
light tailed function if $\alpha>1$. Thus, Weibull can have advantegous when
compared with Gamma for the light tailed empirical distributions. In order to
have the different forms of tails at the positive axis on the real line,
$q$-deformation can be used as an alternative solution to produce heavy-tailed
functions Bercher12a . It should be noted that $\log_{q}(f)={f^{1-q}-1\over
1-q}$, $f\in[0,1]$, $q\in\mathbb{R}\backslash\\{1\\}$ as an objective function
produces heavy-tailed function. It is not easy to know which parametric model
will be true one to represent accurately a population in the questionnaire,
which makes a challenge task for us to choose the true parametric model for
the population. This is an outstanding topic in where it is difficult to
finalize the dispute about choosing the true model. If p.d. function is
defined on the real line, then location can be used to derive the bimodality
around location Cankaya2018  and references therein. In order to derive a
bimodal function defined on the positive axis of the real line, two p.d.
functions can be tried to mix. In this case, the theoretical property of such
mixed function should be examined and must be provided for existence of
moments, entropies, tail behaviour, etc in order to apply it into modelling
issue. Further, the computational complexity can arise due to the increased
number of parameters which will be responsible to model subgroups assumed to
be a member of subpopulation at a phenomena in questionnaire yangcomp . The
managing bimodaliy for distribution defined on the positive axis is an another
challenging task for us if data set shows a bimodality. Another task is about
the applying a p.d. function which has important properties such as existence
of moments, entropies, etc for conducting an accurate modelling on the data
sets. The quadratic transformation for normal distribution is proposed by
bimodal  to produce bimodality and also asymmetry. Using the parameters whose
role is defined exactly can guarantee the smoothness property (absolutely
continuous function which is differantiable according to Radon Nikodym
derivative knapp ) of the derived function as well. In order to derive a
bimodal distribution from a unimodal distribution, we can use quadratic
transformation technique which gives a polynomial movement owing the fact that
the used quadratic function can make a fluctation on the unimodal function. We
will propose a new bimodal distribution generated by quadratic transformation
technique with bimodality parameter $\delta$ and unimodal Weibull distribution
with parameters shape $\alpha$ and scale $\beta$. Thus, we can observe how the
bimodality can be derived. After applying this technique, we can get a
normalizing constant to produce a p.d. function. The property of a newly
generated function should be extensively examined in order to pass the test
about finiteness of moments, existence of entropies, etc. Thus, the proposed
distribution can be used to model a data set.

It can be generally observed that a phenomena in the universe can show a
bimodality. There can be a contamination or an irregularity into underlying or
majority of the data. If the replicated values at an interval on the real line
are observed, then bimodality occurs. In other words, a data set can be a
combination of different parametric models or different values of same
parametric model $f(x;\boldsymbol{\theta})$ to represent bimodality mix1 ;
mix2 ; yangcomp . The light tailed property, tractability of analytical
expression, bimodality kernel $1+(1-\delta x)^{2}$ used to derive only one
mode property occurred due to the polynomial degree $2$ bimodal , and also for
confining yourself for modelling data set having two modes at the different
degree high of peakedness can require to apply bimodality generator for
Weibull distribution.

The organization of this paper is as following forms. Bimodal Weibull
distribution is proposed by Section II. Several properties of the proposed
parametric model are examined and introduced by Sections III and IV. The
section V introduces that the model parameters are estimated by maximum $\log$
and $\log_{q}$ likelihood methods. Section VI gives numerial assessments of
the proposed model and also provides the comparison between recently proposed
Bimodal Gamma (BGamma) distribution Vila2020  and the model in this paper. The
conclusion is given by Conclusions.

##  II Bimodal Weibull distribution

We say that a non-negative random variable $X$ has a bimodal Weibull
distribution with vector parameter
$\boldsymbol{\theta}=(\alpha,\beta,\delta)$, denoted by
$X\sim\text{BWeibull}(\boldsymbol{\theta})$, if its p.d. function is given by

$\displaystyle f(x;\boldsymbol{\theta})={\alpha\over\beta
Z_{\boldsymbol{\theta}}}\,\big{[}1+(1-\delta
x)^{2}\big{]}\,\bigg{(}{x\over\beta}\bigg{)}^{\alpha-1}\exp\Bigg{[}-\bigg{(}{x\over\beta}\bigg{)}^{\alpha}\Bigg{]},\quad
x\geqslant 0;\ \alpha,\beta>0,\ \delta\in\mathbb{R},$ (1)

where $Z_{\boldsymbol{\theta}}$ is the normalization constant given by

$\displaystyle
Z_{\boldsymbol{\theta}}=2+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)},$
(2)

which depends only on the parameter vector $\boldsymbol{\theta}$, $\alpha$ is
the shape parameter, $\beta$ is the scale parameter and $\delta$ is the
parameter that controls the uni- or bimodality of the distribution. Here,
$\Gamma(z)$ denotes the complete gamma function.

The behavior of $f(x;\boldsymbol{\theta})$ with $x\to 0$ or $x\to\infty$ is as
follows:

$\displaystyle\lim_{x\to
0}f(x;\boldsymbol{\theta})=\begin{cases}\infty&\text{for}\quad 0<\alpha<1,\\\
\dfrac{2}{\beta Z_{\boldsymbol{\theta}}}&\text{for}\quad\alpha=1,\\\
0&\text{for}\quad\alpha>1,\end{cases}$
$\displaystyle\lim_{x\to\infty}f(x;\boldsymbol{\theta})=0\quad\forall\alpha>0.$

It is verified that the cumulated distribution (c.d.) function of
$X\sim\text{BWeibull}(\boldsymbol{\theta})$ is given by (see Item 2 of
Proposition IV.7 for more details)

$\displaystyle F(x;\boldsymbol{\theta})=\dfrac{\displaystyle
2-\big{[}1+(1-\delta
x)^{2}\big{]}\,\exp\Bigg{[}-\bigg{(}{x\over\beta}\bigg{)}^{\alpha}\Bigg{]}-{2\delta\beta\over\alpha}\,\bigg{[}\gamma\bigg{(}{1\over\alpha},{x^{\alpha}\over\beta^{\alpha}}\bigg{)}-{\delta\beta}\,\gamma\bigg{(}{2\over\alpha},{x^{\alpha}\over\beta^{\alpha}}\bigg{)}\bigg{]}}{\displaystyle
2+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)}},\quad
x\geqslant 0,$

where $\Gamma(s,x)$, $s>0$, is the upper incomplete gamma function, and
$\gamma(s,x)$, $s>0$, is the lower incomplete gamma function. Related
quantities such as reliability, hazard rate and the mean residual life can be
found in Subsection IV.3.

The p.d and c.d. functions (PDF and CDF) in Figures 1 and 2 are drawn for
different values of parameters.

(a) PDF of BWeibull for $\delta$

(b) PDF of BWeibull for $\alpha$ and $\delta$

Figure 1: PDF of BWeibull with different values of parameters $\alpha,\delta$
and $\beta=2$.

(a) CDF of BWeibull for $\delta$

(b) CDF of BWeibull for $\alpha$ and $\delta$

Figure 2: CDF of BWeibull with different values of parameters $\alpha,\delta$
and $\beta=2$.

##  III Unimodality and Bimodality properties

######  Proposition III.1.

A point $x$ is a mode of the BWeibull density (1), if and only if it is a
solution of the following equation

$\displaystyle\alpha\delta^{2}x^{\alpha+2}-2\alpha\delta x^{\alpha+1}+2\alpha
x^{\alpha}-(\alpha+1)\beta^{\alpha}\delta^{2}x^{2}+2\alpha\beta^{\alpha}\delta
x-2(\alpha-1)\beta^{\alpha}=0.$ (3)

###### Proof.

The proof is trivial and omitted. ∎

######  Proposition III.2.

If $\alpha>2$ is a natural number, the following it hold:

  * 1)

If $\delta<0$ then the equation (3) has an unique positive root;

  * 2)

If $\delta>0$ then the equation (3) has five, or three, or one positive roots.

###### Proof.

The proof follows directly by applying Descartes’ rule of signs (see, e.g. Xue
2012 Xue2012 ). ∎

######  Proposition III.3 (Unimodality and monotonicity).

Let $X\sim\text{BWeibull}(\boldsymbol{\theta})$ with $\alpha=1$. The following
hold:

  * 1)

If $\delta<-1/\beta$ then the BWeibull density is unimodal;

  * 2)

If $\delta>0$ then the BWeibull density is unimodal or strictly decreasing.

###### Proof.

Let $\alpha=1$. In this case note that the equation (3) can be written as

$\displaystyle
p_{2}(x)=\delta^{2}x^{2}-2\delta(1+\beta\delta)x+2(1+\beta\delta)=0.$

Firstly we assume that $\delta<-1/\beta$. By Descartes’ rule of signs the
above equation has a unique positive root, denoted by $x_{0}$. Since
$\lim_{x\to 0}f(x;\boldsymbol{\theta})=2/(\beta Z_{\boldsymbol{\theta}})$ and
$\lim_{x\to\infty}f(x;\boldsymbol{\theta})=0$ we have that $x_{0}$ is the
unique maximum point with maximum value $f(x_{0};\boldsymbol{\theta})>2/(\beta
Z_{\boldsymbol{\theta}})$. This proves the statement in the first item.

In the second case, $\delta>0$, Descartes’ rule of signs guarantees that the
equation $p_{2}(x)=0$ has two roots, or zero roots. Since $\lim_{x\to
0}f(x;\boldsymbol{\theta})=2/(\beta Z_{\boldsymbol{\theta}})$ and
$\lim_{x\to\infty}f(x;\boldsymbol{\theta})=0$, the proof of second item
follows. ∎

Let us define

$\displaystyle\Delta$
$\displaystyle=256a^{3}e^{3}-192a^{2}bde^{2}-128a^{2}c^{2}e^{2}+144a^{2}cd^{2}e-27a^{2}d^{4}$
$\displaystyle+144ab^{2}ce^{2}-6ab^{2}d^{2}e-80abc^{2}de+18abcd^{3}+16ac^{4}e$
(4)
$\displaystyle-4ac^{3}d^{2}-27b^{4}e^{2}+18b^{3}cde-4b^{3}d^{3}-4b^{2}c^{3}e+b^{2}c^{2}d^{2}$

the discriminant of the quartic polynomial $ax^{4}+bx^{3}+cx^{2}+dx+e$, where

$\displaystyle a=2\delta^{2},\quad b=-4\delta,\quad
c=4-3\beta^{2}\delta^{2},\quad d=4\beta^{2}\delta,\quad e=-2\beta^{2}.$ (5)

######  Theorem III.4 (Bimodality and unimodality).

Let $X\sim\text{BWeibull}(\boldsymbol{\theta})$ with $\alpha=2$ and
$\delta>0$. The following hold:

  * 1)

If $\Delta>0$ and $13/12<\beta^{2}\delta^{2}<4/3$, then the BWeibull density
is bimodal;

  * 2)

If $\Delta<0$ and $\beta^{2}\delta^{2}\leqslant 4/{3}$, then the BWeibull
density is unimodal;

  * 3)

If $\Delta=0$ and $\beta^{2}\delta^{2}=4/{3}$, then the BWeibull density is
unimodal;

where $\Delta$ is defined by (III) and (5).

###### Proof.

Let $x$ be a mode of the BWeibull density. By Proposition III.1 this mode
satisfies the equation (3), with $\alpha=2$, or equivalently this one is
solution of the quartic equation

$\displaystyle p_{4}(x)=ax^{4}+bx^{3}+cx^{2}+dx+e=0$

with real coefficients defined in (5). By using that $\delta>0$, in both
cases, $\beta^{2}\delta^{2}<4/{3}$ or $\beta^{2}\delta^{2}=4/{3}$ (in this
case, $c=0$), Descartes’ rule of signs guarantees that:

$\displaystyle\text{The polynomial}\ p_{4}(x)\ \text{has three, or one
positive roots and one negative root.}$ (6)

To check Items 1), 2) and 3), we first define the following quantities:

$\displaystyle P=8ac-3b^{2}=4\delta^{2}(13-12\beta^{2}\delta^{2});$
$\displaystyle
D=64a^{3}e-16a^{2}c^{2}+16ab^{2}c-16a^{2}bd-3b^{4}=64(4-9\beta^{4}\delta^{4});$
$\displaystyle\Delta_{0}=c^{2}-3bd+12ae=16+9\beta^{4}\delta^{4}-24\beta^{2}\delta^{2}.$

1) If $\Delta>0$ then either the equation’s four roots $p_{4}(x)=0$ are all
real or none is. Since $P<0$ and $D<0$ whenever $\beta^{2}\delta^{2}>13/12$
and $\beta^{2}\delta^{2}>2/3$, respectively, then all four roots of $p_{4}(x)$
are real and distinct. Then, from affirmation (6) it follows that the equation
$p_{4}(x)=0$ has three distinct positive roots, denoted by
$x_{1},x_{2},x_{3}$, and one negative root. Let’s assume that
$x_{1}<x_{2}<x_{3}$. Since $\lim_{x\to
0}f(x;\boldsymbol{\theta})=\lim_{x\to\infty}f(x;\boldsymbol{\theta})=0$ we
have that $x_{1}$ and $x_{3}$ are two maximum points and $x_{1}$ is the unique
minimum point. This proves the bimodality property in the first item.

2) If $\Delta<0$ then the equation $p_{4}(x)=0$ has two distinct real roots
and two complex conjugate non-real roots. Then, from (6) it follows that the
equation $p_{4}(x)=0$ has one positive root, denoted by $x_{0}$, and one
negative root. Since $\lim_{x\to
0}f(x;\boldsymbol{\theta})=\lim_{x\to\infty}f(x;\boldsymbol{\theta})=0$, it
follows that $x_{0}$ is the unique maximum point. This proves the second item.

3) If $\Delta=0$ then the polynomial $p_{4}(x)$ has a multiple root. Since
$\Delta_{0}=0$ and $D=-768\neq 0$ when $\beta^{2}\delta^{2}=4/{3}$, there are
a triple root and a simple root, all real. Then, by affirmation (6), the
equation $p_{4}(x)=0$ has a positive triple root, denoted by $x_{0}$, and one
negative root. Since $\lim_{x\to
0}f(x;\boldsymbol{\theta})=\lim_{x\to\infty}f(x;\boldsymbol{\theta})=0$ we
have that $x_{0}$ is the unique maximum point. This completes the proof of
third item. ∎

##  IV Characteristics of Probability Denstiy Function

###  IV.1 Real moments

######  Proposition IV.1.

If $X\sim\text{BWeibull}(\boldsymbol{\theta})$ then

$\displaystyle\mathbb{E}(X^{r})=\beta^{r}\,\dfrac{\displaystyle
2\Gamma\bigg{(}1+{r\over\alpha}\bigg{)}+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{r+2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{r+1\over\alpha}\bigg{)}}{\displaystyle
2+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)}},\quad
r>-\alpha,$

where $\Gamma(z)$ is the complete gamma function.

###### Proof.

A simple observation shows that

$\displaystyle\mathbb{E}(X^{r})={1\over
Z_{\boldsymbol{\theta}}}\,\left[2\mathbb{E}(Y^{r})+\delta^{2}\mathbb{E}(Y^{r+2})-2\delta\mathbb{E}(Y^{r+1})\right],$

where $Y\sim\text{Weibull}(\alpha,\beta)$ and $Z_{\boldsymbol{\theta}}$ is as
in (2). By using the following formula of Item (6) in Cankaya2018 :

$\displaystyle\Gamma\bigg{(}s+{1\over\alpha}\bigg{)}=\alpha p^{\alpha
s+1}\int_{0}^{\infty}y^{\alpha s}\,\exp\big{[}-(yp)^{\alpha}\big{]}\,{\rm
d}y,$

we have

$\displaystyle\mathbb{E}(Y^{r})=\beta^{r}\Gamma\bigg{(}1+{r\over\alpha}\bigg{)},\quad
r>-\alpha.$

By combining the above identities, the proof follows. ∎

As a consequence of the above proposition, the closed expressions for the
moments, variance, skewness and kurtosis of random variable $X$ are easily
obtained. The Figures 3(a) and 3(b) give the skewness and kurtosis for the
different values of bimodality parameter $\delta$, respectively. For the tried
values of $\alpha$ and $\beta$, we observe that when the value of $\delta$ are
not large, we get a large scale of skewness and kurtosis. The Figure 3
represents that the parameter $\delta$ increases the flexiblity of function.

(a) Skewness of BWeibull for $\alpha$

(b) Kurtosis of BWeibull for $\alpha$

Figure 3: Skewness and kurtosis of BWeibull with different values of $\alpha$
and $\beta=2$.

###  IV.2 Moment generating function

######  Theorem IV.2.

If $X\sim\text{BWeibull}(\boldsymbol{\theta})$ then the moment generating
function $M_{X}(t)=\mathbb{E}[\exp({tX})]$ can be expressed as

$\displaystyle M_{X}(t)=\begin{cases}\displaystyle\sum_{n=0}^{\infty}{(\beta
t)^{n}}\,\dfrac{1+\dfrac{1}{2}\,{\delta^{2}\beta^{2}}(n+2)(n+1)-\delta\beta(n+1)}{1+\delta^{2}\beta^{2}-\delta\beta}&\text{for}\
\alpha=1\ \text{and}\ |t|<1/\beta,\\\\[11.38092pt]
\displaystyle\sum_{n=0}^{\infty}{(\beta t)^{n}\over n!}\,\dfrac{\displaystyle
2\Gamma\bigg{(}1+{n\over\alpha}\bigg{)}+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{n+2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{n+1\over\alpha}\bigg{)}}{\displaystyle
2+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)}}&\text{for}\
\alpha>1\ \text{and}\ t\in\mathbb{R}.\end{cases}$

###### Proof.

Using the series expansion of $\exp({tX})$ we have

$\displaystyle M_{X}(t)=\mathbb{E}\Bigg{[}\sum_{n=0}^{\infty}{(tX)^{n}\over
n}\Bigg{]}.$

Assuming the validity of the following identity

$\displaystyle\mathbb{E}\Bigg{[}\sum_{n=0}^{\infty}{(tX)^{n}\over
n}\Bigg{]}=\sum_{n=0}^{\infty}{t^{n}\over
n!}\,\mathbb{E}(X^{n}),\quad\alpha\geqslant 1,$ (7)

note that by using Proposition IV.1 and the identity $\Gamma(n)=(n-1)!$, for
$n=1,2,3,\ldots$ (in the case $\alpha=1$); the proof of the proposition
follows.

In what remains of the proof we prove the identity (7). Indeed, from Monotone
Convergence Theorem knapp  we have

$\displaystyle\mathbb{E}\left[\sum_{n=0}^{\infty}\bigg{|}{(tX)^{n}\over
n}\bigg{|}\right]=\sum_{n=0}^{\infty}\mathbb{E}\Bigg{[}\bigg{|}{(tX)^{n}\over
n}\bigg{|}\Bigg{]}=\sum_{n=0}^{\infty}{|t|^{n}\over n!}\,\mathbb{E}(X^{n}),$
(8)

where, by Proposition IV.1, the series $\sum_{n=0}^{\infty}{|t|^{n}\over
n!}\,\mathbb{E}(X^{n})$ is written as

$\displaystyle\sum_{n=0}^{\infty}{|t|^{n}\over
n!}\,\mathbb{E}(X^{n})=\sum_{n=0}^{\infty}{(\beta|t|)^{n}\over
n!}\,\dfrac{\displaystyle
2\Gamma\bigg{(}1+{n\over\alpha}\bigg{)}+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{n+2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{n+1\over\alpha}\bigg{)}}{\displaystyle
2+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)}}.$
(9)

Note that the series $u(k)$ defined by

$\displaystyle
u(k)=\sum_{n=0}^{\infty}a_{n}=\sum_{n=0}^{\infty}{(\beta|t|)^{n}\over
n!}\,\Gamma\bigg{(}1+{n+k\over\alpha}\bigg{)},\quad\quad k=0,1,2,$

converges when $\alpha\geqslant 1$, because by the ratio test the limit

$\displaystyle L$ $\displaystyle=\lim_{n\to\infty}\left|{a_{n+1}\over
a_{n}}\right|$ $\displaystyle={\beta|t|}\lim_{n\to\infty}\bigg{(}{1\over
n+k}\bigg{)}\,\bigg{(}1+{k\over
n+1}\bigg{)}\,{\displaystyle\Gamma\bigg{(}{n+1+k\over\alpha}\bigg{)}\over\displaystyle\Gamma\bigg{(}{n+k\over\alpha}\bigg{)}},\quad
k=0,1,2,$ $\displaystyle=\begin{cases}0&\text{for}\quad\alpha>1,\\\
{\beta|t|}&\text{for}\quad\alpha=1,\\\ \infty&\text{for}\quad
0<\alpha<1,\end{cases}$

is less than 1 when $\alpha>1$ $\forall t\in\mathbb{R}$; and $\alpha=1$
$\forall|t|<1/\beta$.

Therefore, the series (9)

$\displaystyle\sum_{n=0}^{\infty}{|t|^{n}\over n!}\,\mathbb{E}(X^{n})$
$\displaystyle=\dfrac{\displaystyle 2\sum_{n=0}^{\infty}{(\beta|t|)^{n}\over
n!}\,\Gamma\bigg{(}1+{n\over\alpha}\bigg{)}+\delta^{2}\beta^{2}\sum_{n=0}^{\infty}{(\beta|t|)^{n}\over
n!}\,\Gamma\bigg{(}1+{n+2\over\alpha}\bigg{)}-2\delta\beta\sum_{n=0}^{\infty}{(\beta|t|)^{n}\over
n!}\,\Gamma\bigg{(}1+{n+1\over\alpha}\bigg{)}}{\displaystyle
2+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)}}$
$\displaystyle=\dfrac{\displaystyle 2u(0)+\delta^{2}\beta^{2}u(2)-2\delta\beta
u(1)}{\displaystyle
2+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)}}$

converges for $\alpha\geqslant 1$, and then

$\displaystyle\mathbb{E}\left[\sum_{n=0}^{\infty}\bigg{|}{(tX)^{n}\over
n}\bigg{|}\right]\stackrel{{\scriptstyle\eqref{first-
relation}}}{{=}}\sum_{n=0}^{\infty}{|t|^{n}\over
n!}\,\mathbb{E}(X^{n})<\infty,$

when $\alpha\geqslant 1$. So, applying Fubini’s Theorem knapp  we get

$\displaystyle M_{X}(t)=\mathbb{E}\Bigg{[}\sum_{n=0}^{\infty}{(tX)^{n}\over
n}\Bigg{]}=\sum_{n=0}^{\infty}{t^{n}\over n!}\,\mathbb{E}(X^{n}),$

whenever $\alpha\geqslant 1$. This proves (7). We thus complete the proof. ∎

######  Corollary IV.3 (Light-tailed distribution).

If $X\sim\text{BWeibull}(\boldsymbol{\theta})$ and $\alpha\geqslant 1$, then
there exists $t_{0}>0$ such that $\mathbb{P}(X>x)\leqslant\exp(-t_{0}x)$ for
$x$ large enough.

###### Proof.

In the first case $\alpha=1$, by Theorem IV.2, there exists $|t_{0}|<1/\beta$
such that $M_{X}(t_{0})<\infty$, $X\sim\text{BWeibull}(\boldsymbol{\theta})$.
In the second case $\alpha>1$, again, by Theorem IV.2, there exists
$t_{0}\in\mathbb{R}$ such that $M_{X}(t_{0})<\infty$. Then the proof follows.
∎

######  Remark IV.4.

Let $X$ a continuous random variable with density function $f_{X}(x)$.
Following the reference Klugman1998 , the rate of a random variable is given
by

$\displaystyle\tau_{X}=-\lim_{x\to\infty}{{\rm
d}\log\big{[}f_{X}(x)\big{]}\over{\rm d}x}.$

A simple computation shows that

$\displaystyle\tau_{\rm
BWeibull(\boldsymbol{\theta})}=\lim_{x\to\infty}\left[\dfrac{2\delta(1-\delta
x)}{1+(1-\delta x)^{2}}-(\alpha-1)\,{1\over
x}+{\alpha\over\beta^{\alpha}}\,x^{\alpha-1}\right]=\begin{cases}{1\over\beta}&{\rm
for}\quad\alpha=1,\\\ \infty&{\rm for}\quad\alpha>1,\\\ 0&{\rm
for}\quad\alpha<1.\end{cases}$

Then, far enough out in the tail, every BWeibull distribution looks like an
exponential distribution when $\alpha=1$ and a Normal distribution when
$\alpha>1$. In addition, we have some comparisons between the rates of random
variables with known distributions: Inverse-gamma, Log-normal, Generalized-
Pareto, BWeibull, BGamma Vila2020 , exponential and Normal;

$\displaystyle\tau_{{\rm InvGamma}(\alpha,\beta)}=\tau_{{\rm
LogNorm}(\mu,\sigma^{2})}$ $\displaystyle=\tau_{{\rm
GenPareto}(\alpha,\beta,\xi)}=\tau_{\rm BWeibull(\alpha<1,\beta,\delta)}=0$
$\displaystyle<\tau_{\rm BWeibull(\alpha=1,\beta,\delta)}=\tau_{{\rm
BGamma}(\alpha,1/\beta,\delta)}=\tau_{{\rm exp}(1/\beta)}=1/\beta$
$\displaystyle<\tau_{\rm BWeibull(\alpha>1,\beta,\delta)}=\tau_{{\rm
Normal}(\mu,\sigma^{2})}=\infty.$

By using the very well-known relation $M_{X}(t)=\phi_{X}(-it)$ between moment
generating function and characteristic function of
$X\sim\text{BWeibull}(\boldsymbol{\theta})$, the following result follows
immediately.

######  Proposition IV.5.

If $X\sim\text{BWeibull}(\boldsymbol{\theta})$ then the characteristic
function $\phi_{X}(t)=\mathbb{E}[\exp({itX})]$ can be expressed as

$\displaystyle\phi_{X}(t)=\sum_{n=0}^{\infty}{(i\beta t)^{n}\over
n!}\,\dfrac{\displaystyle
2\Gamma\bigg{(}1+{n\over\alpha}\bigg{)}+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{n+2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{n+1\over\alpha}\bigg{)}}{\displaystyle
2+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)}},\quad
t\in\mathbb{R}.$

###  IV.3 Reliability, hazard rate and the mean residual life

For any $t\geqslant 0$, the reliability, the hazard rate and the mean residual
life functions, associated with a random variable
$X\sim\text{BWeibull}(\boldsymbol{\theta})$, are defined as follows

$\displaystyle R(t;$
$\displaystyle\boldsymbol{\theta})=\int_{t}^{\infty}f(x;\boldsymbol{\theta})\,\textrm{d}x,\quad
H(t;\boldsymbol{\theta})={f(t;\boldsymbol{\theta})\over
R(t;\boldsymbol{\theta})},\quad\textrm{MRL}(t;\boldsymbol{\theta})={1\over
R(t;\boldsymbol{\theta})}\int_{t}^{\infty}R(x;\boldsymbol{\theta})\,\textrm{d}x,$

respectively.

Let $Y\sim\text{Weibull}(\alpha,\beta)$ be a random variable with Weibull
distribution. By using the integration by parts formula we have

$\displaystyle\mathbb{E}\big{[}{1}_{\\{Y\geqslant
t\\}}Y^{s}\big{]}=t^{s}\exp\Bigg{[}-\bigg{(}{t\over\beta}\bigg{)}^{\alpha}\Bigg{]}+s\int_{t}^{\infty}y^{s-1}\exp\Bigg{[}-\bigg{(}{y\over\beta}\bigg{)}^{\alpha}\Bigg{]}\,\textrm{d}y,\quad
s\geqslant 0.$

By taking the change of variable $x=\big{(}{u\over\beta}\big{)}^{\alpha}$ the
right-hand expression above is

$\displaystyle=t^{s}\exp\Bigg{[}-\bigg{(}{t\over\beta}\bigg{)}^{\alpha}\Bigg{]}+{s\beta^{s}\over\alpha}\int_{(t/\beta)^{\alpha}}^{\infty}x^{{s\over\alpha}-1}\exp(-x)\,\textrm{d}x$
$\displaystyle=t^{s}\exp\Bigg{[}-\bigg{(}{t\over\beta}\bigg{)}^{\alpha}\Bigg{]}+{s\beta^{s}\over\alpha}\,\Gamma\bigg{(}{s\over\alpha},{t^{\alpha}\over\beta^{\alpha}}\bigg{)},$

where $\Gamma(s,x)$, $s>0$, is the upper incomplete gamma function, and where
we are adopting the notation $\Gamma(0,x)=0$. Therefore,

$\displaystyle\mathbb{E}\big{[}{1}_{\\{Y\geqslant
t\\}}Y^{s}\big{]}=t^{s}\exp\Bigg{[}-\bigg{(}{t\over\beta}\bigg{)}^{\alpha}\Bigg{]}+{s\beta^{s}\over\alpha}\,\Gamma\bigg{(}{s\over\alpha},{t^{\alpha}\over\beta^{\alpha}}\bigg{)},\quad
s\geqslant 0.$ (10)

Similarly we find that

$\displaystyle\mathbb{E}\big{[}{1}_{\\{Y\leqslant
t\\}}Y^{s}\big{]}=\begin{cases}\displaystyle
1-\exp\Bigg{[}-\bigg{(}{t\over\beta}\bigg{)}^{\alpha}\Bigg{]}&\text{for}\quad
s=0,\\\\[14.22636pt]
\displaystyle-t^{s}\exp\Bigg{[}-\bigg{(}{t\over\beta}\bigg{)}^{\alpha}\Bigg{]}+{s\beta^{s}\over\alpha}\,\gamma\bigg{(}{s\over\alpha},{t^{\alpha}\over\beta^{\alpha}}\bigg{)}&\text{for}\quad
s>0,\end{cases}$ (11)

where $\gamma(s,x)$, $s>0$, is the lower incomplete gamma function.

######  Remark IV.6.

By combining (10) and (11) we have

$\displaystyle\mathbb{E}(X^{s})$
$\displaystyle=\mathbb{E}\big{[}{1}_{\\{Y\leqslant
t\\}}Y^{s}\big{]}+\mathbb{E}\big{[}{1}_{\\{Y\geqslant t\\}}Y^{s}\big{]},\quad
s>0,$
$\displaystyle={s\beta^{s}\over\alpha}\,\bigg{[}\gamma\bigg{(}{s\over\alpha},{t^{\alpha}\over\beta^{\alpha}}\bigg{)}+\Gamma\bigg{(}{s\over\alpha},{t^{\alpha}\over\beta^{\alpha}}\bigg{)}\bigg{]}={s\beta^{s}\over\alpha}\,\Gamma\bigg{(}{s\over\alpha}\bigg{)}={\beta^{s}}\,\Gamma\bigg{(}1+{s\over\alpha}\bigg{)},$

where in the third and fourth equalities we use the identities
$\gamma(s,x)+\Gamma(s,x)=\Gamma(s)$ and $\Gamma(1+s)=s\Gamma(s)$,
respectively. This confirms the validity of Proposition IV.1.

######  Proposition IV.7.

If $X\sim\text{BWeibull}(\boldsymbol{\theta})$ then

  1. 1)

$\displaystyle R(t;\boldsymbol{\theta})=\dfrac{\displaystyle\big{[}1+(1-\delta
t)^{2}\big{]}\,\exp\Bigg{[}-\bigg{(}{t\over\beta}\bigg{)}^{\alpha}\Bigg{]}-{2\delta\beta\over\alpha}\,\bigg{[}\Gamma\bigg{(}{1\over\alpha},{t^{\alpha}\over\beta^{\alpha}}\bigg{)}-{\delta\beta}\,\Gamma\bigg{(}{2\over\alpha},{t^{\alpha}\over\beta^{\alpha}}\bigg{)}\bigg{]}}{\displaystyle
2+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)}};$

  2. 2)

$F(t;\boldsymbol{\theta})=\dfrac{\displaystyle 2-\big{[}1+(1-\delta
t)^{2}\big{]}\,\exp\Bigg{[}-\bigg{(}{t\over\beta}\bigg{)}^{\alpha}\Bigg{]}-{2\delta\beta\over\alpha}\,\bigg{[}\gamma\bigg{(}{1\over\alpha},{t^{\alpha}\over\beta^{\alpha}}\bigg{)}-{\delta\beta}\,\gamma\bigg{(}{2\over\alpha},{t^{\alpha}\over\beta^{\alpha}}\bigg{)}\bigg{]}}{\displaystyle
2+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)}};$

  3. 3)

$H(t;\boldsymbol{\theta})=\dfrac{\displaystyle{\alpha\over\beta}\,\big{[}1+(1-\delta
t)^{2}\big{]}\,\bigg{(}{t\over\beta}\bigg{)}^{\alpha-1}\exp\Bigg{[}-\bigg{(}{t\over\beta}\bigg{)}^{\alpha}\Bigg{]}}{\displaystyle\big{[}1+(1-\delta
t)^{2}\big{]}\,\exp\Bigg{[}-\bigg{(}{t\over\beta}\bigg{)}^{\alpha}\Bigg{]}-{2\delta\beta\over\alpha}\,\bigg{[}\Gamma\bigg{(}{1\over\alpha},{t^{\alpha}\over\beta^{\alpha}}\bigg{)}-{\delta\beta}\,\Gamma\bigg{(}{2\over\alpha},{t^{\alpha}\over\beta^{\alpha}}\bigg{)}\bigg{]}}.$

  4. 4)

$\displaystyle\lim_{t\to 0}H(t;\boldsymbol{\theta})=\begin{cases}0&{\rm for}\
\alpha>1,\\\\[2.84544pt]
\displaystyle{1\over\beta(1-\delta\beta+\delta^{2}\beta^{2})}&{\rm for}\
\alpha=1,\\\\[5.69046pt] +\infty&{\rm for}\
\alpha<1,\end{cases}\quad\lim_{t\to+\infty}H(t;\boldsymbol{\theta})=\begin{cases}+\infty&{\rm
for}\ \alpha>1,\\\\[2.84544pt] \displaystyle{1\over\beta}&{\rm for}\
\alpha=1,\\\\[5.69046pt] 0&{\rm for}\ \alpha<1.\end{cases}$

The hazard function $H$ in Figure 4 is drawn for different values of
parameters.

(a) Hazard of BWeibull for $\delta$

(b) Hazard of BWeibull for $\alpha$ and $\delta$

Figure 4: Hazard of BWeibull with different values of parameters
$\alpha,\delta$ and $\beta=2$.

###### Proof.

A simple algebraic manipulation in the definition of BWeibull density in
equation (1) shows that

$\displaystyle\begin{array}[]{lllll}&\displaystyle
R(t;\boldsymbol{\theta})=\dfrac{1}{Z_{\boldsymbol{\theta}}}\,\left\\{2\,\mathbb{E}\big{[}{1}_{\\{Y\geqslant
t\\}}\big{]}-2\delta\,\mathbb{E}\big{[}{1}_{\\{Y\geqslant
t\\}}Y\big{]}+\delta^{2}\,\mathbb{E}\big{[}{1}_{\\{Y\geqslant
t\\}}Y^{2}\big{]}\right\\},\\\\[11.38092pt] &\displaystyle
F(t;\boldsymbol{\theta})=\dfrac{1}{Z_{\boldsymbol{\theta}}}\,\left\\{2\,\mathbb{E}\big{[}{1}_{\\{Y\leqslant
t\\}}\big{]}-2\delta\,\mathbb{E}\big{[}{1}_{\\{Y\leqslant
t\\}}Y\big{]}+\delta^{2}\,\mathbb{E}\big{[}{1}_{\\{Y\leqslant
t\\}}Y^{2}\big{]}\right\\},\end{array}\quad
Y\sim\text{Weibull}(\alpha,\beta),$

where $Z_{\boldsymbol{\theta}}$ is as in (2). By taking $s=0,1,2$ in (10), we
get the formula of Item 1). By taking $s=0,1,2$ in (11), the formula of Item
2) follows. The proof of Item 3) follows immediately by combining the
definition of $H(t;\boldsymbol{\theta})$ with Item 1). The proof of Item 4)
follows directly by analyzing the form of the hazard function
$H(t;\boldsymbol{\theta})$ in Item 3) through the known limits
$\lim_{x\to+\infty}\Gamma(s,x)=0$, $\lim_{x\to 0}\Gamma(s,x)=\Gamma(s)$, and
standard limit calculations. ∎

######  Proposition IV.8.

If $X\sim\text{BWeibull}(\boldsymbol{\theta})$ then

$\displaystyle\mathbb{E}\big{[}{1}_{\\{X\geqslant
t\\}}X\big{]}=\dfrac{\displaystyle t\big{[}1+(1-\delta
t)^{2}\big{]}\,\exp\Bigg{[}-\bigg{(}{t\over\beta}\bigg{)}^{\alpha}\Bigg{]}+{\beta\over\alpha}\,\bigg{[}2\Gamma\bigg{(}{1\over\alpha},{t^{\alpha}\over\beta^{\alpha}}\bigg{)}-4{\delta\beta}\,\Gamma\bigg{(}{2\over\alpha},{t^{\alpha}\over\beta^{\alpha}}\bigg{)}+3{\delta^{2}\beta^{2}}\,\Gamma\bigg{(}{3\over\alpha},{t^{\alpha}\over\beta^{\alpha}}\bigg{)}\bigg{]}}{\displaystyle
2+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)}}.$

###### Proof.

By a similar decomposition to that of Proposition IV.7, we have

$\displaystyle\mathbb{E}\big{[}{1}_{\\{X\geqslant
t\\}}X\big{]}=\dfrac{1}{Z_{\boldsymbol{\theta}}}\,\left\\{2\,\mathbb{E}\big{[}{1}_{\\{Y\geqslant
t\\}}Y\big{]}-2\delta\,\mathbb{E}\big{[}{1}_{\\{Y\geqslant
t\\}}Y^{2}\big{]}+\delta^{2}\,\mathbb{E}\big{[}{1}_{\\{Y\geqslant
t\\}}Y^{3}\big{]}\right\\},\quad Y\sim\text{Weibull}(\alpha,\beta),$

where $Z_{\boldsymbol{\theta}}$ is as in (2). Hence, by taking $s=1,2,3$ in
(10), the proof follows. ∎

######  Proposition IV.9 (Mean residual life function).

If $X\sim\text{BWeibull}(\boldsymbol{\theta})$ then

$\displaystyle\textrm{MRL}(t;\boldsymbol{\theta})=\dfrac{\displaystyle{2\beta(1+\delta
t)}\,\Gamma\bigg{(}{1\over\alpha},{t^{\alpha}\over\beta^{\alpha}}\bigg{)}-{2\delta\beta^{2}(2+\delta
t)}\,\Gamma\bigg{(}{2\over\alpha},{t^{\alpha}\over\beta^{\alpha}}\bigg{)}+3{\delta^{2}\beta^{2}}\,\Gamma\bigg{(}{3\over\alpha},{t^{\alpha}\over\beta^{\alpha}}\bigg{)}}{\displaystyle\alpha\big{[}1+(1-\delta
t)^{2}\big{]}\,\exp\Bigg{[}-\bigg{(}{t\over\beta}\bigg{)}^{\alpha}\Bigg{]}-{2\delta\beta}\,\bigg{[}\Gamma\bigg{(}{1\over\alpha},{t^{\alpha}\over\beta^{\alpha}}\bigg{)}-{\delta\beta}\,\Gamma\bigg{(}{2\over\alpha},{t^{\alpha}\over\beta^{\alpha}}\bigg{)}\bigg{]}}.$

###### Proof.

Integration by parts gives

$\mathbb{E}\big{[}{1}_{\\{X\geqslant
t\\}}X\big{]}=tR(t;\boldsymbol{\theta})+\int_{t}^{\infty}R(x;\boldsymbol{\theta})\,{\rm
d}x,$

because $xR(x;\boldsymbol{\theta})\to 0$ as $x\to\infty$. Then

$\textrm{MRL}(t;\boldsymbol{\theta})={1\over
R(t;\boldsymbol{\theta})}\,\mathbb{E}\big{[}{1}_{\\{X\geqslant
t\\}}X\big{]}-t,$

where $R(t;\boldsymbol{\theta})$ and $\mathbb{E}\big{[}{1}_{\\{X\geqslant
t\\}}X\big{]}$ are given in Propositions IV.7 and IV.8, respectively. ∎

###  IV.4 Entropies for continuous measure

The Tsallis Tsallis1988 , Quadratic Rao2010  and Shannon Shannon1948
entropies associated with a non-negative random variable $X$ are defined by

$\displaystyle
S_{q}(X)=\dfrac{1}{q-1}\,\bigg{[}1-\int_{0}^{\infty}f^{q}(x;\boldsymbol{\theta})\,{\rm
d}x\bigg{]},\quad q\in\mathbb{R},$ (12) $\displaystyle
H_{2}(X)=-\log\bigg{[}\int_{0}^{\infty}f^{2}(x;\boldsymbol{\theta})\,{\rm
d}x\bigg{]},$ (13) $\displaystyle
H_{1}(X)=-\int_{0}^{\infty}f(x;\boldsymbol{\theta})\log\big{[}f(x;\boldsymbol{\theta})\big{]}\,{\rm
d}x,$ (14)

respectively.

######  Theorem IV.10 (Tsallis entropy).

Let $X\sim\text{BWeibull}(\boldsymbol{\theta})$ and

$\displaystyle\sum_{k=0}^{\infty}\bigg{|}\binom{q}{k}\bigg{|}\int_{0}^{\infty}(1-\delta
x)^{2k}\,\bigg{(}{x\over\beta}\bigg{)}^{q(\alpha-1)}\exp\Bigg{[}-q\bigg{(}{x\over\beta}\bigg{)}^{\alpha}\Bigg{]}\,{\rm
d}x<\infty,$

where $\binom{q}{k}$ is the generalized binomial coefficient and $q\neq 1$. If
$\delta<0$ and $q(\alpha-1)>-1$, then the Tsallis entropy is given by

$\displaystyle
S_{q}(X)=\dfrac{1}{q-1}-\dfrac{\displaystyle\alpha^{q-1}\sum_{k=0}^{\infty}\sum_{l=0}^{2k}\binom{q}{k}\binom{2k}{l}\dfrac{(-\delta)^{l}\beta^{l-1}}{q^{q+(l-q+1)/\alpha}}\,\Gamma\bigg{(}q+{l-q+1\over\alpha}\bigg{)}}{\displaystyle(q-1)\beta^{q}\bigg{[}2+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)}\bigg{]}^{q}}.$

###### Proof.

By definition of Tsallis entropy (12) note that it is enough to prove that

$\displaystyle\int_{0}^{\infty}f^{q}(x;\boldsymbol{\theta})\,{\rm
d}x=\dfrac{\alpha^{q}}{(\beta
Z_{\boldsymbol{\theta}})^{q}}\,\sum_{k=0}^{\infty}\sum_{l=0}^{2k}\binom{q}{k}\binom{2k}{l}\dfrac{(-\delta)^{l}\beta^{l-1}}{\alpha
q^{q+(l-q+1)/\alpha}}\,\Gamma\bigg{(}q+{l-q+1\over\alpha}\bigg{)},$ (15)

whenever $q(\alpha-1)>-1$, where $Z_{\boldsymbol{\theta}}$ is the
normalization constant in equation (2).

Indeed, since

$\displaystyle\int_{0}^{\infty}f^{q}(x;\boldsymbol{\theta})\,{\rm
d}x=\int_{0}^{\infty}{\alpha^{q}\over(\beta
Z_{\boldsymbol{\theta}})^{q}}\,\big{[}1+(1-\delta
x)^{2}\big{]}^{q}\,\bigg{(}{x\over\beta}\bigg{)}^{q(\alpha-1)}\exp\Bigg{[}-q\bigg{(}{x\over\beta}\bigg{)}^{\alpha}\Bigg{]}\,{\rm
d}x,$

by using the Newton’s generalized binomial theorem the above integral is

$\displaystyle={\alpha^{q}\over(\beta
Z_{\boldsymbol{\theta}})^{q}}\,\int_{0}^{\infty}\sum_{k=0}^{\infty}\binom{q}{k}(1-\delta
x)^{2k}\,\bigg{(}{x\over\beta}\bigg{)}^{q(\alpha-1)}\exp\Bigg{[}-q\bigg{(}{x\over\beta}\bigg{)}^{\alpha}\Bigg{]}\,{\rm
d}x$

where the condition $\delta<0$ guarantees that $x\neq\delta/2$, and therefore
the convergence of the above series. Seeing that
$\sum_{k=0}^{\infty}\big{|}\binom{q}{k}\big{|}\int_{0}^{\infty}(1-\delta
x)^{2k}\,\big{(}{x\over\beta}\big{)}^{q(\alpha-1)}\exp\big{[}-q\big{(}{x\over\beta}\big{)}^{\alpha}\big{]}\,{\rm
d}x<\infty$ (hypothesis), by Fubini’s theorem knapp  the above integral is

$\displaystyle={\alpha^{q}\over(\beta
Z_{\boldsymbol{\theta}})^{q}}\,\sum_{k=0}^{\infty}\binom{q}{k}\int_{0}^{\infty}(1-\delta
x)^{2k}\,\bigg{(}{x\over\beta}\bigg{)}^{q(\alpha-1)}\exp\Bigg{[}-q\bigg{(}{x\over\beta}\bigg{)}^{\alpha}\Bigg{]}\,{\rm
d}x$ $\displaystyle={\alpha^{q}\over(\beta
Z_{\boldsymbol{\theta}})^{q}}\,\sum_{k=0}^{\infty}\sum_{l=0}^{2k}\binom{q}{k}\binom{2k}{l}(-\delta)^{l}\int_{0}^{\infty}x^{l}\,\bigg{(}{x\over\beta}\bigg{)}^{q(\alpha-1)}\exp\Bigg{[}-q\bigg{(}{x\over\beta}\bigg{)}^{\alpha}\Bigg{]}\,{\rm
d}x,$ (16)

where in the second equality we have used the classic binomial expansion. By
using the formula of Item (6) in Cankaya2018 :

$\displaystyle\Gamma\bigg{(}s+{1\over\alpha}\bigg{)}=\alpha p^{\alpha
s+1}\int_{0}^{\infty}y^{\alpha s}\,\exp\big{[}-(yp)^{\alpha}\big{]}\,{\rm
d}y,$

the expression in (IV.4) is rewritten as

$\displaystyle=\dfrac{\alpha^{q}}{(\beta
Z_{\boldsymbol{\theta}})^{q}}\,\sum_{k=0}^{\infty}\sum_{l=0}^{2k}\binom{q}{k}\binom{2k}{l}\dfrac{(-\delta)^{l}\beta^{l-1}}{\alpha
q^{q+(l-q+1)/\alpha}}\,\Gamma\bigg{(}q+{l-q+1\over\alpha}\bigg{)},$

whenever $q(\alpha-1)>-1$. Therefore, combining the above identities we have
proved that the statement in (15) is valid. This completes the proof of
theorem. ∎

######  Theorem IV.11 (Quadratic entropy).

Let $X\sim\text{BWeibull}(\boldsymbol{\theta})$. If $\alpha>1/2$ then the
quadratic entropy can be written as

$\displaystyle H_{2}(X)$
$\displaystyle=2\log(\beta)-\log(\alpha)+2\log{\displaystyle\left[2+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)}\right]}$
$\displaystyle-\log\left[{\displaystyle{2^{1/\alpha}\over\beta}\,\Gamma\bigg{(}2-{1\over\alpha}\bigg{)}+{2\delta^{2}\beta\over
2^{1/\alpha}}\,\Gamma\bigg{(}2+{1\over\alpha}\bigg{)}-{\delta^{3}\beta^{2}\over
2^{2/\alpha}}\,\Gamma\bigg{(}2+{2\over\alpha}\bigg{)}+{\delta^{4}\beta^{3}\over
2^{2+3/\alpha}}\,\Gamma\bigg{(}2+{3\over\alpha}\bigg{)}-2\delta}\right].$

###### Proof.

A simple algebraic manipulation shows that

$\displaystyle\int_{0}^{\infty}f^{2}(x;\boldsymbol{\theta})\,{\rm d}x$
$\displaystyle=\int_{0}^{\infty}{\alpha^{2}\over\beta^{2}Z_{\boldsymbol{\theta}}^{2}}\,\big{[}1+(1-\delta
x)^{2}\big{]}^{2}\,\bigg{(}{x\over\beta}\bigg{)}^{2(\alpha-1)}\exp\Bigg{[}-2\bigg{(}{x\over\beta}\bigg{)}^{\alpha}\Bigg{]}\,{\rm
d}x$
$\displaystyle={\alpha^{2}\over\beta^{2}Z_{\boldsymbol{\theta}}^{2}}\,\int_{0}^{\infty}\big{(}4-8\delta
x+8\delta^{2}x^{2}-4\delta^{3}x^{3}+\delta^{4}x^{4}\big{)}\bigg{(}{x\over\beta}\bigg{)}^{2(\alpha-1)}\exp\Bigg{[}-\bigg{(}{x\over\beta}\,2^{1/\alpha}\bigg{)}^{\alpha}\Bigg{]}\,{\rm
d}x,$

where $Z_{\boldsymbol{\theta}}$ is as in (2). By using the formula of Item (6)
in reference Cankaya2018 :

$\displaystyle\Gamma\bigg{(}s+{1\over\alpha}\bigg{)}=\alpha p^{\alpha
s+1}\int_{0}^{\infty}y^{\alpha s}\,\exp\big{[}-(yp)^{\alpha}\big{]}\,{\rm
d}y,$

we have

$\displaystyle{\alpha^{2}\over\beta^{2}Z_{\boldsymbol{\theta}}^{2}}\,\int_{0}^{\infty}\big{(}4-8\delta
x+8\delta^{2}x^{2}-4\delta^{3}x^{3}+\delta^{4}x^{4}\big{)}\bigg{(}{x\over\beta}\bigg{)}^{2(\alpha-1)}\exp\Bigg{[}-\bigg{(}{x\over\beta}\,2^{1/\alpha}\bigg{)}^{\alpha}\Bigg{]}\,{\rm
d}x$
$\displaystyle={\alpha^{2}\over\beta^{2}Z_{\boldsymbol{\theta}}^{2}}\,\left[{\displaystyle
4v_{\boldsymbol{\theta}}(0)-8\delta
v_{\boldsymbol{\theta}}(1)+8\delta^{2}v_{\boldsymbol{\theta}}(2)-4\delta^{3}v_{\boldsymbol{\theta}}(3)+\delta^{4}v_{\boldsymbol{\theta}}(4)}\right],$

where

$\displaystyle v_{\boldsymbol{\theta}}(s)$
$\displaystyle=\int_{0}^{\infty}x^{s}\bigg{(}{x\over\beta}\bigg{)}^{2(\alpha-1)}\exp\Bigg{[}-\bigg{(}{x\over\beta}\,2^{1/\alpha}\bigg{)}^{\alpha}\Bigg{]}\,{\rm
d}x$ $\displaystyle=\dfrac{\beta^{s-1}}{\alpha
2^{2+(s-1)/\alpha}}\,\Gamma\bigg{(}2+{s-1\over\alpha}\bigg{)},\quad
s>1-2\alpha.$

Therefore,

$\displaystyle\int_{0}^{\infty}f^{2}(x;\boldsymbol{\theta})\,{\rm
d}x={\alpha^{2}\over\beta^{2}Z^{2}}\,\left[{\displaystyle
4v_{\boldsymbol{\theta}}(0)-8\delta
v_{\boldsymbol{\theta}}(1)+8\delta^{2}v_{\boldsymbol{\theta}}(2)-4\delta^{3}v_{\boldsymbol{\theta}}(3)+\delta^{4}v_{\boldsymbol{\theta}}(4)}\right].$

Finally, taking logarithms to both sides of the above equation and then
multiplying by $-1$, by definition (13) of quadratic entropy and by standard
algebraic manipulations, we complete the proof.

∎

Let $Y\sim\text{Weibull}(\alpha,\beta)$. By taking the change of variable
$x=\big{(}{y\over\beta}\big{)}^{\alpha}$, for each $s>-\alpha$, we have

$\displaystyle\mathbb{E}[Y^{s}\log(Y)]$
$\displaystyle={\alpha\over\beta}\,\int_{0}^{\infty}y^{s}\log(y)\,\bigg{(}{y\over\beta}\bigg{)}^{\alpha-1}\exp\Bigg{[}-\bigg{(}{y\over\beta}\bigg{)}^{\alpha}\Bigg{]}\,{\rm
d}y$
$\displaystyle=\beta^{s}\log(\beta)\int_{0}^{\infty}x^{s\over\alpha}\exp(-x)\,{\rm
d}x+{\beta^{s}\over\alpha}\,\int_{0}^{\infty}\log(x)x^{s\over\alpha}\exp(-x)\,{\rm
d}x.$

Since

$\displaystyle\int_{0}^{\infty}x^{s\over\alpha}\exp(-x)\,{\rm
d}x=\Gamma\bigg{(}1+\dfrac{s}{\alpha}\bigg{)},$
$\displaystyle\int_{0}^{\infty}\log(x)x^{s\over\alpha}\exp(-x)\,{\rm
d}x=\Gamma\bigg{(}1+\dfrac{s}{\alpha}\bigg{)}\Psi^{(0)}\bigg{(}1+\dfrac{s}{\alpha}\bigg{)},$

the expectation $\mathbb{E}[Y^{s}\log(Y)]$ can be written as follows

$\displaystyle\mathbb{E}[Y^{s}\log(Y)]=\beta^{s}\Gamma\bigg{(}1+\dfrac{s}{\alpha}\bigg{)}\bigg{[}\log(\beta)+{1\over\alpha}\,\Psi^{(0)}\bigg{(}1+\dfrac{s}{\alpha}\bigg{)}\bigg{]},$
(17)

where $\Psi^{(m)}(z)={{\rm d}^{m+1}\log[\Gamma(z)]\over{\rm d}z^{m+1}}$ is the
the polygamma function of order $m$.

######  Proposition IV.12.

If $X\sim\text{BWeibull}(\boldsymbol{\theta})$ then

$\displaystyle\mathbb{E}[\log(X)]$ $\displaystyle=\dfrac{\displaystyle
2\bigg{[}\log(\beta)-{\gamma\over\alpha}\bigg{]}}{\displaystyle
2+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)}}$
$\displaystyle-\dfrac{\displaystyle
2\delta\beta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)}\bigg{[}\log(\beta)+{1\over\alpha}\,\Psi^{(0)}\bigg{(}1+{1\over\alpha}\bigg{)}\bigg{]}-\delta^{2}\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}\bigg{[}\log(\beta)+{1\over\alpha}\,\Psi^{(0)}\bigg{(}1+{2\over\alpha}\bigg{)}\bigg{]}}{\displaystyle
2+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)}},$

where $\gamma=-{{\rm d}\Gamma(x)\over{\rm d}x}\big{|}_{x=1}\approx 0.57721$ is
the Euler-Mascheroni constant.

###### Proof.

Through a simple calculation we have

$\displaystyle\mathbb{E}[\log(X)]=\dfrac{1}{Z_{\boldsymbol{\theta}}}\,\left\\{2\,\mathbb{E}\big{[}\log(Y)\big{]}-2\delta\,\mathbb{E}\big{[}Y\log(Y)\big{]}+\delta^{2}\,\mathbb{E}\big{[}Y^{2}\log(Y)\big{]}\right\\},\quad
Y\sim\text{Weibull}(\alpha,\beta).$

Then, by taking $s=0,1,2$ in (17), from the above identity, the statement of
proposition follows. ∎

######  Theorem IV.13 (Shannon entropy).

Let $X\sim\text{BWeibull}(\boldsymbol{\theta})$, $\mathbb{E}(X)=\mu$,
$\mathrm{Var}(X)=\sigma^{2}$, $G(x)=1+(1-\delta x)^{2}$ for $x\geqslant 0$,
and $\mu_{G}=\mathbb{E}[G(X)]=1+(1-\delta\mu)^{2}+\delta^{2}\sigma^{2}$. If

$\displaystyle\sum_{n=1}^{\infty}\dfrac{1}{n\mu_{G}^{n}}\,\Big{|}\mathbb{E}\big{[}\big{(}G(X)-\mu_{G}\big{)}^{n}\big{]}\Big{|}<\infty,$

then the Shannon entropy can be written as

$\displaystyle
H_{1}(X)=\alpha\log(\beta)-\log(\alpha)-\log\big{[}1+(1-\delta\mu)^{2}+\delta^{2}\sigma^{2}\big{]}$
$\displaystyle+\log\bigg{[}2+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)}\bigg{]}+\dfrac{\displaystyle
2+\delta^{2}\beta^{2}\Gamma\bigg{(}2+{2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}2+{1\over\alpha}\bigg{)}}{\displaystyle
2+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)}}$
$\displaystyle-\sum_{n=1}^{\infty}\sum_{k=0}^{n}\sum_{j=0}^{k}\sum_{l=0}^{2j}\binom{n}{k}\binom{k}{j}\binom{2j}{l}\,\dfrac{(-1)^{l-k+1}(\beta\delta)^{l}}{n\big{[}1+(1-\delta\mu)^{2}+\delta^{2}\sigma^{2}\big{]}^{k}}\,\dfrac{\displaystyle
2\Gamma\bigg{(}1+{l\over\alpha}\bigg{)}+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{l+2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{l+1\over\alpha}\bigg{)}}{\displaystyle
2+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)}}$
$\displaystyle-(\alpha-1)\,\dfrac{\displaystyle
2\bigg{[}\log(\beta)-{\gamma\over\alpha}\bigg{]}-2\delta\beta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)}\bigg{[}\log(\beta)+{1\over\alpha}\,\Psi^{(0)}\bigg{(}1+{1\over\alpha}\bigg{)}\bigg{]}+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}\bigg{[}\log(\beta)+{1\over\alpha}\,\Psi^{(0)}\bigg{(}1+{2\over\alpha}\bigg{)}\bigg{]}}{\displaystyle
2+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)}}.$

###### Proof.

The Shannon entropy formula for the random variable
$X\sim\text{BWeibull}(\boldsymbol{\theta})$, as defined in 14, is given by

$H_{1}(X)=-\mathbb{E}\big{[}\log f(X;\boldsymbol{\theta})\big{]}.$

Developing the logarithm $\log[f(X;\boldsymbol{\theta})]$ in the above
identity, we have

$\displaystyle
H_{1}(X)=\log(Z_{\boldsymbol{\theta}})+\alpha\log(\beta)-\log(\alpha)-\mathbb{E}\big{[}\log
G(X)\big{]}-(\alpha-1)\mathbb{E}[\log(X)]+{1\over\beta^{\alpha}}\,\mathbb{E}(X^{\alpha}),$
(18)

where $\mathbb{E}[\log(X)]$ and $\mathbb{E}(X^{\alpha})$ are given in
Propositions IV.12 and IV.1, respectively.

It is worth mentioning that the expectation $\mathbb{E}[\log G(X)]$ cannot be
expressed through known mathematical functions. To obtain an expression for
$\mathbb{E}[\log G(X)]$ we consider Taylor’s expansion of the function
$\log[G(X)]$ around the point
$\mu_{G}=\mathbb{E}[G(X)]=1+(1-\delta\mu)^{2}+\delta^{2}\sigma^{2}$. Hence

$\displaystyle\log[G(X)]=\log(\mu_{G})+\sum_{n=1}^{\infty}\dfrac{(-1)^{n+1}}{n\mu_{G}^{n}}\,\big{(}G(X)-\mu_{G}\big{)}^{n}.$
(19)

Taking expectation on both sides of (19) we get

$\displaystyle\mathbb{E}[\log
G(X)]=\log(\mu_{G})+\mathbb{E}\left[\sum_{n=1}^{\infty}\dfrac{(-1)^{n+1}}{n\mu_{G}^{n}}\,\big{(}G(X)-\mu_{G}\big{)}^{n}\right].$
(20)

Since $0\leqslant\mathbb{E}[\log G(X)]\leqslant\mathbb{E}\big{[}(1-\delta
X)^{2}\big{]}<\infty$ and
$\sum_{n=1}^{\infty}\frac{1}{n\mu_{G}^{n}}\,\big{|}\mathbb{E}\big{[}\big{(}G(X)-\mu_{G}\big{)}^{n}\big{]}\big{|}<\infty$
(hypothesis), by applying Fubini’s Theorem knapp , the expression on the
right-hand side of (20) is

$\displaystyle=\log(\mu_{G})+\sum_{n=1}^{\infty}\dfrac{(-1)^{n+1}}{n\mu_{G}^{n}}\,\mathbb{E}\big{[}\big{(}G(X)-\mu_{G}\big{)}^{n}\big{]}.$

By using binomial expansion repeatedly in the above expression, this is

$\displaystyle=\log(\mu_{G})+\sum_{n=1}^{\infty}\sum_{k=0}^{n}\sum_{j=0}^{k}\sum_{l=0}^{2j}\binom{n}{k}\binom{k}{j}\binom{2j}{l}\,\dfrac{(-1)^{l-k+1}\delta^{l}}{n\mu_{G}^{k}}\,\mathbb{E}(X^{l}).$

Therefore, we get

$\displaystyle\mathbb{E}[\log
G(X)]=\log(\mu_{G})+\sum_{n=1}^{\infty}\sum_{k=0}^{n}\sum_{j=0}^{k}\sum_{l=0}^{2j}\binom{n}{k}\binom{k}{j}\binom{2j}{l}\,\dfrac{(-1)^{l-k+1}\delta^{l}}{n\mu_{G}^{k}}\,\mathbb{E}(X^{l}).$
(21)

Finally, by combining (21), Proposition IV.1 and (18), the proof of theorem
follows. ∎

##  V Inference: Estimation Methods and Fisher information

###  V.1 Maximum likelihood estimation method

Maximum likelihood estimation (MLE) method is the standard approach for
estimations of parameters of $f(x;\boldsymbol{\theta})$, mainly due to the
desirable asymptotic properties of consistency, efficiency and asymptotic
normality Fer196 ; Fer296 .

The log likelihood function for $\boldsymbol{\theta}$ is given by

$\displaystyle l(\boldsymbol{\theta};\boldsymbol{x})$
$\displaystyle=\sum_{i=1}^{n}\left\\{\log\Big{(}\frac{\alpha}{\beta
Z_{\boldsymbol{\theta}}}\Big{)}+\log\left[1+(1-\delta
x_{i})^{2}\right]+(\alpha-1)\log\Big{(}\frac{x_{i}}{\beta}\Big{)}-\Big{(}{x_{i}\over\beta}\Big{)}^{\alpha}\right\\}$
(22)
$\displaystyle=-n\log(Z_{\boldsymbol{\theta}})+\sum_{i=1}^{n}\log\left[1+(1-\delta
x_{i})^{2}\right]+l_{Y}(\alpha,\beta;\boldsymbol{x}),$

where
$Z_{\boldsymbol{\theta}}=2+\delta^{2}\beta^{2}\Gamma\big{(}1+{2\over\alpha}\big{)}-2\delta\beta\Gamma\big{(}1+{1\over\alpha}\big{)}$
is as in (2) and $l_{Y}(\alpha,\beta;\boldsymbol{x})$,
$Y\sim\text{Weibull}(\alpha,\beta)$, is the log likelihood function for the
parameter vector $(\alpha,\beta)$.

A standard calculation shows that the first-order partial derivatives of
$l(\boldsymbol{\theta};\boldsymbol{x})$ are

$\displaystyle\begin{array}[]{llllll}\displaystyle{\partial
l(\boldsymbol{\theta};\boldsymbol{x})\over\partial\alpha}&=\displaystyle-{n\over
Z_{\boldsymbol{\theta}}}\,{\partial
Z_{\boldsymbol{\theta}}\over\partial\alpha}+{\partial
l_{Y}(\alpha,\beta;\boldsymbol{x})\over\partial\alpha},\\\\[14.22636pt]
\displaystyle{\partial
l(\boldsymbol{\theta};\boldsymbol{x})\over\partial\beta}&=\displaystyle-{n\over
Z_{\boldsymbol{\theta}}}\,{\partial
Z_{\boldsymbol{\theta}}\over\partial\beta}+{\partial
l_{Y}(\alpha,\beta;\boldsymbol{x})\over\partial\beta},\\\\[14.22636pt]
\displaystyle{\partial
l(\boldsymbol{\theta};\boldsymbol{x})\over\partial\delta}&=\displaystyle-{n\over
Z_{\boldsymbol{\theta}}}\,{\partial
Z_{\boldsymbol{\theta}}\over\partial\delta}-\sum_{i=1}^{n}{2\delta(1-\delta
x_{i})\over 1+(1-\delta x_{i})^{2}},\end{array}$ (26)

where

$\displaystyle{\partial
l_{Y}(\alpha,\beta;\boldsymbol{x})\over\partial\alpha}$
$\displaystyle={n\over\alpha}-n\log(\beta)+\sum_{i=1}^{n}\log(x_{i})-{1\over\beta^{\alpha}}\sum_{i=1}^{n}x_{i}^{\alpha}\log(x_{i}),$
$\displaystyle{\partial l_{Y}(\alpha,\beta;\boldsymbol{x})\over\partial\beta}$
$\displaystyle=-{n\alpha\over\beta}+{\alpha\over\beta^{\alpha+1}}\sum_{i=1}^{n}x_{i}^{\alpha},$

and

$\displaystyle{\partial Z_{\boldsymbol{\theta}}\over\partial\alpha}$
$\displaystyle=-{2\delta^{2}\beta^{2}\over\alpha^{2}}\,\Psi^{(0)}\bigg{(}1+{2\over\alpha}\bigg{)}+{2\delta\beta\over\alpha^{2}}\,\Psi^{(0)}\bigg{(}1+{1\over\alpha}\bigg{)},$
$\displaystyle{\partial Z_{\boldsymbol{\theta}}\over\partial\beta}$
$\displaystyle=2\delta^{2}\beta\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}-2\delta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)},$
$\displaystyle{\partial Z_{\boldsymbol{\theta}}\over\partial\delta}$
$\displaystyle=2\delta\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}-2\beta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)}.$

Here, $\Psi^{(m)}(z)={{\rm d}^{m+1}\log[\Gamma(z)]\over{\rm d}z^{m+1}}$ is the
the polygamma function of order $m$.

The second-order partial derivatives of
$l(\boldsymbol{\theta};\boldsymbol{x})$ can be written as

$\displaystyle{\partial^{2}l(\boldsymbol{\theta};\boldsymbol{x})\over\partial\alpha^{2}}$
$\displaystyle={n\over Z_{\boldsymbol{\theta}}^{2}}\,\bigg{(}{\partial
Z_{\boldsymbol{\theta}}\over\partial\alpha}\bigg{)}^{2}-{n\over
Z_{\boldsymbol{\theta}}}\,{\partial^{2}Z_{\boldsymbol{\theta}}\over\partial\alpha^{2}}+{\partial^{2}l_{Y}(\alpha,\beta;\boldsymbol{x})\over\partial\alpha^{2}},$
$\displaystyle{\partial^{2}l(\boldsymbol{\theta};\boldsymbol{x})\over\partial\beta^{2}}$
$\displaystyle={n\over Z_{\boldsymbol{\theta}}^{2}}\,\bigg{(}{\partial
Z_{\boldsymbol{\theta}}\over\partial\beta}\bigg{)}^{2}-{n\over
Z_{\boldsymbol{\theta}}}\,{\partial^{2}Z_{\boldsymbol{\theta}}\over\partial\beta^{2}}+{\partial^{2}l_{Y}(\alpha,\beta;\boldsymbol{x})\over\partial\beta^{2}},$
$\displaystyle{\partial^{2}l(\boldsymbol{\theta};\boldsymbol{x})\over\partial\delta^{2}}$
$\displaystyle={n\over Z_{\boldsymbol{\theta}}^{2}}\,\bigg{(}{\partial
Z_{\boldsymbol{\theta}}\over\partial\delta}\bigg{)}^{2}-{n\over
Z_{\boldsymbol{\theta}}}\,{\partial^{2}Z_{\boldsymbol{\theta}}\over\partial\delta^{2}}-2\sum_{i=1}^{n}{(1-2\delta
x_{i})[1+(1-\delta x_{i})^{2}]+2\delta^{2}(1-\delta
x_{i})^{2}\over[1+(1-\delta x_{i})^{2}]^{2}},$
$\displaystyle{\partial^{2}l(\boldsymbol{\theta};\boldsymbol{x})\over\partial\alpha\partial\beta}$
$\displaystyle={\partial^{2}l(\boldsymbol{\theta};\boldsymbol{x})\over\partial\beta\partial\alpha}={n\over
Z_{\boldsymbol{\theta}}^{2}}\,{\partial
Z_{\boldsymbol{\theta}}\over\partial\alpha}\,{\partial
Z_{\boldsymbol{\theta}}\over\partial\beta}-{n\over
Z_{\boldsymbol{\theta}}}\,{\partial^{2}Z_{\boldsymbol{\theta}}\over\partial\alpha\partial\beta}+{\partial^{2}l_{Y}(\alpha,\beta;\boldsymbol{x})\over\partial\alpha\partial\beta},$
$\displaystyle{\partial^{2}l(\boldsymbol{\theta};\boldsymbol{x})\over\partial\alpha\partial\delta}$
$\displaystyle={\partial^{2}l(\boldsymbol{\theta};\boldsymbol{x})\over\partial\delta\partial\alpha}={n\over
Z_{\boldsymbol{\theta}}^{2}}\,{\partial
Z_{\boldsymbol{\theta}}\over\partial\alpha}\,{\partial
Z_{\boldsymbol{\theta}}\over\partial\delta}-{n\over
Z_{\boldsymbol{\theta}}}\,{\partial^{2}Z_{\boldsymbol{\theta}}\over\partial\alpha\partial\delta},$
$\displaystyle{\partial^{2}l(\boldsymbol{\theta};\boldsymbol{x})\over\partial\beta\partial\delta}$
$\displaystyle={\partial^{2}l(\boldsymbol{\theta};\boldsymbol{x})\over\partial\delta\partial\beta}={n\over
Z_{\boldsymbol{\theta}}^{2}}\,{\partial
Z_{\boldsymbol{\theta}}\over\partial\beta}\,{\partial
Z_{\boldsymbol{\theta}}\over\partial\delta}-{n\over
Z_{\boldsymbol{\theta}}}\,{\partial^{2}Z_{\boldsymbol{\theta}}\over\partial\beta\partial\delta},$

where

$\displaystyle{\partial^{2}l_{Y}(\alpha,\beta;\boldsymbol{x})\over\partial\alpha^{2}}$
$\displaystyle=-{n\over\alpha^{2}}+{\log(\beta)\over\beta^{\alpha}}\sum_{i=1}^{n}x_{i}^{\alpha}\log(x_{i})-{1\over\beta^{\alpha}}\sum_{i=1}^{n}x_{i}^{\alpha}\log^{2}(x_{i}),$
$\displaystyle{\partial^{2}l_{Y}(\alpha,\beta;\boldsymbol{x})\over\partial\beta^{2}}$
$\displaystyle={n\alpha\over\beta^{2}}-{\alpha(\alpha+1)\over\beta^{\alpha+2}}\sum_{i=1}^{n}x_{i}^{\alpha},$
$\displaystyle{\partial^{2}l_{Y}(\alpha,\beta;\boldsymbol{x})\over\partial\alpha\partial\beta}$
$\displaystyle={\partial^{2}l_{Y}(\alpha,\beta;\boldsymbol{x})\over\partial\beta\partial\alpha}=-{n\over\beta}+{\alpha\over\beta^{\alpha+1}}\sum_{i=1}^{n}x_{i}^{\alpha}\log(x_{i}),$

and

$\displaystyle\begin{array}[]{lllll}\displaystyle{\partial^{2}Z_{\boldsymbol{\theta}}\over\partial\alpha^{2}}&=\displaystyle{4\delta^{2}\beta^{2}\over\alpha^{3}}\,\Psi^{(0)}\bigg{(}1+{2\over\alpha}\bigg{)}+{4\delta^{2}\beta^{2}\over\alpha^{4}}\,\Psi^{(1)}\bigg{(}1+{2\over\alpha}\bigg{)}-{4\delta\beta\over\alpha^{3}}\,\Psi^{(0)}\bigg{(}1+{1\over\alpha}\bigg{)}-{2\delta\beta\over\alpha^{4}}\,\Psi^{(1)}\bigg{(}1+{1\over\alpha}\bigg{)},\\\\[14.22636pt]
\displaystyle{\partial^{2}Z_{\boldsymbol{\theta}}\over\partial\beta^{2}}&=\displaystyle
2\delta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)},\\\\[14.22636pt]
\displaystyle{\partial^{2}Z_{\boldsymbol{\theta}}\over\partial\delta^{2}}&=\displaystyle
2\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)},\\\\[14.22636pt]
\displaystyle{\partial^{2}Z_{\boldsymbol{\theta}}\over\partial\alpha\partial\beta}&=\displaystyle{\partial^{2}Z_{\boldsymbol{\theta}}\over\partial\beta\partial\alpha}=-{4\delta^{2}\beta\over\alpha^{2}}\,\Psi^{(0)}\bigg{(}1+{2\over\alpha}\bigg{)}+{2\delta\over\alpha^{2}}\,\Psi^{(0)}\bigg{(}1+{1\over\alpha}\bigg{)},\\\\[14.22636pt]
{\partial^{2}Z_{\boldsymbol{\theta}}\over\partial\alpha\partial\delta}&=\displaystyle{\partial^{2}Z_{\boldsymbol{\theta}}\over\partial\delta\partial\alpha}=-{4\delta\beta^{2}\over\alpha^{2}}\,\Psi^{(0)}\bigg{(}1+{2\over\alpha}\bigg{)}+{2\beta\over\alpha^{2}}\,\Psi^{(0)}\bigg{(}1+{1\over\alpha}\bigg{)},\\\\[14.22636pt]
\displaystyle{\partial^{2}Z_{\boldsymbol{\theta}}\over\partial\beta\partial\delta}&=\displaystyle{\partial^{2}Z_{\boldsymbol{\theta}}\over\partial\delta\partial\beta}=4\delta\beta\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}-2\Gamma\bigg{(}1+{1\over\alpha}\bigg{)}.\end{array}$
(33)

###  V.2 Maximum $\log_{q}$ likelihood estimation method

The generalization of MLE with $\log$ form is defined as maximum $\log_{q}$
likelihood estimation (MLqE). $\log_{q}$ is $q$-deformed logarithm of natural
logarithm ($\log$) in MLE Tsallisbook09 . The concavity property of $\log_{q}$
guarantees to replace $\log$ in MLE by $\log_{q}$ Lindsay94 ; CanKor18 ; Jan1
; Jan2 .

$\displaystyle
l_{q}(\boldsymbol{\theta};\boldsymbol{x})=\sum_{i=1}^{n}\log_{q}\big{[}f(x_{i};\boldsymbol{\theta})\big{]},$
(34)

where $\log_{q}(f)={f^{1-q}-1\over 1-q},q\in\mathbb{R}\backslash\\{1\\}$. If
$q\rightarrow 1$, then we have $\log$ and so equation (22) is dropped. When
$\log_{q}$ likelihood is compared with $\log$ likelihood, the deformation
gives an advantage to manage the modelling capability of a parametric model.
$\log_{q}(f)$ as an objective or loss function in MLqE in the M-estimation is
used to manage efficiency and robustness for contamination. Note that
$\log_{q}(f)$ is M-function MLqEBio ; CanKor18 . This estimation method is
applied to BWeibull$(\boldsymbol{\theta})$.

The $l_{q}$ function in equation (34) is optimized according to parameters
$\boldsymbol{\theta}$ in order to get the estimators
$\widehat{\boldsymbol{\theta}}$ of $\boldsymbol{\theta}$. Alternatively, a
system of estimating equations can be solved simultaneously according to the
corresponding estimating equation for the parameters ${\alpha}$, ${\beta}$ and
${\delta}$. We omit to rewrite the system
$\sum_{i=1}^{n}f(x_{i};\boldsymbol{\theta})^{1-q}\,{\partial\log[f(x_{i};\boldsymbol{\theta})]\over\partial\boldsymbol{\theta}}=\boldsymbol{0}$.
Note that
${\partial\log[f(x_{i};\boldsymbol{\theta})]\over\partial\boldsymbol{\theta}}$
is given by equation (26) (with $n=1$) for $\alpha,\beta$ and $\delta$.

###  V.3 Fisher information based on $\log$ and $\log_{q}$

Fisher information matrix (FI) is given by

$\mathbb{E}\left[\bigg{(}{\partial
l\over\partial\boldsymbol{\theta}}\bigg{)}\bigg{(}{\partial
l\over\partial\boldsymbol{\theta}}\bigg{)}^{\intercal}\,\right]=\mathbb{E}\left[{\partial^{2}l\over\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{\intercal}}\right]$
(35)

and it is a tool to provide the variances of $\widehat{\boldsymbol{\theta}}$,
i.e., Var($\widehat{\boldsymbol{\theta}}$). As it is well-known, the inverse
of FI gives Var($\widehat{\boldsymbol{\theta}}$) if the inverse of FI matrix
exists Fisher25 ; LehmannCas98 . If the FI matrix is singular, then the
generalized inverse techniques are used to get the inverse of FI even though
we have information loss due to the used generalized inverse LiYeh12gen .

FI matrix based on $\log_{q}$ is given by the following form:

$F=n\begin{bmatrix}E_{\alpha\alpha}&E_{\alpha\beta}&E_{\alpha\delta}\\\
E_{\beta\alpha}&E_{\beta\beta}&E_{\beta\delta}\\\
E_{\delta\alpha}&E_{\delta\beta}&E_{\delta\delta}\end{bmatrix},$ (36)

where $n$ is sample size. $E$ is integral for partial derivatives of $\log(L)$
according to parameters and it is taken over probability density function
$f(x;\boldsymbol{\theta})$. The subscript in $E$ represents second-order
partial derivatives of $\log(L)$ according to parameters $\alpha,\beta$ and
$\delta$. In other words, if $X\sim\text{BWeibull}(\boldsymbol{\theta})$ then

$\displaystyle E_{\alpha\alpha}$ $\displaystyle={1\over
Z_{\boldsymbol{\theta}}^{2}}\,\bigg{(}{\partial
Z_{\boldsymbol{\theta}}\over\partial\alpha}\bigg{)}^{2}-{1\over
Z_{\boldsymbol{\theta}}}\,{\partial^{2}Z_{\boldsymbol{\theta}}\over\partial\alpha^{2}}+-{1\over\alpha^{2}}+{\log(\beta)\over\beta^{\alpha}}\,\mathbb{E}\big{[}X^{\alpha}\log(X)\big{]}-{1\over\beta^{\alpha}}\,\mathbb{E}\big{[}X^{\alpha}\log^{2}(X)\big{]},$
$\displaystyle E_{\beta\beta}$ $\displaystyle={1\over
Z_{\boldsymbol{\theta}}^{2}}\,\bigg{(}{\partial
Z_{\boldsymbol{\theta}}\over\partial\beta}\bigg{)}^{2}-{1\over
Z_{\boldsymbol{\theta}}}\,{\partial^{2}Z_{\boldsymbol{\theta}}\over\partial\beta^{2}}+{\alpha\over\beta^{2}}-{\alpha(\alpha+1)\over\beta^{\alpha+2}}\,\mathbb{E}\big{(}X^{\alpha}\big{)},$
$\displaystyle E_{\delta\delta}$ $\displaystyle={1\over
Z_{\boldsymbol{\theta}}^{2}}\,\bigg{(}{\partial
Z_{\boldsymbol{\theta}}\over\partial\delta}\bigg{)}^{2}-{1\over
Z_{\boldsymbol{\theta}}}\,{\partial^{2}Z_{\boldsymbol{\theta}}\over\partial\delta^{2}}-2\mathbb{E}\left[{(1-2\delta
X)\big{(}1+(1-\delta X)^{2}\big{)}+2\delta^{2}(1-\delta
X)^{2}\over\big{(}1+(1-\delta X)^{2}\big{)}^{2}}\right],$ $\displaystyle
E_{\alpha\beta}$ $\displaystyle=E_{\beta\alpha}={1\over
Z_{\boldsymbol{\theta}}^{2}}\,{\partial
Z_{\boldsymbol{\theta}}\over\partial\alpha}\,{\partial
Z_{\boldsymbol{\theta}}\over\partial\beta}-{1\over
Z_{\boldsymbol{\theta}}}\,{\partial^{2}Z_{\boldsymbol{\theta}}\over\partial\alpha\partial\beta}-{1\over\beta}+{\alpha\over\beta^{\alpha+1}}\,\mathbb{E}\big{[}X^{\alpha}\log(X)\big{]},$
$\displaystyle E_{\alpha\delta}$ $\displaystyle=E_{\delta\alpha}={1\over
Z_{\boldsymbol{\theta}}^{2}}\,{\partial
Z_{\boldsymbol{\theta}}\over\partial\alpha}\,{\partial
Z_{\boldsymbol{\theta}}\over\partial\delta}-{1\over
Z_{\boldsymbol{\theta}}}\,{\partial^{2}Z_{\boldsymbol{\theta}}\over\partial\alpha\partial\delta},$
$\displaystyle E_{\beta\delta}$ $\displaystyle=E_{\delta\beta}={1\over
Z_{\boldsymbol{\theta}}^{2}}\,{\partial
Z_{\boldsymbol{\theta}}\over\partial\beta}\,{\partial
Z_{\boldsymbol{\theta}}\over\partial\delta}-{1\over
Z_{\boldsymbol{\theta}}}\,{\partial^{2}Z_{\boldsymbol{\theta}}\over\partial\beta\partial\delta},$

where
${\partial^{2}Z_{\boldsymbol{\theta}}\over\partial\theta\partial\theta^{\prime}}$,
$\theta,\theta^{\prime}\in\\{\alpha,\beta,\delta\\}$, are given in Item (33)
and

$\displaystyle\mathbb{E}\big{[}X^{\alpha}\log(X)\big{]}=\dfrac{\beta^{\alpha}\big{[}2-2\gamma+2\alpha\log(\beta)\big{]}}{\displaystyle{\alpha\over
2}\bigg{[}2+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)}\bigg{]}-\alpha}$
(37)
$\displaystyle+\dfrac{\displaystyle\beta^{\alpha+1}\delta\bigg{\\{}-2\Gamma\bigg{(}2+{1\over\alpha}\bigg{)}\bigg{[}\alpha\log(\beta)+\psi^{(0)}\bigg{(}2+{1\over\alpha}\bigg{)}\bigg{]}+\beta\delta\Gamma\bigg{(}2+{2\over\alpha}\bigg{)}\bigg{[}\alpha\log(\beta)+\psi^{(0)}\bigg{(}2+{1\over\alpha}\bigg{)}\bigg{]}\bigg{\\}}}{\displaystyle{\alpha\over
2}\bigg{[}2+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)}\bigg{]}-\alpha},$

$\displaystyle\mathbb{E}\big{[}X^{\alpha}\log^{2}(X)\big{]}=\dfrac{\displaystyle\beta^{\alpha}\big{[}-12\gamma+6\gamma^{2}+\pi^{2}+12\alpha\log(\beta)-12\gamma\alpha\log(\beta)+6\alpha^{2}\log(\beta)^{2}\big{]}}{\displaystyle
3\alpha^{2}\bigg{[}2+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)}\bigg{]}}$
(38) $\displaystyle-\dfrac{\displaystyle
6\beta^{\alpha+1}\delta\Gamma\bigg{(}2+{1\over\alpha}\bigg{)}\bigg{[}\alpha^{2}\log^{2}(\beta)+2\alpha\log(\beta)\psi^{(0)}\bigg{(}2+{1\over\alpha}\bigg{)}+\psi^{(0)}\bigg{(}2+{1\over\alpha}\bigg{)}^{2}+\psi^{(1)}\bigg{(}2+{1\over\alpha}\bigg{)}\bigg{]}}{\displaystyle
3\alpha^{2}\bigg{[}2+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)}\bigg{]}}$
$\displaystyle+\dfrac{\displaystyle
3\beta^{\alpha+2}\delta^{2}\Gamma\bigg{(}2+{2\over\alpha}\bigg{)}\bigg{[}\alpha^{2}\log(\beta)^{2}+2\alpha\log(\beta)\psi^{(0)}\bigg{(}2+{2\over\alpha}\bigg{)}+\psi^{(0)}\bigg{(}2+{2\over\alpha}\bigg{)}^{2}+\psi^{(1)}\bigg{(}2+{2\over\alpha}\bigg{)}\bigg{]}}{\displaystyle
3\alpha^{2}\bigg{[}2+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)}\bigg{]}},$

where $\gamma=-{{\rm d}\Gamma(x)\over{\rm d}x}\big{|}_{x=1}\approx 0.57721$ is
the Euler-Mascheroni constant, and

$\displaystyle\mathbb{E}\big{(}X^{\alpha}\big{)}=\beta^{\alpha}\,\dfrac{2+\beta\delta\bigg{[}\beta\delta\Gamma\bigg{(}2+\dfrac{2}{\alpha}\bigg{)}-2\Gamma\bigg{(}2+\dfrac{1}{\alpha}\bigg{)}\bigg{]}}{\displaystyle
2+\delta^{2}\beta^{2}\Gamma\bigg{(}1+{2\over\alpha}\bigg{)}-2\delta\beta\Gamma\bigg{(}1+{1\over\alpha}\bigg{)}}.$
(39)

The proofs of equations (37)-(39) are followed by the integral kernel given by
item (6) Cankaya2018  and taking derivative of both sides of integral or the
software Mathematica 12.0$\copyright$ can be used.

The definition of Fisher information based on $\log_{q}$ is given by, see
Plastinoetal97 ; CanKor18 ,

$\displaystyle_{q}F(\boldsymbol{\theta})\Big{|}_{\boldsymbol{\theta}\coloneqq\widehat{\boldsymbol{\theta}}}=_{q}\mathbb{E}\left[{\partial
l\over\partial\boldsymbol{\theta}}\left({\partial
l\over\partial\boldsymbol{\theta}}\right)^{\intercal}f^{1-q}(X;\boldsymbol{\theta})\right]=\int_{0}^{\infty}{\partial
l\over\partial\boldsymbol{\theta}}\left({\partial\l\over\partial\boldsymbol{\theta}}\right)^{\intercal}f^{2-q}(x;\boldsymbol{\theta})\,\mathrm{d}x,$
(40)

where $l=\sum_{i=1}^{n}\log[f(x_{i};\boldsymbol{\theta})]$ and
$f(x;\boldsymbol{\theta})$ is a parametric model. The calculation of integral
is performed for one variable case $x_{i}$ from observations
$x_{1},x_{2},\dots,x_{n}$. Since it is replicated for $n$ case, FI matrix can
be rewritten as the form in equation (36). If $q=1$ in equation (40), then
${}_{q}F$ drops to $F$ in equation (36). The connection between Fisher
information and Tsallis $q$-entropy is proved by CanKor18 . Since the
analytical tractability of integral calculation cannot be easy to follow, the
numerical integration technique in Mathematica 12.0$\copyright$ is used to get
the values of elements in matrix in equation (40).

##  VI Application on real data sets

The real data sets are modelled by using BWeibull and BGamma distributions
which have bimodality property. The parameter $\delta$ that gives an advantage
to model the bimodal data which can occur due to contamination or irregularity
into underlying or main distribution as an indicator for abnormalization in
the working principle of a phenomena in the universe is used in Weibull and
Gamma distributions. As it is clearly observed, the parameter $\delta$ with
$1+(1-\delta x)^{2}$ produces a modality due to fact that the function is
polynomial with order $2$. We confine ourself the real data sets which can be
modelled by function having the bimodality property.

MLE and MLqE methods are used to estimate the model parameters. If a data set
has two peaks, then it is expected that BWeibull and BGamma distributions
capture the peaks and so they are flexible when they are compared with Weibull
and Gamma distributions. Since BWeibull and BGamma distributions have
analytical expression for their CDF, it is advantegous to use them for getting
statistics from Kolmogorov-Smirnov (KS) and Cramér–von Mises (CVM). Since we
use $\log$ and $\log_{q}$ in MLE and MLqE methods for getting the estimates of
parameters, using information criteria (IC) such as Akaike and Bayesian, etc
is not appropriate in view of the fact that $\log$ and $\log_{q}$ are not
comparable. As it is well-known Hamparsum87 , the lack of fit part of IC
depends on $\log$. There should be an equilibrium among the values of IC
applied to different probability density functions. Since there does not exist
have the equilibrium between $\log$ and $\log_{q}$ to make a comparison for
lack of fit part of IC, it is not appropriate to use IC and we consult
goodness of fit tests which are KS and CVM in order to test the modelling
competence of the used objective functions produced by $\log(f)$ and
$\log_{q}(f)$ God60 ; Hub64  as objective functions from MLE and MLqE,
respectively. Note that MLE and MLqE are methods used to estimate the
parameters of $f(x;\boldsymbol{\theta})$. However, $\log(f)$ are alternative
step for calculation in order to reach the estimators
$\widehat{\boldsymbol{\theta}}$ from MLE. This alternative situation is a way
to generalize MLE method as M-estimation method God60 ; Hub64 . $q$ in
$\log_{q}$ is a parameter to derive the different forms of
$f(x;\boldsymbol{\theta})$ Bercher12a . For this reason, $\log_{q}$ likelihood
and its special form ($\log$ likelihood) are applied to estimate the
parameters $\boldsymbol{\theta}=(\alpha,\beta,\delta)$.

The functions in equations (22)-(34) are nonlinear and the estimating
equations for the parameters $\alpha,\beta$ and $\delta$ are also nonlinear to
get the estimators $\widehat{\alpha},\widehat{\beta}$ and $\widehat{\delta}$
analytically. The stochasticity and metaheuristic algorithms are tools used to
optimize the nonlinear functions according to the parameters to get the
estimates of $\widehat{\alpha},\widehat{\beta}$ and $\widehat{\delta}$. The
optimization of functions in equations (22)-(34) is performed by a
’metaheuristicOpt’ package with ’metaOpt’ function in free statistical
software $R$ version $4.0.2$. The chosen algorithm in ’metaOpt’ is Harmony
Search (HS) harmony  due to fact that we have the highest p-values of KS and
CVM, which can show that BWeibull and BGamma with $\log$ and $\log_{q}$ should
be optimized by ’metaOpt’ with HS. After applying the ’metaOpt’ with HS to
functions in equations (22) and (34), we have the estimated values of
$\widehat{\boldsymbol{\theta}}$ from $\log$ and $\log_{q}$ likelihood
estimaton methods. The parameter $q$ in $\log_{q}$ likelihood estimation
method is chosen according the biggest values of probability (p-value) for KS
and CVM test statistics. Note that the objective function $\log_{q}(f)$ takes
different forms for each value of $q$. Thus, for only one parametric model
$f(x;\boldsymbol{\theta})$, we can have neighborhoods of
$f(x;\boldsymbol{\theta})$ for different value of $q$, which helps us to
manage the efficiency and robustness while performing the modelling on data.

BWeibull and BGamma distributions with their $\log$, $\log_{q}$ forms, i.e.,
$\log(f)$ and $\log_{q}(f)$ as objective functions are applied at $5$ real
data sets in order to test the modelling competence on the data sets. The c.d.
and p.d. functions abbreviated as CDF and PDF are superimposed on the
empirical c.d. function and histogram of data to depict the best modelling
chosen according to the p-values of KS and CVM (see Figures 5-9). In Figures
5(a)-9(a), CDF is more accurate illustration to depict the modelling of
function with estimates and also the illustrations of CDF and PDF from Figures
5-9) show that the modelling on data sets is accomplished well. Using the
probability values of KS and CVM is necessary to provide the homogenity for
comparison with their corresponding test statistics of KS and CVM.

Example 1 : The data are the breaking stress of carbon fibers. First 10
subgroups in control which consists of $50$ observations are modelled by using
Weibull distribution carbonfibers . Table 1 represents the values of estimates
for the carbon fibers data modelled by BWeibull and BGamma with $\log$ and
$\log_{q}$ as objective functions. The role of MLqE is observed. That is,
there is an important increasing on the modelling competence when
$\log_{q}(f)$ from MLqE is used. Note that it is observed that there is a good
accommodation between this data as an empricial distribution and $\log_{q}(f)$
from MLqE, which can lead to have the highest p-value of KS. The parameter $q$
as tuning constant produces the neighborhoods of a parametric model, i.e.,
objective function, in order to increase the modelling competence of a
parametric model Bercher12a .

Table 1: Inference values for parameters, KS and CVM values for carbon fibers
data

Model | $\hat{\alpha}$ | $\hat{\beta}$ | $\hat{\delta}$ | KS | p-value(KS) | CVM | p-value(CVM)
---|---|---|---|---|---|---|---
BWeibull($\hat{{\theta}}_{\text{MLE}}$) | 3.6961(0.0807) | 2.7482(0.0306) | 2.3073(1.2630) | 0.14 | 0.7112 | 0.0690 | 0.1556
BWeibull($\hat{{\theta}}_{\text{MLqE}}$), $q=0.8$ | 4.8707(0.1450) | 2.8434(0.0378) | 1.8986(2.4058) | 0.06 |  $\approx$ 1 | 0.0158 | 0.1641
BGamma($\hat{{\theta}}_{\text{MLE}}$) | 14.9970(0.6521) | 5.8832(0.2449) | 1.5823(1.0677) | 0.12 | 0.8643 | 0.0734 | 0.1549
BGamma($\hat{{\theta}}_{\text{MLqE}}$), $q=0.75$ | 14.9994(0.7332) | 5.9217(0.3195) | 1.7915(3.0210) | 0.10 | 0.9639 | 0.0706 | 0.1553

KS: Kolmogorov-Smirnov, CVM: Cramér–von Mises, Bold represents the best
fitting.

(a) (Empirical) CDF of BWeibull($\widehat{{\theta}}_{\text{MLqE}}$)

(b) Frequency and PDF of BWeibull($\widehat{{\theta}}_{\text{MLqE}}$)

Figure 5: The fitted values from MLqE of parameters of bimodal Weibull
distribution for carbon fibers data

Example 2 : The data are maximum ozone concentrations labeled as ”o3max” in
”goft” package in free statistical software $R$ version $4.0.2$. 24
observations are included by o3max. Table 2 shows MLE and MLqE with BWeibull
are better than that of BGamma for the p-values of KS and CVM. Note that
BWeibull with $\log$ and $\log_{q}$ has an important fitting competence when
compared with that of BGamma.

Table 2: Inference values for parameters, KS and CVM values for o3max data

Model | $\hat{\alpha}$ | $\hat{\beta}$ | $\hat{\delta}$ | KS | p-value(KS) | CVM | p-value(CVM)
---|---|---|---|---|---|---|---
BWeibull($\hat{{\theta}}_{\text{MLE}}$) | 10.0851(0.7039) | 0.1728(0.0018) | 14.9983(10.9941) | 0.1667 | 0.8928 | 0.0556 | 0.1577
BWeibull($\hat{{\theta}}_{\text{MLqE}}$), $q=0.99$ | 10.1609(0.7079) | 0.1728(0.0018) | 14.9995(11.1710) | 0.1667 | 0.8928 | 0.0556 | 0.1577
BGamma($\hat{{\theta}}_{\text{MLE}}$) | 3.2524(0.2537) | 14.9963(1.0863) | 1.6134(1.0448) | 0.4583 | 0.0129 | 0.6667 | 0.0856
BGamma($\hat{{\theta}}_{\text{MLqE}}$), $q=0.99$ | 3.2513(0.2542) | 14.9935(1.0867) | 1.6130(1.0585) | 0.4583 | 0.0129 | 0.6667 | 0.0856

(a) (Empirical) CDF of BWeibull($\widehat{{\theta}}_{\text{MLqE}}$)

(b) Frequency and PDF of BWeibull($\widehat{{\theta}}_{\text{MLqE}}$)

Figure 6: The fitted values from MLqE of parameters of bimodal Weibull
distribution for o3max data

Example 3 : Growth hormone data was also modelled by Alizadeh . The number of
observation is 35. MLqEs with BWeibull and BGamma distributions outperform
when compared with parametric model in Alizadeh .

Table 3:  Inference values for parameters, KS and CVM values for growth
hormone data

Model | $\hat{\alpha}$ | $\hat{\beta}$ | $\hat{\delta}$ | KS | p-value(KS) | CVM | p-value(CVM)
---|---|---|---|---|---|---|---
BWeibull($\hat{{\theta}}_{\text{MLE}}$) | 1.1344(0.0370) | 2.1303(0.1213) | 3.8210(1.0365) | 0.1429 | 0.8674 | 0.0651 | 0.1562
BWeibull($\hat{{\theta}}_{\text{MLqE}}$), $q=0.85$ | 1.2489(0.0548) | 2.2989(0.1541) | 3.5132(1.4327) | 0.0857 | 0.9995 | 0.0357 | 0.1608
BGamma($\hat{{\theta}}_{\text{MLE}}$) | 1.9136(0.1685) | 0.7620(0.0319) | 3.6132(1.7301) | 0.1429 | 0.8674 | 0.0765 | 0.1544
BGamma($\hat{{\theta}}_{\text{MLqE}}$), $q=0.87$ | 2.4648(0.2891) | 0.9224(0.0513) | 3.2896(2.9941) | 0.0857 | 0.9995 | 0.0349 | 0.1610

Italic represents comparison for the best fitting for MLE of BWeibull and
BGamma.

Table 3 represents the values of estimates for growth hormone data modelled by
objective functions. Italic represents comparison between CVM values when MLE
method is used. Note that MLE method reflects directly parametric models for
which BWeibull and BGamma are used. In other words, deformation on the
parametric model does not exist when we compare with $\log_{q}$. When MLEs of
parameters of BWeibull and BGamma are compared, we observe that MLE of
parameters of BWeibull is a little ahead for the performance of the fitting
competence according to the values of CVM and its corresponding p-values. KS
statistic and the corresponding p-values are not capable to assess the
modelling competence. According to the values of p-value(CVM), MLqE with
BGamma give the best performance on the modelling data when compared with MLqE
with BWeibull.

(a) (Empirical) CDF of BWeibull($\widehat{{\theta}}_{\text{MLE}}$)

(b) Frequency and PDF of BWeibull($\widehat{{\theta}}_{\text{MLE}}$)

Figure 7: The fitted values from MLE of parameters of bimodal Weibull
distribution for growth hormone data

Example 4 : Wheaton River data are analyzed by Vila2020  and references
therein. The number of observation is 72. MLqE with BWeibull give
outperforming when we compare by MLqE with BGamma. BGamma from Vila2020  has a
good performance for fitting among the distributions Ortega11  and references
therein. The advantage of using MLqE is observed when we look at the p-values
of KS and CVM in Table 4.

Table 4: Inference values for parameters, KS and CVM values for Wheaton River
data

Model | $\hat{\alpha}$ | $\hat{\beta}$ | $\hat{\delta}$ | KS | p-value(KS) | CVM | p-value(CVM)
---|---|---|---|---|---|---|---
BWeibull($\hat{{\theta}}_{\text{MLE}}$) | 0.9721(0.0102) | 5.4258(0.1397) | 0.1924(0.0050) | 0.0694 | 0.9951 | 0.0260 | 0.1624
BWeibull($\hat{{\theta}}_{\text{MLqE}}$), $q=0.99$ | 0.9770(0.0103) | 5.4536(0.1411) | 0.1910(0.0050) | 0.0694 | 0.9951 | 0.0206 | 0.1633
BGamma($\hat{{\theta}}_{\text{MLE}}$) | 1.0591(0.0182) | 0.1768(0.0023) | 0.1776(0.0039) | 0.0833 | 0.9639 | 0.0336 | 0.1612
BGamma($\hat{{\theta}}_{\text{MLqE}}$), $q=0.99$ | 1.0568(0.0184) | 0.1780(0.0024) | 0.1783(0.0040) | 0.0833 | 0.9639 | 0.0287 | 0.1620

(a) (Empirical) CDF of BWeibull($\widehat{{\theta}}_{\text{MLqE}}$)

(b) Frequency and PDF of BWeibull($\widehat{{\theta}}_{\text{MLqE}}$)

Figure 8: The fitted values from MLqE of parameters of bimodal Weibull
distribution for Wheaton River data

Example 5: The death occurring due to cancer increases. Gastric cancer data
consisting of 125 observations are the censored data because of the difficulty
of observering the patients who can die in the progress of medical treatment
of patients. If data are censored, then there exists a replication for some
observations, which leads to bimodaliy for frequency. For this reason, we
prefer to use this data set which will be modelled by BWeibull and BGamma with
$\log$ and $\log_{q}$ as objective functions. Let us note that censoring
designs censoring ; Ngcensoring  can be originally regarded as an objective
function Hampeletal86 ; MLqEBio ; PhDthesis ; CanKor18 ; Canorder20  used to
fit a data set. Thus, using $\log(f)$ and $\log_{q}(f)$ is reasonable and $f$
can be chosen as BWeibull and BGamma. Table 5 shows that MLqE with BWeibull
should be considered as good fitting performance assessed by p-value of KS
even if p-value of CVM of BGamma with $\log$ as a little ahead from others is
the best one among them.

Table 5: Inference values for parameters, KS and CVM values for gastric cancer
data

Model | $\hat{\alpha}$ | $\hat{\beta}$ | $\hat{\delta}$ | KS | p-value(KS) | CVM | p-value(CVM)
---|---|---|---|---|---|---|---
BWeibull($\hat{{\theta}}_{\text{MLE}}$) | 1.1181(0.0092) | 8.4517(0.1373) | 0.1764(0.0031) | 0.104 | 0.5085 | 0.0962 | 0.1514
BWeibull($\hat{{\theta}}_{\text{MLqE}}$), $q=0.9$ | 1.0092(0.0097) | 6.7659(0.1403) | 0.2106(0.0043) | 0.096 | 0.6121 | 0.1390 | 0.1450
BGamma($\hat{{\theta}}_{\text{MLE}}$) | 0.9531(0.0133) | 0.1433(0.0030) | 0.2193(0.0031) | 0.112 | 0.4131 | 0.0849 | 0.1531
BGamma($\hat{{\theta}}_{\text{MLqE}}$), $q=0.9$ | 0.9439(0.0147) | 0.1481(0.0012) | 0.2225(0.0036) | 0.104 | 0.5085 | 0.1422 | 0.1446

(a) (Empirical) CDF of BWeibull($\widehat{{\theta}}_{\text{MLqE}}$)

(b) Frequency and PDF of BWeibull($\widehat{{\theta}}_{\text{MLqE}}$)

Figure 9: The fitted values from MLqE of parameters of bimodal Weibull
distribution for gastric cancer data

## Conclusions

A bimodal form of Weibull distribution has been derived. The properties such
as unimodality, bimodality, moment, moment generation functions, entropies,
etc. have been obtained. In order to get the estimators
$\widehat{\boldsymbol{\theta}}$ of parameters, MLE and MLqE methods have been
used. The concavity property of $\log$ and $\log_{q}$ are very important to
apply for estimation of parameters Lindsay94 . Further, it should be noted
that the existences of Shannon and Tsallis $q$-entropies can also be important
to apply MLE and MLqE methods. The estimates of
$\widehat{\boldsymbol{\theta}}$ from MLE and MLqE are obtained by using of the
heuristic algorithm because of nonlinearity of $\log(f)$ and $\log_{q}(f)$.
The p.d. function $f$ has been chosen as BWeibull and BGamma distributions.
Instead of using a parametric model directly for modelling a data set via
$\log$, it should be preferred to apply $\log_{q}$ which can derive the
different forms of a parametric model for the different values of $q$ from
MLqE method. Thus, we have obtained the results from numerical experiments
showing that using MLqE method for parameters of BWeibull distribution is
suggested. The anaytical expressions for the elements of FI based on $\log$
have been obtained. The square root of variance values of estimators
$\widehat{\boldsymbol{\theta}}$ from MLE and MLqE are given in the numerical
experiments.

The application of the proposed distribution for case of the censored data
will be studied and comparative study for different forms of the censoring
designs censoring ; Ngcensoring  will be applied to test performance of the
censoring designs when we have the empirical distribution which has
bimodality. The properties such as existence, uniqueness of roots for $\log$
and $\log_{q}$ likelihood functions of the censoring designs which can be
regarded as objective function Hampeletal86  will be studied.

## Acknowledgments

This study was financed in part by the Coordenação de Aperfeiçoamento de
Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001.

## Appendix: Main buildings of codes used for computation of parameters and
goodness of fit test statistics

The PDF is computed by

    biwblpdf<- function(x,a,b,d){
    Ψf=dweibull(x, a, b, log = FALSE);
    ΨZ = 2 + ((2 *b *d *(-gamma(1/a) + b *d* gamma(2/a)))/a);
    Ψg=(1 + (1 - d*x)^2)*f/Z;
    }

The optimization of $l_{q}$ in equation (34) is performed by following
functions:

    ΨlogqLikFun <- function(p) {
    ΨΨa <- p[1]; b <- p[2]; d <- p[3];
    ΨΨ-sum((biwblpdf(x,a,b,d)^(1-q)-1)/(1-q))
    Ψ}
    Ψmlqe <- metaOpt(logqLikFun, optimType = "MIN", algorithm = "HS",
    ΨnumVar,rangeVar,control = list(), seed = NULL)

The CDF arranged by manipulation and adopted for the computation of R platform
is computed by

    ΨCDFbiWei<- function(x,a,b,d){
    ΨΨx=sort(x);
    ΨΨG=(-2 * b * d * gamma(1/a) + a * (2 - 2 * exp(-(x/b)^a) +
    ΨΨb^2 * d^2 * gamma(2/a+1) +   2 * b * d * (gamma(1 + 1/a)
    ΨΨ* pgamma((x/b)^a, 1 + 1/a, 1, lower = FALSE))
    ΨΨ- b^2 * d^2 * (gamma(2/a+1) * pgamma((x/b)^a, 2/a+1, 1, lower = FALSE))))
    ΨΨ/(a * (2 + (2 * b * d * (-gamma(1/a) + b * d * gamma(2/a)))/a))
    Ψ}

The values of goodness of fit test statistics are obtained by r1 ; r2

    Ψinstall.packages("CDFt")
    Ψlibrary("CDFt")
    Ψval_CDFbiWeiq=CDFbiWei(sort(x),mleBiWeiabq[1],mleBiWeiabq[2],mleBiWeiabq[3])
    Ψfun.ecdf=ecdf(x);my.ecdf <- fun.ecdf(sort(x));
    Ψks.test(my.ecdf, val_CDFbiWeiq,alternative = c("two.sided", "less", "greater"),
    Ψexact = NULL, tol=1e-8,simulate.p.value=FALSE,B=2000)
    Ψresq = CramerVonMisesTwoSamples(my.ecdf,val_CDFbiWeiq);
    ΨpvalueCVMq = 1/6*exp(-resq);

## References

## References

  * (1) Ahsanullah, M., Hamedani, G. G. (2010). Exponential distribution: theory and methods. Nova Science Publ.
  * (2) Castellares, F., Lemonte, A. J. (2015). A new generalized Weibull distribution generated by gamma random variables. Journal of the Egyptian Mathematical Society, 23(2), 382-390.
  * (3) Xie, M., Tang, Y., Goh, T. N. (2002). A modified Weibull extension with bathtub-shaped failure rate function. Reliability Engineering & System Safety, 76(3), 279-285.
  * (4) Martinez, E. Z., Achcar, J. A., Jácome, A. A., Santos, J. S. (2013). Mixture and non-mixture cure fraction models based on the generalized modified Weibull distribution with an application to gastric cancer data. Computer methods and programs in biomedicine, 112(3), 343-355.
  * (5) Yang, G. (2007). Life cycle reliability engineering. John Wiley & Sons.
  * (6) Silva, G. O., Ortega, E. M., Cordeiro, G. M. (2010). The beta modified Weibull distribution. Lifetime data analysis, 16(3), 409-430.
  * (7) Rinne, H. (2008). The Weibull distribution: a handbook. CRC press.
  * (8) Elal-Olivero, D. (2010). Alpha-skew-normal distribution. Proyecciones (Antofagasta), 29(3), 224-240.
  * (9) Çankaya, M. N., Yalçınkaya, A., Altındaǧ, Ö., Arslan, O. (2019). On the robustness of an epsilon skew extension for Burr III distribution on the real line. Computational Statistics, 34(3), 1247-1273.
  * (10) Çankaya M.N. Asymmetric bimodal exponential power distribution on the real line, Entropy, 2018, 20(1), 23.
  * (11) Knapp, A. W. (2005). Basic real analysis. Springer Science & Business Media.
  * (12) S. Klugman; H. Panjer; G. Willmot. Loss models: From data to decisions, Wiley, New York, 1998.
  * (13) Rao C. R. Quadratic entropy and analysis of diversity, Sankhya A, 2010, 72, 70-80.
  * (14) Shannon C.E. A mathematical theory of communication, Bell System Technical Journal, 1948, 27, 379-423, 623-656.
  * (15) Tsallis C. Possible generalization of Boltzmann-Gibbs statistics, Journal of Statistical Physics, 1988, 52, 479-487.
  * (16) Vila R.; Ferreira L.; Saulo H.; Prataviera F.; Ortega E. A bimodal gamma distribution: properties, regression model and applications, Statistics, 2020, 54(3), 469-493.
  * (17) J. Xue. Loop Tiling for Parallelism, The Springer International Series in Engineering and Computer Science, Springer US, 2012.
  * (18) Lindsay, B. G. (1994). Efficiency versus robustness: the case for minimum Hellinger distance and related methods. The annals of statistics, 22(2), 1081-1114.
  * (19) C. Tsallis, Introduction to Nonextensive Statistical Mechanics: Approaching a Complex World, Springer, New York, 2009.
  * (20) Ferguson, T. S. (1996). Strong Consistency of Maximum-Likelihood Estimates. In A Course in Large Sample Theory (pp. 112-118). Springer US.
  * (21) Ferguson, T. S. (1996). Asymptotic Normality of the Maximum-Likelihood Estimate. In A Course in Large Sample Theory (pp. 119-125). Springer US.
  * (22) Plastino, A., Plastino, A. R., Miller, H. G. (1997). Tsallis nonextensive thermostatistics and Fisher’s information measure. Physica A: Statistical Mechanics and its Applications, 235(3-4), 577-588.
  * (23) Çankaya, M. N., Korbel, J. (2018). Least informative distributions in maximum q-log-likelihood estimation. Physica A: Statistical Mechanics and its Applications, 509, 140-150.
  * (24) Geem, Z. W. (Ed.). (2009). Music-inspired harmony search algorithm: theory and applications (Vol. 191). Springer.
  * (25) Sultan, K. S., Ismail, M. A., Al-Moisheer, A. S. (2007). Mixture of two inverse Weibull distributions: Properties and estimation. Computational Statistics Data Analysis, 51(11), 5377-5387.
  * (26) Sultan, K. S., Al-Moisheer, A. S. (2013). Estimation of a discriminant function from a mixture of two inverse Weibull distributions. Journal of Statistical Computation and Simulation, 83(3), 405-416.
  * (27) Hasegawa, Y., Arita, M. (2009). Properties of the maximum q-likelihood estimator for independent random variables. Physica A: Statistical Mechanics and its Applications, 388(17), 3399-3412.
  * (28) Nichols, M. D., Padgett, W. J. (2006). A bootstrap control chart for Weibull percentiles. Quality and reliability engineering international, 22(2), 141-151.
  * (29) Bozdogan, H. (1987). Model selection and Akaike’s information criterion (AIC): The general theory and its analytical extensions. Psychometrika, 52(3), 345-370.
  * (30) Godambe, V. P. (1960). An optimum property of regular maximum likelihood estimation. The Annals of Mathematical Statistics, 31(4), 1208-1211.
  * (31) Huber, P. J. (1964). Robust estimation of a location parameter: Annals Mathematics Statistics, 35.
  * (32) Fisher, R. A. (1925, July). Theory of statistical estimation. In Mathematical Proceedings of the Cambridge Philosophical Society (Vol. 22, No. 5, pp. 700-725). Cambridge University Press.
  * (33) E.L. Lehmann, G. Casella, Theory of point estimation, Wadsworth & Brooks/Cole. Pacific Grove, CA, 589, USA, 1998.
  * (34) L. Pardo, Statistical inference based on divergence measures, CRC Press, Taylor & Francis Group, 2005.
  * (35) Korbel, J., Hanel, R., Thurner, S. (2019). Information geometric duality of $\phi$-deformed exponential families. Entropy, 21(2), 112.
  * (36) Korbel, J., Hanel, R., Thurner, S. (2020). Information geometry of scaling expansions of non-exponentially growing configuration spaces. The European Physical Journal Special Topics, 229(5), 787-807.
  * (37) Li, Y. H., Yeh, P. C. (2012). An interpretation of the Moore-Penrose generalized inverse of a singular Fisher Information Matrix. IEEE Transactions on Signal Processing, 60(10), 5532-5536.
  * (38) Bercher, J. F. (2012). A simple probabilistic construction yielding generalized entropies and divergences, escort distributions and q-Gaussians. Physica A: Statistical Mechanics and its Applications, 391(19), 4460-4469.
  * (39) Alizadeh, M., Bagheri, S. F., Samani, E. B., Ghobadi, S., Nadarajah, S. (2018). Exponentiated power Lindley power series class of distributions: Theory and applications. Communications in Statistics-Simulation and Computation, 47(9), 2499-2531.
  * (40) De Pascoa, M. A., Ortega, E. M., Cordeiro, G. M. (2011). The Kumaraswamy generalized gamma distribution with application in survival analysis. Statistical methodology, 8(5), 411-433.
  * (41) F.R. Hampel, E. M. Ronchetti, P. J. Rousseeuw, W. A. Stahel, Robust statistics: The approach based on influence functions, Wiley Series in Probability and Statistics, New York, 1986.
  * (42) Çankaya, M. N. 2015. M-Estimators with asymmetric influence function: Properties and their applications. PhD diss., University of Ankara.
  * (43) Çankaya, M. N. (2020). M-Estimations of Shape and Scale Parameters by Order Statistics in Least Informative Distributions on q-deformed logarithm. Iǧdır Üniversitesi Fen Bilimleri Enstitüsü Dergisi, 10(3), 1984-1996.
  * (44) McCool, J. I. (2012). Using the Weibull distribution: reliability, modeling, and inference (Vol. 950). John Wiley & Sons.
  * (45) Lio, Y., Ng, H. K. T., Tsai, T. R., Chen, D. G. (Eds.). (2019). Statistical Quality Technologies: Theory and Practice. Springer.
  * (46) https://www.r-bloggers.com/2015/01/goodness-of-fit-test-in-r/, 09/17/2020
  * (47) https://www.rdocumentation.org/packages/CDFt/versions/1.0.1/topics/CramerVonMisesTwoSamples, 09/17/2020