# A General Class of Trimodal Distributions: Properties and Inference

Roberto Vila<EMAIL_ADDRESS>Department of Statistics, University of
Brasília, Brasília, Brazil  Victor Serra<EMAIL_ADDRESS>Department of
Statistics, University of Brasília, Brasília, Brazil  Mehmet N. Çankaya
<EMAIL_ADDRESS>Department of International Trading and Finance,
Faculty of Applied Sciences, Uşak University, Uşak, Turkey;  Department of
Statistics, Faculty of Art and Sciences, Uşak University, Uşak, Turkey  Felipe
Quintino<EMAIL_ADDRESS>Department of Mathematics and Statistics,
Federal University of Rondônia, Paraná, Brazil

###### Abstract

The modality is important topic for modelling. Using parametric models is an
efficient way when real data set shows trimodality. In this paper we propose a
new class of trimodal probability distributions, that is, probability
distributions that have up to three modes. Trimodality itself is achieved by
applying a proper transformation to density function of certain continuous
probability distributions. At first, we obtain preliminary results for an
arbitratry density function $g(x)$ and, next, we focus on the Gaussian case,
studying trimodal Gaussian model more deeply. The Gaussian distribution is
applied to produce the trimodal form of Gaussian known as normal distribution.
The tractability of analytical expression of normal distribution, and
properties of the trimodal normal distribution are important reasons why we
choose normal distribution. Furthermore, the existing distributions should be
improved to be capable of modelling efficiently when there exists a trimodal
form in a data set. After new density function is proposed, estimating its
parameters is important. Since Mathematica 12.0 software has optimization
tools and important modelling techniques, computational steps are performed by
using this software. The bootstrapped form of real data sets are applied to
show the modelling ability of the proposed distribution when real data sets
show trimodality.

Keywords. Class of distributions $\cdot$ unimodality $\cdot$ bimodality
$\cdot$ trimodality $\cdot$ inference.
Mathematics Subject Classification (2010). MSC 60E05 $\cdot$ MSC 62Exx $\cdot$
MSC 62Fxx.

##  1 Introduction

The modality is an important topic when the nature of phenomena can be
modelled by using the function which can be capable to have different forms of
peaks. The modality can occur when there exists an irregularity in the output
of an experiment. In the statistical view point, the random variables are
nonidentically distributed. In other words, there can be a mixing of some
populations even if same experiment is conducted while getting outputs of the
corresponding experiment [1, 4, 5, 6]. In such a situation, location and scale
parameters of the mixed populations are important to get the central tendency
and disperson (statistics) of the mixed form of the population in the
experiment after shape, scale and bimodality parameters of function in the
statistical theory are necessary components of a function which is used for
conducting an efficient modelling [2, 3, 4, 8, 9, 12, 17].

There are many well-known distributions used at the statistical inference in
which regression and its counter parts such as time series, design of
experiments, structural equation modelling in the applications from social
science. There are different techniques to produce a probability density
function [7]. The analytical tractability and properties of the proposed
distribution are important when the new distribution is used for modelling.
For example, the existence of moments and entropy function are important when
the modelling on the real data sets are performed [24, 34]. The well-known
normal distribution which is also known as Gaussian distribution should be
transformed into a trimodal form when data sets show trimodality. The
advantage of using a trimodal distribution is that a real data set can be a
combination of two, three or more normal distributions with different
parameter values for location and scale. The apperance of such mixed
distributions can be in a trimodal form. Especially, since the working
principle of a phenomena depends on many factors, it is reasonable to expect
that a trimodal form for a real data set can occur. In other words, it is
assumed that the random variables are identically distributed. However, the
identicality is a restrictive assumptation for modelling a data set. The
parametric models are necessary to perform an efficient modelling when
trimodal representation in a data set exists due to the structure of non-
identicality (hetero data) or the mixing form of distributions. On the other
side, if observations $x_{1},x_{2},\cdots,x_{n}$ are distributed as a
parametric model such as trimodal form, then they are identical, that is, it
depends on where and how you look at the results of experiments, because we do
already have finite sample size. While managing an efficient modelling on real
data sets, it is reasonable to consult parametric models which can be capable
of dealing with different forms of modality. The bimodality parameters $\rho$
and $\delta$ with shape parameter $\alpha$ of Maxwell distribution in [18] are
tools for us to generate different forms of modality. There are different
degrees of mixing via Maxwell distribution and the expression which can help
us to have different forms of modality of function for generating modality via
compounding distributions [9, 12]. The bimodal form on the positive part of
the real line produces a trimodal form when it is reflected to the negative
part of the real line via mirror imaging [1, 13, 14, 15]. The different degree
of bimodality, i.e. the length of periodicity of modality, on the positive
part of the real line can be constructed by producing a new objective function
based on the deformation [9, 36], because deformation which can be regarded as
a kind of rescaling can make the different length for periodicity [2, 10, 11].

The smooth kernel distribution will be used to fit data sets; because, the
strict and soft forms of trimodal normal distribution should be compared with
smooth kernel technique to perform a comparison among them. Three normal
distributions, represented by functions $g(x;\mu_{1},\sigma_{1})$,
$g(x;\mu_{2},\sigma_{2})$ and $g(x;\mu_{3},\sigma_{3})$, can be mixed to get
the trimodal form of normal distribution $g$. Thus, the modelling performance
of the trimodal distribution constructed for the normal distribution $g$ can
be tested for the mixed data sets. In fact, the mixed form of two or three
function with one mode and symmetric can show the symmetric form with trimodal
representation. In other words, the groups around location parameter can be
divided into two forms [1, 4, 5, 6].

The main aim in this paper is to propose a distribution with trimodal form on
the real line via using a technique, as is given by [9, 5, 6, 16]. We keep to
follow the symmetric case and our aim is to focus on the trimodal form on the
real line. When there exists a trimodal form in a data set, the location and
scale parameters should be estimated efficiently. In other words, each group
coming from groups $g(x;\mu_{1},\sigma_{1})$, $g(x;\mu_{2},\sigma_{2})$ and
$g(x;\mu_{3},\sigma_{3})$ has its values for location and scale parameters. In
our case, we try to estimate one location parameter and one scale parameter
when the mixed data sets for two or three groups are used. It is important to
note that the apperance of trimodality can occur via mixing the different
values of parameters of functions. Since the true model for the mixed three
normal distributions is chosen to estimate the location and scale parameters
precisely, the performance of modelling will be increased greatly when the
probability density functions having modality property are taken into account.
For example, the mixed distribution has parameters which are the mixing
proportion $w_{1},w_{2}$ and $w_{3}=1-w_{1}-w_{2}$ for three groups which are
necessary to estimate. For the mixed normal distributions, we have
$\mu_{1},\sigma_{1},\mu_{2},\sigma_{2},\mu_{3},\sigma_{3},w_{1}$ and $w_{2}$
(see Section 5). In totally, there are 8 parameters which have to be
estimated. However, in our case we have three main parameters and also $\mu$
and $\sigma$ of the distribution. Thus, there will be 5 parameters which will
be estimated. The optimization of the log-likelihood function according to
these 5 parameters of a function can be easy to reach the global point of the
$\log(f)$ when compared with a function including 8 parameters. In addition,
we have only one location and scale parameter which can be free from the
mixing proportion $w_{1},w_{2}$ and $w_{3}=1-w_{1}-w_{2}$. The numerical
computation while conducting the optimization of log-likelihood according to
parameters can include the less numerical errors. Note that the numerical
errors in the function with 8 parameters can be bigger than that of 5
parameters. Thus, the more precise evaluation can be achieved for the
numerically precise evaluation of estimating the parameters $\mu$ and
$\sigma$, which is why we prefer to consider proposing such a trimodality
while conducting an efficient fitting on the real data sets. On the other side
of the modelling perspective, the structure of grouping cannot be determined
only for the used mixing proportions $w_{1},w_{2}$ and $w_{3}=1-w_{1}-w_{2}$.
There can be irrational proportions such as 1/6, 1/9, etc. for the mixing in a
data set. The precise evaluation in computation for the true value of
irrational proportion cannot be performed accurately, which makes a
disadvantage for us when we use the estimated values of these kinds of the
proportions such as irrational ones.

This paper is organized as follows. In Section 2, we define the class of
trimodal probabilistic models. In Section 3, some structural properties of the
proposed model are examined. We provide a formal proof for the trimodality of
a class of symmetric kernel densities, present a stochastic representation,
and we provide the closed formulas for the moments, entropies and stochastic
representation. The existence of these expressions is important to use the
proposed distrbution for fitting the data sets. Section 4 is divided to
describe the proposed model when the normal (Gaussian) is applied (see Table
1). For this case, some properties such as modality, moments, Shannon entropy
among others also are discussed. In Sections 5, we introduce the mixing form.
Section 6, represents the $\log_{q}$ likelihood function used as a method for
parameter estimation. In Section 7, the real data sets are applied. Section 8
is for conclusion and future works. In Appendices as a supplementary material,
we provide proofs and codes of Mathematica software.

##  2 A class of continuous probability distributions

Let $g:D\subset\mathbbm{R}\to[0,\infty)$, $D={\rm supp}(g)\neq\emptyset$, be a
kernel density with corresponding cumulative distribution function (CDF)
denoted by $G$. The function $g$ can be associated (or not) with an additional
parameter $\xi$ (or vector $\bm{\xi}$). For a random variable $X$ we define
the following probability density function (PDF)

$\displaystyle f(x;\bm{\theta})={\sigma\over
Z_{\bm{\theta}}}\,\biggl{[}\rho+\delta
T\biggl{(}{x-\mu\over\sigma};\alpha,p\biggr{)}\biggr{]}\,g\biggl{(}{x-\mu\over\sigma}\biggr{)},\quad{x-\mu\over\sigma}\in
D,$ (2.1)

where $\bm{\theta}=(\mu,\sigma,\alpha,\rho,\delta)$ is a parameter vector such
that $\mu\in\mathbbm{R}$ is a location parameter, $\sigma>0$ is a scale
parameter, $\alpha>0$ is a shape parameter and $\rho\geqslant 0$ and
$\delta\geqslant 0$ are parameters that control different forms of modality of
distribution. Note that $\rho$ and $\delta$ cannot be zero simultaneously. The
function $Z_{\bm{\theta}}$ appearing in the definition of $f$ is a normalizing
factor and the function $T:D\subset\mathbbm{R}\to(0,1)$ is given by

$\displaystyle
T(x;\alpha;p)={\gamma(p,x^{2}/\alpha^{2})\over\Gamma(p)}\quad\text{with}\ p>0\
\text{known}.$ (2.2)

Here $\gamma(p,u)=\int_{0}^{u}w^{p-1}{\rm e}^{-w}{\rm d}w$ is the incomplete
gamma function and $\Gamma(p)$ is the gamma function. For $p=3/2$ and
$D=(0,\infty)$, the function $T(x;\alpha,p)$ on $D$ defines the Maxwell
distribution with scale parameter $\alpha$.

Hereafter, we will denote $X\sim{\rm TD}(\bm{\theta})$ for a random variable
$X$ that follows the trimodal distribution (2.1).

When $\delta=0$ and $\rho\neq 0$ fixed in (2.1), the original density $g$ is
recovered.

Since $0<T(x;\alpha,p)<1$ for almost all $x\in D$, we have
$0<Z_{\bm{\theta}}\leqslant(\rho+\delta)\sigma$. A simple calculation shows
that (see Corollary 3.6)

$\displaystyle Z_{\bm{\theta}}$
$\displaystyle=(\rho+\delta)\sigma+\delta\sigma\left\\{{\mathbbm{E}}\big{[}G(-\alpha\sqrt{Y})\big{]}-{\mathbbm{E}}\big{[}G(\alpha\sqrt{Y})\big{]}\right\\},$
(2.3)

where $Y\sim{\rm Gamma}(p,1)$ and $G$ is the corresponding CDF of $g$.
Furthermore, the CDF of $X\sim{\rm TD}(\bm{\theta})$, denoted by
$F(x;\bm{\theta})$, is written as

$\displaystyle F(x;\bm{\theta})$ $\displaystyle={\rho\sigma\over
Z_{\bm{\theta}}}\,G\biggl{(}{x-\mu\over\sigma}\biggr{)}+{\delta\sigma\over
Z_{\bm{\theta}}}\left\\{\mathbbm{E}\big{[}G(-\alpha\sqrt{Y})\big{]}+T\biggl{(}{x-\mu\over\sigma};\alpha,p\biggr{)}G\biggl{(}{x-\mu\over\sigma}\biggr{)}\right\\}$
$\displaystyle-{\delta\sigma\over
Z_{\bm{\theta}}}\left\\{\mathbbm{E}\Big{[}\mathds{1}_{\\{Y\leqslant({x-\mu\over\sigma\alpha})^{2}\\}}G(-\alpha\sqrt{Y})\Big{]}\mathds{1}_{\\{x<\mu\\}}+\mathbbm{E}\Big{[}\mathds{1}_{\\{Y\leqslant({x-\mu\over\sigma\alpha})^{2}\\}}G(\alpha\sqrt{Y})\Big{]}\mathds{1}_{\\{x\geqslant\mu\\}}\right\\},$
(2.4)

for each $x\in\mathbbm{R}$. For more details, see Corollary 3.6.

Taking $x=\mu$ in (2), we get
$F(\mu;\bm{\theta})=\sigma\\{{\rho}G(0)+{\delta}\mathbbm{E}[G(-\alpha\sqrt{Y})]\\}/Z_{\bm{\theta}}.$
Let $W$ be a random variable with corresponding CDF $G$. If the distribution
$G$ of $W$ is symmetric about zero, then
$Z_{\bm{\theta}}=\rho\sigma+2\delta\sigma{\mathbbm{E}}[G(-\alpha\sqrt{Y})]$,
$G(0)=1/2$, and then $F(\mu;\bm{\theta})=1/2$. So, in this case, $\mu$ is
location parameter of $X$. Moreover, by letting $x\to\infty$ in (2), a simple
observation shows that $F(x;\bm{\theta})$ tends to 1, showing that the
parametric function in (2.1) is in fact a PDF.

Some natural examples of kernel densities $g$ to be plugged into (2.1), with
$p$ given, where trimodality shape is observed, are presented in Table 1.

Table 1: Some kernel densities $(g)$ that generate multimodality in the model (2.1). Distribution | $g$ | $G$ | $\bm{\xi}$ | $D$
---|---|---|---|---
Trimodal Gumbel | ${\rm e}^{-x-{\rm e}^{-x}}$ | ${\rm e}^{-{\rm e}^{-x}}$ | $-$ | $\mathbbm{R}$
Trimodal Laplace | ${1\over 2}{\rm e}^{-|x|}$ | $1+{1\over 2}{\rm e}^{-x}[\mathds{1}_{\\{x\leqslant 0\\}}-\mathds{1}_{\\{x\geqslant 0\\}}]$ | $-$ | $\mathbbm{R}$
Trimodal Logistic | ${{\rm e}^{-x}\over(1+{\rm e}^{-x})^{2}}$ | ${1\over 1+{\rm e}^{-x}}$ | $-$ | $\mathbbm{R}$
Trimodal Cauchy | ${1\over\pi(1+x^{2})}$ | ${1\over\pi}{\rm arctan}(x)+{1\over 2}$ | $-$ | $\mathbbm{R}$
Trimodal Student-$t$ | ${\Gamma({\nu+1\over 2})\over\sqrt{\nu\pi}\Gamma({\nu\over 2})}(1+{x^{2}\over\nu})^{-(\nu+1)/2}$ | ${1\over 2}+\Gamma({\nu+1\over 2}){{}_{2}F_{1}({1\over 2},{\nu+1\over 2};{3\over 2};-{x^{2}\over\nu})\over\sqrt{\nu\pi}\Gamma({\nu\over 2})}$ | $\nu>0$ | $\mathbbm{R}$
Trimodal Normal | $\phi(x)={1\over\sqrt{2\pi}}{\rm e}^{-x^{2}/2}$ | $\Phi(x)=\int_{-\infty}^{x}\phi(t){\rm d}t$ | $-$ | $\mathbbm{R}$

Here ${}_{2}F_{1}$ is the hypergeometric function and ${\rm
erf}(x)={2}\int_{0}^{x}{\rm e}^{-t^{2}}{\rm d}t/\sqrt{\pi}$ is the error
function (also called the Gauss error function).

##  3 Structural properties

In this section, some basic properties such as trimodality for symmetric
kernels, moments and truncated moments, and entropies for $X\sim{\rm
TD}(\bm{\theta})$ are discussed in detail.

###  3.1 Trimodality for a class of symmetric kernel densities

In this subsection, we suposse that the kernel density $g$ in (2.1) has the
following form

$\displaystyle g(x)=g(0)\,{\rm
e}^{{}^{\textstyle-\int_{0}^{x}t\mathfrak{h}(t^{2})\,{\rm d}t}},\quad
x\in\mathbbm{R},$ (3.1)

for some positive real function $\mathfrak{h}$ such that the integral
$\int_{0}^{x}t\mathfrak{h}(t^{2})\,{\rm d}t$ exists. Notice that (3.1) is
equivalent to

$\displaystyle g^{\prime}(x)=-x\mathfrak{h}(x^{2})g(x),\quad x\in\mathbbm{R}.$
(3.2)

It is immediate to verify that $g$, as defined in (3.1), is symmetric about
zero, that is, $g(x)=g(-x)$ on the real line $D=\mathbbm{R}$. For example, in
the Laplace, Cauchy, Student-$t$ and Normal kernel densities (see Table 1) we
have $\mathfrak{h}(y)=1/\sqrt{y}$, $\mathfrak{h}(y)=2/(1+y)$,
$\mathfrak{h}(y)=(\nu+1)/(\nu+{y})$ and $\mathfrak{h}(y)=1$, $\forall y>0$,
respectively.

Moreover, we assume that $\mathfrak{h}$ has the following form

$\displaystyle\mathfrak{h}(y)={C\over(\alpha^{2}A+y)^{\beta-p}}\quad\text{for
}\,C\geqslant 1,A\geqslant 0\ \text{and}\ \beta>p.$ (3.3)

That is, $\mathfrak{h}(y)$, $y>0$, decays polynomially. The function
$\mathfrak{h}(y)=1$ corresponding to the Normal distribution is not of the
form (3.3), then the next result cannot be applied and a separate study must
be carried out. In this paper, the Gaussian case will be studied in detail in
Section 4.

For a formal proof of following lemma, see Section A of Appendix.

######  Lemma 3.1.

Let $\mathfrak{h}$ be as in (3.3). For some $\rho>0$ the function
$\mathfrak{R}$, defined by

$\displaystyle\mathfrak{R}(y)={2\delta}\biggl{[}{(1/\alpha^{2})^{p}\over\Gamma(p)}\,y^{p-1}{\rm
e}^{-y/\alpha^{2}}\biggr{]}-\biggl{[}\rho+\delta\,{\gamma(p,y/\alpha^{2})\over\Gamma(p)}\biggr{]}\mathfrak{h}(y),\quad
y>0,$

has at most two real roots.

######  Proposition 3.2.

Let $g$ be a kernel density as in (3.1), with $\mathfrak{h}$ as in (3.3). A
point $x\in\mathbbm{R}$ is a critical point of density (2.1) if $x=\mu$ or
$\mathfrak{R}[(x-\mu)^{2}/\sigma^{2}]=0,$ where $\mathfrak{R}$ is as in Lemma
3.1.

###### Proof.

The proof is immediate since, by using Equation (3.2), the first-order
derivative of $f(x;\bm{\theta})$, with respect to $x$, is given by
$f^{\prime}(x;\bm{\theta})=[(x-\mu)/\sigma]g((x-\mu)/\sigma)\mathfrak{R}[(x-\mu)^{2}/\sigma^{2}]/(\sigma
Z_{\bm{\theta}})$. ∎

######  Theorem 3.3 (Uni- bi- or trimodality).

If $X\sim{\rm TD}(\bm{\theta})$ then the following hold:

  * (1)

If $\mathfrak{R}$ has no real roots then $f(x;\bm{\theta})$ is unimodal with
mode $x=\mu$.

  * (2)

If $\mathfrak{R}$ has one real root then $f(x;\bm{\theta})$ is bimodal with
minimum point $x=\mu$.

  * (3)

If $\mathfrak{R}$ has two distinct real roots then $f(x;\bm{\theta})$ is
trimodal where $x=\mu$ is one of the modes.

###### Proof.

It is clear that if $\mathfrak{R}$ has no real roots, by Proposition 3.2,
$x=\mu$ is the only critical point of the density $f$. Since
$\lim_{x\to\pm\infty}f(x;\bm{\theta})=0$, the point $x=\mu$ is a mode. This
proves Item (1).

In order to prove Item (2), we suppose that $\mathfrak{R}$ has one real root,
denoted by $a$. By Proposition 3.2, it follows that $x=\mu$ and
$x=\mu\pm\sigma\sqrt{a}$ are three critical points of $f$. Since
$\lim_{x\to\pm\infty}f(x;\bm{\theta})=0$, the point $x=\mu$ is a minimum and
$x=\mu\pm\sigma\sqrt{a}$ are two symmetrical modes. This proves the second
item.

Now, we assume that $\mathfrak{R}$ has two distinct real roots, denoted by $a$
and $b$. Without loss of generality, we can assume that $a<b$. Again, by
Proposition 3.2 we have that $x=\mu$, $x=\mu\pm\sigma\sqrt{a}$ and
$x=\mu\pm\sigma\sqrt{b}$ are five critical points of $f$. Since
$\lim_{x\to\pm\infty}f(x;\bm{\theta})=0$ and $a<b$, the critical points
$x=\mu$ and $x=\mu\pm\sigma\sqrt{b}$ are modes and $x=\mu\pm\sigma\sqrt{a}$
are minimum points. Hence, the proof of Item (3) follows. ∎

###  3.2 Moments

######  Theorem 3.4.

Let $X\sim{\rm TD}(\bm{\theta})$ and $L:\mathbbm{R}\to\mathbbm{R}$ be a Borel-
measurable function. Then, the expectation of random variable $L(X)$ with
$X\leqslant b$ and $b\in\mathbbm{R}$, is given by

$\displaystyle\mathbbm{E}\big{[}\mathds{1}_{\\{X\leqslant b\\}}L(X)\big{]}$
$\displaystyle={\rho\sigma\over
Z_{\bm{\theta}}}\,\mathbbm{E}\big{[}\mathds{1}_{\\{W\leqslant{b-\mu\over\sigma}\\}}L(W_{\mu,\sigma})\big{]}$
$\displaystyle+{\delta\sigma\over
Z_{\bm{\theta}}}\left\\{\mathbbm{E}\left[\mathds{1}_{\\{Y\leqslant({b-\mu\over\alpha\sigma})^{2},\,W\leqslant{b-\mu\over\sigma}\\}}L(W_{\mu,\sigma})\right]+\mathbbm{E}\left[\mathds{1}_{\\{Y\geqslant({b-\mu\over\alpha\sigma})^{2},\,W\leqslant-\alpha\sqrt{Y}\\}}L(W_{\mu,\sigma})\right]\right\\}\mathds{1}_{\\{b<\mu\\}}$
$\displaystyle+{\delta\sigma\over
Z_{\bm{\theta}}}\left\\{\mathbbm{E}\left[\mathds{1}_{\\{W\leqslant-\alpha\sqrt{Y}\\}}L(W_{\mu,\sigma})\right]+\mathbbm{E}\left[\mathds{1}_{\\{Y\leqslant({b-\mu\over\alpha\sigma})^{2},\,\alpha\sqrt{Y}\leqslant
W\leqslant{b-\mu\over\sigma}\\}}L(W_{\mu,\sigma})\right]\right\\}\mathds{1}_{\\{b\geqslant\mu\\}},$

where $W_{\mu,\sigma}=\sigma W+\mu$, $W$ is a continuous random variable with
CDF $G$ (that for brevity we write $W\stackrel{{\scriptstyle d}}{{=}}G$),
$Y\sim{\rm Gamma}(p,1)$, and $W$ and $Y$ are independent.

###### Proof.

By using the definition of expectation and by taking the change of variables
$w=(x-\mu)/\sigma$ and ${\rm d}x=\sigma{\rm d}w$, we have

$\displaystyle\mathbbm{E}\big{[}\mathds{1}_{\\{X\leqslant b\\}}L(X)\big{]}$
$\displaystyle={\rho\over Z_{\bm{\theta}}}\,\int_{\sigma
D+\mu}\mathds{1}_{\\{x\leqslant
b\\}}L(x)g\biggl{(}{x-\mu\over\sigma}\biggr{)}\,{\rm d}x$
$\displaystyle+{\delta\over Z_{\bm{\theta}}}\,\int_{\sigma
D+\mu}\mathds{1}_{\\{x\leqslant
b\\}}L(x)T\biggl{(}{x-\mu\over\sigma};\alpha,p\biggr{)}g\biggl{(}{x-\mu\over\sigma}\biggr{)}\,{\rm
d}x$ $\displaystyle={\rho\sigma\over
Z_{\bm{\theta}}}\int_{D}\mathds{1}_{\\{w\leqslant{b-\mu\over\sigma}\\}}L(\sigma
w+\mu)g(w){\rm d}w+{\delta\sigma\over Z_{\bm{\theta}}}\
\iint\limits_{\begin{subarray}{c}\scriptscriptstyle
w\leqslant{(b-\mu)/\sigma}\\\ \scriptscriptstyle 0<y\leqslant
w^{2}/\alpha^{2}\end{subarray}}\mathds{1}_{D}(\omega)\tau(w,y)\,{\rm d}y{\rm
d}w,$ (3.4)

where, for notational simplicity, we denote

$\displaystyle\tau(w,y)=L(\sigma w+\mu)g(w)\biggl{[}{y^{p-1}{\rm
e}^{-y}\over\Gamma(p)}\biggr{]}.$

There are two cases to consider according to whether
$\xi\coloneqq(b-\mu)/\sigma<0$ or $\xi\coloneqq(b-\mu)/\sigma\geqslant 0$; see
Figure 1 (a) and (b).

Figure 1: (a) $\\{w\leqslant\xi,0<y\leqslant w^{2}/\alpha^{2},\xi<0\\}$; (b)
$\\{w\leqslant\xi,0<y\leqslant w^{2}/\alpha^{2},\xi\geqslant 0\\}$.

In the former case,

$\displaystyle\iint\limits_{\begin{subarray}{c}\scriptscriptstyle
w\leqslant{(b-\mu)/\sigma}\\\ \scriptscriptstyle 0<y\leqslant
w^{2}/\alpha^{2}\end{subarray}}\mathds{1}_{D}(\omega)\tau(w,y)\,{\rm d}y{\rm
d}w$
$\displaystyle=\int_{0}^{({b-\mu\over\alpha\sigma})^{2}}\Biggl{(}\int_{-\infty}^{{b-\mu\over\sigma}}\mathds{1}_{D}(\omega)\tau(w,y)\,{\rm
d}w\Biggr{)}{\rm d}y$
$\displaystyle+\int_{({b-\mu\over\alpha\sigma})^{2}}^{\infty}\Biggl{(}\int_{-\infty}^{-\alpha\sqrt{y}}\mathds{1}_{D}(\omega)\tau(w,y)\,{\rm
d}w\Biggr{)}{\rm d}y$

and in the latter case,

$\displaystyle\iint\limits_{\begin{subarray}{c}\scriptscriptstyle
w\leqslant{(b-\mu)/\sigma}\\\ \scriptscriptstyle 0<y\leqslant
w^{2}/\alpha^{2}\end{subarray}}\mathds{1}_{D}(\omega)\tau(w,y)\,{\rm d}y{\rm
d}w$
$\displaystyle=\int_{0}^{\infty}\Biggl{(}\int_{-\infty}^{-\alpha\sqrt{y}}\mathds{1}_{D}(\omega)\tau(w,y)\,{\rm
d}w\Biggr{)}{\rm d}y$
$\displaystyle+\int_{0}^{({b-\mu\over\alpha\sigma})^{2}}\Biggl{(}\int_{\alpha\sqrt{y}}^{{b-\mu\over\sigma}}\mathds{1}_{D}(\omega)\tau(w,y)\,{\rm
d}w\Biggr{)}{\rm d}y.$

Hence, by combining the last two integral identities with (3.2), when $b<\mu$,
we have

$\displaystyle\mathbbm{E}\big{[}\mathds{1}_{\\{X\leqslant b\\}}L(X)\big{]}$
$\displaystyle={\rho\sigma\over
Z_{\bm{\theta}}}\int_{D}\mathds{1}_{\\{w\leqslant{b-\mu\over\sigma}\\}}L(\sigma
w+\mu)g(w){\rm d}w+{\delta\sigma\over
Z_{\bm{\theta}}}\int_{0}^{({b-\mu\over\alpha\sigma})^{2}}\Biggl{(}\int_{-\infty}^{{b-\mu\over\sigma}}\mathds{1}_{D}(\omega)\tau(w,y)\,{\rm
d}w\Biggr{)}\,{\rm d}y$ $\displaystyle+{\delta\sigma\over
Z_{\bm{\theta}}}\int_{({b-\mu\over\alpha\sigma})^{2}}^{\infty}\Biggl{(}\int_{-\infty}^{-\alpha\sqrt{y}}\mathds{1}_{D}(\omega)\tau(w,y)\,{\rm
d}w\Biggr{)}\,{\rm d}y$

and, when $b\geqslant\mu$,

$\displaystyle\mathbbm{E}\big{[}\mathds{1}_{\\{X\leqslant b\\}}L(X)\big{]}$
$\displaystyle={\rho\sigma\over
Z_{\bm{\theta}}}\int_{D}\mathds{1}_{\\{w\leqslant{b-\mu\over\sigma}\\}}L(\sigma
w+\mu)g(w){\rm d}w+{\delta\sigma\over
Z_{\bm{\theta}}}\int_{0}^{\infty}\Biggl{(}\int_{-\infty}^{-\alpha\sqrt{y}}\mathds{1}_{D}(\omega)\tau(w,y)\,{\rm
d}w\Biggr{)}\,{\rm d}y$ $\displaystyle+{\delta\sigma\over
Z_{\bm{\theta}}}\int_{0}^{({b-\mu\over\alpha\sigma})^{2}}\Biggl{(}\int_{\alpha\sqrt{y}}^{{b-\mu\over\sigma}}\mathds{1}_{D}(\omega)\tau(w,y)\,{\rm
d}w\Biggr{)}\,{\rm d}y.$

Then there are $W\stackrel{{\scriptstyle d}}{{=}}G$ and $Y\sim{\rm
Gamma}(p,1)$ independent so that, for $b<\mu$,

$\displaystyle\mathbbm{E}\big{[}\mathds{1}_{\\{X\leqslant b\\}}L(X)\big{]}$
$\displaystyle={\rho\sigma\over
Z_{\bm{\theta}}}\,\mathbbm{E}\big{[}\mathds{1}_{\\{W\leqslant{b-\mu\over\sigma}\\}}L(W_{\mu,\sigma})\big{]}$
$\displaystyle+{\delta\sigma\over
Z_{\bm{\theta}}}\left\\{\mathbbm{E}\left[\mathds{1}_{\\{Y\leqslant({b-\mu\over\alpha\sigma})^{2},\,W\leqslant{b-\mu\over\sigma}\\}}L(W_{\mu,\sigma})\right]+\mathbbm{E}\left[\mathds{1}_{\\{Y\geqslant({b-\mu\over\alpha\sigma})^{2},\,W\leqslant-\alpha\sqrt{Y}\\}}L(W_{\mu,\sigma})\right]\right\\}$

and, for $b\geqslant\mu$,

$\displaystyle\mathbbm{E}\big{[}\mathds{1}_{\\{X\leqslant b\\}}L(X)\big{]}$
$\displaystyle={\rho\sigma\over
Z_{\bm{\theta}}}\,\mathbbm{E}\big{[}\mathds{1}_{\\{W\leqslant{b-\mu\over\sigma}\\}}L(W_{\mu,\sigma})\big{]}$
$\displaystyle+{\delta\sigma\over
Z_{\bm{\theta}}}\left\\{\mathbbm{E}\left[\mathds{1}_{\\{W\leqslant-\alpha\sqrt{Y}\\}}L(W_{\mu,\sigma})\right]+\mathbbm{E}\left[\mathds{1}_{\\{Y\leqslant({b-\mu\over\alpha\sigma})^{2},\,\alpha\sqrt{Y}\leqslant
W\leqslant{b-\mu\over\sigma}\\}}L(W_{\mu,\sigma})\right]\right\\}.$

This completes the proof. ∎

By taking $b\to\infty$ in Theorem 3.4, with ${(X-\mu)/\sigma}$ instead of $X$,
we get

######  Corollary 3.5.

Under the hypotheses of Theorem 3.4,

$\displaystyle\mathbbm{E}\biggl{[}L\biggl{(}{X-\mu\over\sigma}\biggr{)}\biggr{]}={(\rho+\delta)\sigma\over
Z_{\bm{\theta}}}\,\mathbbm{E}[L(W)]+{\delta\sigma\over
Z_{\bm{\theta}}}\left\\{\mathbbm{E}\left[\mathds{1}_{\\{W\leqslant-\alpha\sqrt{Y}\\}}L(W)\right]-\mathbbm{E}\left[\mathds{1}_{\\{W\leqslant\alpha\sqrt{Y}\\}}L(W)\right]\right\\}.$

######  Corollary 3.6.

Under the hypotheses of Theorem 3.4, if

  * (a)

$b\to\infty$ and $L(x)=1$, $\forall x\in\sigma D+\mu$, then the formula (2.3)
for the normalizing factor $Z_{\bm{\theta}}$ is obtained.

  * (b)

$b=x\in\mathbbm{R}$ fixed and $L(x)=1$, $\forall x\in\sigma D+\mu$, then the
formula (2) of the CDF $F(x;\bm{\theta})$ is obtained.

By taking $b\to\infty$ in Theorem 3.4, with $L(x)=x^{n}$, $\forall x\in\sigma
D+\mu$ and $n\geqslant 1$ integer, by a binomial expansion, we have the next
formula for the moments.

######  Corollary 3.7.

The $n$-th moment of $X\sim{\rm TD}(\bm{\theta})$ is written as

$\displaystyle\mathbbm{E}(X^{n})=\sum_{k=0}^{n}\binom{n}{k}\mu^{n-k}\sigma^{k}\left\\{{{({{\rho}+\delta})\sigma\over
Z_{\bm{\theta}}}\,\mathbbm{E}(W^{k})+{\delta\sigma\over
Z_{\bm{\theta}}}\left[\mathbbm{E}\left(\mathds{1}_{\\{W\leqslant-\alpha\sqrt{Y}\\}}W^{k}\right)-\mathbbm{E}\left(\mathds{1}_{\\{W\leqslant\alpha\sqrt{Y}\\}}W^{k}\right)\right]}\right\\}.$

The above formula informs that the moments of $X$ (whenever they exist)
depends on the existence of moments of
$\mathds{1}_{\\{W\leqslant\pm\alpha\sqrt{Y}\\}}W$, with
$W\stackrel{{\scriptstyle d}}{{=}}G$ and $Y\sim{\rm Gamma}(p,1)$ independent.

###  3.3 Entropies

The Tsallis [31] entropy associated with the random variable $X\sim{\rm
TD}(\bm{\theta})$ is defined as

$\displaystyle S_{q}(X)=\begin{cases}\displaystyle-\int_{\sigma
D+\mu}f^{q}(x;\bm{\theta})\log_{q}f(x;\bm{\theta})\,{\rm d}x,&\text{ if}\
q\neq 1,\\\\[17.07182pt] \displaystyle-\int_{\sigma D+\mu}f(x;\bm{\theta})\log
f(x;\bm{\theta})\,{\rm d}x,&\text{ if}\ q=1,\end{cases}$

where, for $x>0$,

$\displaystyle\log_{q}(x)=\begin{cases}\displaystyle{x^{1-q}-1\over
1-q},&\text{ if}\ q\neq 1,\\\\[5.69046pt] \log(x),&\text{ if}\
q=1,\end{cases}$ (3.5)

represents a Box-Cox transformation in statistics (often called deformed
logarithm [20, 32]). Since $\log_{q}(x)\to\log(x)$ when $q\to 1$, we have
$S_{q}(X)\to S_{1}(X)$ when $q\to 1$. That is, when $q\to 1$ the usual
definition of Shannon’s entropy $S_{1}(X)$ [30] is recovered.

######  Proposition 3.8.

Under the hypotheses of Theorem 3.4,

$\displaystyle\mathbbm{E}\big{[}f^{q-1}(X;\bm{\theta})\big{]}$
$\displaystyle={({\rho}+\delta)\sigma\over
Z_{\bm{\theta}}}\,\mathbbm{E}\big{[}f^{q-1}(\sigma W+\mu;\bm{\theta})\big{]}$
$\displaystyle+{{\delta}\sigma\over
Z_{\bm{\theta}}}\,\left\\{\mathbbm{E}\left[\mathds{1}_{\\{W\leqslant-\alpha\sqrt{Y}\\}}f^{q-1}(\sigma
W+\mu;\bm{\theta})\right]-\mathbbm{E}\left[\mathds{1}_{\\{W\leqslant\alpha\sqrt{Y}\\}}f^{q-1}(\sigma
W+\mu;\bm{\theta})\right]\right\\},$

where $W\stackrel{{\scriptstyle d}}{{=}}G$, $Y\sim{\rm Gamma}(p,1)$, and $W$
and $Y$ are independent. Consequently, the Tsallis entropy of $X$ (whenever it
exists) depends on the existence of truncated moments of
$f^{q-1}(W;\bm{\theta})$.

###### Proof.

By taking $b\to\infty$ in Theorem 3.4, with $L(x)=f^{q-1}(x;\bm{\theta})$,
$\forall x\in\sigma D+\mu$, the proof follows. ∎

For a formal proof of next result, see Section A of Appendix.

######  Proposition 3.9.

If $S_{q}(W)$ exists and $q>0,$ then $S_{q}(X)$, with $X\sim{\rm
TD}(\bm{\theta})$, also exists.

######  Remark 3.10.

As consequence of Proposition 3.9, by letting $q\to 1$, the Shannon entropy
$S_{1}(X)$ exists whenever $S_{1}(W)$ also exists.

For a rigorous proof of next result, see Section A of Appendix.

######  Proposition 3.11.

Under the hypotheses of Theorem 3.4, the Shannon entropy of $X\sim{\rm
TD}(\bm{\theta})$ is written as

$\displaystyle S_{1}(X)$
$\displaystyle=\log(Z_{\bm{\theta}})-{{({\rho}+\delta)\sigma\over
Z_{\bm{\theta}}}\,\mathbbm{E}\big{[}\log(\rho+\delta T(W;\alpha,p))\big{]}}$
$\displaystyle-{{{\delta}\sigma\over
Z_{\bm{\theta}}}\,\left\\{\mathbbm{E}\left[\mathds{1}_{\\{W\leqslant-\alpha\sqrt{Y}\\}}\log(\rho+\delta
T(W;\alpha,p))\right]-\mathbbm{E}\left[\mathds{1}_{\\{W\leqslant\alpha\sqrt{Y}\\}}\log(\rho+\delta
T(W;\alpha,p))\right]\right\\}}$ $\displaystyle+{{({\rho}+\delta)\sigma\over
Z_{\bm{\theta}}}\,S_{1}(W)+{{\delta}\sigma\over
Z_{\bm{\theta}}}\,\left\\{{\mathbbm{E}}\left[S_{1}\big{(}\mathds{1}_{\\{W\leqslant-\alpha\sqrt{Y}\\}}W\big{)}\right]-{\mathbbm{E}}\left[S_{1}\big{(}\mathds{1}_{\\{W\leqslant\alpha\sqrt{Y}\\}}W\big{)}\right]\right\\}},$

where $T$ was defined in (2.2).

###  3.4 Stochastic representation

Let $h(u)$, $0<u<1$, be a PDF with corresponding CDF $H$. Let
$\mathcal{S}:D\to(0,1)$ be an injective and increasing transformation, where
$D$ is a non-empty set of $\mathbbm{R}$. We consider the following CDF:

$\displaystyle F(z)=\int_{0}^{\mathcal{S}(z)}h(u){\rm
d}u=H(\mathcal{S}(z)),\quad z\in D.$ (3.6)

We also define by $f$ to the corresponding PDF of $F$. That is,
$F^{\prime}(z)=f(z)=h(\mathcal{S}(z))\mathcal{S}^{\prime}(z)$ for almost all
$z\in D$.

We define the PDF $h$ as follows

$\displaystyle h(u)={\sigma\over Z_{\bm{\theta}}}\,{\big{[}\rho+\delta
T\big{(}G^{-1}(u);\alpha,p\big{)}\big{]}},\quad 0<u<1,\rho\geqslant
0,\delta\geqslant 0,\alpha>0,\sigma>0,$ (3.7)

where $G$ and $G^{-1}$, respectively, are the CDF defined in Section 2 and its
inverse function, and $T$ is as in (2.2). When $\delta=0$, $h$ reduces to the
continuous uniform distribution on the interval $(0,1)$. The CDF $H$ of $h$ is
given by

$\displaystyle H(u)$ $\displaystyle={\rho\over
Z_{\bm{\theta}}}\,u+{\delta\over
Z_{\bm{\theta}}}\left\\{\mathbbm{E}\big{[}G(-\alpha\sqrt{Y})\big{]}+\mathcal{S}\big{(}G^{-1}(u);\alpha,p\big{)}u\right\\}$
$\displaystyle-{\delta\over
Z_{\bm{\theta}}}\left\\{\mathbbm{E}\Big{[}\mathbbm{1}_{\\{Y\leqslant[{G^{-1}(u)\over\alpha}]^{2}\\}}G(-\alpha\sqrt{Y})\Big{]}\mathbbm{1}_{\\{G^{-1}(u)<0\\}}+\mathbbm{E}\Big{[}\mathbbm{1}_{\\{Y\leqslant[{G^{-1}(u)\over\alpha}]^{2}\\}}G(\alpha\sqrt{Y})\Big{]}\mathbbm{1}_{\\{G^{-1}(u)\geqslant
0\\}}\right\\},$

where $Y\sim{\rm Gamma}(p,1)$.

If $\mathcal{S}:D\to(0,1)$ is defined by $\mathcal{S}(z)=G(z)$, $\forall z\in
D$, by (3.6), the family of trimodal distributions in (2.1) is obtained. That
is,

$\displaystyle
F(x;\bm{\theta})=F\biggl{(}{x-\mu\over\sigma}\biggr{)}=H\biggl{(}G\Big{(}{x-\mu\over\sigma}\Big{)}\biggr{)}$
(3.8)

and $f(x;\bm{\theta})={1\over\sigma}\,f((x-\mu)/\sigma)$.

If $U$ is distributed according to (3.7) and $X\sim{\rm TD}(\bm{\theta})$,
then, by (3.8), the random variable $X$ admits the following stochastic
representation:

$\displaystyle X=\mu+\sigma G^{-1}(U).$

##  4 The Gaussian case

In this section, the standard normal kernel density $g(x)=\phi(x)$,
$x\in\mathbbm{R}$, is plugged into (2.1). Some structural properties as
modality, moments, entropies and rate of the distribution are discussed.

The CDF corresponding to the normal kernel density $g$ is $G(x)=\Phi(x)$, with

$\displaystyle\Phi(x)=\int_{-\infty}^{x}\phi(t)\,{\rm d}t={1\over
2}\biggl{[}1\pm{\rm erf}\biggl{(}\pm{x\over\sqrt{2}}\biggl{)}\biggr{]},$ (4.1)

where ${\rm erf}(x)$ is the error function. So, in this section, we consider
the following PDF

$\displaystyle f(x;\bm{\theta})={1\over Z_{\bm{\theta}}}\,\biggl{[}\rho+\delta
T\biggl{(}{x-\mu\over\sigma};\alpha,p\biggr{)}\biggr{]}\,\phi\biggl{(}{x-\mu\over\sigma}\biggr{)},\quad
x\in\mathbbm{R},$ (4.2)

where $T$ is as in (2.2) and $Z_{\bm{\theta}}$ is the normalizing factor
(2.3), which is given by

$\displaystyle
Z_{\bm{\theta}}=(\rho+\delta)\sigma+\delta\sigma\left\\{{\mathbbm{E}}\big{[}\Phi(-\alpha\sqrt{Y})\big{]}-{\mathbbm{E}}\big{[}\Phi(\alpha\sqrt{Y})\big{]}\right\\}.$
(4.3)

We will denote $X\sim{\rm TD}_{\Phi}(\bm{\theta})$ for a random variable $X$
that follows (4.2). Plots of the TDΦ density, where trimodality is observed,
are given in Figures 2(a) and 2(b).

((a)) The strict form of modality for trimodal normal distribution.

((b)) The soft form of modality for trimodal normal distribution.

Figure 2: PDFs of trimodal normal distribution.

###  4.1 A study on the modality

For a mathematical proof of next result, see Section A of Appendix.

######  Lemma 4.1.

The function $R$, defined by

$\displaystyle{R}(y)={2\delta}\,{{(1/\alpha^{2})^{p-1}{\rm
e}^{-y/\alpha^{2}}\over\Gamma(p)}}-\biggl{[}\rho+\delta\,{\gamma(p,y/\alpha^{2})\over\Gamma(p)}\biggr{]},\quad
y>0,$

has at most two real roots.

######  Remark 4.2.

Notice that when $\rho$ is sufficiently large, the function $R$ has no roots.
For sufficiently small $\rho$, $R$ have one, or two roots depending on whether
$p\leqslant 1$ or $p>1$. Furthermore, when $\delta=0$ (with $\rho>0$) or
$\alpha\to\infty$, the function $R$ has no roots.

######  Proposition 4.3.

A point $x\in\mathbbm{R}$ is a critical point of density (4.2) if $x=\mu$ or
$R[(x-\mu)^{2}/\sigma^{2}]=0$, where $R$ is as in Lemma 4.1.

###### Proof.

The proof follows from identity
$f^{\prime}(x;\bm{\theta})=[(x-\mu)/\sigma]g((x-\mu)/\sigma){R}[(x-\mu)^{2}/\sigma^{2}]/(\sigma
Z_{\bm{\theta}})$. ∎

######  Theorem 4.4 (Uni- bi- or trimodality).

Let $X\sim{\rm TD}_{\Phi}(\bm{\theta})$ and $R$ be the function defined in
Lemma 4.1. The following hold:

  * (1)

If $R$ has no real roots then $f(x;\bm{\theta})$ is unimodal with mode
$x=\mu$.

  * (2)

If $R$ has one real root then $f(x;\bm{\theta})$ is bimodal with minimum point
$x=\mu$.

  * (3)

If $R$ has two distinct real roots then $f(x;\bm{\theta})$ is trimodal where
$x=\mu$ is one of the modes.

###### Proof.

The proof follows the same steps of proof of Theorem 3.3 by taking $R$ instead
$\mathfrak{R}$, and by using Proposition 4.3 instead of Proposition 3.2. ∎

###  4.2 The normalizing factor

In what follows we find a closed expression for the normalizing factor
$Z_{\bm{\theta}}$ in (4.2). By (4.3) it is necessary to calculate
${\mathbbm{E}}[\Phi(\pm\alpha\sqrt{Y})]$, where $Y\sim{\rm Gamma}(p,1)$.
Indeed, using the formula (4.1) and then taking the change of variables
$z=\sqrt{y}$ and ${\rm d}y=2z{\rm d}z$, we obtain

$\displaystyle{\mathbbm{E}}\big{[}\Phi(\pm\alpha\sqrt{Y})\big{]}=\int_{0}^{\infty}\Phi(\pm\alpha\sqrt{y})\,{y^{p-1}{\rm
e}^{-y}\over\Gamma(p)}\,{\rm d}y$ $\displaystyle={1\over 2}+{1\over
2}\,\int_{0}^{\infty}{\rm
erf}\biggl{(}\pm{\alpha\sqrt{y}\over\sqrt{2}}\biggr{)}\,{y^{p-1}{\rm
e}^{-y}\over\Gamma(p)}\,{\rm d}y$ (4.4) $\displaystyle={1\over
2}+\int_{0}^{\infty}{\rm erf}\biggl{(}\pm{\alpha
z\over\sqrt{2}}\biggr{)}\,{z^{2p-1}{\rm e}^{-z^{2}}\over\Gamma(p)}\,{\rm d}z.$
(4.5)

By using the formula (see Item 8 in Subsection 4.3 of reference [25]):

$\displaystyle\int_{0}^{\infty}{\rm erf}(ax)x^{p}{\rm e}^{-b^{2}x^{2}}\,{\rm
d}x={a\over\sqrt{\pi}}\,b^{-p-2}\Gamma\Big{(}{p\over
2}+1\Big{)}\,_{2}F_{1}\biggl{(}{1\over 2},{p\over 2}+1;{3\over 2};-{a^{2}\over
b^{2}}\biggr{)},\quad b^{2}>0,p>-2,$

where ${}_{2}F_{1}(a_{1},a_{2};b_{1};x)$ is the generalized Hypergeometric
function, we have

$\displaystyle\int_{0}^{\infty}{\rm erf}\biggl{(}\pm{\alpha
z\over\sqrt{2}}\biggr{)}\,{z^{2p-1}{\rm e}^{-z^{2}}\over\Gamma(p)}\,{\rm
d}z=\pm{\alpha\Gamma(p+{1\over
2})\over\sqrt{2\pi}\,\Gamma(p)}\,_{2}F_{1}\biggl{(}{1\over 2},p+{1\over
2};{3\over 2};-{\alpha^{2}\over 2}\biggr{)}.$

Therefore,

$\displaystyle{\mathbbm{E}}\big{[}\Phi(\pm\alpha\sqrt{Y})\big{]}={1\over
2}\pm{\alpha\Gamma(p+{1\over
2})\over\sqrt{2\pi}\,\Gamma(p)}\,_{2}F_{1}\biggl{(}{1\over 2},p+{1\over
2};{3\over 2};-{\alpha^{2}\over 2}\biggr{)}.$ (4.6)

Replacing (4.6) in (4.3), we obtain the following closed expression for the
normalizing factor $Z_{\bm{\theta}}$:

$\displaystyle
Z_{\bm{\theta}}=(\rho+\delta)\sigma-{{2}\delta\sigma\alpha\Gamma(p+{1\over
2})\over\sqrt{2\pi}\,\Gamma(p)}\,_{2}F_{1}\biggl{(}{1\over 2},p+{1\over
2};{3\over 2};-{\alpha^{2}\over 2}\biggr{)}.$ (4.7)

The flexibility of the TD model is shown. Note that the TD model can be
unimodal, bimodal or trimodal. Figures 2(a) and 2(b) show how the TD density
function is influenced by parameters $\alpha$, $\rho$ and $\delta$.

###  4.3 Cumulative distribution function

To determine the CDF of $X\sim{\rm TD}_{\Phi}(\bm{\theta})$ we use formula
(2). So, by (2), it is essential to calculate the expectation
$\mathbbm{E}[\mathds{1}_{\\{Y\leqslant({x-\mu\over\sigma\alpha})^{2}\\}}\Phi(\pm\alpha\sqrt{Y})],$
$x\in\mathbbm{R},$ where $Y\sim{\rm Gamma}(p,1)$. Indeed, similarly to the one
done in (4.5),

$\displaystyle{\mathbbm{E}}\Big{[}\mathds{1}_{\\{Y\leqslant({x-\mu\over\alpha\sigma})^{2}\\}}\Phi(\pm\alpha\sqrt{Y})\Big{]}={1\over
2}\,T\biggl{(}{x-\mu\over\sigma};\alpha,p\biggr{)}+{1\over\Gamma(p)}\,I\biggl{(}{x-\mu\over\alpha\sigma};\pm{\alpha\over\sqrt{2}},2p-1,1\biggr{)},$
(4.8)

where $I(u;a,b,c)=\int_{0}^{u}{\rm erf}(ax)x^{b}{\rm e}^{-c^{2}x^{2}}\,{\rm
d}x,\ u>0,a\in\mathbbm{R},b>0,c>0.$ Since $|{\rm erf}(x)|\leqslant 1$, we get
$|I(u;a,b,c)|\leqslant I(u;0,b,c)<\infty$, and therefore, $I(u;a,b,c)$ always
exists. In general, closed form solutions for the definite integral
$I(u;a,b,c)$ are not available in terms of commonly used functions.

Replacing the formulas (4.6) and (4.8) in (2), we obtain the following closed
formula for the CDF of $X\sim{\rm TD}_{\Phi}(\bm{\theta})$:

$\displaystyle F(x;\bm{\theta})=\begin{cases}\displaystyle{\sigma\over
Z_{\bm{\theta}}}\displaystyle\biggl{[}\rho+\delta
T\biggl{(}{x-\mu\over\sigma};\alpha,p\biggr{)}\biggr{]}\Phi\biggl{(}{x-\mu\over\sigma}\biggr{)}+{\delta\sigma\over
Z_{\bm{\theta}}}\biggl{[}{1\over 2}-{\alpha\Gamma(p+{1\over
2})\over\sqrt{2\pi}\,\Gamma(p)}\,_{2}F_{1}\biggl{(}{1\over 2},p+{1\over
2};{3\over 2};-{\alpha^{2}\over 2}\biggr{)}\biggr{]}&\\\\[14.22636pt]
\displaystyle\hskip 116.65646pt-{\delta\sigma\over
Z_{\bm{\theta}}}\,\left[{1\over
2}\,T\biggl{(}{x-\mu\over\sigma};\alpha,p\biggr{)}-{1\over\Gamma(p)}\,I\biggl{(}{x-\mu\over\alpha\sigma};{\alpha\over\sqrt{2}},2p-1,1\biggr{)}\right],&{\rm
if}\,x<\mu,\\\\[28.45274pt] \displaystyle{\sigma\over
Z_{\bm{\theta}}}\displaystyle\biggl{[}\rho+\delta
T\biggl{(}{x-\mu\over\sigma};\alpha,p\biggr{)}\biggr{]}\Phi\biggl{(}{x-\mu\over\sigma}\biggr{)}+{\delta\sigma\over
Z_{\bm{\theta}}}\biggl{[}{1\over 2}-{\alpha\Gamma(p+{1\over
2})\over\sqrt{2\pi}\,\Gamma(p)}\,_{2}F_{1}\biggl{(}{1\over 2},p+{1\over
2};{3\over 2};-{\alpha^{2}\over 2}\biggr{)}\biggr{]}&\\\\[14.22636pt]
\displaystyle\hskip 116.65646pt-{\delta\sigma\over
Z_{\bm{\theta}}}\,\left[{1\over
2}\,T\biggl{(}{x-\mu\over\sigma};\alpha,p\biggr{)}+{1\over\Gamma(p)}\,I\biggl{(}{x-\mu\over\alpha\sigma};{\alpha\over\sqrt{2}},2p-1,1\biggr{)}\right],&{\rm
if}\,x\geqslant\mu,\end{cases}$

where $Z_{\bm{\theta}}$ and $T$ are as in (4.7) and (2.2), respectively.

######  Remark 4.5.

As expected, since $\Phi$ has a symmetric distribution around $0$,
$F(\mu;\bm{\theta})={1/2}$ and then, $Q_{2}=\mu$ is the median and the mean
for $X\sim{\rm TD}_{\Phi}(\bm{\theta})$.

######  Remark 4.6.

Some examples where
$I({x-\mu\over\alpha\sigma};\pm{\alpha\over\sqrt{2}},2p-1,1)$ in (4.8) admits
a closed form are:

  1. 1.

By taking $\alpha=\sqrt{2}$ and by using the following formula (see Item 4 in
Subsection 1.5.3, pp. 31, of [29]):

$\displaystyle\int_{0}^{u}{\rm erf}(ax)x^{\lambda}{\rm e}^{-a^{2}x^{2}}\,{\rm
d}x={2a\over\sqrt{\pi}(\lambda+2)}\,u^{\lambda+2}\,_{2}F_{2}\biggl{(}{1},{\lambda\over
2}+1;{3\over 2},{\lambda\over 2}+2;{a^{2}u^{2}}\biggr{)},\quad\lambda>-2,$

where ${}_{2}F_{2}(a_{1},a_{2};b_{1},b_{2};x)$ is the generalized
Hypergeometric, we have

$\displaystyle
I\biggl{(}{x-\mu\over\alpha\sigma};\pm{\alpha\over\sqrt{2}},2p-1,1\biggr{)}=\pm{2\over\sqrt{\pi}(2p+1)}\,\biggl{(}{x-\mu\over\alpha\sigma}\biggr{)}^{2p+1}\,_{2}F_{2}\biggl{(}{1},p+{1\over
2};{3\over 2},p+{3\over
2};\Big{(}{\mu-x\over\alpha\sigma}\Big{)}^{2}\biggr{)}.$

  2. 2.

By taking $\alpha=\sqrt{2}$ and $p=3/2$,

$\displaystyle
I\biggl{(}{x-\mu\over\alpha\sigma};\pm{\alpha\over\sqrt{2}},2p-1,1\biggr{)}=\mp{1\over
2}\,\biggl{(}{x-\mu\over\sigma\sqrt{2}}\biggr{)}{\rm
erf}\biggl{(}{x-\mu\over\sigma\sqrt{2}}\biggr{)}\exp\biggl{[}{-\biggl{(}{x-\mu\over\sigma\sqrt{2}}\biggr{)}^{2}}\biggr{]}$
$\displaystyle\pm{1\over
4\sqrt{\pi}}\,\biggl{\\{}1-\exp\biggl{[}{-2\biggl{(}{x-\mu\over\sigma\sqrt{2}}\biggr{)}^{2}}\biggr{]}\biggr{\\}}\pm{\sqrt{\pi}\over
8}\,\biggl{[}{\rm
erf}\biggl{(}{x-\mu\over\sigma\sqrt{2}}\biggr{)}\biggr{]}^{2}.$

###  4.4 Moments and Shannon entropy

A routine calculation shows that the moments of $X\sim{\rm
TD}_{\Phi}(\bm{\theta})$ can be written as (for more details, see Section B in
the Appendix)

$\displaystyle\mathbbm{E}(X^{n})$ $\displaystyle={({{\rho}+\delta})\sigma\over
Z_{\bm{\theta}}}\sum_{k=0}^{n}\binom{n}{k}\mu^{n-k}\sigma^{k}\,{2^{-k/2}k!\over(k/2)!}\,\mathds{1}_{\\{k\,{\rm
even}\\}}$ $\displaystyle+{\delta\sigma\over
Z_{\bm{\theta}}\sqrt{2\pi}\Gamma(p)}\sum_{k=0}^{n}\binom{n}{k}\mu^{n-k}\sigma^{k}{\big{[}1-(-1)^{k-1}\big{]}}\sum_{m=0}^{\lfloor
k/2\rfloor}\sum_{j=0}^{\lfloor
k_{m}/2\rfloor}{c_{k,m}d_{k_{m},j}\,\alpha^{k_{m}-2j}\over(1+{\alpha^{2}\over
2})^{{k_{m}-2j\over 2}+p}}\,\Gamma\biggl{(}{k_{m}-2j\over 2}+p\biggr{)},$

where $p\geqslant 1$ is an integer, $k_{m}=k-2m-1$, $Z_{\bm{\theta}}$ is as in
(4.7), $c_{n,m}$ and $d_{n_{m},j}$ are as in (B.2). In particular,
$\mathbbm{E}(X)=\mu$ and ${\rm Var}(X)=\sigma^{2}$. As a consequence, from
Corollary 3.5 the closed expressions for skewness and kurtosis of random
variable $X\sim{\rm TD}_{\Phi}(\bm{\theta})$ are easily obtained.

On the other hand, a formula for the Shannon entropy is given by (for more
details, see Section C in the Appendix)

$\displaystyle S_{1}(X)$
$\displaystyle=\log(Z_{\bm{\theta}})-{{2({\rho}+\delta)\sigma\over
Z_{\bm{\theta}}}\,\sum_{k=0}^{\infty}{1\over(2k+1)}\biggl{[}1+{1\over(\rho+1)^{2k+1}}\sum_{i=p}^{\infty}\widetilde{c}_{i,k}\,{2^{-i}(2i)!\over
i!}\biggr{]}}$ $\displaystyle+{2{\delta}\sigma\over
Z_{\bm{\theta}}}\sum_{k=0}^{\infty}{1\over(2k+1)}\biggl{\\{}2\mathbbm{E}\big{[}\Phi(\alpha\sqrt{Y})\big{]}-1+{1\over(\rho+1)^{2k+1}}\sum_{i=p}^{\infty}\widetilde{c}_{i,k}\mathbbm{E}\big{(}\mathds{1}_{\\{|Z|\leqslant\alpha\sqrt{Y}\\}}{Z^{2i}}\big{)}\biggr{\\}}$
$\displaystyle+{({\rho}+\delta)\sigma\over
Z_{\bm{\theta}}}\,\biggl{[}\log(\sqrt{2\pi})+{1\over
2}\biggr{]}+{{\delta}\sigma\over
Z_{\bm{\theta}}}\left\\{{\alpha\Gamma(p+{1\over
2})\over\sqrt{2\pi}\Gamma(p)}\,{1\over(1+{\alpha^{2}\over 2})^{p+{1\over
2}}}-\biggl{[}\log(\sqrt{2\pi})+{1\over
2}\biggr{]}\big{(}2{\mathbbm{E}}\big{[}\Phi(\alpha\sqrt{Y})\big{]}-1\big{)}\right\\},$

where $Z_{\bm{\theta}}$ is as in (4.7), the coefficients
$\widetilde{c}_{i,k}$’s are determined by relation (C.6),
${\mathbbm{E}}\big{[}\Phi(\alpha\sqrt{Y})\big{]}$ is given in (4.6), and
$\mathbbm{E}(\mathds{1}_{\\{|Z|\leqslant\alpha\sqrt{Y}\\}}{Z^{2i}})=\mathbbm{E}(\mathds{1}_{\\{Z\leqslant\alpha\sqrt{Y}\\}}{Z^{2i}})-\mathbbm{E}(\mathds{1}_{\\{Z\leqslant-\alpha\sqrt{Y}\\}}{Z^{2i}})$
is determined by (B.3) in the Appendix.

###  4.5 Rate of the random variable $X\sim{\rm TD}_{\Phi}(\bm{\theta})$

Following the reference [23], the rate of a continuous random variable $X$ is
given by

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

where $f_{X}(x)$ denotes its respective PDF.

A simple computation shows that, the rate of $X\sim{\rm
TD}_{\Phi}(\bm{\theta})$ is

$\displaystyle\tau_{{\rm
TD}_{\Phi}(\bm{\theta})}={1\over\sigma}\lim_{x\to\infty}\left[{x-\mu\over\sigma}-{2\over\alpha\Gamma(p)}\,{({x-\mu\over\alpha\sigma})^{2p-1}{\rm
e}^{-({x-\mu\over\alpha\sigma})^{2}}\over\rho+T({x-\mu\over\sigma};\alpha,p)}\right]=\infty.$

Then, far enough out in the tail, the distribution of $X\sim{{\rm
TD}_{\Phi}(\bm{\theta})}$ looks like a Normal distribution, as expected. In
addition, we have some comparisons between the rate of $X\sim{{\rm
TD}_{\Phi}(\bm{\theta})}$ with the rates of some random variables with known
distributions in the literature: Inverse-gamma, Log-normal, Generalized-
Pareto, BGumbel [26], BWeibull [9], BGamma [8], 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 BGumbel(\mu,\beta,\delta)}=\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
TD}_{\Phi}(\bm{\theta})}=\tau_{{\rm Normal}(\mu,\sigma^{2})}=\infty.$

In other words, the tail of the normal distribution, of $X\sim{\rm
TD}_{\Phi}(\bm{\theta})$ and of ${\rm BWeibull}(\alpha>1,\beta,\delta)$ are
lighter than the tail of the other distributions specified above.

##  5 The Model arising from a Generalized Mixture

The weighted distributions of random variable $(W-\mu)/\sigma$, with weight
function $w_{k}$, have its PDF defined by

$\displaystyle f_{k}(x;\mu,\sigma)={1\over\int_{D}w_{k}(y)g(y){\rm
d}y}\,w_{k}{\biggl{(}{x-\mu\over\sigma}\biggr{)}}g\biggl{(}{x-\mu\over\sigma}\biggr{)},\quad{x-\mu\over\sigma}\in
D,\ k=0,1,\ldots.$

By using the power series expansion of the incomplete gamma function in (C.4),
note that the PDF $f(x;\bm{\theta})$ in (2.1) interprets as an infinite
(generalized) mixture of weighted distributions of $(W-\mu)/\sigma$ with the
weight functions $w_{k}(y)=y^{2k}$, $k=0,1,\ldots$, and with same parameter
vector $(\mu,\sigma)$. That is,

$\displaystyle
f(x;\bm{\theta})=c_{0}f_{0}(x;\mu,\sigma)+\sum_{k=p}^{\infty}c_{k}f_{k}(x;\mu,\sigma),$
(5.1)

where the constants $c_{0}$, $c_{k}$, with $k=p,p+1,\ldots$, that depends only
on $\bm{\theta}$, respectively, are given by

$\displaystyle c_{0}={\rho\sigma\over Z_{\bm{\theta}}};\quad
c_{k}={\delta\sigma\over\Gamma(p)Z_{\bm{\theta}}}\,{(-1)^{k-p}\mathbbm{E}(W^{2k})\over
k(k-p)!\alpha^{2k}},\quad k=p,p+1,\ldots.$

It is straightforward to verify that $c_{0}+\sum_{k=p}^{\infty}c_{k}=1$. Note
that the weights $c_{k}$ can take on negative values.

If $X\sim{\rm TD}_{\Phi}(\bm{\theta})$ and $L:\mathbbm{R}\to\mathbbm{R}$ is a
Borel-measurable function, then, by representation (5.1), the expectation of
truncated random variable $\mathds{1}_{\\{X\leqslant b\\}}L(X)$, with
$b\in\mathbbm{R}$, is given by

$\displaystyle\mathbbm{E}\big{[}\mathds{1}_{\\{X\leqslant
b\\}}L(X)\big{]}=c_{0}\mathbbm{E}_{0}\big{[}\mathds{1}_{\\{X\leqslant
b\\}}L(X)\big{]}+\sum_{k=p}^{\infty}c_{k}\mathbbm{E}_{k}\big{[}\mathds{1}_{\\{X\leqslant
b\\}}L(X)\big{]}.$

But this representation as infinite sum is not as informative as that
representation of Theorem 3.4. For this reason in the paper, we had been
concerned with providing mathematical representations that allow us to find
the closed expressions for some characteristics of the distribution (such as
normalizing factor, CDF, moments and entropies) as a function of the known
mathematical functions.

##  6 Inference

The maximum logq-likelihood estimation (MLqE) and Fisher information matrix
for this estimation method are introduced. MLqE drops to the maximum
likelihood estimation (MLE) if $q\to 1$ [19, 20, 21, 32].

###  6.1 The Logq-Likelihood Function

For a random sample $X_{1},\ldots,X_{n}$ from the random variable $X\sim{\rm
TD}(\bm{\theta})$, with parameter vector
$\bm{\theta}=(\mu,\sigma,\alpha,\rho,\delta)$, let us suppose that
$x_{1},\ldots,x_{n}$ are the observed values of $X_{1},\ldots,X_{n}$. The MLqE
of $\bm{\theta}$ is defined as

$\displaystyle
l_{q}(\bm{\theta};\bm{x})=\sum_{i=1}^{n}\log_{q}f(x_{i};\bm{\theta}),\quad\bm{x}=(x_{1},\ldots,x_{n})\in\mathbbm{R}^{n},$
(6.1)

where $\log_{q}(x)$ is the deformed logarithm defined in (3.5) [20, 32].

By taking partial derivatives in (6.1), with respect to
$\theta\in\\{\mu,\sigma,\alpha,\rho,\delta\\}$, we have

$\displaystyle{\partial
l_{q}(\bm{\theta};\bm{x})\over\partial\theta}=\sum_{i=1}^{n}f^{1-q}(x_{i};\bm{\theta})\,{\partial\log
f(x_{i};\bm{\theta})\over\partial\theta},$

with

$\displaystyle\log
f(x;\bm{\theta})=-\log(Z_{\bm{\theta}})+\log\biggl{[}\rho+\delta
T\biggl{(}{x-\mu\over\sigma};\alpha,p\biggr{)}\biggr{]}+\log
g\biggl{(}{x-\mu\over\sigma}\biggr{)}$ (6.2)

and $Z_{\bm{\theta}}$, $T$ are as in (2.3) and (2.2), respectively. In
general, the estimating equations for the parameters have the form

$\displaystyle\sum_{i=1}^{n}f^{1-q}(x_{i};\bm{\theta})\,{\partial\log
f(x_{i};\bm{\theta})\over\partial\theta}=0,\quad\theta\in\\{\mu,\sigma,\alpha,\rho,\delta\\},$
(6.3)

where the first-order partial derivatives of the function $\log
f(x;\bm{\theta})$ are given in Section D of the Appendix. A solution of the
system of equations in (6.3) is called MLqE estimator. It is not possible to
derive analytical solution for the MLqE $\widehat{\bm{\theta}}$.

###  6.2 Fisher Information based on $\log_{q}$

The definition of Fisher information (FI) based on $\log_{q}$ is given in [19,
28]. The elements of FI matrix are defined by

$\displaystyle[\mathcal{I}(\bm{\theta})]_{j,k}=\sum_{i=1}^{n}\mathbbm{E}\bigg{[}f^{1-q}(X_{i};\bm{\theta})\,{\partial\log
f(X_{i};\bm{\theta})\over\partial\theta_{j}}\,{\partial\log
f(X_{i};\bm{\theta})\over\partial\theta_{k}}\bigg{]},\quad\theta_{j},\theta_{k}\in\\{\mu,\sigma,\alpha,\rho,\delta\\},$

where $\log f(x;\bm{\theta})$ is given in (6.2) and $f(x;\bm{\theta})$ is the
parametric model in (2.1). When the inverse of FI matrix exists, it is well-
known that the diagonal elements of inverse of FI give ${\rm
Var}(\widehat{\bm{\theta}})$ based on logq [10, 19]. In general, there is no
closed form expression for the FI matrix (see Appendices D and F.6). It is
clear that when $q=1$ in the above identity, under standard regularity
conditions, we obtain the classical FI matrix [21, 24].

##  7 Application on real data sets

We present applications to illustrate the performance of the trimodal normal
model compared with smooth kernel distribution as semiparametric distribution
(SmoothKernelDistribution) and estimation of distribution which is performed
by using a function named as FindDistribution embedded into Mathematica 12.0
software to find an appropriate distribution for the data set. Since
Mathematica software are capable for performing the optimization, bootstrap
and also includes the numerical evaluation of hypergeometric function
(${}_{2}F_{1}$) while conducting the modelling the data sets, practitioners
can use these codes for their aims in researches. The supplementary metarials
provide the building codes for practitioners [37, 38, 39].

The Examples 1 and 2 represent the real data called as heterodatatrain$V5 and
heterodata$V4, respectively in the "Rmixmod" package at R software with
version 4.1.3. The numbers of sample size $n$ are $300$ and $200$ for Examples
1 and 2, respectively. The package "multilevel" includes data called as
"bh1996". The columns 11 and 13 of bh1996 data represent the modality. The
Examples 3 and 4 are for data sets with sample sizes $n=7382$. The $q$ values
for Examples 1-4 are chosen as $0.98,0.95,0.98,0.99$, respectively, according
to the GOFTs.

Table 2: The models, the estimates of parameters, statistics and information
criteria for assesment of models

Example 1
---
| Estimates | Goodness of fit tests (GOFTs) | log(L) & Information criteria
PM | $\hat{\mu}$ | $\hat{\sigma}$ | KS | CVM | AD | log(L) | AIC | BIC
qTDΦ | -2.25117 | 1.32990 | 0.0778466 | 0.551431 | 620.933 | -570.632 | 1155.25 | 1172.67
SK | -2.24486 | 1.62363 | 0.0680419 | 0.402510 | 593.949 | -563.837 | 1130.25 | 1137.49
EstD | -2.22717 | 1.61703 | 0.0706750 | 0.410858 | 597.038 | -563.777 | 1134.51 | 1141.45
TDΦ | -2.25803 | 1.38998 | 0.0764954 | 0.559258 | 615.978 | -570.669 | 1147.46 | 1166.25
Example 2
PM | $\hat{\mu}$ | $\hat{\sigma}$ | KS | CVM | AD | log(L) | AIC | BIC
qTDΦ | -1.94546 | 1.29762 | 0.120948 | 0.988454 | 471.754 | -391.988 | 794.657 | 814.271
SK | -1.98323 | 1.63448 | 0.0862193 | 0.420766 | 400.882 | -365.779 | 740.391 | 744.768
EstD | -1.96097 | 1.65964 | 0.0873551 | 0.43558 | 401.620 | -370.828 | 743.939 | 752.477
TDΦ | -1.91852 | 1.47189 | 0.120597 | 0.927997 | 430.148 | -377.753 | 772.283 | 787.288
Example 3
PM | $\hat{\mu}$ | $\hat{\sigma}$ | KS | CVM | AD | log(L) | AIC | BIC
qTDΦ | 2.77324 | 0.856429 | 0.0350506 | 1.66183 | 14510.4 | -9762.83 | 19535.7 | 19570.2
SK | 2.78055 | 0.915343 | 0.0150676 | 0.293907 | 14661.0 | -9659.05 | 19322.1 | 19335.9
EstD | 2.78047 | 0.909687 | 0.0326443 | 1.31802 | 14397.0 | -9764.47 | 19532.9 | 19546.8
TDΦ | 2.77505 | 0.882003 | 0.0333040 | 1.52912 | 14414.5 | -9761.59 | 19533.2 | 19567.7
Example 4
PM | $\hat{\mu}$ | $\hat{\sigma}$ | KS | CVM | AD | log(L) | AIC | BIC
qTDΦ | -0.001341870 | 0.865658 | 0.989796 | 43.8267 | 677.270 | -9550.79 | 19111.6 | 19146.1
SK | -0.000056456 | 0.889187 | 0.981952 | 43.0978 | 600.388 | -9472.98 | 18950.0 | 18963.8
EstD | 0.000581656 | 0.884030 | 0.988856 | 43.7301 | 663.308 | -9553.08 | 19110.2 | 19124.0
TDΦ | -0.002632980 | 0.869382 | 0.989468 | 43.7947 | 672.560 | -9550.61 | 19111.2 | 19145.8

  * •

PM: (Semi)-parametric models

  * •

qTDΦ: Objective function $\log_{q}$ from trimodal normal distribution.

  * •

SK: The smooth kernel distribution used Gaussian (normal) distribution
(semiparametric model).

  * •

EstD: The automatically chosen function by "FindDistribution" in Mathematica
software.

  * •

TDΦ: Objective function $\log$ from trimodal normal distribution.

  * •

Italic represents closeness to the values produced by SK and EstD or almost
best ones

The location ($\mu$) and scale ($\sigma$) are important parameters to
summarize the data set. The efficient estimations of these parameters depend
on the chosen function used for modelling. Table 2 provides them and other
statistics for testing the modelling competence of the used functions. When
all of PM (qTDΦ, SK, EstD and TDΦ ) in Table 2 are compared, SK and EstD are
rival ones among models. On the other side, the statistics and information
criteria of qTDΦ and TDΦ can be near to SK and EstD as alternative models (see
also discussion in Appendix F.3).

Table 3: Standard errors of estimators and estimates of shape and bimodality
parameters in qTDΦ and TDΦ

Example 1
---
PM | $\sqrt{{\rm Var}(\hat{\mu})}$ | $\sqrt{{\rm Var}(\hat{\sigma})}$ |  $\hat{\alpha}$ ($\sqrt{{\rm Var}(\hat{\alpha})}$) |  $\hat{\rho}$ ($\sqrt{{\rm Var}(\hat{\rho})}$) |  $\hat{\delta}$ ($\sqrt{{\rm Var}(\hat{\delta})}$)
qTDΦ | 0.152203 | 0.124745 | 1.22607(0.253034) | 1.06657(7.41619) | 0.233824(1.63256)
TDΦ | 0.167649 | 0.143923 | 1.21113(0.257971) | 1.08073(7.41113) | 0.230591(1.58755)
Example 2
PM | $\sqrt{{\rm Var}(\hat{\mu})}$ | $\sqrt{{\rm Var}(\hat{\sigma})}$ |  $\hat{\alpha}$ ($\sqrt{{\rm Var}(\hat{\alpha})}$) |  $\hat{\rho}$ ($\sqrt{{\rm Var}(\hat{\rho})}$) |  $\hat{\delta}$ ($\sqrt{{\rm Var}(\hat{\delta})}$)
qTDΦ | 0.071482 | 0.103308 | 37.1884(205.348) | 5.53004(3048.8) | 50.6147(27963.2)
TDΦ | 0.0985512 | 0.140166 | 31.9374(217.858) | 11.6716(5270.27) | 48.0121(21714.5)
Example 3
PM | $\sqrt{{\rm Var}(\hat{\mu})}$ | $\sqrt{{\rm Var}(\hat{\sigma})}$ |  $\hat{\alpha}$ ($\sqrt{{\rm Var}(\hat{\alpha})}$) |  $\hat{\rho}$ ($\sqrt{{\rm Var}(\hat{\rho})}$) |  $\hat{\delta}$ ($\sqrt{{\rm Var}(\hat{\delta})}$)
qTDΦ | 0.000889860 | 0.000992885 | 1.161980(0.00549161) | 1.39357(0.0514929) | 0.400877(0.0145018)
TDΦ | 0.000932255 | 0.000732616 | 0.888664(0.00562133) | 1.17315(0.0706782) | 0.182165(0.0109832)
Example 4
PM | $\sqrt{{\rm Var}(\hat{\mu})}$ | $\sqrt{{\rm Var}(\hat{\sigma})}$ |  $\hat{\alpha}$ ($\sqrt{{\rm Var}(\hat{\alpha})}$) |  $\hat{\rho}$ ($\sqrt{{\rm Var}(\hat{\rho})}$) |  $\hat{\delta}$ ($\sqrt{{\rm Var}(\hat{\delta})}$)
qTDΦ | 0.0002215930 | 0.000690205 | 0.499788(0.00379787) | 1.405910(0.00379787) | 0.12482400(0.00240774)
TDΦ | 0.0000789721 | 0.000723072 | 0.509288(0.00275330) | 0.346905(0.00799665) | 0.0312284(0.000700832)

Table 3 shows the square root of variance of estimators and also the estimates
of $\hat{\alpha}$, $\hat{\rho}$ and $\hat{\delta}$ for the examples analyzed.

Table 4: The estimates for $\hat{\mu}$ and $\hat{\sigma}$ from normal
distribution and robust form

|  | Example 1 | Example 2 | Example 3 | Example 4
---|---|---|---|---|---
N | $\hat{\mu}$ | -2.24483 | -1.96589 | 2.77989 | 1.24627 × 10-12
$\hat{\sigma}$ | 1.58461 | 1.58535 | 0.908344 | 0.882701
R | Median | -2.34344 | -2.23693 | 2.83333 | 0.0309632
MAD | 1.08095 | 0.981103 | 0.611111 | 0.607600

  * •

N: Normal distribution

  * •

R: Robust statistics

  * •

M: Median, MAD: Median absolute deviation (Median(|x-Median(x)|))

Table 4 introduces the basic statistics to see the role of distribution with
one mode property and trimodality.

The estimates of $\hat{\mu}$ and $\hat{\sigma}$ from different PM show that we
can have a clue to imply that the existence of modality can be observed,
because the estimates of $\hat{\mu}$ and $\hat{\sigma}$ from SK as a smooth
kernel technique based on working on the data-adaptive approach (which is
capable to fit the modality whether or not it exists in the reality–see also
discussions on Appendix E.2) instead of parametric approach for modelling can
be close to the estimates of $\hat{\mu}$ and $\hat{\sigma}$ from qTDΦ and TDΦ.
The kernel estimation method as a smoothing tehnique is the best one
generally. Even if 1000 replication for the different design of samplings
constructed by use of bootstrap technique is applied, the numerical error(s)
in computation for optimization can be tricker to consider and make an
accurate judgement among the modelling performance of the used four models.
For example, CDF and PDF of TDΦ depend on the Hypergeometric2F1 in
Mathematica. Soft forms of PDF of TDΦ in Figure 2(b) and smooth kernel
technique in Ref. [40] can be alternative to each other when the estimates of
$\hat{\mu}$ and $\hat{\sigma}$, the statistics from GOFTs, the values of
log(L) and IC are taken into account. On the other side, it is very difficult
to know which function will be the best one for modelling when the data sets
in the finite sample size are tried to be fitted by the functions. Even if the
sample sizes of Examples 3 and 4 are 7382, eventually we have finite sample
size whatever it is. The population in reality will not known exactly.
Consequently, an alternative function can be necessary for driving modality
via parameters $\rho$ and $\delta$. This is the reason why we make a
comparison between the SK and other parametric models to observe what and how
the estimates of $\hat{\mu}$ and $\hat{\sigma}$ will be changed if PDF are
changed [33, 2]. Note that the modelling and numerical error(s) are topics
which can affect each other.

##  8 Conclusions

Since recent times show that an increasing popularity has been observed in the
modelling for data sets having modality, producing the trimodal form of any
PDF has been proposed. The trimodal form is constructed by using the technique
which includes the cumulative function of Maxwell distribution, the existing
unimodal distribution and the corresponding normalizing constant of the
proposed distribution. The properties of trimodality have been examined. The
application of producing the trimodality has been conducted for the normal
distribution which is symmetric and unimodal form with two parameters which
are location and scale. The properties of trimodal normal distribution have
been examined. Thus, the applicability of this distribution have been tested.
The trimodal normal distribution can have different forms such as strict and
soft modalities to perform a precise fitting when there exists three modes in
the empirical distribution of the data sets.

A comparison among trimodal normal, the kernel type estimation method, the
probable parametric distribution driven by Mathematica software has been
performed in order to make applications for numerical evaluation of TDΦ. The
$\log_{q}$-likelihood estimation method and its special form with $q\to 1$
have been used to estimate of parameters of trimodal normal distribution. The
proofs, properties of TDΦ distribution and codes used for application have
been given by appendices if the researchers perform to model the data sets by
use of TDΦ distribution.

The future will be an application on the different areas of statistics such as
regression modelling, the tools in the multivariate statistics and other tools
based on the distribution theory. The order statistic form of TDΦ in the least
informative distribution will be studied for the trimodal forms of the
existing distributions in the applied field of science. Additionally, the
precise modelling for inliers into data sets can also be performed by use of
trimodality, the generalized logarithms, entropy functions, order statistic
and different estimation methods all together [35, 36]. A package in R
software will be prepared for practitioners after the special function in R
software are improved.

#### Acknowledgements

We acknowledge the anonymous referees for their helpful comments, suggestions
and references provided in their reports. R. V. thanks A. V. Medino, J. Roldan
and E. M. M. Ortega for partial discussions of Theorem 3.4 and for general
paper questions.

#### Disclosure statement

There are no conflicts of interest to disclose.

#### Funding

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

#### ORCID

Roberto Vila https://orcid.org/0000-0003-1073-0114
Victor Serra https://orcid.org/0000-0003-3061-2094
Mehmet N. Çankaya https://orcid.org/0000-0002-2933-857X
Felipe Quintino https://orcid.org/0000-0003-0286-0541

## References

  * [1] M. N. Çankaya, Y. M. Bulut, F. Z. Doğru, O. Arslan, _A bimodal extension of the generalized gamma distribution_ , Revista Colombiana de Estadística (2015), 38:371–384.
  * [2] Çankaya, M. N., _On the estimating equations and objective functions for parameters of exponential power distribution: Application for disorder_ , arXiv e-prints, (2021), arXiv-2102.
  * [3] F. Z. Doğru, Y. M. Bulut, O. Arslan, _Doubly reweighted estimators for the parameters of the multivariate t-distribution_ , Communications in Statistics-Theory and Methods, 47(19), (2018), pp. 4751-4771.
  * [4] B. Everitt, (2013). Finite mixture distributions. Springer Science & Business Media.
  * [5] E. Gómez-Déniz, J. M. Sarabia, E. Calderín-Ojeda, _Bimodal normal distribution: Extensions and applications_ , Journal of Computational and Applied Mathematics, 388, (2021), 113292.
  * [6] R. Vila, H. Saulo, J. Roldan, _On some properties of the bimodal normal distribution and its bivariate version_ , Chilean Journal of Statistics, 12(2), (2021), pp. 125-144.
  * [7] C. Lee, F. Famoye, A. Y. Alzaatreh Methods for generating families of univariate continuous distributions in the recent decades. Wiley Interdisciplinary Reviews: Computational Statistics, 5(3), (2013), pp. 219-238.
  * [8] R. Vila, L. Ferreira, H. Saulo, F. Prataviera, E. M. M. Ortega, _A bimodal gamma distribution: Properties, regression model and applications_ , Statistics, 54 (2020), pp. 469-493.
  * [9] R. Vila, M. N. Çankaya, _A Bimodal Weibull Distribution: Properties and Inference_ , Journal of Applied Statistics, (2021), pp. 1-19.
  * [10] J. F. Bercher, _Some properties of generalized Fisher information in the context of nonextensive thermostatistics_ , Physica A: Statistical Mechanics and its Applications, 392(15), (2013), pp. 3140-3154.
  * [11] J. F. Bercher, _A simple probabilistic construction yielding generalized entropies and divergences, escort distributions and q-Gaussians_ , Physica A: Statistical Mechanics and its Applications, 391(19), (2012), pp. 4460-4469.
  * [12] F. Domma, B.V. Popović, and S. Nadarajah, _An extension of Azzalini’s method_ , J. Comput. Appl. Math. 278 (2015), pp. 37-47.
  * [13] M.N. Çankaya, _Asymmetric bimodal exponential power distribution on the real line_ , Entropy, 20(1), (2018), 23.
  * [14] F. Domma, F. Condino, and B.V. Popović, _A new generalized weighted Weibull distribution with decreasing, increasing, upside-down bathtub, N-shape and M-shape hazard rate_ , J. Appl. Stat. 44 (2017), pp. 2978-2993.
  * [15] M. Rahman, B. Al-Zahrani, M.Q. Shahbaz, Cubic transmuted Pareto distribution. Annals of Data Science, 7(1), (2020), pp. 91-108.
  * [16] D. Elal-Olivero, _Alpha-skew-normal distribution_ , Proyecciones (Antofagasta) 29 (2010), pp. 224-240.
  * [17] R. Vila, L. Alfaia, A. F. Menezes, M.N. Çankaya, M. Bourguignon, _A Model for Bimodal Rates and Proportions_ , (2021), arXiv preprint arXiv:2108.07934.
  * [18] R.C. Dunbar, _Deriving the Maxwell distribution_ , Journal of Chemical Education, 59(1), (1982), 22.
  * [19] M.N. Çankaya, J. Korbel, _Least informative distributions in maximum q-log-likelihood estimation_ , Physica A 509 (2018), pp. 140-150.
  * [20] D. Ferrari, Y. Yang, _Maximum Lq-likelihood estimation_ , Ann. Stat. 38 (2010), pp. 753-783.
  * [21] R.A. Fisher, Theory of Statistical Estimation. _In Mathematical Proceedings of the Cambridge Philosophical Society_ , (Vol. 22, No. 5, pp. 700-725). Cambridge University Press, 1925.
  * [22] I.S. Gradshteyn, I.M. Ryzhik, Table of Integrals, Series and Products, Academic Press, San Diego, 2000.
  * [23] S. Klugman, H. Panjer, G. Willmot, _Loss models: From data to decisions_ , Wiley, New York, 1998.
  * [24] E. L. Lehmann, G. Casella, _Theory of Point Estimation, Wadsworth & Brooks/Cole_, Pacific Grove, CA, 589, USA, 1998.
  * [25] E. W. Ng, M. Geller, _A table of integrals of the Error functions_ , Journal of Research of the National Bureau of Standards Section B Mathematical Sciences, 73B(1), 1 - January 1969.
  * [26] C. E. G. Otiniano, R. Vila, P. C. Brom, M. Bourguignon, _On the bimodal Gumbel model with application to environmental data_ , To appear in Austrian Journal of Statistics, (2021).
  * [27] J. Pender, _The truncated normal distribution: Applications to queues with impatient customers_ , Operations Research Letters, 43 (2015), pp. 40-45.
  * [28] A. Plastino, A. R. Plastino, H.G. Miller, _Tsallis nonextensive thermostatistics and Fisher’s information measure_ , Physica A 235 (1997), pp. 577-588.
  * [29] A. P. Prudnikov, IU. A. Brychkov, O. I. Marichev, _Integrals and series. Vol 2, Special functions_ , London: Taylor & Francis, 2002.
  * [30] C. E. Shannon, _A mathematical theory of communication_ , Bell Labs. Tech. J. 27 (1948), pp. 623-656.
  * [31] C. Tsallis, _Possible generalization of Boltzmann-Gibbs statistics_ , J. Stat. Phys. 52 (1988), pp. 479-487.
  * [32] C. Tsallis, _Introduction to Nonextensive Statistical Mechanics: Approaching a Complex World_ , Springer, New York, 2009.
  * [33] P.J. Huber, _Robust estimation of a location parameter_ , Ann. Math. Stat. (1964), 35.
  * [34] Çankaya, M.N., Arslan, O., _On the robustness properties for maximum likelihood estimators of parameters in exponential power and generalized T distributions_ , Communications in Statistics-Theory and Methods, 49(3), (2020), pp. 607-630.
  * [35] M.N. Çankaya, _M-estimations of shape and scale parameters by order statistics in least informative distributions on q-deformed logarithm_ , Journal of the Institute of Science and Technology, (2020), 10(3), 1984-1996.
  * [36] M.N. Çankaya, _Derivatives by ratio principle for $q$-sets on the time scale calculus_ , Fractals, 29(8) (2021), 2140040, DOI: 10.1142/S0218348X21400405.
  * [37] León, J. G. S. (2017). Mathematica$\circledR$ Beyond Mathematics: The Wolfram Language in the Real World. Chapman and Hall/CRC.
  * [38] Baglivo, J. A. (2005). Mathematica laboratories for mathematical statistics: Emphasizing simulation and computer intensive methods. Society for industrial and applied mathematics.
  * [39] Feagin, J. F., Quantum methods with Mathematica$\circledR$, (2002), Springer Science & Business Media.
  * [40] Härdle, W., Müller, M., Sperlich, S., Werwatz, A. (2004). Nonparametric and semiparametric models (Vol. 1). Berlin: Springer.
  * [41] Härdle, W. K. (1991). Smoothing techniques: with implementation in S. Springer Science & Business Media.
  * [42] Somani, A. K., Deka, G. C. (Eds.). Big data analytics: Tools and technology for effective planning. (2017), CRC Press.

##  Appendix A Proof of some results of Sections 3 and 4

###### Proof of Lemma 3.1.

Let $Y$ be a random variable with ${\rm Gamma}(p,1/\alpha^{2})$ distribution.
The corresponding PDF and CDF of $Y$, respectively, are given by
$f(y;\alpha,p)={(1/\alpha^{2})^{p}}\,y^{p-1}{\rm e}^{-y/\alpha^{2}}/\Gamma(p)$
and $F(y;\alpha,p)={\gamma(p,y/\alpha^{2})/\Gamma(p)}$. Therefore,

$\displaystyle\mathfrak{R}(y)={2\delta}f(y;\alpha,p)-[\rho+\delta
F(y;\alpha,p)]\mathfrak{h}(y),\quad y>0.$ (A.1)

Suppose $y\geqslant\alpha^{2}(1-A)$. Then

$\displaystyle{\rm
e}^{(\alpha^{2}A+y)/\alpha^{2}}=\sum_{k=0}^{\infty}{(\alpha^{2}A+y)^{k}/\alpha^{2k}\over
k!}\geqslant{(\alpha^{2}A+y)^{n+1}/\alpha^{2(n+1)}\over(n+1)!},$

for each $n\in\mathbbm{N}$, such that

$\displaystyle{\rm e}^{-y/\alpha^{2}}\leqslant{(n+1)!\alpha^{2(n+1)}{\rm
e}^{A}\over(\alpha^{2}A+y)^{n+1}}.$

Hence,

$\displaystyle{2\delta}f(y;\alpha,p)$
$\displaystyle\leqslant{2\delta}\biggl{[}{(1/\alpha^{2})^{p}\over\Gamma(p)}\,y^{p-1}\,{(n+1)!\alpha^{2(n+1)}{\rm
e}^{A}\over(\alpha^{2}A+y)^{n+1}}\biggr{]}$
$\displaystyle\leqslant{2\delta}\biggl{[}{(1/\alpha^{2})^{p}\over\Gamma(p)}\,{(n+1)!\alpha^{2(n+1)}{\rm
e}^{A}\over(\alpha^{2}A+y)^{n+1-p}}\biggr{]}\leqslant{\rho
C\over(\alpha^{2}A+y)^{\beta-p}}=\rho\mathfrak{h}(y),$ (A.2)

for $\beta<n+1$, $C\geqslant 1$ and some $\rho>0$, so that

$\displaystyle\rho\geqslant{2\delta}\biggl{[}{(1/\alpha^{2})^{p}\over\Gamma(p)}\,{(n+1)!\alpha^{2(n+1)}{\rm
e}^{A}}\biggr{]}.$ (A.3)

Since $F(y;\alpha,p)$ is a CDF, from (A.2) it follows that

$\displaystyle{2\delta}f(y;\alpha,p)\leqslant[\rho+\delta
F(y;\alpha,p)]\mathfrak{h}(y)\leqslant(\rho+\delta)\mathfrak{h}(y),\quad\forall
y\geqslant\alpha^{2}(1-A).$ (A.4)

In other words, for every $y\geqslant\alpha^{2}(1-A)$ and $\rho$ large enough,
the tail of Gamma distribution $2\delta f(y;\alpha,p)$ is lighter than the
tail of $[\rho+\delta F(y;\alpha,p)]\mathfrak{h}(y)$ which decays
polynomially.

On the other hand, it is well-known that: when $p<1$, the Gamma distribution
is exponentially shaped and asymptotic to both the vertical and horizontal
axes; the Gamma distribution with shape parameter $p=1$ and scale parameter
$1/\alpha^{2}$ is the same as an exponential distribution of scale parameter
(or mean) $1$; when $p$ is greater than one, the Gamma distribution assumes a
maximum value (unimodal), but skewed shape. The skewness reduces as the value
of $p$ increases.

Based on the shapes of the Gamma distribution and on the inequality (A.4), in
Figure 3, we graphically sketch the functions $2\delta f(y;\alpha,p)$ and
$[\rho+\delta F(y;\alpha,p)]\mathfrak{h}(y)$, and consider all possible cases,
by varying the parameters $\alpha$, $\rho$ as in (A.3), $\delta$ and $p$
(known), in which the graphs of $2\delta f(y;\alpha,p)$ and $[\rho+\delta
F(y;\alpha,p)]\mathfrak{h}(y)$ intersect (or not).

Figure 3: The graphs of $2\delta f(y;\alpha,p)$ and $[\rho+\delta
F(y;\alpha,p)]\mathfrak{h}(y)$ have at most two common points.

Then, in both cases $p\leqslant 1$ or $p>1$, the function $\mathfrak{R}$ in
(A.1) has at most two real zeros. This completes the proof. ∎

###### Proof of Proposition 3.9.

A simple algebric manipulation shows that

$\displaystyle S_{q}(X)=\dfrac{1}{q-1}\,\bigg{[}1-\int_{\sigma
D+\mu}f^{q}(x;\bm{\theta})\,{\rm d}x\bigg{]},\quad q\neq 1.$

Since $0<T(x;\alpha,p)<1$ for almost all $x\in D$, we have

$\displaystyle 0\leqslant f(x;\bm{\theta})\leqslant{(\rho+\delta)\over
Z_{\bm{\theta}}}\,g\biggl{(}{x-\mu\over\sigma}\biggr{)}.$

Then, by using that the function $x\mapsto x^{q}$ is increasing, for $x>0$ and
$q>0$, we have

$\displaystyle 0\leqslant f^{q}(x;\bm{\theta})\leqslant{(\rho+\delta)^{q}\over
Z^{q}_{\bm{\theta}}}\,g^{q}\biggl{(}{x-\mu\over\sigma}\biggr{)},$

for all $x\in\sigma D+\mu$. Consequently,

$\displaystyle\int_{\sigma D+\mu}f^{q}(x;\bm{\theta})\,{\rm
d}x\leqslant{(\rho+\delta)^{q}\over
Z^{q}_{\bm{\theta}}}\,\int_{D}g^{q}\biggl{(}{x-\mu\over\sigma}\biggr{)}\,{\rm
d}x={(\rho+\delta)^{q}\sigma\over
Z^{q}_{\bm{\theta}}}\,\big{[}1-(q-1)S_{q}(W)\big{]},\quad q\neq 1.$

Hence, if $S_{q}(W)$ exists then of $S_{q}(X)$ also exists. ∎

###### Proof of Proposition 3.11.

By taking the logarithm of each side of (2.1), by definition of Shannon
entropy, we have

$\displaystyle
S_{1}(X)=\log(Z_{\bm{\theta}})+\mathbbm{E}\biggl{[}-\log\biggr{(}\rho+\delta
T\biggl{(}{X-\mu\over\sigma};\alpha,p\biggr{)}\biggr{)}\biggr{]}+\mathbbm{E}\biggl{[}-\log
g\biggl{(}{X-\mu\over\sigma}\biggr{)}\biggr{]}.$ (A.5)

By taking $L(x)=\log(\rho+\delta T(x;\alpha,p))$ and $L(x)=-\log g(x)$,
$\forall x\in\sigma D+\mu$, in Corollary 3.5, respectively, we have

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\mathbbm{E}\biggl{[}-\log\biggr{(}\rho$
$\displaystyle+\delta
T\biggl{(}{X-\mu\over\sigma};\alpha,p\biggr{)}\biggr{)}\biggr{]}=-{{({\rho}+\delta)\sigma\over
Z_{\bm{\theta}}}\,\mathbbm{E}\big{[}\log(\rho+\delta T(W;\alpha,p))\big{]}}$
$\displaystyle-{{{\delta}\sigma\over
Z_{\bm{\theta}}}\,\left\\{\mathbbm{E}\left[\mathds{1}_{\\{W\leqslant-\alpha\sqrt{Y}\\}}\log(\rho+\delta
T(W;\alpha,p))\right]-\mathbbm{E}\left[\mathds{1}_{\\{W\leqslant\alpha\sqrt{Y}\\}}\log(\rho+\delta
T(W;\alpha,p))\right]\right\\}}$ (A.6)

and

$\displaystyle\\!\\!\\!\\!\mathbbm{E}\biggl{[}-\log
g\biggl{(}{X-\mu\over\sigma}\biggr{)}\biggr{]}={{({\rho}+\delta)\sigma\over
Z_{\bm{\theta}}}\,S_{1}(W)+{{\delta}\sigma\over
Z_{\bm{\theta}}}\,\left\\{{\mathbbm{E}}\left[S_{1}\big{(}\mathds{1}_{\\{W\leqslant-\alpha\sqrt{Y}\\}}W\big{)}\right]-{\mathbbm{E}}\left[S_{1}\big{(}\mathds{1}_{\\{W\leqslant\alpha\sqrt{Y}\\}}W\big{)}\right]\right\\}}.$
(A.7)

By replacing the identities (A) and (A.7) in (A.5), the proof follows. ∎

###### Proof of Lemma 4.1.

The PDF and CDF of $Y\sim{\rm Gamma}(p,1/\alpha^{2})$, respectively, are given
by $f(y;\alpha,p)={{(1/\alpha^{2})^{p-1}{\rm e}^{-y/\alpha^{2}}/\Gamma(p)}}$
and $F(y;\alpha,p)={\gamma(p,y/\alpha^{2})/\Gamma(p)}$. Then, $R$ takes on the
following form

$\displaystyle{R}(y)={2\delta}f(y;\alpha,p)-[\rho+\delta F(y;\alpha,p)].$
(A.8)

Figure 4:  The graphs of $2\delta f(y;\alpha,p)$ and $\rho+\delta
F(y;\alpha,p)$ have at most two points of intersection.

Based on the shapes of the Gamma law, in Figure 4, we graphically sketch the
functions $2\delta f(y;\alpha,p)$ and $\rho+\delta F(y;\alpha,p)$, and
consider all possible cases, by varying the parameters $\alpha$, $\rho$,
$\delta$ and $p$ (known), in which the graphs of $2\delta f(y;\alpha,p)$ and
$\rho+\delta F(y;\alpha,p)$ intersect (or not).

So, for the cases when $p\leqslant 1$ or $p>1$, the function $R$ in (A.8) has
at most two real zeros. This completes the proof. ∎

##  Appendix B Moments in the Gaussian case

In this section, for simplicity, we consider $p\geqslant 1$ an integer. To get
the moments of $X\sim{\rm TD}_{\Phi}(\bm{\theta})$, by Corollary 3.7, it is
necessary to determine
$\mathbbm{E}(\mathds{1}_{\\{Z\leqslant\pm\alpha\sqrt{Y}\\}}Z^{n})$,
$n\geqslant 1$, where $Z\sim N(0,1)$ and $Y\sim{\rm Gamma}(p,1)$ are
independent. Indeed, since $Z$ and $Y$ are independent, by definition of
expectation,

$\displaystyle\mathbbm{E}\big{(}\mathds{1}_{\\{Z\leqslant\pm\alpha\sqrt{Y}\\}}Z^{n}\big{)}=\int_{0}^{\infty}\mathbbm{E}\big{(}\mathds{1}_{\\{Z\leqslant\pm\alpha\sqrt{y}\\}}Z^{n}\big{)}\,{y^{p-1}{\rm
e}^{-y}\over\Gamma(p)}\,{\rm d}y.$ (B.1)

Let us define $n_{m}=n-2m-1$,

$\displaystyle c_{n,m}={n!2^{-m}n_{m}!\over m!(n_{m}+1)!}\quad\text{and}\quad
d_{n_{m},j}={(-1)^{j}\over j!(n_{m}-2j)!2^{j}}.$ (B.2)

By using the following known formula (see Theorem 2.3. of [27]): for each real
numbers $a,b$ such that $a<b$,

$\displaystyle\mathbbm{E}\big{(}\mathds{1}_{\\{a\leqslant Z\leqslant
b\\}}Z^{n}\big{)}=\sum_{m=0}^{\lfloor n/2\rfloor}\sum_{j=0}^{\lfloor
n_{m}/2\rfloor}c_{n,m}d_{n_{m},j}\,\big{[}a^{n_{m}-2j}\phi(a)-b^{n_{m}-2j}\phi(b)\big{]},$

where ${\lfloor x\rfloor}$ is the greatest integer less than or equal to $x$,
we have

$\displaystyle\mathbbm{E}\big{(}\mathds{1}_{\\{Z\leqslant\pm\alpha\sqrt{y}\\}}Z^{n}\big{)}=-\sum_{m=0}^{\lfloor
n/2\rfloor}\sum_{j=0}^{\lfloor
n_{m}/2\rfloor}c_{n,m}d_{n_{m},j}\,(\pm\alpha\sqrt{y})^{n_{m}-2j}\phi(\pm\alpha\sqrt{y}).$

Replacing the above identity in (B.1),

$\displaystyle\mathbbm{E}\big{(}\mathds{1}_{\\{Z\leqslant\pm\alpha\sqrt{Y}\\}}Z^{n}\big{)}=-\sum_{m=0}^{\lfloor
n/2\rfloor}\sum_{j=0}^{\lfloor
n_{m}/2\rfloor}c_{n,m}d_{n_{m},j}\,\int_{0}^{\infty}(\pm\alpha\sqrt{y})^{n_{m}-2j}\phi(\pm\alpha\sqrt{y})\,{y^{p-1}{\rm
e}^{-y}\over\Gamma(p)}\,{\rm d}y.$

But, by definition of gamma distribution,

$\displaystyle\int_{0}^{\infty}(\pm\alpha\sqrt{y})^{n_{m}-2j}\phi(\pm\alpha\sqrt{y})\,{y^{p-1}{\rm
e}^{-y}\over\Gamma(p)}\,{\rm
d}y={(\pm\alpha)^{n_{m}-2j}\over\sqrt{2\pi}\Gamma(p)(1+{\alpha^{2}\over
2})^{{n_{m}-2j\over 2}+p}}\,\Gamma\biggl{(}{n_{m}-2j\over 2}+p\biggr{)}.$

Hence

$\displaystyle\mathbbm{E}\big{(}\mathds{1}_{\\{Z\leqslant\pm\alpha\sqrt{Y}\\}}Z^{n}\big{)}=-{1\over\sqrt{2\pi}\Gamma(p)}\sum_{m=0}^{\lfloor
n/2\rfloor}\sum_{j=0}^{\lfloor
n_{m}/2\rfloor}{c_{n,m}d_{n_{m},j}\,(\pm\alpha)^{n_{m}-2j}\over(1+{\alpha^{2}\over
2})^{{n_{m}-2j\over 2}+p}}\,\Gamma\biggl{(}{n_{m}-2j\over 2}+p\biggr{)}.$
(B.3)

By using (B.3) and the known formula

$\displaystyle\mathbbm{E}(Z^{n})={2^{-n/2}n!\over(n/2)!}\,\mathds{1}_{\\{n\,{\rm
even}\\}},$

in Corollary 3.7, we get the formula of moments for $X\sim{\rm
TD}_{\Phi}(\bm{\theta})$ given in Subsection 4.4.

##  Appendix C Entropy in the Gaussian case

By Remark 3.10, the Shannon entropy $S_{1}(X)$ of $X\sim{\rm
TD}_{\Phi}(\bm{\theta})$ exists whenever $S_{1}(Z)$ of $Z\sim N(0,1)$ also
exists. Since $S_{1}(Z)=\log(\sqrt{2\pi})+{1/2}$, the existence of $S_{1}(X)$
is guaranteed. So it makes sense to find a closed expression for $S_{1}(X)$.
To find this expression, by Proposition 3.11, it is enough to determine the
expectations:
$\mathbbm{E}[\mathds{1}_{\\{Z\leqslant\pm\alpha\sqrt{Y}\\}}\log(\rho+\delta
T(Z;\alpha,p))]$ and
$\mathbbm{E}[S_{1}(\mathds{1}_{\\{Z\leqslant\pm\alpha\sqrt{Y}\\}}Z)]$, where
$T$ is as in (2.2) and $Z\sim N(0,1)$ and $Y\sim{\rm Gamma}(p,1)$ are
independent. As in the previous subsection, for simplicity, in this subsection
we consider $p\geqslant 1$ an integer.

A simple calculation shows that, for each real numbers $a,b$ such that $a<b$,

$\displaystyle S_{1}(\mathds{1}_{\\{a\leqslant Z\leqslant b\\}}Z)$
$\displaystyle=-\int_{a}^{b}\phi(z)\log\phi(z)\,{\rm d}z$
$\displaystyle={1\over 2}\,\big{[}a\phi(a)-b\phi(b)\big{]}+{1\over
2}\,\biggl{[}\log(\sqrt{2\pi})+{1\over 2}\biggr{]}\left[{\rm
erf}\biggl{(}{b\over\sqrt{2}}\biggr{)}-{\rm
erf}\biggl{(}{a\over\sqrt{2}}\biggr{)}\right].$

Then

$\displaystyle\mathbbm{E}\big{[}S_{1}(\mathds{1}_{\\{Z\leqslant\pm\alpha\sqrt{Y}\\}}Z)\big{]}$
$\displaystyle=\mp{1\over
2}\int_{0}^{\infty}\alpha\sqrt{y}\phi(\pm\alpha\sqrt{y})\,{y^{p-1}{\rm
e}^{-y}\over\Gamma(p)}\,{\rm d}y$ $\displaystyle+{1\over
2}\biggl{[}\log(\sqrt{2\pi})+{1\over 2}\biggr{]}\int_{0}^{\infty}\left[{\rm
erf}\biggl{(}{\pm\alpha\sqrt{y}\over\sqrt{2}}\biggr{)}+1\right]{y^{p-1}{\rm
e}^{-y}\over\Gamma(p)}\,{\rm d}y.$

But,

  * •

by definition of Gamma distribution,

$\displaystyle\int_{0}^{\infty}\alpha\sqrt{y}\phi(\pm\alpha\sqrt{y})\,{y^{p-1}{\rm
e}^{-y}\over\Gamma(p)}\,{\rm d}y={\alpha\Gamma(p+{1\over
2})\over\sqrt{2\pi}\Gamma(p)}\,{1\over(1+{\alpha^{2}\over 2})^{p+{1\over
2}}},$

  * •

and by using (4.4),

$\displaystyle\int_{0}^{\infty}\left[{\rm
erf}\biggl{(}{\pm\alpha\sqrt{y}\over\sqrt{2}}\biggr{)}+1\right]{y^{p-1}{\rm
e}^{-y}\over\Gamma(p)}\,{\rm
d}y=2{\mathbbm{E}}\big{[}\Phi(\pm\alpha\sqrt{Y})\big{]}.$

Hence

$\displaystyle\mathbbm{E}\big{[}S_{1}(\mathds{1}_{\\{Z\leqslant\pm\alpha\sqrt{Y}\\}}Z)\big{]}=\mp{\alpha\Gamma(p+{1\over
2})\over 2\sqrt{2\pi}\Gamma(p)}\,{1\over(1+{\alpha^{2}\over 2})^{p+{1\over
2}}}+\biggl{[}\log(\sqrt{2\pi})+{1\over
2}\biggr{]}{\mathbbm{E}}\big{[}\Phi(\pm\alpha\sqrt{Y})\big{]}.$ (C.1)

The proof of the following result is technical.

######  Proposition C.1.

If $Z\sim N(0,1)$ the following hold

$\displaystyle\mathbbm{E}\big{[}\log(\rho+\delta
T(Z;\alpha,p))\big{]}=2\sum_{k=0}^{\infty}{1\over(2k+1)}\biggl{[}1+{1\over(\rho+1)^{2k+1}}\sum_{i=p}^{\infty}\widetilde{c}_{i,k}\,{2^{-i}(2i)!\over
i!}\biggr{]};$ (C.2)
$\displaystyle\mathbbm{E}\big{[}\mathds{1}_{\\{Z\leqslant\pm\alpha\sqrt{Y}\\}}\log(\rho+\delta
T(Z;\alpha,p))\big{]}$ $\displaystyle\hskip
76.82234pt=2\sum_{k=0}^{\infty}{1\over(2k+1)}\biggl{\\{}\mathbbm{E}\big{[}\Phi(\pm\alpha\sqrt{Y})\big{]}+{1\over(\rho+1)^{2k+1}}\sum_{i=p}^{\infty}\widetilde{c}_{i,k}\mathbbm{E}\big{(}\mathds{1}_{\\{Z\leqslant\pm\alpha\sqrt{Y}\\}}{Z^{2i}}\big{)}\biggr{\\}}.$
(C.3)

###### Proof.

We consider the power series expansions of $\log(z)$ and of the incomplete
gamma function:

$\displaystyle\log(z)=2\sum_{k=0}^{\infty}{1\over 2k+1}\biggl{(}{z-1\over
z+1}\biggr{)}^{2k+1},\quad\gamma(p,z)=\sum_{k=0}^{\infty}{(-1)^{k}\over
k!}\,{z^{p+k}\over
p+k}=\sum_{k=p}^{\infty}{(-1)^{k-p}\over(k-p)!}\,{z^{k}\over k},$ (C.4)

respectively. Using these expansions, we have

$\displaystyle\log(\rho+\delta T(z;\alpha,p))$
$\displaystyle=2\sum_{k=0}^{\infty}{1\over 2k+1}\biggl{[}{\rho-1+\delta
T(z;\alpha,p)\over\rho+1+\delta T(z;\alpha,p)}\biggr{]}^{2k+1}$
$\displaystyle=2\sum_{k=0}^{\infty}{1\over
2k+1}\Biggl{[}{\rho-1+{\delta\over\Gamma(p)}\sum_{i=p}^{\infty}{(-1)^{i-p}\over
2i(i-p)!\alpha^{2i}}\,{z^{2i}}\over\rho+1+{\delta\over\Gamma(p)}\sum_{j=p}^{\infty}{(-1)^{j-p}\over
2j(j-p)!\alpha^{2j}}\,{z^{2j}}}\Biggr{]}^{2k+1}.$

Letting $b_{0}=\rho-1$ and
$b_{i}=[{\delta(-1)^{i-p}/2i(i-p)!\alpha^{2i}\Gamma(p)}]\mathds{1}_{\\{i\geqslant
p\\}}$; and similarly, defining $a_{0}=\rho+1$ and
$a_{j}=[{\delta(-1)^{j-p}/2j(j-p)!\alpha^{2j}}\Gamma(p)]\mathds{1}_{\\{j\geqslant
p\\}}$; the last series can be expressed as

$\displaystyle=2\sum_{k=0}^{\infty}{1\over
2k+1}\Biggl{(}{\sum_{i=0}^{\infty}b_{i}{z^{2i}}\over\sum_{j=0}^{\infty}a_{j}{z^{2j}}}\Biggr{)}^{2k+1}.$

By application of the equation in Section 0.313 of [22] for division of power
series, the above series is written as

$\displaystyle=2\sum_{k=0}^{\infty}{1\over
2k+1}\Biggl{(}{1\over\rho+1}\sum_{i=0}^{\infty}c_{i}{z^{2i}}\Biggr{)}^{2k+1},$
(C.5)

where the coefficients $c_{i}$’s are determined from the recurrence equation

$\displaystyle
c_{n}=\biggl{(}b_{n}-{1\over\rho+1}\sum_{i=p}^{n}c_{n-i}a_{i}\biggr{)}\mathds{1}_{\\{n\geqslant
p\\}}.$

By application of the equation in Section 0.314 of [22] for power series
raised to powers, the series (C.5) is equal to

$\displaystyle
2\sum_{k=0}^{\infty}{1\over(2k+1)(\rho+1)^{2k+1}}\sum_{i=0}^{\infty}\widetilde{c}_{i,k}{z^{2i}}=2\sum_{k=0}^{\infty}{1\over(2k+1)}\biggl{[}1+{1\over(\rho+1)^{2k+1}}\sum_{i=p}^{\infty}\widetilde{c}_{i,k}{z^{2i}}\biggr{]},$

where the coefficients $\widetilde{c}_{i,k}$’s are determined by
$\widetilde{c}_{0,k}=(\rho+1)^{2k+1}$ and from the recurrence relation

$\displaystyle\widetilde{c}_{m,k}={1\over
m(\rho+1)}\sum_{i=p}^{m}\big{[}2i(k+1)-m\big{]}a_{i}\widetilde{c}_{m-i,k}\mathds{1}_{\\{m\geqslant
p\\}}.$ (C.6)

In short, we have

$\displaystyle\log(\rho+\delta
T(z;\alpha,p))=2\sum_{k=0}^{\infty}{1\over(2k+1)}\biggl{[}1+{1\over(\rho+1)^{2k+1}}\sum_{i=p}^{\infty}\widetilde{c}_{i,k}{z^{2i}}\biggr{]}.$

By using the above expansion and assuming we can change the series with the
integral sign, we have

$\displaystyle\mathbbm{E}\big{[}\mathds{1}_{\\{a\leqslant Z\leqslant
b\\}}\log(\rho+\delta T(Z;\alpha,p))\big{]}$
$\displaystyle=\int_{a}^{b}\log(\rho+\delta T(z;\alpha,p))\phi(z)\,{\rm d}z$
$\displaystyle=2\sum_{k=0}^{\infty}{1\over(2k+1)}\biggl{[}\Phi(b)-\Phi(a)+{1\over(\rho+1)^{2k+1}}\sum_{i=p}^{\infty}\widetilde{c}_{i,k}\mathbbm{E}\big{(}\mathds{1}_{\\{a\leqslant
Z\leqslant b\\}}{Z^{2i}}\big{)}\biggr{]}.$

From the above identity the identities in (C.2) and (C.3) follow. ∎

Since $S_{1}(Z)=\log(\sqrt{2\pi})+{1/2}$, by substituting identities (C.2),
(C.3) and (C.1) in Proposition 3.11, we get the formula of the Shannon entropy
for $X\sim{\rm TD}_{\Phi}(\bm{\theta})$ given in Subsection 4.4.

##  Appendix D Derivatives of the function $\log f(x;\bm{\theta})$

Let $Y\sim{\rm Gamma}(p,1)$ and let $T$ be as in (2.2). Then

$\displaystyle\begin{array}[]{lllll}\displaystyle{\partial\log
f(x;\bm{\theta})\over\partial\mu}=-\frac{2(\frac{x-\mu}{\sigma\alpha})^{2(p-1)}\,{\rm
e}^{-(\frac{x-\mu}{\alpha\sigma})^{2}}}{\sigma^{2}\alpha^{2}\Gamma(p)}\,\log\left[\rho+\delta
T\left(\frac{x-\mu}{\sigma};\alpha,p\right)\right]-{1\over\sigma}\frac{g^{\prime}\left(\frac{x-\mu}{\sigma}\right)}{g\left(\frac{x-\mu}{\sigma}\right)};\\\\[19.91684pt]
\displaystyle{\partial\log
f(x;\bm{\theta})\over\partial\sigma}=-\frac{\rho+\delta}{Z_{\bm{\theta}}}-\frac{\delta}{Z_{\bm{\theta}}}\big{\\{}\mathbbm{E}[G(-\alpha\sqrt{Y})]-\mathbbm{E}[G(\alpha\sqrt{Y})]\big{\\}}\\\\[17.07182pt]
\displaystyle\hskip
244.97824pt-\frac{2\delta(\frac{x-\mu}{\sigma\alpha})^{2p-1}{\rm
e}^{-(\frac{x-\mu}{\alpha\sigma})^{2}}}{\sigma^{3}\left[\rho+\delta
T(\frac{x-\mu}{\sigma};\alpha;p)\right]\Gamma(p)}-\frac{(x-\mu)}{\sigma^{2}}\frac{g^{\prime}\left(\frac{x-\mu}{\sigma}\right)}{g\left(\frac{x-\mu}{\sigma}\right)};\\\\[19.91684pt]
\displaystyle{\partial\log
f(x;\bm{\theta})\over\partial\alpha}=-\frac{2\sigma\delta\Gamma(p+1/2)}{\sqrt{2\pi}Z_{\bm{\theta}}\Gamma(p)}\,{{}_{2}F_{1}}\biggl{(}\frac{1}{2};p+\frac{1}{2};\frac{3}{2};\frac{-\alpha^{2}}{2}\biggr{)}\\\\[17.07182pt]
\displaystyle\hskip
65.44142pt-\frac{2\sigma\Gamma(p+1/2)}{\sqrt{2\pi}Z_{\bm{\theta}}\Gamma(p)}\biggl{[}\biggl{(}\frac{\alpha^{2}}{2}+1\biggr{)}^{-(p+{1\over
2})}\\!\\!\\!-{{}_{2}F_{1}}\biggl{(}\frac{1}{2};p+\frac{1}{2};\frac{3}{2};\frac{-\alpha^{2}}{2}\biggr{)}\biggr{]}\\!-\\!\frac{2\delta(\frac{x-\mu}{\sigma\alpha})^{2p}\,{\rm
e}^{-(\frac{x-\mu}{\alpha\sigma})^{2}}}{\alpha\left[\rho+\delta
T(\frac{x-\mu}{\sigma};\alpha,p)\right]};\\\\[19.91684pt]
\displaystyle{\partial\log
f(x;\bm{\theta})\over\partial\rho}=-\frac{\sigma}{Z_{\bm{\theta}}}+\frac{1}{\rho+\delta
T(\frac{x-\mu}{\sigma};\alpha;p)};\\\\[19.91684pt] \displaystyle{\partial\log
f(x;\bm{\theta})\over\partial\delta}=-\frac{\delta}{Z_{\bm{\theta}}}\big{\\{}1+\mathbbm{E}[G(-\alpha\sqrt{Y})]-\mathbbm{E}[G(\alpha\sqrt{Y})]\big{\\}}+\frac{T\left(\frac{x-\mu}{\sigma};\alpha,p\right)}{\rho+\delta
T\left(\frac{x-\mu}{\sigma};\alpha,p\right)}.\end{array}$

##  Appendix E Inference, Optimization, Bootstrap and Applications

###  E.1 Rao-Cramer lower bounds of estimators from MLqE and MLE

Since we provide the square roots of Rao-Cramer lower bounds of estimators
from MLqE and MLE, the confidence interval of estimators can be constructed as
the following form:

$\hat{\theta}\mp z_{\tau/2}\sqrt{\text{Var}_{q}(\hat{\theta})},$

where $\hat{\theta}$ is chosen as
$(\hat{\mu},\hat{\sigma},\hat{\alpha},\hat{\rho},\hat{\delta})$. $z_{\tau/2}$
is the critical value from the standard normal distribution when the
significance level $\tau$ is chosen by researcher [33, 24, 20].
$\text{Var}_{q}(\hat{\theta})$ is the theoretical variance of estimator
evaluated numerically by the codes in Section F.6. Thus, one can observe the
limit values of estimates of parameters to generate the different kind of
artificial data sets if it is necessary to do so at an experiment.

The lower limit values of confidence intervals for the estimates must be as
follows: $\hat{\sigma}_{L}>0$, $\hat{\alpha}_{L}>0$, $\hat{\rho}_{L}\geqslant
0$ and $\hat{\delta}_{L}\geqslant 0$. In other words, these values should be
greater than zero due to the defined values of these parameters.

###  E.2 Benefit of bootstrapped data when semiparametric and parametric
models are used for fitting data sets

The general tendency for modelling is performed by the existing techniques
which are the estimated distribution including parametric models, kernel
smoothing. As is expected, the kernel smoothing technique as a semiparametric
approach for modelling can show better performance for finite sample size when
compared with parametric models [40]. However, we can need to assume that a
data set is member of a parametric model. In this case, if an artificial data
set is necessary to model the probable results in the future or if conducting
an experiment is expensive for researchers, the necessity of parametric models
such as bimodal [5] and trimodal normal distribution can be inevitable;
because, the smooth kernel distribution used for fitting on the data set is
data-adaptive or semiparametric approach [41]. In this case, the parametric
model puts a restriction in where we have situations which give the results
mimicing the parametric model perfectly because of the probabilities coming
from PDF. If we have such results, then we need to make a comparison between
kernel smoothing and the trimodal forms of a function which can be chosen as
normal, Student $t$, etc. [38].

Since the real data set includes results taken at the moment or sometimes
experiments have to be conducted at only one time due to the cost of running
of experiment, different collections of the same data should be generated by
use of the bootstrap technique. Any mixed models as a hetero form can be
modelled by using TDΦ distribution if the trimodality exists. Since hetero-
mixing forms a non-identically distributed data set, trimodal distribution can
be necessary to perform an efficient fitting on the data sets observed from an
experiment.

Since bootstrap technique performs a random choice from a data set of
experiment, we need to make a comparison among the semiparametric and
parametric models whether or not which one will be performable in fitting on
the data set [38, 41]. The soft and strict forms of trimodal normal
distributions are capable to perform an efficient fitting on data set and is
an alternative approach when compared with smoothing technique if trimodality
exists. Because, it is concluded that a parametric model must be necessary to
generate the random numbers from corresponding probabilities of parametric
model, i.e PDF, for research instead of doing an another experiment again. If
we generate artificial numbers by use of smoothing technique, these numbers
will be generated according to the used smoothing technique which has its
probabilities corresponding empirical probability function representing the
probabilities coming from each event in the finite sample size (see Subsection
F.1) [41]. When we come across big data analytics, the population starts to
evolve and take what its real form is as. In this case, a parametric model can
be necessary [42].

###  E.3 Driving the bootstrap and the optimization of $\log_{q}$ likelihood
function for estimations of parameters

GOFTs which perform a testing for assesment of the used (semi)-parametric
model can show different performance for each model, which is why the
different GOFTs have been applied to test the fitting performance of models.
The bootstrap technique in Mathematica is applied to the real data set (see
Appendix F.3). The replication number of bootstrap is $1000$. Thus, the
probable mistakings in convergence of the optimization and the different
scenarios of real data sets are tried to be clarified. It is also noted that
the data observed after the experiment is conducted can include the
measurement error and the measurements can depend on many factors which are
known as the randomness, mistakes, unhidden factors, etc. in the experiment
[38, 39].

The bootstrapped form of the real data is not only beneficial to jump on these
kind of problems in the process of measurement but also the replicated
optimization should be useful for us to make the convergence around global
point (GP) or reach the real GP where the optimization can reach if the
replication is performed. Further, SK technique depends on the chosen function
for kernel and it is a semiparametric technique [40, 41]. As is shown by Table
2, SK has performance on the modelling. However, the performance of trimodal
normal distribution (TDΦ) can have superior performance when the closeness to
statistics generated by SK and the estimated distribution (EstD) which is
parametric method is taken into account. EstD integrated into Mathematica
software tries to find the best function (mixed form of the parametric models
or other parametric models) while conducting the fitting on the real data set.
For example, a data set can be a mixing form of two or three normal
distributions with the corresponding parameters from populations $1$ $2$ or
$3$, respectively. If we have a mixed form for the real data set, it can also
be member of trimodal distribution which is eventually not known by
researchers, which is an important gap in conducting a research. For this
reason, TDΦ has been proposed and also we make a comparison among them to
observe the performance of parametric model (see also discussion in
introduction Section 1).

The initial values for starting the optimization performed by FindMaximum with
constraints in the parameters $\alpha$, $\rho$ and $\delta$ searched at the
interval $(0,100]$ are generated by uniform distribution with $[0,1]$ [38,
39].

##  Appendix F Mathematica 12.0 codes

###  F.1 Generating procedure for random numbers if SmoothKernelDistribution
is used

The codes for generating random numbers from the smooth kernel distribution in
Mathematica is given by the following lines:

    ΨSK = SmoothKernelDistribution[x, {"Adaptive", Automatic, .1},PerformanceGoal
    Ψ-> "Quality"];
    ΨF[m_, s_] := CDF[SK, x];
    Ψm := Moment[SK, 1]; s := Sqrt[Moment[SK, 2] - Moment[SK, 1]^2];
    ΨFor[i = 0, i < n, i++,
    Ψz = Table[
    Ψx /. FindRoot[F[m, s] - RandomReal[], {x, Lower point, Upper point}], {i, n}]
    Ψ]

One can get the random number from trimodal normal distribution by use of
subsection 3.4 or CDF of TDΦ in Subsection 4.3.

###  F.2 Optimization and statistics

    iv := RandomReal[];
    FM = FindMaximum[{Total[Log[f]],
    Ψa< \[Alpha] <= b && a <= \[Rho] <= b &&
    Ψa <= \[Delta] <= b}, {{\[Alpha], iv}, {\[Rho] ,
    ΨΨiv}, {\[Delta], iv}, {\[Mu], Median[x]}, {\[Sigma],
    ΨΨMedian[Abs[x - Median[x]]]}}];
    ΨEstD = FindDistribution[x];
    ΨSK = SmoothKernelDistribution[x, {"Adaptive", Automatic, .1},PerformanceGoal
    Ψ -> "Quality"];
     statistics =
     {
     ΨMoment[EstD, 1], Moment[SK, 1],  Mean[x], Median[x],
    Ψ  Sqrt[Moment[EstD, 2] - Moment[EstD, 1]^2],
    Ψ  Sqrt[Moment[SK, 2] - Moment[SK, 1]^2],
    Ψ  StandardDeviation[x], Median[Abs[x - Median[x]]]
     };
    Ψ

###  F.3 The codes for bootstrap

A part for bootstrap is given by the following form [38]:

    ΨFor[i = 1, i <= replication, i++,
    Ψ     Ψx := Table[RandomChoice[data, Length[data]], {i,1,Length[data]}]
    Ψ   ]

###  F.4 CDF of TNΦ

    p := 3/2; ns := n;
    x = Sort[x];
    nsn = 0; nsp = 0;
    For[i = 1, i <= ns, i++,
    If[x[[i]] < \[Mu], nsn = nsn + 1]
    ]
    For[i = 1, i <= nsn, i++,
    If[x[[i]] < \[Mu], xn[i] := x[[i]]]
    ]
    For[i = nsn + 1, i <= ns, i++,
    If[x[[i]] >= \[Mu], xp[i] := x[[i]]]
    ]
    For[i = 1, i <= nsn, i++,
    FN[i] = (\[Sigma]/
    Z)*(\[Rho] + \[Delta]*
    T[xn[i]/\[Sigma] - \[Mu]/\[Sigma], \[Alpha], p])*
    CDF[NormalDistribution[\[Mu], \[Sigma]],
    xn[i]] + ((\[Delta]*\[Sigma])/
    Z)*(1/2 - (\[Alpha]*Gamma[p + 1/2]*
    Hypergeometric2F1[1/2, p + 1/2, 3/2, -\[Alpha]^2/2])/(Sqrt[
    2*Pi]*Gamma[p])) - ((\[Delta]*\[Sigma])/
    Z)*((1/2)*T[xn[i]/\[Sigma] - \[Mu]/\[Sigma], \[Alpha], p] -
    NIntegrate[
    Erf[(\[Alpha]/Sqrt[2])*z]*z^(2*p - 1)*Exp[-1^2*z^2], {z,
    Ψ0, (xn[i] - \[Mu])/(\[Alpha]*\[Sigma])},
    Method -> {"LobattoKronrodRule"}]/Gamma[p]);
    ]
    For[i = nsn + 1, i <= ns, i++,
    FP[i] = (\[Sigma]/
    Z)*(\[Rho] + \[Delta]*
    T[xp[i]/\[Sigma] - \[Mu]/\[Sigma], \[Alpha], p])*
    CDF[NormalDistribution[\[Mu], \[Sigma]],
    xp[i]] + ((\[Delta]*\[Sigma])/
    Z)*(1/2 - (\[Alpha]*Gamma[p + 1/2]*
    Hypergeometric2F1[1/2, p + 1/2, 3/2, -\[Alpha]^2/2])/(Sqrt[
    2*Pi]*Gamma[p])) - ((\[Delta]*\[Sigma])/
    Z)*((1/2)*T[xp[i]/\[Sigma] - \[Mu]/\[Sigma], \[Alpha], p] +
    NIntegrate[
    Erf[(\[Alpha]/Sqrt[2])*z]*z^(2*p - 1)*Exp[-1^2*z^2], {z,
    Ψ0, (xp[i] - \[Mu])/(\[Alpha]*\[Sigma])},
    Method -> {"LobattoKronrodRule"}]/Gamma[p]);
    ]
    FNFP := Join[{Table[FN[j], {j, nsn}], Table[FP[j], {j, nsn + 1, ns}]}];
    JoinFNP := Join[FNFP[[1]], FNFP[[2]]];

###  F.5 Goodness of fit tests for the proposed and used distributions

The goodness of fit tests such as Kolmogorov-Smirnov (KS), Cramér–von Mises
and Anderson-Darling are evaluated by using the the following codes for TNΦ,
EstD and SK.

    For[i = 1, i <= ns, i++,
    CVM[i] = (JoinFNP[[i]] - (2*i - 1)/(2*ns))^2 + 1/(12*ns);
    AD[i] = (-1/
    ns)*((2*i - 1)*(Log[JoinFNP[[i]]] + Log[1 - JoinFNP[[i]]]));
    DP[i] = i/ns - JoinFNP[[i]]; DN[i] = JoinFNP[[i]] - (i - 1)/ns;
    CVMED[i] = (CDF[estimated\[ScriptCapitalD], x][[
    i]] - (2*i - 1)/(2*ns))^2 + 1/(12*ns);
    ADED[i] = (-1/
    ns)*((2*i - 1)*(Log[CDF[EstD], x][[i]]] +
    Log[1 - CDF[EstD], x][[i]]]));
    DPED[i] = i/ns - CDF[EstD], x][[i]];
    DNED[i] = CDF[estimated\[ScriptCapitalD], x][[i]] - (i - 1)/ns;
    CVMSK[i] = (CDF[SK, x][[i]] - (2*i - 1)/(2*ns))^2 + 1/(12*ns);
    ADSK[i] = (-1/ns)*((2*i - 1)*(Log[CDF[SK, x][[i]]] + Log[1 - CDF[SK, x][[i]]]));
    DPSK[i] = i/ns - CDF[SK, x][[i]];
    DNSK[i] = CDF[SK, x][[i]] - (i - 1)/ns;
    ]
    DPmax = Max[Table[DP[j], {j, ns}]]; DNmax = Max[Table[DN[j], {j, ns}]];
    KS = Max[DPmax, DNmax];StaCVM = Total[Table[CVM[j], {j, ns}]];
    StaAD = Total[Table[AD[j], {j, ns}]];{KS, StaCVM, StaAD}

In order to get information criteria (IC) such as Akaike and Bayesian, CDF of
Mathematica must be replaced with its corresponding PDF and the formulae of IC
are used. For example,

    Ψ -2*Total[Log[PDF[SK, x]]] + 2*p

###  F.6 The codes for variances of estimators

The PDF of TDΦ is given by

    f[\[Mu]_, \[Sigma]_, \[Alpha]_, \[Rho]_, \[Delta]_] := \
    \[Sigma]*((1/((\[Rho] + \[Delta])*\[Sigma] - ((Sqrt[
    2]*\[Delta]*\[Sigma]*\[Alpha]*
    Gamma[p + 1/2])/(Sqrt[Pi]*Gamma[p]))*
    Hypergeometric2F1[1/2, p + 1/2,
    3/2, -\[Alpha]^2/2]))*(\[Rho] + \[Delta]*(1 -
    Gamma[p, (x - \[Mu])^2/(\[Sigma]^2* \[Alpha]^2)]/
    Gamma[p]))*(exp^(-((x - \[Mu])^2/(2 \[Sigma]^2)))/(
    Sqrt[2 \[Pi]] \[Sigma])));
    g[\[Mu]_, \[Sigma]_] := exp^(-((x - \[Mu])^2/(2 \[Sigma]^2)))/(
    Sqrt[2 \[Pi]] \[Sigma]);
    T[p_, \[Mu]_, \[Sigma]_, \[Alpha]_] :=
    1 - Gamma[p, (x - \[Mu])^2/(\[Sigma]^2* \[Alpha]^2)]/Gamma[p];

    Dgm2[\[Mu], \[Sigma]] := D[g[\[Mu], \[Sigma]], \[Mu]];

    Z := (\[Rho] + \[Delta])*\[Sigma] - ((Sqrt[
    2]*\[Delta]*\[Sigma]*\[Alpha]*Gamma[p + 1/2])/(Sqrt[Pi]*
    Gamma[p]))*
    Hypergeometric2F1[1/2, p + 1/2, 3/2, -\[Alpha]^2/2];

The score ${\partial\log f(x;\bm{\theta})\over\partial\mu}$ for $\mu$ is given
by

    em := ((-2*((x - \[Mu])/(\[Sigma]*\[Alpha]))^(2*p - 2)*
    Exp[-((x - \[Mu])/(\[Sigma]*\[Alpha]))^2])/(\[Sigma]^2*\
    \[Alpha]^2*Gamma[p]))*
    Log[\[Rho] + \[Delta]*
    T[p, \[Mu], \[Sigma], \[Alpha]]] - (1/\[Sigma])*
    Dgm2[\[Mu], \[Sigma]];

The score ${\partial\log f(x;\bm{\theta})\over\partial\sigma}$ for $\sigma$ is
given by

    es := -\[Rho]/Z - \[Delta]/
    Z - (\[Delta]/
    Z)*(NIntegrate[
    CDF[NormalDistribution[\[Mu], \[Sigma]], -\[Alpha]*y^0.5]*
    PDF[GammaDistribution[p, 1], y], {y, 0, Infinity},
    Method -> {"LobattoKronrodRule"}] -
    NIntegrate[
    CDF[NormalDistribution[\[Mu], \[Sigma]], \[Alpha]*y^0.5]*
    PDF[GammaDistribution[p, 1], y], {y, 0, Infinity},
    Method -> {"LobattoKronrodRule"}]) - ((2*\[Delta]*((x - \
    \[Mu])/(\[Sigma]*\[Alpha]))^(2*p - 1)*
    Exp[-((x - \[Mu])/(\[Sigma]*\[Alpha]))^2])/(\[Sigma]^3*(\[Rho] \
    + \[Delta]*T[p, \[Mu], \[Sigma], \[Alpha]])*
    Gamma[p])) - (((x - \[Mu])/\[Sigma]^2)*Dgs2[\[Mu], \[Sigma]]);

The score ${\partial\log f(x;\bm{\theta})\over\partial\alpha}$ for $\alpha$ is
given by

    ea := -(2*\[Sigma]*\[Delta]*Gamma[p + 0.5])/(Sqrt[2*Pi]*Z*Gamma[p])*
    Hypergeometric2F1[1/2, p + 1/2,
    3/2, -\[Alpha]^2/
    2] - (2*\[Sigma]*Gamma[p + 1/2])/(Sqrt[2*Pi]*Z*
    Gamma[p])*((\[Alpha]^2/2 + 1)^(-p - 1/2) -
    Hypergeometric2F1[1/2, p + 1/2,
    3/2, -\[Alpha]^2/
    2]) - (2*\[Delta]*((x - \[Mu])/(\[Sigma]*\[Alpha]))^(2*p)*
    Exp[-((x - \[Mu])/(\[Sigma]*\[Alpha]))^2])/(\[Alpha]*(\[Rho] + \
    \[Delta]*T[p, \[Mu], \[Sigma], \[Alpha]]));
    NIntegrate[
    em^2*f[\[Mu], \[Sigma], \[Alpha], \[Rho], \[Delta]]^(2 -
    q), {x, -Infinity, Infinity}, Method -> {"LobattoKronrodRule"}];

The score ${\partial\log f(x;\bm{\theta})\over\partial\rho}$ for $\rho$ is
given by

    er := -\[Sigma]/Z + 1/(\[Rho] + \[Delta]*T[p, \[Mu], \[Sigma], \[Alpha]]);

The score ${\partial\log f(x;\bm{\theta})\over\partial\delta}$ for $\delta$ is
given by

    ed := -\[Delta]/Z*(1 + NIntegrate[CDF[NormalDistribution[\[Mu], \[Sigma]],
    -\[Alpha]*y^0.5]*PDF[GammaDistribution[p, 1], y], {y, 0, Infinity},
    Method -> {"LobattoKronrodRule"}] -
    NIntegrate[CDF[NormalDistribution[\[Mu], \[Sigma]], \[Alpha]*y^0.5]*
    PDF[GammaDistribution[p, 1], y], {y, 0, Infinity},
    Method -> {"LobattoKronrodRule"}]) +
    T[p, \[Mu], \[Sigma], \[Alpha]] /
    (\[Rho] + \[Delta]*T[p, \[Mu], \[Sigma], \[Alpha]]);

An element for ${\partial\log f(x;\bm{\theta})\over\partial\mu}$ of Fisher
information matrix based on $\log_{q}$ is computed by the expression and its
corresponding elements of matrix [19]:

    Imm := NIntegrate[em^2*f[\[Mu], \[Sigma], \[Alpha], \[Rho], \[Delta]]^(2 - q),
    {x, -Infinity, Infinity}, Method -> {"LobattoKronrodRule"}]

The Fisher matrix M and inverse of M are given by

    M = n*{{Imm, Ims, Ima, Imr, Imd}, {Ims, Iss, Isa, Isr, Isd}, {Ima,
    ΨΨIsa, Iaa, Iar, Iad}, {Imr, Isr, Iar, Irr, Ird}, {Imd, Isd, Iad,
    ΨΨIrd, Idd}};IM := Inverse[M]/n