††thanks: These authors contribute equally to this work††thanks: These authors
contribute equally to this work††thanks: These authors contribute equally to
this work

# Quantum simulation of the two-dimensional Weyl equation in a magnetic field

Y. Jiang  Center for Quantum Information, Institute for Interdisciplinary
Information Sciences, Tsinghua University, Beijing, 100084, PR China  M.-L.
Cai  Center for Quantum Information, Institute for Interdisciplinary
Information Sciences, Tsinghua University, Beijing, 100084, PR China  HYQ Co.,
Ltd., Beijing, 100176, P. R. China  Y.-K. Wu  Center for Quantum Information,
Institute for Interdisciplinary Information Sciences, Tsinghua University,
Beijing, 100084, PR China  Q.-X. Mei  Center for Quantum Information,
Institute for Interdisciplinary Information Sciences, Tsinghua University,
Beijing, 100084, PR China  W.-D. Zhao  Center for Quantum Information,
Institute for Interdisciplinary Information Sciences, Tsinghua University,
Beijing, 100084, PR China  X.-Y. Chang  Center for Quantum Information,
Institute for Interdisciplinary Information Sciences, Tsinghua University,
Beijing, 100084, PR China  L. Yao  Center for Quantum Information, Institute
for Interdisciplinary Information Sciences, Tsinghua University, Beijing,
100084, PR China  HYQ Co., Ltd., Beijing, 100176, P. R. China  L. He  Center
for Quantum Information, Institute for Interdisciplinary Information Sciences,
Tsinghua University, Beijing, 100084, PR China  Z.-C. Zhou  Center for Quantum
Information, Institute for Interdisciplinary Information Sciences, Tsinghua
University, Beijing, 100084, PR China  L.-M. Duan<EMAIL_ADDRESS>Center for Quantum Information, Institute for Interdisciplinary Information
Sciences, Tsinghua University, Beijing, 100084, PR China

###### Abstract

Quantum simulation of 1D relativistic quantum mechanics has been achieved in
well-controlled systems like trapped ions, but properties like spin dynamics
and response to external magnetic fields that appear only in higher dimensions
remain unexplored. Here we simulate the dynamics of a 2D Weyl particle. We
show the linear dispersion relation of the free particle and the discrete
Landau levels in a magnetic field, and we explicitly measure the spatial and
spin dynamics from which the conservation of helicity and properties of anti-
particles can be verified. Our work extends the application of an ion trap
quantum simulator in particle physics with the additional spatial and spin
degrees of freedom.

Relativistic quantum mechanics Greiner (2000); Peskin and Schroeder (2018)
combines two most important theories of modern physics: special relativity and
quantum mechanics. It possesses negative-energy solutions which naturally lead
to the prediction of antimatter, and it also explains the half-integer value
of spins and their magnetic moments. Weyl equation Weyl (1929) is one of the
simplest relativistic quantum mechanical equations. Its solution, a Weyl
fermion, is a spin-$1/2$ particle with zero mass and a featured linear
dispersion relation, and can be derived from the famous Dirac equation Greiner
(2000); Peskin and Schroeder (2018) in the massless limit. As one of the
fundamental particles allowed by relativistic quantum mechanics Greiner
(2000); Peskin and Schroeder (2018), the Weyl fermion was long believed to
describe neutrinos, which are however extremely difficult to detect and are
now known to have nonzero masses Fukuda  _et al._ (1998); Ahmad  _et al._
(2001). This leaves many properties of the Weyl particles only to be analyzed
theoretically, or through the idea of quantum simulation Georgescu  _et al._
(2014); Cirac and Zoller (2012) using other well-controlled quantum systems.
Indeed, quantum simulation of relativistic quantum mechanical systems has been
proposed Freedman  _et al._ (2002); Lamata  _et al._ (2007); Lan  _et al._
(2011); Jordan  _et al._ (2012); Hauke  _et al._ (2013); Li  _et al._ (2019)
and performed Gerritsma  _et al._ (2010, 2011); Zhang  _et al._ (2018); Kokail
_et al._ (2019) in various physical systems like trapped ions Leibfried  _et
al._ (2003); Blatt and Roos (2012); Bruzewicz  _et al._ (2019).

To date, Weyl fermions have been realized in photonic crystals Lu  _et al._
(2015) and in condensed matter systems Xu  _et al._ (2015); Lv  _et al._
(2015), but in these systems only the spectral or the transport properties can
be measured Armitage  _et al._ (2018), while direct study of the Weyl particle
dynamics is still lacking. On the other hand, massive Dirac particles have
been simulated in ion trap Gerritsma  _et al._ (2010, 2011), which can reduce
to the massless Weyl particles by tuning the experimental parameters.
Nevertheless, to minimize the required degrees of freedom to be controlled,
the experiments so far are restricted to dynamics in 1D, where interactions
with external magnetic fields and evolution of spin states become trivial. In
this work, we report the quantum simulation of a 2D Weyl fermion, which allows
us to explore much richer spatial and spin dynamics.

Let us start from the 3D Dirac equation for a charged particle in a magnetic
field Greiner (2000) (we have chosen the natural units by setting $\hbar=1$
and $c=1$ for simplicity)

$i\frac{\partial\psi}{\partial
t}=\hat{H}\psi=[\boldsymbol{\hat{\alpha}}\cdot(\boldsymbol{\hat{p}}-e\boldsymbol{\hat{A}})+m\hat{\beta}]\psi,$
(1)

where $m$ is the mass of the particle, $e$ the electric charge,
$\boldsymbol{\hat{A}}$ the vector potential of the magnetic field, and
$\boldsymbol{\hat{p}}$ the momentum operator. $\hat{\alpha}_{j}$ ($j=x,y,z$)
and $\hat{\beta}$ are Dirac matrices satisfying
$\\{\hat{\alpha}_{i},\hat{\alpha}_{j}\\}=2\hat{I}\delta_{ij}$,
$\hat{\beta}^{2}=\hat{I}$ and $\\{\hat{\alpha}_{j},\hat{\beta}\\}=0$. In 3D
space, we thus need four matrices anti-commuting with each other, and
therefore the Dirac matrices, as well as the Dirac spinor $\psi$, need to have
a dimension of at least four Greiner (2000).

Figure 1: Schematic of simulating a 2D Weyl fermion using ion trap. We apply a
pair of Raman laser beams on a trapped ion, with different frequency
components resonant to the red and the blue sidebands of both the $x$ and $y$
oscillation modes, to couple the internal states of the ion to its spatial
oscillation modes. By further tuning the phase of these frequency components
in the Raman laser, we can choose to couple to $\hat{\sigma}_{x}$ or
$\hat{\sigma}_{y}$ of the internal states, and to different quadratures
$\hat{x}$ or $\hat{p}$ of each mode, which finally gives us the Hamiltonian in
Eq. (2).  Figure 2: Linear dispersion relation of a free Weyl particle and its
spectrum in a magnetic field. (a) Spin dynamics
$\langle\hat{\sigma}_{y}(t)\rangle$ from an initial state
$|+z\rangle|\alpha_{x}=ip/\sqrt{2}\rangle|\alpha_{y}=0\rangle$ with $B=0$ in
Eq. (2). Each data point consists of 500 experimental shots and the error bars
are estimated from the standard deviation of 10 repetitions. The solid curves
are theoretical predictions under the same parameters. (b) From the slope of
the early-time spin dynamics (inset), we can extract the average energy $E(p)$
for a free Weyl particle with momentum $p$. A linear dispersion relation is
observed as the black fitting line. The error bars are standard deviations of
5 repetitions. (c) Spin dynamics $\langle\hat{\sigma}_{z}(t)\rangle$ in a
magnetic field with $eB=1$ from an initial state
$|+z\rangle|\alpha_{x}=i\rangle|\alpha_{y}=0\rangle$. Due to the larger number
of data points and the longer evolution time, here we only repeat each point
by 200 times and the error bars are one standard deviation of the average
value. The red curve is the theoretical prediction with the motional
decoherence included. (d) The energy spectrum of a Weyl fermion in a magnetic
field can be obtained by a Fourier transform of the spin dynamics. The first
three peaks at $2E_{0}=0$, $2E_{1}=2\pi\times 8.3\,$kHz and $2E_{2}=2\pi\times
11.7\,$kHz agree well with the theoretical $E_{n}\propto\sqrt{n}$ scaling,
while higher peaks are more difficult to distinguish due to the frequency
resolution limited by the total evolution time of $600\,\mu$s. Inset is the
ideal result for an evolution time of $5\,$ms without decoherence.

Now if we go to the massless limit $m\to 0$, the matrix $\hat{\beta}$ is
removed and hence we only require three anti-commuting matrices. In this
situation we can simply set $\hat{\alpha}_{j}=\hat{\sigma}_{j}$ ($j=x,y,z$)
where $\hat{\sigma}_{j}$’s are the Pauli operators. This gives us the 3D Weyl
equation Greiner (2000) in a magnetic field. We can further define a spin
operator $\boldsymbol{\hat{S}}=\boldsymbol{\hat{\sigma}}/2$ such that the
total angular momentum
$\boldsymbol{\hat{J}}=\boldsymbol{\hat{L}}+\boldsymbol{\hat{S}}$ is conserved
($[\boldsymbol{\hat{J}},\hat{H}]=0$) when the magnetic field is zero
($\boldsymbol{\hat{A}}=0$), where
$\hat{L}_{i}=\epsilon_{ijk}\hat{x}_{j}\hat{p}_{k}$ is the orbital angular
momentum and $\epsilon_{ijk}$ the Levi-Civita symbol Greiner (2000).

Finally, for a minimal model to demonstrate nontrivial spatial and spin
dynamics in a magnetic field, we consider a uniform field along the $z$ axis
such that $\hat{A}_{y}=B\hat{x}$ and $\hat{A}_{x}=\hat{A}_{z}=0$. The momentum
in the $z$ direction is conserved so that we can restrict to the 2D case,
which gives us the Hamiltonian to be simulated

$\hat{H}=\hat{\sigma}_{x}\hat{p}_{x}+\hat{\sigma}_{y}(\hat{p}_{y}-eB\hat{x}).$
(2)

Figure 3: Spin and spatial dynamics and conservation of helicity. (a), (b)
Spin dynamics in a magnetic field $eB=1$ from the initial state
$|+x\rangle|\alpha_{x}=i\rangle|\alpha_{y}=0\rangle$. Each point is measured
by 500 times and the error bar is estimated as one standard deviation of the
average. Solid curves are theoretical results under the same parameters. (c),
(d) Evolution of the kinetic momentum
$\boldsymbol{\hat{\pi}}=\boldsymbol{\hat{p}}-e\boldsymbol{\hat{A}}$ under the
same condition. Each quadrature is measured by applying an additional spin-
dependent force and fitting the early-time evolution (see Supplementary
Materials supplementary ), with error bars estimated from one standard
deviation of the fitting. The error bar of
$\langle\hat{\pi}_{y}\rangle=\langle\hat{p}_{y}\rangle-\langle\hat{x}\rangle$
is further computed from those of $\langle\hat{p}_{y}\rangle$ and
$\langle\hat{x}\rangle$ as shown in the inset. (e) We plot the ratio between
the $x$ and $y$ components of the spin (red) and the kinetic momentum (blue).
The two values are close to each other and follow the same tendency when the
ratios by themselves change orders of magnitudes, although small difference
exists. Solid curves are the ideal theoretical results. (f) We further compute
the azimuthal angles of the spin and the kinetic momentum in the $x$-$y$
plane. The data points distribute around the diagonal, which indicates that
the spin and the kinetic momentum roughly align with each other during the
time evolution in a magnetic field, namely the conservation of helicity.

Our experimental scheme is sketched in Fig. 1. A pair of counter-propagating
$355\,$nm laser beams is shined on a trapped ${}^{171}\mathrm{Yb}^{+}$ ion to
create a spin-dependent force Choi  _et al._ (2014); Cai  _et al._ (2021). The
Raman laser has an angle of $45^{\circ}$ to both the $x$ and the $y$
directions with trap frequencies $\omega_{x}=2\pi\times 2.35\,$MHz and
$\omega_{y}=2\pi\times 1.98\,$MHz. By tuning the Raman laser resonant to the
red (blue) sideband of the mode $j$ ($j=x,\,y$), we get the Hamiltonian
$\hat{H}^{r(b)}_{j}(\Omega,\,\phi)=\Omega[\hat{\sigma}_{-(+)}\hat{a}_{j}^{{\dagger}}e^{i\phi}+h.c.]/2$
where $\hat{\sigma}_{+(-)}$ is the raising (lowering) operator of the qubit,
$\hat{a}_{j}(\hat{a}_{j}^{\dagger})$ the annihilation (creation) operator of
the corresponding mode, $\Omega$ the sideband Rabi frequency, and $\phi$
controlled by the phase of the Raman laser. Now we apply four frequency
components, driving both sidebands of the two oscillation modes simultaneously
as
$\hat{H}=\hat{H}_{x}^{r}((1-r)\Omega,\,\pi/2)+\hat{H}_{x}^{b}((1+r)\Omega,\,\pi/2)+\hat{H}_{y}^{r}(\Omega,\,\pi)+\hat{H}_{y}^{b}(\Omega,\,0)=(\Omega/\sqrt{2})[\hat{\sigma}_{x}\hat{p}_{x}+\hat{\sigma}_{y}(\hat{p}_{y}-r\hat{x})]$
where we define $\hat{x}=(\hat{a}_{x}+\hat{a}_{x}^{\dagger})/\sqrt{2}$ and
$\hat{p}_{x}=i(\hat{a}_{x}^{\dagger}-\hat{a}_{x})/\sqrt{2}$ and similarly for
the $y$ mode. This gives us the desired Hamiltonian with $r=eB$.

A characteristic property of the Weyl fermion is its linear dispersion
relation in free space. With the above spin-dependent force, we can initialize
the phonon states into
$|\psi_{m}\rangle=|\alpha_{x}=ip_{x}/\sqrt{2}\rangle|\alpha_{y}=ip_{y}/\sqrt{2}\rangle$
with expectation values of the momentum as $\langle\hat{p}_{x}\rangle=p_{x}$
and $\langle\hat{p}_{y}\rangle=p_{y}$. Without squeezing, such a coherent
state has a continuous distribution of the momentum, hence we also expect a
continuous distribution of the energy of the Weyl particle and will focus on
its expectation value. Note that, for a momentum on the $x$-$y$ plane with an
angle $\tan\theta=p_{y}/p_{x}$, we have the positive-energy and negative-
energy states when the spin is parallel or anti-parallel to it, namely the two
eigenstates of
$\hat{\sigma}_{\theta}\equiv\hat{\sigma}_{x}\cos\theta+\hat{\sigma}_{y}\sin\theta$.
Now if we initialize the momentum state $|\psi_{m}\rangle$ with spin
$|+z\rangle$, evolve the system under the free Weyl fermion Hamiltonian, and
measure the early-time dynamics for a spin perpendicular to
$\hat{\sigma}_{\theta}$ as
$\hat{\sigma}_{\theta}^{\perp}\equiv-\hat{\sigma}_{x}\sin\theta+\hat{\sigma}_{y}\cos\theta$,
we get
$\langle\hat{\sigma}_{\theta}^{\perp}(t)\rangle\approx-t[(\langle+\sigma_{\theta}|\langle\psi_{m}|)\hat{H}(|+\sigma_{\theta}\rangle|\psi_{m}\rangle)-(\langle-\sigma_{\theta}|\langle\psi_{m}|)\hat{H}(|-\sigma_{\theta}\rangle|\psi_{m}\rangle)]=-2E(p)t$
(see Supplementary Materials supplementary ). For example, in Fig. 2(a) we
plot the spin dynamics $\langle\sigma_{y}(t)\rangle$ from the initial state
$|+z\rangle|\alpha_{x}=ip/\sqrt{2}\rangle|\alpha_{y}=0\rangle$ under various
values of $p$, from which we extract a linear dispersion relation $E(p)\propto
p$ as shown in Fig. 2(b).

With the existence of a magnetic field, the continuous spectrum collapses into
discrete Landau levels Landau (1930). For the massless Weyl particle, these
energy levels have a unique scaling $E_{n}^{\mathrm{Weyl}}=\sqrt{2neB}$
compared with the Dirac particle $E_{n}^{\mathrm{Dirac}}=\sqrt{m^{2}+2neB}$
Lamata  _et al._ (2011) which in the nonrelativistic limit gives us the well-
known result of $E_{n}^{\mathrm{NR}}=neB/m$ Landau (1930). In Fig. 2(c), we
plot the spin dynamics $\langle\sigma_{z}(t)\rangle$ in a magnetic field
$eB=1$, which is oscillating at the frequency difference between different
energy levels. Through a Fourier transform, we can identify discrete peaks in
the spectrum as shown in Fig. 2(d). Theoretically, these peaks locate at
$2E_{n}$ (inset, see Supplementary Materials supplementary ) where the factor
of $2$ comes from the positive and negative energy states. In this experiment,
the evolution time up to $600\,\mu$s (which should be much shorter than the
decoherence time of the motional states) restricts our frequency resolution so
that high energy levels cannot be distinguished, but the first three peaks at
$n=0,\,1,\,2$ already show good agreement with the predicted $\sqrt{n}$
scaling.

Next we examine the spatial and spin dynamics of a Weyl particle in the
magnetic field. As we can see from Fig. 3(a)-(d), in general, the components
of the spin or the momentum are not conserved (except for $p_{y}$ in the inset
which is related to our gauge choice). However, as shown in Fig. 3(e) and (f),
the directions of the spin $\boldsymbol{S}$ and the kinetic momentum
$\boldsymbol{\pi}\equiv\boldsymbol{p}-e\boldsymbol{A}$ stay roughly the same
during the evolution. This is known as the conservation of helicity
$\hat{h}\equiv\boldsymbol{\hat{\sigma}}\cdot\boldsymbol{\hat{\pi}}/|\boldsymbol{\hat{\pi}}|$
Greiner (2000) in a magnetic field where $|\boldsymbol{\hat{\pi}}|$ represents
the magnitude of the kinetic momentum. Strictly speaking, this conservation is
proved for the scattering problem where the incoming and the outgoing
particles do not feel the magnetic field Goldhaber (1977), and from the
theoretical curves we do observe small fluctuation between the orientations of
these two vectors even under the ideal evolution. On the other hand, what we
simulate here is more like the Larmor procession of a spin in a uniform
magnetic field. By showing the spin procession speed to be matched by the
spatial rotation rate $\omega=eB/m$ of the particle, this naturally explains
the Lande $g$-factor of 2 for the spin-1/2 particles, which had to be added
into the theory by hand before the relativistic quantum mechanics Greiner
(2000).

Figure 4: Trajectories of 2D Weyl particles with opposite helicities. Using
the method of Fig. 3, we can plot the trajectory in a magnetic field $eB=1$
from the initial state $|+x\rangle|\alpha_{x}=i\rangle|\alpha_{y}=0\rangle$
(red) and $|-x\rangle|\alpha_{x}=i\rangle|\alpha_{y}=0\rangle$ (blue). These
correspond to a Weyl particle with positive helicity and an anti-particle with
negative helicity, respectively. Theoretically, the particles circles in real
space with 2D Zitterbewegung as shown by the solid curves. Here our
measurement precision cannot resolve the Zitterbewegung, but it can clearly be
seen that the particle and the anti-particle have opposite charges as they
rotate in opposite directions in a magnetic field.

Following the same method, we can measure the spatial dynamics of a Weyl
fermion with different initial spin states and we plot the 2D trajectories in
Fig. 4. For an initial state
$|+x\rangle|\alpha_{x}=i\rangle|\alpha_{y}=0\rangle$ (red), namely a helicity
of $+1$, the trajectory starts to the right and goes clockwise. On the other
hand, for an initial state
$|-x\rangle|\alpha_{x}=i\rangle|\alpha_{y}=0\rangle$ (blue) with a helicity of
$-1$, the trajectory points to the opposite direction and goes
counterclockwise. This is because the two situations correspond to a particle
and its antiparticle, respectively, with opposite masses and charges. By
initializing them with the same momentum, we thus get opposite initial
velocities, and they bend into different directions due to the different signs
in their charges.

To sum up, we have simulated a 2D Weyl particle with characteristic spectral
properties using a single trapped ion and two of its spatial oscillation
modes. Different from a neutrino, the simulated particle has a nonzero charge
and demonstrates nontrivial spatial and spin dynamics in a magnetic field.
Therefore our scheme provides a direct approach to study the dynamics of a
class of elementary particles allowed by relativistic quantum mechanics but
have not been discovered in nature. Theoretically, there will also be 2D
Zitterbewegung on the trajectories (solid curves in Fig. 4), whose amplitudes
are however comparable to the error bars in our measurement results and
therefore require future efforts to improve the experimental precision. By
further setting the trap frequencies in the three spatial directions to be
comparable and by pointing the Raman laser to be at a nonzero angle to all
these axes, we can generalize this scheme to a 3D Weyl particle with even
richer dynamics. Besides, our scheme of using spatial oscillation modes in
higher dimensions can also be generalized to multi-ion cases such as the
quantum simulation of the spin model Monroe  _et al._ (2021) and the spin-
phonon coupled system Mei  _et al._ (2021). Our work thus demonstrates the
trapped ion system as a powerful quantum simulator, and significantly extends
its application in particle physics by providing more spatial and spin degrees
of freedom beyond 1D.

###### Acknowledgements.

This work was supported by Tsinghua University Initiative Scientific Research
Program, Beijing Academy of Quantum Information Sciences, and Frontier Science
Center for Quantum Information of the Ministry of Education of China. Y.-K. W.
acknowledges support from the start-up fund from Tsinghua University.

## References

  * Greiner (2000) Walter Greiner,  _Relativistic Quantum Mechanics: Wave Equations_ , 3rd ed. (Springer, 2000).
  * Peskin and Schroeder (2018) Michael Peskin and Daniel V. Schroeder,  _An introduction to quantum field theory_ (CRC press, 2018).
  * Weyl (1929) Hermann Weyl, “Elektron und gravitation. i,” Zeitschrift für Physik 56, 330–352 (1929).
  * Fukuda  _et al._ (1998) Y. Fukuda, T. Hayakawa, E. Ichihara, K. Inoue, K. Ishihara, H. Ishino, Y. Itow, T. Kajita, J. Kameda, S. Kasuga, K. Kobayashi, Y. Kobayashi, Y. Koshio, M. Miura, M. Nakahata, S. Nakayama, A. Okada, K. Okumura, N. Sakurai, M. Shiozawa, Y. Suzuki, Y. Takeuchi, Y. Totsuka, S. Yamada, M. Earl, A. Habig, E. Kearns, M. D. Messier, K. Scholberg, J. L. Stone, L. R. Sulak, C. W. Walter, M. Goldhaber, T. Barszczxak, D. Casper, W. Gajewski, P. G. Halverson, J. Hsu, W. R. Kropp, L. R. Price, F. Reines, M. Smy, H. W. Sobel, M. R. Vagins, K. S. Ganezer, W. E. Keig, R. W. Ellsworth, S. Tasaka, J. W. Flanagan, A. Kibayashi, J. G. Learned, S. Matsuno, V. J. Stenger, D. Takemori, T. Ishii, J. Kanzaki, T. Kobayashi, S. Mine, K. Nakamura, K. Nishikawa, Y. Oyama, A. Sakai, M. Sakuda, O. Sasaki, S. Echigo, M. Kohama, A. T. Suzuki, T. J. Haines, E. Blaufuss, B. K. Kim, R. Sanford, R. Svoboda, M. L. Chen, Z. Conner, J. A. Goodman, G. W. Sullivan, J. Hill, C. K. Jung, K. Martens, C. Mauger, C. McGrew, E. Sharkey, B. Viren, C. Yanagisawa, W. Doki, K. Miyano, H. Okazawa, C. Saji, M. Takahata, Y. Nagashima, M. Takita, T. Yamaguchi, M. Yoshida, S. B. Kim, M. Etoh, K. Fujita, A. Hasegawa, T. Hasegawa, S. Hatakeyama, T. Iwamoto, M. Koga, T. Maruyama, H. Ogawa, J. Shirai, A. Suzuki, F. Tsushima, M. Koshiba, M. Nemoto, K. Nishijima, T. Futagami, Y. Hayato, Y. Kanaya, K. Kaneyuki, Y. Watanabe, D. Kielczewska, R. A. Doyle, J. S. George, A. L. Stachyra, L. L. Wai, R. J. Wilkes, and K. K. Young (Super-Kamiokande Collaboration), “Evidence for oscillation of atmospheric neutrinos,” Phys. Rev. Lett. 81, 1562–1567 (1998).
  * Ahmad  _et al._ (2001) Q. R. Ahmad, R. C. Allen, T. C. Andersen, J. D. Anglin, G. Bühler, J. C. Barton, E. W. Beier, M. Bercovitch, J. Bigu, S. Biller, R. A. Black, I. Blevis, R. J. Boardman, J. Boger, E. Bonvin, M. G. Boulay, M. G. Bowler, T. J. Bowles, S. J. Brice, M. C. Browne, T. V. Bullard, T. H. Burritt, K. Cameron, J. Cameron, Y. D. Chan, M. Chen, H. H. Chen, X. Chen, M. C. Chon, B. T. Cleveland, E. T. H. Clifford, J. H. M. Cowan, D. F. Cowen, G. A. Cox, Y. Dai, X. Dai, F. Dalnoki-Veress, W. F. Davidson, P. J. Doe, G. Doucas, M. R. Dragowsky, C. A. Duba, F. A. Duncan, J. Dunmore, E. D. Earle, S. R. Elliott, H. C. Evans, G. T. Ewan, J. Farine, H. Fergani, A. P. Ferraris, R. J. Ford, M. M. Fowler, K. Frame, E. D. Frank, W. Frati, J. V. Germani, S. Gil, A. Goldschmidt, D. R. Grant, R. L. Hahn, A. L. Hallin, E. D. Hallman, A. Hamer, A. A. Hamian, R. U. Haq, C. K. Hargrove, P. J. Harvey, R. Hazama, R. Heaton, K. M. Heeger, W. J. Heintzelman, J. Heise, R. L. Helmer, J. D. Hepburn, H. Heron, J. Hewett, A. Hime, M. Howe, J. G. Hykawy, M. C. P. Isaac, P. Jagam, N. A. Jelley, C. Jillings, G. Jonkmans, J. Karn, P. T. Keener, K. Kirch, J. R. Klein, A. B. Knox, R. J. Komar, R. Kouzes, T. Kutter, C. C. M. Kyba, J. Law, I. T. Lawson, M. Lay, H. W. Lee, K. T. Lesko, J. R. Leslie, I. Levine, W. Locke, M. M. Lowry, S. Luoma, J. Lyon, S. Majerus, H. B. Mak, A. D. Marino, N. McCauley, A. B. McDonald, D. S. McDonald, K. McFarlane, G. McGregor, W. McLatchie, R. Meijer Drees, H. Mes, C. Mifflin, G. G. Miller, G. Milton, B. A. Moffat, M. Moorhead, C. W. Nally, M. S. Neubauer, F. M. Newcomer, H. S. Ng, A. J. Noble, E. B. Norman, V. M. Novikov, M. O’Neill, C. E. Okada, R. W. Ollerhead, M. Omori, J. L. Orrell, S. M. Oser, A. W. P. Poon, T. J. Radcliffe, A. Roberge, B. C. Robertson, R. G. H. Robertson, J. K. Rowley, V. L. Rusu, E. Saettler, K. K. Schaffer, A. Schuelke, M. H. Schwendener, H. Seifert, M. Shatkay, J. J. Simpson, D. Sinclair, P. Skensved, A. R. Smith, M. W. E. Smith, N. Starinsky, T. D. Steiger, R. G. Stokstad, R. S. Storey, B. Sur, R. Tafirout, N. Tagg, N. W. Tanner, R. K. Taplin, M. Thorman, P. Thornewell, P. T. Trent, Y. I. Tserkovnyak, R. Van Berg, R. G. Van de Water, C. J. Virtue, C. E. Waltham, J.-X. Wang, D. L. Wark, N. West, J. B. Wilhelmy, J. F. Wilkerson, J. Wilson, P. Wittich, J. M. Wouters, and M. Yeh (SNO Collaboration), “Measurement of the rate of ${\nu}_{e}+\mathit{d}\rightarrow\mathit{p}+\mathit{p}+{\mathit{e}}^{-}$ interactions produced by ${}^{8}\mathrm{B}$ solar neutrinos at the Sudbury Neutrino Observatory,” Phys. Rev. Lett. 87, 071301 (2001).
  * Georgescu  _et al._ (2014) I. M. Georgescu, S. Ashhab, and Franco Nori, “Quantum simulation,” Rev. Mod. Phys. 86, 153–185 (2014).
  * Cirac and Zoller (2012) J. Ignacio Cirac and Peter Zoller, “Goals and opportunities in quantum simulation,” Nat. Phys. 8, 264–266 (2012).
  * Freedman  _et al._ (2002) Michael H. Freedman, Alexei Kitaev, and Zhenghan Wang, “Simulation of topological field theories by quantum computers,” Communications in Mathematical Physics 227, 587–603 (2002).
  * Lamata  _et al._ (2007) L. Lamata, J. León, T. Schätz, and E. Solano, “Dirac equation and quantum relativistic effects in a single trapped ion,” Phys. Rev. Lett. 98, 253005 (2007).
  * Lan  _et al._ (2011) Z. Lan, N. Goldman, A. Bermudez, W. Lu, and P. Öhberg, “Dirac-Weyl fermions with arbitrary spin in two-dimensional optical superlattices,” Phys. Rev. B 84, 165115 (2011).
  * Jordan  _et al._ (2012) Stephen P. Jordan, Keith S. M. Lee, and John Preskill, “Quantum algorithms for quantum field theories,” Science 336, 1130–1133 (2012).
  * Hauke  _et al._ (2013) P. Hauke, D. Marcos, M. Dalmonte, and P. Zoller, “Quantum simulation of a lattice Schwinger model in a chain of trapped ions,” Phys. Rev. X 3, 041018 (2013).
  * Li  _et al._ (2019) De-Sheng Li, Chun-Wang Wu, Lin-Ze He, Wei Wu, and Ping-Xing Chen, “Quantum simulation of the Weyl equation with a trapped ion,” Quantum Information Processing 18, 151 (2019).
  * Gerritsma  _et al._ (2010) R. Gerritsma, G. Kirchmair, F. Zähringer, E. Solano, R. Blatt, and C. F. Roos, “Quantum simulation of the Dirac equation,” Nature 463, 68–71 (2010).
  * Gerritsma  _et al._ (2011) R. Gerritsma, B. P. Lanyon, G. Kirchmair, F. Zähringer, C. Hempel, J. Casanova, J. J. García-Ripoll, E. Solano, R. Blatt, and C. F. Roos, “Quantum simulation of the Klein paradox with trapped ions,” Phys. Rev. Lett. 106, 060503 (2011).
  * Zhang  _et al._ (2018) Xiang Zhang, Kuan Zhang, Yangchao Shen, Shuaining Zhang, Jing-Ning Zhang, Man-Hong Yung, Jorge Casanova, Julen S. Pedernales, Lucas Lamata, Enrique Solano, and Kihwan Kim, “Experimental quantum simulation of fermion-antifermion scattering via boson exchange in a trapped ion,” Nature Communications 9, 195 (2018).
  * Kokail  _et al._ (2019) C. Kokail, C. Maier, R. van Bijnen, T. Brydges, M. K. Joshi, P. Jurcevic, C. A. Muschik, P. Silvi, R. Blatt, C. F. Roos, and P. Zoller, “Self-verifying variational quantum simulation of lattice models,” Nature 569, 355–360 (2019).
  * Leibfried  _et al._ (2003) D. Leibfried, R. Blatt, C. Monroe, and D. Wineland, “Quantum dynamics of single trapped ions,” Rev. Mod. Phys. 75, 281–324 (2003).
  * Blatt and Roos (2012) R. Blatt and C. F. Roos, “Quantum simulations with trapped ions,” Nature Physics 8, 277–284 (2012).
  * Bruzewicz  _et al._ (2019) Colin D. Bruzewicz, John Chiaverini, Robert McConnell, and Jeremy M. Sage, “Trapped-ion quantum computing: Progress and challenges,” Applied Physics Reviews 6, 021314 (2019).
  * Lu  _et al._ (2015) Ling Lu, Zhiyu Wang, Dexin Ye, Lixin Ran, Liang Fu, John D. Joannopoulos, and Marin Soljačić, “Experimental observation of Weyl points,” Science 349, 622–624 (2015).
  * Xu  _et al._ (2015) Su-Yang Xu, Ilya Belopolski, Nasser Alidoust, Madhab Neupane, Guang Bian, Chenglong Zhang, Raman Sankar, Guoqing Chang, Zhujun Yuan, Chi-Cheng Lee, Shin-Ming Huang, Hao Zheng, Jie Ma, Daniel S. Sanchez, BaoKai Wang, Arun Bansil, Fangcheng Chou, Pavel P. Shibayev, Hsin Lin, Shuang Jia, and M. Zahid Hasan, “Discovery of a Weyl fermion semimetal and topological Fermi arcs,” Science 349, 613–617 (2015).
  * Lv  _et al._ (2015) B. Q. Lv, H. M. Weng, B. B. Fu, X. P. Wang, H. Miao, J. Ma, P. Richard, X. C. Huang, L. X. Zhao, G. F. Chen, Z. Fang, X. Dai, T. Qian, and H. Ding, “Experimental discovery of Weyl semimetal taas,” Phys. Rev. X 5, 031013 (2015).
  * Armitage  _et al._ (2018) N. P. Armitage, E. J. Mele, and Ashvin Vishwanath, “Weyl and Dirac semimetals in three-dimensional solids,” Rev. Mod. Phys. 90, 015001 (2018).
  * Choi  _et al._ (2014) T. Choi, S. Debnath, T. A. Manning, C. Figgatt, Z.-X. Gong, L.-M. Duan, and C. Monroe, “Optimal quantum control of multimode couplings between trapped ion qubits for scalable entanglement,” Phys. Rev. Lett. 112, 190502 (2014).
  * Cai  _et al._ (2021) M.-L. Cai, Z.-D. Liu, W.-D. Zhao, Y.-K. Wu, Q.-X. Mei, Y. Jiang, L. He, X. Zhang, Z.-C. Zhou, and L.-M. Duan, “Observation of a quantum phase transition in the quantum Rabi model with a single trapped ion,” Nature Communications 12, 1126 (2021).
  * (27) See Supplementary Materials for details about the experimental setup, measurement and numerical simulation methods, which includes Ref. [28].
  * Johansson  _et al._ (2013) J.R. Johansson, P.D. Nation, and Franco Nori, “Qutip 2: A python framework for the dynamics of open quantum systems,” Computer Physics Communications 184, 1234–1240 (2013).
  * Landau (1930) L. Landau, “Diamagnetismus der metalle,” Zeitschrift für Physik 64, 629–637 (1930).
  * Lamata  _et al._ (2011) L Lamata, J Casanova, R Gerritsma, C F Roos, J J García-Ripoll, and E Solano, “Relativistic quantum mechanics with trapped ions,” New Journal of Physics 13, 095003 (2011).
  * Goldhaber (1977) Alfred S. Goldhaber, “Dirac particle in a magnetic field: Symmetries and their breaking by monopole singularities,” Phys. Rev. D 16, 1815–1827 (1977).
  * Monroe  _et al._ (2021) C. Monroe, W. C. Campbell, L.-M. Duan, Z.-X. Gong, A. V. Gorshkov, P. W. Hess, R. Islam, K. Kim, N. M. Linke, G. Pagano, P. Richerme, C. Senko, and N. Y. Yao, “Programmable quantum simulations of spin systems with trapped ions,” Rev. Mod. Phys. 93, 025001 (2021).
  * Mei  _et al._ (2021) Q. Mei, B. Li, Y. Wu, M. Cai, Y. Wang, L. Yao, Z. Zhou, and L. Duan, “Experimental realization of Rabi-Hubbard model with trapped ions,” (2021), arXiv:2110.03227 [quant-ph] .


# On the bimodal Gumbel model with application to environmental data

Cira E. G. Otiniano<EMAIL_ADDRESS>Departamento de Estatística,
Universidade de Brasília, 70910-900, Brasília, Brazil  Roberto Vila
<EMAIL_ADDRESS>Departamento de Estatística, Universidade de Brasília,
70910-900, Brasília, Brazil  Pedro C. Brom<EMAIL_ADDRESS>Departamento de
Estatística, Universidade de Brasília, 70910-900, Brasília, Brazil  Marcelo
Bourguignon∗<EMAIL_ADDRESS>Departamento de Estatística,
Universidade Federal do Rio Grande do Norte, 59078-970, Natal/RN, Brazil

###### Abstract

The Gumbel model is a very popular statistical model due to its wide
applicability for instance in the course of certain survival, environmental,
financial or reliability studies. In this work, we have introduced a bimodal
generalization of the Gumbel distribution that can be an alternative to model
bimodal data. We derive the analytical shapes of the corresponding probability
density function and the hazard rate function and provide graphical
illustrations. Furthermore, We have discussed the properties of this density
such as mode, bimodality, moment generating function and moments. Our results
were verified using the Markov chain Monte Carlo simulation method. The
maximum likelihood method is used for parameters estimation. Finally, we also
carry out an application to real data that demonstrates the usefulness of the
proposed distribution.

Keywords. Gumbel distribution $\cdot$ BG $\cdot$ MCMC.
Mathematics Subject Classification (2010). MSC 60E05 $\cdot$ MSC 62Exx $\cdot$
MSC 62Fxx.

###### Contents

  1. 1 Introduction
  2. 2 The bimodal Gumbel model
  3. 3 Some properties of the BG distribution
    1. 3.1 Cumulative distribution function
    2. 3.2 Stochastic representation
    3. 3.3 Rate of a random variable with BG distribution
    4. 3.4 Bimodality
    5. 3.5 The hazard function
    6. 3.6 Moment-Generating Function
    7. 3.7 Moments
  4. 4 Graphical illustrations and simulation results
    1. 4.1 Graphical illustration
    2. 4.2 Simulation results
  5. 5 Maximum likelihood estimation
  6. 6 Real-world data analysis
  7. 7 Conclusions

##  1 Introduction

Asymmetrical models for a real-valued random variable such as the Gumbel and
generalized extreme value distributions have been extensively utilized for
modeling various random phenomena encountered for instance in the course of
certain environmental, financial or reliability studies.

Let $Y_{1},Y_{2},\dots,Y_{n}$ be a series of independent random variables with
common distribution function $F$, and
$M_{n}=\max\\{Y_{1},Y_{2},\dots,Y_{n}\\}$. The Gumbel distribution is one of
the extreme value distributions, characterized by Fisher and Tippet (1928) [4]
as the limit distribution for maxima. That is, if exist the normalization
sequences $a_{n}$ and $b_{n}>0$ such that $P[(M_{n}-a_{n})/b_{n}]$ converges
to a non-degenerate distribution $G$, then $G$ is an extreme value
distribution. The $G$ distribution must be one of three types: Fréchet, Gumbel
or negative Weibull. Among these three distributions the Gumbel distribution
(see Gumbel (1958) [5] ) is the only one with a light tail. For this reason it
can be considered as a good alternative for modeling extreme data whose tails
are not heavy.

A random variable $Y$ have a Gumbel distribution with location parameter
$\mu\in\mathbb{R}$ and scale paramer $\sigma>0$, denoted by $Y\sim
F_{G}(\cdot;\mu,\sigma)$, if its probability density function (PDF) and
cumulated distribution function (CDF) are given, respectively, by

$\displaystyle
f_{G}(y;\mu,\sigma)=\frac{1}{\sigma}\exp\biggl{\\{}-\Big{(}\frac{y-\mu}{\sigma}\Big{)}-\exp\Big{[}-\Big{(}\frac{y-\mu}{\sigma}\Big{)}\Big{]}\biggr{\\}}\quad$
(1.1)

and

$\displaystyle
F_{G}(y;\mu,\sigma)=\displaystyle\exp\biggl{\\{}-\exp\Big{[}-\Big{(}\frac{y-\mu}{\sigma}\Big{)}\Big{]}\biggr{\\}},\quad
y\in\mathbb{R}.$ (1.2)

Generalizations of the Gumbel distribution have been proposed by several
authors. Pinheiro and Ferrari (2015) [13] carried out a comprehensive review
of the generalization of Gumbel distribution and after comparing them they
conclude that some distributions suffer from overparameterization. Another
generalization of Gumbel is the Exponentiated Gumbel Type-2 by Okorie (2016)
[11] and the references in [12] and [13].

Despite its broad applicability in many fields the Gumbel is not suitable to
model bimodal data. Furthermore, all the models cited above are not suitable
for capturing this. In this context, in this paper, we propose the bimodal
Gumbel (BG) distribution as an alternative model of extreme data with more
than one mode. Our approach consists of introduce bimodality in (1.1) as Elal-
Olivero (2010) [3]. Results involving bimodality in other related
probabilistic models can be found, for example, in Martinez et al. (2013) [9],
Çankaya et al. (2015) [2]; and more recently, in Vila et al. (2020-2021) [16,
17, 18, 19]. The advantage of our model in comparison to other generalizations
of the Gumbel distribution is the number of parameters and the fact that it
can be used to model extreme data with one or two modes.

This paper is organized as follows. In Section 2, we define the BG
distribution. In Section 3, we provide general properties of the BG
distribution including the cumulative distribution function, hazard rate
function, mode, bimodality, moment generating function and moments, and
stochastic representation. In Section 4, we provide graphical illustrations.
Estimation of the parameters by maximum likelihood is investigated in Section
5. In Section 6, we discuss an application to real data. Some conclusions are
addressed in Section 7.

##  2 The bimodal Gumbel model

We said that a real-valued random variable $X$ has a bimodal Gumbel (BG)
distribution with parameters $\mu\in\mathbb{R},\sigma>0$ and
$\delta\in\mathbb{R}$, denoted by $X\sim F_{\rm BG}(\cdot;\mu,\sigma,\delta)$,
if its PDF is given by

$\displaystyle f_{\rm BG}(x;\mu,\sigma,\delta)={\displaystyle\big{[}(1-\delta
x)^{2}+1\big{]}\exp\biggl{\\{}-\Big{(}\frac{x-\mu}{\sigma}\Big{)}-\exp\Big{[}-\Big{(}\frac{x-\mu}{\sigma}\Big{)}\Big{]}\biggr{\\}}\over\textstyle\sigma\big{[}1+\delta^{2}\sigma^{2}{\pi^{2}\over
6}+(\delta\mu+\delta\sigma\gamma-1)^{2}\big{]}},\quad x\in\mathbb{R},$ (2.1)

where $\gamma$ is the Euler’s constant.

Let

$\displaystyle Z_{\delta}=1+\delta^{2}\sigma^{2}{\pi^{2}\over
6}+(\delta\mu+\delta\sigma\gamma-1)^{2}$ (2.2)

be the normalization constant of the BG distribution (2.1). Using this
notation, note that $f_{\rm BG}(x;\mu,\sigma,\delta)={1\over
Z_{\delta}}\,[(1-\delta x)^{2}+1]\,f_{\rm G}(x;\mu,\sigma)$ and that $f_{\rm
BG}(x;\mu,\sigma,0)=f_{\rm G}(x;\mu,\sigma)$. In other words, we introduce
bimodality in the Gumbel distribution (1.1) through of a quadratic
transformation $(1-\delta x)^{2}+1$.

The BG distribution function (2.1) is well defined, because

$\displaystyle\int_{-\infty}^{\infty}f_{\rm BG}(x;\mu,\sigma,\delta)\,{\rm
d}x$ $\displaystyle=$ $\displaystyle{1\over
Z_{\delta}}\,[1+\mathbb{E}(1-\delta Y)^{2}],\quad Y\sim
F_{G}(\cdot;\mu,\sigma)$ $\displaystyle=$ $\displaystyle{1\over
Z_{\delta}}\,\big{\\{}1+\delta^{2}\mathrm{Var}(Y)+[\delta\,\mathbb{E}(Y)-1]^{2}\big{\\}}=1,$

where $\mathbb{E}(Y)=\mu+\sigma\gamma$ and
$\mathrm{Var}(Y)=\sigma^{2}{\pi^{2}\over 6}$.

The CDF of a BG random variable $X\sim F_{\rm BG}(\cdot;\mu,\sigma,\delta)$,
defined for $x\in\mathbb{R}$, is given by (see Subsection 3.1 for more
details)

$\displaystyle F_{\rm BG}(x;\mu,\sigma,\delta)=$
$\displaystyle{\displaystyle\big{[}2-\delta\mu(2-\delta\mu)\big{]}\exp\biggl{\\{}-\exp\Big{[}-\Big{(}\frac{x-\mu}{\sigma}\Big{)}\Big{]}\biggr{\\}}+\delta^{2}\sigma^{2}I\biggl{(}2;\exp\Big{[}-\Big{(}\frac{x-\mu}{\sigma}\Big{)}\Big{]},+\infty\biggr{)}\over
1+\delta^{2}\sigma^{2}{\pi^{2}\over 6}+(\delta\mu+\delta\sigma\gamma-1)^{2}}$
$\displaystyle+{\displaystyle
2\delta(1-\delta\mu)\Biggl{\\{}{(x-\mu)}\,\exp\biggl{\\{}-\exp\Big{[}-\Big{(}\frac{x-\mu}{\sigma}\Big{)}\Big{]}\biggr{\\}}-\sigma\Gamma\biggl{(}0,\exp\Big{[}-\Big{(}\frac{x-\mu}{\sigma}\Big{)}\Big{]}\biggr{)}\Biggr{\\}}\over
1+\delta^{2}\sigma^{2}{\pi^{2}\over 6}+(\delta\mu+\delta\sigma\gamma-1)^{2}},$
(2.3)

where $\Gamma(a,b)$ is the upper incomplete gamma function and

$\displaystyle I(k;a,b)=(-1)^{k}\int_{a}^{b}\ln^{k}(v)\exp(-v)\,{\rm d}v,\quad
k\in\mathbb{N}\cup\\{0\\},\ 0\leqslant a<b\leqslant+\infty$ (2.4)

is the incomplete moments of the random variable $V\sim F_{G}(\cdot;0,1)$. For
$k>1$, closed form solutions for the definite integral $I(k;a,b)$ are not
available in terms of commonly used functions.

From the formula in (2), when $\delta=0$, $F_{\rm
BG}(x;\mu,\sigma,0)=\exp\\{-\exp[-(\frac{x-\mu}{\sigma})]\\}=F_{G}(x;\mu,\sigma)$,
$x\in\mathbb{R}$, what was known in (1.2).

######  Remark 1.

Let $X\sim F_{\rm BG}(\cdot;\mu,\sigma,\delta)$ and let $g(\cdot)$ be a real-
valued Borel measurable function. From definition of expectation and by using
the PDF of the BG distribution (2.1), we have

$\displaystyle\mathbb{E}[g(X)]$ $\displaystyle=$ $\displaystyle{1\over
Z_{\delta}}\,\mathbb{E}\big{[}g(Y)(1-\delta Y)^{2}+g(Y)\big{]}$
$\displaystyle=$ $\displaystyle{1\over
Z_{\delta}}\,\big{\\{}2\mathbb{E}[g(Y)]-2\delta\mathbb{E}[Yg(Y)]+\delta^{2}\mathbb{E}[Y^{2}g(Y)]\big{\\}},$

where $Y\sim F_{G}(\cdot;\mu,\sigma)$ and $Z_{\delta}$ is as in (2.2).

##  3 Some properties of the BG distribution

In this section, some mathematical properties as closed expression for the
CDF, rate, modes, bimodality, hazard function and moments of the BG
distribution are discussed.

###  3.1 Cumulative distribution function

In this subsection, we derive in detail the closed expression in (2) for the
CDF of a BG random variable $X\sim F_{\rm BG}(\cdot;\mu,\sigma,\delta)$.

Indeed, from Remark 1 with $g(X)=\mathds{1}_{X\leqslant x}$, we obtain

$\displaystyle F_{\rm
BG}(x;\mu,\sigma,\delta)=\frac{1}{Z_{\delta}}\big{[}2F_{\rm
G}(y;\mu,\sigma)-2\delta\mathbb{E}\left(Y\mathds{1}_{Y\leq
x}\right)+\delta^{2}\mathbb{E}\left(Y^{2}\mathds{1}_{Y\leq x}\right)\big{]},$
(3.1)

where $Y\sim F_{G}(\cdot;\mu,\sigma)$ and $Z_{\delta}$ is as in 2.2.

Taking the following change of variables $z=\exp[-({y-\mu})/{\sigma}],\ {\rm
d}z=-({z/\sigma}){\rm d}y,\ y=\mu-\sigma\ln(z),$ and using a binomial
expansion we have

$\displaystyle\mathbb{E}\big{(}Y^{k}\mathds{1}_{Y\leq x}\big{)}$
$\displaystyle=$
$\displaystyle\int_{\exp[-(\frac{x-\mu}{\sigma})]}^{+\infty}[\mu-\sigma\ln(z)]^{k}\exp(-z)\,{\rm
d}z$ (3.2) $\displaystyle=$
$\displaystyle\sum_{i=0}^{k}\binom{k}{i}(-1)^{i}\mu^{k-i}\sigma^{i}\int_{\exp[-(\frac{x-\mu}{\sigma})]}^{+\infty}\ln^{i}(z)\exp(-z)\,{\rm
d}z$ $\displaystyle=$
$\displaystyle\sum_{i=0}^{k}\binom{k}{i}\mu^{k-i}\sigma^{i}\,I\biggl{(}i;\exp\Big{[}-\Big{(}\frac{x-\mu}{\sigma}\Big{)}\Big{]},+\infty\biggr{)},$

where $k\in\mathbb{N}\cup\\{0\\}$ and $I(i;a,b)$ is as in (2.4). By combining
(3.1) and (3.2), and by using the relations

$\displaystyle I(0;a,+\infty)$ $\displaystyle=$ $\displaystyle\exp(-a),$
$\displaystyle I(1;a,+\infty)$ $\displaystyle=$
$\displaystyle-\exp(-a)\ln(a)-\Gamma(0,a),$

we get formula (2).

###  3.2 Stochastic representation

Suppose $Y_{k}$ has a weighted Gumbel distribution with parameters
$\mu\in\mathbb{R}$ and $\sigma>0$. That is, if $Y\sim F_{\rm
G}(\cdot;\mu,\sigma)$, by (3.2), $Y_{k}$ has CDF given by, for each
$x\in\mathbb{R}$ and $k=0,1,2,\ldots$,

$\displaystyle F_{Y_{k}}(x)={\mathbb{E}\big{(}Y^{k}\mathds{1}_{Y\leq
x}\big{)}\over\mathbb{E}(Y^{k})}=\frac{\displaystyle\sum_{i=0}^{k}\binom{k}{i}\mu^{k-i}\sigma^{i}\,I\biggl{(}i;\exp\Big{[}-\Big{(}\frac{x-\mu}{\sigma}\Big{)}\Big{]},+\infty\biggr{)}}{\displaystyle\sum_{i=0}^{k}\binom{k}{i}\mu^{k-i}\sigma^{i}\,I(i;0,+\infty)},$
(3.3)

where $I(k;a,b)$ is the incomplete moments defined in (2.4). Note that
$F_{Y_{0}}(x)=F_{\rm G}(x;\mu,\sigma)$.

Let $W$ be a discrete random variable, so that $W=1$ or $W=2$ or $W=3$, each
with probability

$\displaystyle p_{1}={2\over Z_{\delta}},\quad
p_{2}=-{2(\mu+\sigma\gamma)\delta\over Z_{\delta}},\quad
p_{3}={[\sigma^{2}{\pi^{2}\over 6}+(\mu+\sigma\gamma)^{2}]\delta^{2}\over
Z_{\delta}},$

respectively, where $[\delta>0\,\text{and}\,\mu+\sigma\gamma<0]$ or
$[\delta<0\,\text{and}\,\mu+\sigma\gamma>0]$, and $Z_{\delta}$ is as in (2.2).
It is straightforward see that $p_{1}+p_{2}+p_{3}=1$.

Assume that

$T=\sum_{k=1}^{3}{1\over k}\,Y_{k-1}\delta_{W,k}\delta_{W,l},\quad l=1,2,3,$

and that $W$ is independent of $Y_{l}$, for each $l=1,2,3$. Here
$\delta_{x,y}$ is the Kronecker delta function, i.e., this function is 1 if
the variables are equal, and 0 otherwise.

######  Proposition 3.1.

The following holds

$\displaystyle X=WT\ \ \text{if and only if}\ \ X\sim F_{\rm
BG}(\cdot;\mu,\sigma,\delta).$ (3.4)

###### Proof.

By Law of total probability and by independence, we get

$\displaystyle\mathbb{P}(X\leqslant x)=\mathbb{P}(WT\leqslant x)$
$\displaystyle=\sum_{l=1}^{3}\mathbb{P}(WT\leqslant x|W=l)\mathbb{P}(W=l)$
$\displaystyle=\sum_{l=1}^{3}\mathbb{P}(Y_{l-1}\leqslant
x|W=l)\mathbb{P}(W=l)=\sum_{l=1}^{3}F_{Y_{l-1}}(x)\,p_{l}.$

By using the CDF of the weighted Gumbel distribution, given in (3.3), and by
definitions of $p_{l}$’s and $Z_{\delta}$, the above expression is

$\displaystyle={2\over
Z_{\delta}}\,F_{Y_{0}}(x)-{2(\mu+\sigma\gamma)\delta\over
Z_{\delta}}\,F_{Y_{1}}(x)+{[\sigma^{2}{\pi^{2}\over
6}+(\mu+\sigma\gamma)^{2}]\delta^{2}\over Z_{\delta}}\,F_{Y_{2}}(x)$
$\displaystyle={2\over Z_{\delta}}\,F_{\rm
G}(x;\mu,\sigma)-{2(\mu+\sigma\gamma)\delta\over
Z_{\delta}}\,\frac{\displaystyle\sum_{i=0}^{1}\binom{1}{i}\mu^{1-i}\sigma^{i}\,I\biggl{(}i;\exp\Big{[}-\Big{(}\frac{x-\mu}{\sigma}\Big{)}\Big{]},+\infty\biggr{)}}{\displaystyle\sum_{i=0}^{1}\binom{1}{i}\mu^{1-i}\sigma^{i}\,I(i;0,+\infty)}$
$\displaystyle\quad+{[\sigma^{2}{\pi^{2}\over
6}+(\mu+\sigma\gamma)^{2}]\delta^{2}\over
Z_{\delta}}\,\frac{\displaystyle\sum_{i=0}^{2}\binom{2}{i}\mu^{2-i}\sigma^{i}\,I\biggl{(}i;\exp\Big{[}-\Big{(}\frac{x-\mu}{\sigma}\Big{)}\Big{]},+\infty\biggr{)}}{\displaystyle\sum_{i=0}^{2}\binom{2}{i}\mu^{2-i}\sigma^{i}\,I(i;0,+\infty)}$
$\displaystyle={\displaystyle\big{[}2-\delta\mu(2-\delta\mu)\big{]}\exp\biggl{\\{}-\exp\Big{[}-\Big{(}\frac{x-\mu}{\sigma}\Big{)}\Big{]}\biggr{\\}}+\delta^{2}\sigma^{2}I\biggl{(}2;\exp\Big{[}-\Big{(}\frac{x-\mu}{\sigma}\Big{)}\Big{]},+\infty\biggr{)}\over
1+\delta^{2}\sigma^{2}{\pi^{2}\over 6}+(\delta\mu+\delta\sigma\gamma-1)^{2}}$
$\displaystyle\quad+{\displaystyle
2\delta(1-\delta\mu)\Biggl{\\{}{(x-\mu)}\,\exp\biggl{\\{}-\exp\Big{[}-\Big{(}\frac{x-\mu}{\sigma}\Big{)}\Big{]}\biggr{\\}}-\sigma\Gamma\biggl{(}0,\exp\Big{[}-\Big{(}\frac{x-\mu}{\sigma}\Big{)}\Big{]}\biggr{)}\Biggr{\\}}\over
1+\delta^{2}\sigma^{2}{\pi^{2}\over 6}+(\delta\mu+\delta\sigma\gamma-1)^{2}}$
$\displaystyle\stackrel{{\scriptstyle\eqref{cdfbgev2-1}}}{{=}}F_{\rm
BG}(x;\mu,\sigma,\delta),$

where $\Gamma(a,b)$ is the upper incomplete gamma function.

Then the statement in (3.4) follows. ∎

###  3.3 Rate of a random variable with BG distribution

Following Klugman et al. (1998) [7], for a continuous random variable $X$ with
density function $f_{X}(x)$, the rate of a random variable is given by

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

A simple computation shows that $\tau_{{\rm
BG}(\mu,\sigma,\delta)}={1/\sigma}.$

In what follows we present some comparisons between the rates of random
variables with known distributions: Inverse-gamma, Log-normal, Generalized-
Pareto, BWeibull (see Vila et al. 2020 [18]), BGamma (see Vila et al. 2020
[17]), BG, exponential and Normal;

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

In other words, far enough out in the tail, every BG distribution looks like
an exponential distribution.

###  3.4 Bimodality

######  Proposition 3.2.

A point $x\in\mathbb{R}$ is a mode of the BG density (2.1) if it is a root of
the following non-polynomial function:

$\displaystyle
g(x)={1\over\sigma}\biggl{\\{}\exp\Big{[}-\Big{(}{x-\mu\over\sigma}\Big{)}\Big{]}-1\biggl{\\}}-{2\delta(1-\delta
x)\over(1-\delta x)^{2}+1}.$ (3.5)

###### Proof.

Taking the derivative of $f_{\rm BG}(x;\mu,\sigma)$ with respect to $x$, we
have

$\displaystyle{f^{\prime}_{\rm BG}(x;\mu,\sigma)}=f_{\rm
BG}(x;\mu,\sigma)g(x).$ (3.6)

Hence, the proof follows. ∎

Let $\mathcal{C}$ be the set formed for all
$(\mu,\sigma,\delta)\in\mathbb{R}\times(0,+\infty)\times\mathbb{R}$ such that
the following hold:

$\displaystyle\delta>\max\biggl{\\{}1,{1\over\sigma}\Big{[}\exp\Big{(}{\mu\over\sigma}\Big{)}-1\Big{]}\biggr{\\}},$
(3.7)
$\displaystyle{2\delta(1+\delta)\over(1+\delta)^{2}+1}<{1\over\sigma}\Big{\\{}\exp\Big{(}{1+\mu\over\sigma}\Big{)}-1\Big{\\}},$
(3.8)
$\displaystyle{2\delta(1-2\delta)\over(1-2\delta)^{2}+1}<{1\over\sigma}\biggl{\\{}\exp\Big{[}-\Big{(}{2-\mu\over\sigma}\Big{)}\Big{]}-1\biggl{\\}},$
(3.9)
$\displaystyle{2\delta(1-3\delta)\over(1-3\delta)^{2}+1}>{1\over\sigma}\biggl{\\{}\exp\Big{[}-\Big{(}{3-\mu\over\sigma}\Big{)}\Big{]}-1\biggl{\\}}.$
(3.10)

######  Remark 2.

By considering $\mu=\sigma=1$ and $\delta>{\rm e}-1$, we have that
$(\mu,\sigma,\delta)\in\mathcal{C}$. That is, the set $\mathcal{C}$ is non-
empty.

######  Lemma 3.3.

If $(\mu,\sigma,\delta)\in\mathcal{C}$ then the function $g(x)$ has at least
three distinct real roots.

###### Proof.

Since $(\mu,\sigma,\delta)\in\mathcal{C}$, a simple observation in the
definition (3.5) of $g(x)$ shows that

  * •

$\displaystyle
g(-1)={1\over\sigma}\Big{\\{}\exp\Big{(}{1+\mu\over\sigma}\Big{)}-1\Big{\\}}-{2\delta(1+\delta)\over(1+\delta)^{2}+1}>0,$
because of condition (3.8);

  * •

$\displaystyle
g(0)={1\over\sigma}\Big{\\{}\exp\Big{(}{\mu\over\sigma}\Big{)}-1\Big{\\}}-\delta<0$,
because of condition (3.7);

  * •

$\displaystyle
g(2)={1\over\sigma}\biggl{\\{}\exp\Big{[}-\Big{(}{2-\mu\over\sigma}\Big{)}\Big{]}-1\biggl{\\}}-{2\delta(1-2\delta)\over(1-2\delta)^{2}+1}>0$,
because of condition (3.9);

  * •

$\displaystyle
g(3)={1\over\sigma}\biggl{\\{}\exp\Big{[}-\Big{(}{3-\mu\over\sigma}\Big{)}\Big{]}-1\biggl{\\}}-{2\delta(1-3\delta)\over(1-3\delta)^{2}+1}<0$,
because of condition (3.10).

Since $g(x)$ is a continuous real-valued function, by Intermediate Value
Theorem, there are some points $r_{1},r_{2},r_{3}$, with $-1<r_{1}<0$,
$0<r_{2}<2$ and $2<r_{3}<3$, so that $g(r_{i})=0$, $i=1,2,3$. ∎

Let $\mathcal{D}$ be the set formed for all $x\in\mathbb{R}$ such that

$\displaystyle 2\delta^{2}\,{(1-\delta x)^{2}-1\over(1-\delta x)^{2}+1}$
$\displaystyle<$
$\displaystyle-{1\over\sigma^{2}}\exp\Big{[}-\Big{(}{x-\mu\over\sigma}\Big{)}\Big{]}\quad\text{where}\
(\mu,\sigma,\delta)\in\mathcal{C}.$

######  Remark 3.

The set $\mathcal{D}$ is non-empty. To see this just take $\mu=\sigma=1$ and
$\delta=2$ $(>{\rm e}-1)$. By Remark 2, $(\mu,\sigma,\delta)\in\mathcal{C}$.
In this case we have $\mathcal{D}=(0.132178,0.937349)$.

######  Lemma 3.4.

If $(\mu,\sigma,\delta)\in\mathcal{C}$ then the function $g(x)$ has no root
outside the interval $(-1,3)$.

Futhermore, if $-1<r_{1}<0$, $0<r_{2}<2$ and $2<r_{3}<3$ are the roots of
$g(x)$ found in Lemma 3.3, so that $r_{2}\in\mathcal{D}\subset(0,2)$, then
these roots are unique.

###### Proof.

Taking the derivative of $g(x)$ with respect to $x$, we have

$\displaystyle
g^{\prime}(x)=-{1\over\sigma^{2}}\exp\Big{[}-\Big{(}{x-\mu\over\sigma}\Big{)}\Big{]}-2\delta^{2}\,{(1-\delta
x)^{2}-1\over(1-\delta x)^{2}+1}.$

Since, for $x\leqslant 0$ or $x\geqslant 2/\delta$, $2\delta^{2}\,{[(1-\delta
x)^{2}-1]/[(1-\delta x)^{2}+1]}\geqslant 0,$ we obtain

$\displaystyle
g^{\prime}(x)\leqslant-{1\over\sigma^{2}}\exp\Big{[}-\Big{(}{x-\mu\over\sigma}\Big{)}\Big{]}<0\quad\text{for}\
x\leqslant 0\ \text{ou}\ x\geqslant 2/\delta.$

In other words, outside of the interval $(-1,3)$ the function $g(x)$ has no
zeros because this one is decreasing on $(-\infty,0]\cup[2/\delta,+\infty)$
with $2/\delta<2$ (by condition (3.7)). This implies that the roots
$-1<r_{1}<0$ and $2<r_{3}<3$, of $g(x)$, are unique in at their respective
intervals. Finally, since $r_{2}\in\mathcal{D}$ we have that the function
$g(x)$ crosses the abscissa axis at the only point $r_{2}$. Therefore, we
conclude that the root $r_{2}$ is also unique in $\mathcal{D}\subset(0,2)$. ∎

######  Remark 4.

Let $\mu=\sigma=1$ and $\delta=2$ $(>{\rm e}-1)$. By Remark 3,
$\mathcal{D}=(0.132178,0.937349)$. Computationally it can be verified that
$-1<r_{1}=-0.0896138<0,0<r_{2}=0.389792<2$ and $2<r_{3}=2.79117<3$ are the
only roots of $g(x)$, with $r_{2}\in\mathcal{D}$.

######  Theorem 1 (Bimodality).

If $r_{2}\in\mathcal{D}$ then the BG distribution (2.1) is bimodal.

###### Proof.

By Lemma 3.4 the function $g(x)$ in (3.5) has exactly three roots
$r_{1},r_{2},r_{3}$ so that $r_{1}<r_{2}<r_{3}$. Since $f_{\rm
BG}(x;\mu,\sigma)\to 0$ as $x\to\pm\infty$, it follows that the BG
distribution (2.1) increases on the intervals $(-\infty,r_{1})$ and
$(r_{2},r_{3})$, and decreases on $(r_{1},r_{2})$ and $(r_{3},+\infty)$. That
is, $r_{1}$ and $r_{3}$ are two maximum points and $r_{2}$ is the unique
minimum point. Therefore, the bimodality of $f_{\rm BG}(x;\mu,\sigma)$ is
guaranteed. ∎

###  3.5 The hazard function

The survival and hazard functions, denoted by SF and HR of the BG is given by
$S_{\rm BG}(x;\mu,\sigma,\delta)=1-F_{\rm BG}(x;\mu,\sigma,\delta)$ and
$H_{\rm BG}(x;\mu,\sigma,\delta)={f_{\rm BG}(x;\mu,\sigma,\delta)/[1-F_{\rm
BG}(x;\mu,\sigma,\delta)]}.$

Considering $r_{1},r_{2},r_{3}$ the three distinct real roots, of function
$g(x)$, found in Lemma 3.3, we state the following result.

######  Proposition 3.5.

If $r_{2}\in\mathcal{D}$ then the HR $H_{\rm BG}(x;\mu,\sigma,\delta)$ of the
BG distribution (2.1) has the following monotonicity properties:

  * 1)

It is increasing for each $x<r_{1}$ or $r_{2}<x<r_{3}$,

  * 2)

It is decreasing for each $x\in\mathcal{D}$.

###### Proof.

As a sub-product of the proof of Theorem 1, note that the density $f_{\rm
BG}(x;\mu,\sigma,\delta)$ is increasing on the intervals $x<r_{1}$ or
$r_{2}<x<r_{3}$. Since $1-F_{\rm BG}(x;\mu,\sigma,\delta)$ is a decreasing
function, we have that the hazard function on $x<r_{1}$ or $r_{2}<x<r_{3}$ is
the product of two increasing and nonnegative functions, then the proof of
first item follows.

By Glaser (1980), for to prove the second item, it is enough to prove that the
function $G_{X}(x)$, defined as

$\displaystyle G_{X}(x)=-{f_{\rm BG}^{\prime}(x;\mu,\sigma,\delta)\over f_{\rm
BG}(x;\mu,\sigma,\delta)}\stackrel{{\scriptstyle\eqref{f-der-
relation}}}{{=}}-g(x),$

is decreasing for each $x\in\mathcal{D}$. But this is immediate since $g(x)$
is increasing on this interval. ∎

###  3.6 Moment-Generating Function

######  Theorem 2.

If $X\sim F_{\rm BG}(\cdot;\mu,\sigma,\delta)$, $t<\min\\{0,-m/\sigma\\}$ and
$m\in\mathbb{N}\cup\\{0\\}$, we have

$\displaystyle\mathbb{E}\big{[}X^{m}\exp(tX)\big{]}$
$\displaystyle=(-1)^{m+2}\exp(t\mu)\,\dfrac{\delta^{2}\sigma^{m+2}\,\Gamma^{(m+2)}(1-\sigma
t)+\delta\sigma^{m+1}\big{[}2-\delta\mu(m+2)\big{]}\,\Gamma^{(m+1)}(1-\sigma
t)}{1+\delta^{2}\sigma^{2}{\pi^{2}\over
6}+(\delta\mu+\delta\sigma\gamma-1)^{2}}$
$\displaystyle+\dfrac{\exp(t\mu)\,\sum_{i=0}^{m}(-1)^{i}\sigma^{i}\mu^{m-i}\big{[}2\binom{m}{i}-2\delta\mu\binom{m+1}{i}+\delta^{2}\mu^{2}\binom{m+2}{i}\big{]}\Gamma^{(i)}(1-\sigma
t)}{1+\delta^{2}\sigma^{2}{\pi^{2}\over
6}+(\delta\mu+\delta\sigma\gamma-1)^{2}},$

where $\Gamma^{(i)}(x)={{\rm d}^{i}\Gamma(x)/{\rm d}x^{i}}$.

###### Proof.

Taking $g(X)=X^{m}\exp(tX)$ in Remark 1,

$\displaystyle\mathbb{E}\big{[}X^{m}\exp(t\,X)\big{]}={1\over
Z_{\delta}}\,\big{\\{}2\,\mathbb{E}\big{[}Y^{m}\exp(t\,Y)\big{]}-2\delta\,\mathbb{E}\big{[}Y^{m+1}\exp(t\,Y)\big{]}+\delta^{2}\,\mathbb{E}\big{[}Y^{m+2}\exp(t\,Y)\big{]}\big{\\}},$
(3.11)

where $Y\sim F_{\rm G}(\cdot;0,\mu,\sigma)$ and $Z_{\delta}$ are as in (1.2)
and (2.2), respectively.

Taking the change of variables $z=\exp[-({y-\mu})/{\sigma}],\ {\rm
d}z=-({z/\sigma}){\rm d}y,\ y=\mu-\sigma\ln(z),$ and using a binomial
expansion we get, for $k=m,m+1,m+2$,

$\displaystyle\mathbb{E}\big{[}Y^{k}\exp(tY)\big{]}$ $\displaystyle=$
$\displaystyle\int_{0}^{+\infty}[\mu-\sigma\ln(z)]^{k}z^{(1-\sigma
t)-1}\exp(-z)\,{\rm d}z,\quad k\in\mathbb{N}\cup\\{0\\}$ (3.12)
$\displaystyle=$ $\displaystyle\exp(\mu
t)\sum_{i=0}^{k}\binom{k}{i}(-1)^{i}\mu^{k-i}\sigma^{i}\int_{0}^{+\infty}\ln^{i}(z)\,z^{(1-\sigma
t)-1}\exp(-z)\,{\rm d}z.$

We claim that

$\displaystyle\int_{0}^{+\infty}\ln^{i}(z)\,z^{(1-\sigma t)-1}\exp(-z)\,{\rm
d}z=\Gamma^{(i)}(1-\sigma t),\quad i=0,1,\ldots,k;\ k=m,m+1,m+2.$ (3.13)

Note that, by combining (3.13), (3.12) and (3.11), the proof of theorem
follows.

Indeed, since
${\partial^{j}z^{\alpha-1}\over\partial\alpha^{j}}=\ln^{j}(z)\,z^{\alpha-1}$,
$j=1,2,\ldots i$,

$\displaystyle\int_{0}^{+\infty}\ln^{i}(z)\,z^{(1-\sigma t)-1}\exp(-z)\,{\rm
d}z=\int_{0}^{+\infty}\dfrac{\partial^{j}}{\partial(1-\sigma
t)^{j}}\big{[}\ln^{i-j}(z)\,z^{(1-\sigma t)-1}\exp(-z)\big{]}\,{\rm d}z.$
(3.14)

Let us note that the following conditions hold true.

  * •

The derivatives

$\displaystyle\dfrac{\partial^{j}}{\partial(1-\sigma
t)^{j}}\big{[}\ln^{i-j}(z)\,z^{(1-\sigma
t)-1}\exp(-z)\big{]}=\ln^{i}(z)\,z^{(1-\sigma t)-1}\exp(-z)$

exists and are continuous in $1-\sigma t$ for all $z$ and all $1-\sigma t$,
for $j=1,2,\ldots i$.

  * •

By using the following inequalities

$\displaystyle|\ln(z)|\leqslant{1\over
z}\,\mathds{1}_{\\{x>0:x<1\\}}(z)+z\,\mathds{1}_{\\{x>0:x\geqslant 1\\}}(z),$
$\displaystyle|\ln^{i}(z)|\leqslant{1\over
z^{i-1}}\,\mathds{1}_{A_{i}}(z)+z^{i-1}\,\mathds{1}_{B_{i}}(z),\quad
i=2,3,\ldots,k,$

where $A_{i}=\big{\\{}x>0:x<{\rm e}^{iW({i-1\over i})/(i-1)}\big{\\}}$,
$B_{i}=\big{\\{}x>0:x\geqslant{\rm e}^{iW({i-1\over i})/(i-1)}\big{\\}}$, and
$W$ is the product logarithm function (or Lambert $W$ function), we get

$\displaystyle\biggl{|}\dfrac{\partial^{j}}{\partial(1-\sigma
t)^{j}}\big{[}\ln^{i-j}(z)\,z^{(1-\sigma
t)-1}\exp(-z)\big{]}\biggr{|}=\big{|}\ln^{i}(z)\big{|}\,z^{(1-\sigma
t)-1}\exp(-z)$ $\displaystyle\leqslant\begin{cases}z^{(1-\sigma
t)-1}\exp(-z)&\text{if}\ i=0,\\\\[5.69046pt] \big{[}z^{(1-\sigma
t)-2}\,\mathds{1}_{\\{x>0:x<1\\}}(z)+z^{(1-\sigma
t)}\,\mathds{1}_{\\{x>0:x\geqslant 1\\}}(z)\big{]}\exp(-z)&\text{if}\
i=1,\\\\[5.69046pt] \big{[}z^{(1-\sigma
t)-i}\,\mathds{1}_{A_{i}}(z)+z^{(1-\sigma
t)+i-2}\,\mathds{1}_{B_{i}}(z)\big{]}\exp(-z)&\text{if}\
i=2,3,\ldots,k,\end{cases}\eqqcolon G(z,1-\sigma t).$

Furthermore, the integral

$\displaystyle\int_{0}^{+\infty}G(z,1-\sigma t)\,{\rm
d}z\leqslant\begin{cases}\Gamma(1-\sigma t)&\text{if}\ i=0,\
t<1/\sigma,\\\\[5.69046pt] \Gamma(-\sigma t)+\Gamma(2-\sigma t)&\text{if}\
i=1,\ t<0,\\\\[5.69046pt] \Gamma(2-\sigma t-i)+\Gamma(i-\sigma t)&\text{if}\
i=2,3,\ldots,k,\ t<(2-k)/\sigma,\end{cases}$

is finite, for $k=m,m+1,m+2$.

  * •

The integral $\int_{0}^{+\infty}\ln^{i-j}(z)\,z^{(1-\sigma t)-1}\exp(-z)\,{\rm
d}z$ exists because of last item above.

Under the above three conditions, by Leibniz integral rule, we can interchange
the derivative with the integral in (3.14). Hence, for $j=1,2,\ldots i$,

$\displaystyle\int_{0}^{+\infty}\ln^{i}(z)\,z^{(1-\sigma t)-1}\exp(-z)\,{\rm
d}z=\dfrac{{\rm d}^{j}}{{\rm d}(1-\sigma
t)^{j}}\int_{0}^{+\infty}\ln^{i-j}(z)\,z^{(1-\sigma t)-1}\exp(-z)\,{\rm d}z.$

Letting $j=i$, we obtain

$\displaystyle\int_{0}^{+\infty}\ln^{i}(z)\,z^{(1-\sigma t)-1}\exp(-z)\,{\rm
d}z$ $\displaystyle=$ $\displaystyle\dfrac{{\rm d}^{i}}{{\rm d}(1-\sigma
t)^{i}}\int_{0}^{+\infty}z^{(1-\sigma t)-1}\exp(-z)\,{\rm d}z$
$\displaystyle=$ $\displaystyle\Gamma^{(i)}(1-\sigma t),\quad i=0,1,\ldots,$

where in the last line we used the definition of gamma function,
$\Gamma(\alpha)=\int_{0}^{+\infty}z^{\alpha-1}\exp(-z)\,{\rm d}z$. Then the
claimed (3.13) follows.

Thus, this complete the proof of theorem. ∎

By taking $m=0$ in Theorem 2 we obtain a clossed expression for the moment-
generating function of BGumbel distribution.

######  Corollary 3.6.

If $X\sim F_{\rm BG}(\cdot;\mu,\sigma,\delta)$ then, the moment-generating
function $M_{X}(t)=\mathbb{E}[\exp(tX)]$ with $t<0$, is given by

$\displaystyle M_{X}(t)={\exp(\mu t)\Gamma(1-\sigma
t)\big{[}2-2\mu\delta+\mu^{2}\delta^{2}+2\sigma\delta(1-\mu\delta)\psi(1-\sigma
t)+{\sigma^{2}\delta^{2}\over\Gamma(1-\sigma t)}\,\Gamma^{(2)}(1-\sigma
t)\big{]}\over 1+\delta^{2}\sigma^{2}{\pi^{2}\over
6}+(\delta\mu+\delta\sigma\gamma-1)^{2}},$

where $\psi(x)=\Gamma^{\prime}(x)/\Gamma(x)$ is the digamma function and
$\Gamma^{(2)}(1-\sigma t)={{\rm d}^{2}\Gamma(x)/{\rm d}x^{2}}$.

######  Remark 5.

Letting $\delta=0$ in Corollary 3.6 we obtain the known formula
$M_{Y}(t)=\exp(\mu t)\Gamma(1-\sigma t)$, where $Y\stackrel{{\scriptstyle
d}}{{=}}X\sim F_{\rm G}(\cdot;0,\mu,\sigma)$. Note that, by a by-product of
proof of Theorem 2, in this case it is sufficient take $t<1/\sigma$ instead
$t<0$.

######  Remark 6.

The characteristic function of $X\sim F_{\rm BG}(\cdot;\mu,\sigma,\delta)$,
denoted by $\phi_{X}(t)$, can be obtained from the moment-generating function
by the relation $M_{X}(t)=\phi_{X}(-it)$.

###  3.7 Moments

######  Theorem 3.

If $X\sim F_{\rm BG}(\cdot;\mu,\sigma,\delta)$ then

$\displaystyle\mathbb{E}(X^{k})={\displaystyle\delta^{2}\sigma^{k+2}I(k+2;0,+\infty)-\delta\sigma^{k+1}\big{[}2-\delta\mu(k+2)\big{]}I(k+1;0,+\infty)\over
1+\delta^{2}\sigma^{2}{\pi^{2}\over 6}+(\delta\mu+\delta\sigma\gamma-1)^{2}}$
$\displaystyle+{\sum_{i=0}^{k}\sigma^{i}\mu^{k-i}\big{[}2\binom{k}{i}-2\delta\mu\binom{k+1}{i}+\delta^{2}\mu^{2}\binom{k+2}{i}\big{]}I(i;0,+\infty)\over
1+\delta^{2}\sigma^{2}{\pi^{2}\over 6}+(\delta\mu+\delta\sigma\gamma-1)^{2}},$

where $I(i;a,b)$ is defined in (2.4).

###### Proof.

From Remark 1 with $g(X)=X^{k}$, we obtain

$\displaystyle\mathbb{E}(X^{k})$ $\displaystyle=$ $\displaystyle{1\over
Z_{\delta}}\,\big{[}2\,\mathbb{E}(Y^{k})-2\delta\,\mathbb{E}(Y^{k+1})+\delta^{2}\,\mathbb{E}(Y^{k+2})\big{]},$
(3.15)

where $Y\sim F_{\rm G}(\cdot;0,\mu,\sigma)$ is as in (1.2).

By taking $x\longrightarrow+\infty$ in (3.2), Lebesgue dominated convergence
theorem gives

$\displaystyle\mathbb{E}(Y^{k})=\lim_{x\to+\infty}\mathbb{E}\big{(}Y^{k}\mathds{1}_{Y\leq
x}\big{)}=\sum_{i=0}^{k}\binom{k}{i}\mu^{k-i}\sigma^{i}\,I(i;0,+\infty).$
(3.16)

By combining (3.15) and (3.16), the proof follows. ∎

######  Corollary 3.7.

If $X\sim F_{\rm BG}(\cdot;\mu,\sigma,\delta)$ then

$\displaystyle\mathbb{E}(X)$
$\displaystyle={\delta^{2}\sigma^{3}\big{[}2\zeta(3)+\gamma^{3}+{\gamma\pi^{2}\over
2}\big{]}-\delta\sigma^{2}(2-3\delta\mu)\big{(}\gamma^{2}+{\pi^{2}\over
6}\big{)}+\mu\big{[}2-\delta\mu(2-\delta\mu)\big{]}+\sigma\big{[}2-\delta\mu(4-3\delta\mu)\big{]}\gamma\over
1+\delta^{2}\sigma^{2}{\pi^{2}\over 6}+(\delta\mu+\delta\sigma\gamma-1)^{2}},$
(3.17) $\displaystyle\mathbb{E}(X^{2})$
$\displaystyle={\delta^{2}\sigma^{4}\big{[}8\gamma\zeta(3)+\gamma^{4}+\gamma^{2}\pi^{2}+{3\pi^{4}\over
20}\big{]}-2\delta\sigma^{3}(1-2\delta\mu)\big{[}2\zeta(3)+\gamma^{3}+{\gamma\pi^{2}\over
2}\big{]}+\mu^{2}\big{[}2-\delta\mu(2-\delta\mu)]\over
1+\delta^{2}\sigma^{2}{\pi^{2}\over 6}+(\delta\mu+\delta\sigma\gamma-1)^{2}}$
$\displaystyle+{2\sigma\mu\big{[}2-\delta\mu(3-2\delta\mu)\big{]}\gamma+2\sigma^{2}\big{[}1-3\delta\mu(1-\delta\mu)\big{]}\big{(}\gamma^{2}+{\pi^{2}\over
6}\big{)}\over 1+\delta^{2}\sigma^{2}{\pi^{2}\over
6}+(\delta\mu+\delta\sigma\gamma-1)^{2}},$ (3.18)
$\displaystyle\mathbb{E}(X^{3})$
$\displaystyle={\delta^{2}\sigma^{5}\big{[}20\gamma^{2}\zeta(3)+{10\pi^{2}\zeta(3)\over
3}+24\zeta(5)+\gamma^{5}+{5\gamma^{3}\pi^{2}\over 3}+{3\gamma\pi^{4}\over
4}\big{]}\over 1+\delta^{2}\sigma^{2}{\pi^{2}\over
6}+(\delta\mu+\delta\sigma\gamma-1)^{2}}$
$\displaystyle+{\mu^{3}\big{[}2-\delta\mu(2-\delta\mu)\big{]}+\sigma\mu^{2}\big{[}6-\delta\mu(8-5\delta\mu)\big{]}\gamma+2\sigma^{2}\mu\big{[}3-\delta\mu(6-5\delta\mu)\big{]}(\gamma^{2}+{\pi^{2}\over
6})\over 1+\delta^{2}\sigma^{2}{\pi^{2}\over
6}+(\delta\mu+\delta\sigma\gamma-1)^{2}}$
$\displaystyle+{2\sigma^{3}\big{[}1-\delta\mu(4-5\delta\mu)\big{]}\big{[}2\zeta(3)+\gamma^{3}+{\gamma\pi^{2}\over
2}\big{]}-\delta\sigma^{4}(2-5\delta\mu)\big{[}8\gamma\zeta(3)+\gamma^{4}+\gamma^{2}\pi^{2}+{3\pi^{4}\over
20}\big{]}\over 1+\delta^{2}\sigma^{2}{\pi^{2}\over
6}+(\delta\mu+\delta\sigma\gamma-1)^{2}},$

where $\gamma$ is the Euler-Mascheroni constant and $\zeta(s)$ is the Riemann
zeta function.

###### Proof.

By combining Theorem 3 with the following values for the improper integral
$I(k;0,+\infty)$ in (2.4), with $k=0,1,2,3,4$ and $5$,

$\displaystyle\begin{array}[]{lllll}I(0;0,+\infty)&=&1,\\\\[2.84544pt]
I(1;0,+\infty)&=&\gamma,\\\\[2.84544pt]
I(2;0,+\infty)&=&\gamma^{2}+{\pi^{2}\over 6},\\\\[2.84544pt]
I(3;0,+\infty)&=&2\zeta(3)+\gamma^{3}+{\gamma\pi^{2}\over 2},\\\\[2.84544pt]
I(4;0,+\infty)&=&8\gamma\zeta(3)+\gamma^{4}+\gamma^{2}\pi^{2}+{3\pi^{4}\over
20},\\\\[2.84544pt]
I(5;0,+\infty)&=&20\gamma^{2}\zeta(3)+{10\pi^{2}\zeta(3)\over
3}+24\zeta(5)+\gamma^{5}+{5\gamma^{3}\pi^{2}\over 3}+{3\gamma\pi^{4}\over
4},\end{array}$ (3.25)

the proof follows. ∎

For the variance of the BGumbel distribution

${\rm Var}(X)=\mathbb{E}(X^{2})-\big{[}\mathbb{E}(X)\big{]}^{2},$ (3.26)

just replace (3.17) and (3.7) in (3.26).

######  Remark 7 (Standardized moments).

Let $X\sim F_{\rm BG}(\cdot;\mu,\sigma,\delta)$ with $\mathbb{E}(X)$ and
$\sqrt{{\rm Var}(X)}>0$, both finite. Newton’s binomial expansion gives

$\displaystyle\mathbb{E}\biggl{[}{X-\mathbb{E}(X)\over\sqrt{{\rm
Var}(X)}}\,\biggr{]}^{n}={(-1)^{n}\big{[}\mathbb{E}(X)\big{]}^{n}\over\big{[}{\rm
Var}(X)\big{]}^{n/2}}\,\sum_{k=0}^{n}\binom{n}{k}(-1)^{-k}\,\big{[}\mathbb{E}(X)\big{]}^{-k}\,\mathbb{E}(X^{k}).$
(3.27)

By combining (3.27), Theorem 3, identities in (3.25) and the following
relation

$\displaystyle\begin{array}[]{lllll}I(6;0,+\infty)=20\gamma(2\gamma^{2}+\pi^{2})\zeta(3)+40[\zeta(3)]^{2}+144\gamma\zeta(5)+\gamma^{6}+{5\gamma^{4}\pi^{2}\over
2}+{9\gamma^{2}\pi^{4}\over 4}+{61\pi^{6}\over 168},\end{array}$

closed formulas for skewness and kurtosis of the random variable $X\sim F_{\rm
BG}(\cdot;\mu,\sigma,\delta)$ it is possible to obtain.

##  4 Graphical illustrations and simulation results

###  4.1 Graphical illustration

Here, the flexibility of the BG model is shown. Note that the BG model can be
unimodal or bimodal. Figures 1 and 2 show how the BG density function is
influenced by the shape parameters $\delta$ and $\mu$, since $\sigma$ is a
scale parameter. Figure 1 shows that the bimodality of the BG model is more
evident the greater the absolute value of $\delta$. The density tends to
become unimodal when the values of $\delta$ and $\mu$ are far apart. In the
bimodal case, the left mode is smaller and the opposite for the right mode as
$\mu$ increases, see Figure 2. All graphics were generated with our bgumbel
package, [1].

Figure 1: BG density graphs for $\mu=1$, $\sigma=2$ and $\delta$ varying as
shown in the caption.  Figure 2: BG density graphs for $\delta=1$, $\sigma=2$
and $\mu$ varying as shown in the caption.

###  4.2 Simulation results

In this section, we compare the sample mean and variance with the population
mean and variance (3.17) and (3.26). In addition, we compared the cumulative
distribution function given in (2) with the empirical distribution. For that,
we present a Markov chain Monte Carlo (MCMC) simulation study with various
sample sizes. In the MCMC simulation method we considered the Metropolis-
Hastings method as pseudo-random numbers (PRN) generator [10],[6]. We chose
this method instead of the inverse transform method, because the BG cumulative
distribution (3.1) does not admit analytical inversion. To generate PRN was
used the MCMCpack package [8] available in the R program [14]. The mean,
variance and the corresponding bias of the sample estimates were computed over
$10^{5}$ iterations of MCMC. These results are shown in Tables 1. The
performance of the MCMC algorithm was tested by varying the number of
iterations, n. Graphically, the convergence of the Markov Chain is shown in
Figure 3, for the first set of parameters in Table 1. For the other parameters
sets , in the Table 1, the performance of the algorithm is illustrated in
Appendix A, Figures 6 \- 14. The value of n varies from 1000 to 100000 for
each parameters set. In all these figures the transition probability density
function (TPDF) ,defined by the MCMCpack package, is in blue and the
estacionary density in red. We can conclude that the algorithm, [1], is
efficient for n greater than 100000.

Table 1:  Comparison between the sample moments and population moments ($n=10^{5}$). $\mu$ | $\sigma$ | $\delta$ | sample mean | $E(X)$ | Bias(mean) | variance sample | $Var(X)$ | Bias (variance)
---|---|---|---|---|---|---|---|---
$-$2 | 1 |  $-$1 |  $-$1.0931 |  $-$1.0640 |  $-$0.0291 | 4.9268 | 5.0126 |  $-$0.0858
$-$1 | 2 |  $-$1 | 3.9504 | 4.0170 |  $-$0.0666 | 17.126 | 18.016 |  $-$0.8902
$-$1 | 2 |  $-$2 | 3.9901 | 3.9909 |  $-$0.0008 | 21.032 | 21.575 |  $-$0.5435
$-$2 | 2 |  $-$1 | 2.0065 | 1.9512 | 0.0553 | 24.104 | 24.592 |  $-$0.4883
Figure 3:  Transition probability density function (TPDF)(left) and transition
cumulative distribution function (TCDF)(right), varying $n$ from $10^{3}$ to
$10^{5}$.

##  5 Maximum likelihood estimation

In this section, we determine the maximum likelihood estimates (MLEs) of the
parameters of the BG distribution. Let $X\sim F_{\rm
BG}(\cdot;\mu,\sigma,\delta)$ be a random variable with PDF $f_{\rm
BG}(x;\Theta)$ defined in (2.1), where $\Theta=(\mu,\sigma,\delta)$. Let
$X_{1},\dots,X_{n}$ be a random sample of $X$ and $x=(x_{1},\dots,x_{n})$ the
corresponding observed values of the random sample $X_{1},\dots,X_{n}$. The
log-likelihood function for $\Theta=(\mu,\sigma,\delta)$ is given by

$\displaystyle\\!\ell(\Theta;x)=-n\ln
Z_{\delta}-n\ln\sigma+\sum_{i=1}^{n}\left\\{\ln\big{[}(1-\delta
x_{i})^{2}+1\big{]}-\Big{(}\frac{x_{i}-\mu}{\sigma}\Big{)}-\exp\Big{[}-\Big{(}\frac{x_{i}-\mu}{\sigma}\Big{)}\Big{]}\right\\}.$
(5.1)

The first-order partial derivatives of $Z_{\delta}$ are given by

$\displaystyle\textstyle\frac{\partial Z_{\delta}}{\partial\mu}=\textstyle
2\delta\big{[}\delta(\mu+\sigma\gamma)-1\big{]};$
$\displaystyle\textstyle\frac{\partial
Z_{\delta}}{\partial\sigma}=\frac{\sigma\delta\pi^{2}}{3}+2\delta\gamma\big{[}\delta(\mu+\sigma\gamma)-1\big{]};$
$\displaystyle\textstyle\frac{\partial
Z_{\delta}}{\partial\delta}=\frac{\sigma\delta\pi^{2}}{3}+2\big{[}\delta(\mu+\sigma\gamma)-1\big{]}(\mu+\sigma\gamma).$

Then, the mle of $\mu,\sigma,\delta$ are the solutions of the following system
of equations

$\displaystyle\frac{\partial\ell(\Theta;x)}{\partial\mu}$
$\displaystyle=-\frac{n}{Z_{\delta}}\frac{\partial
Z_{\delta}}{\partial\mu}+\frac{n}{\sigma}-\frac{1}{\sigma}\sum_{i=1}^{n}\exp\Big{[}-\Big{(}\frac{x_{i}-\mu}{\sigma}\Big{)}\Big{]}=0;$
$\displaystyle\frac{\partial\ell(\Theta;x)}{\partial\sigma}$
$\displaystyle=-\frac{n}{Z_{\delta}}\frac{\partial
Z_{\delta}}{\partial\sigma}-\frac{n}{\sigma}-\frac{1}{\sigma^{2}}\sum_{i=1}^{n}(x_{i}-\mu)\left\\{1-\exp\Big{[}-\Big{(}\frac{x_{i}-\mu}{\sigma}\Big{)}\Big{]}\right\\}=0;$
$\displaystyle\frac{\partial\ell(\Theta;x)}{\partial\delta}$
$\displaystyle=-\frac{n}{Z_{\delta}}\frac{\partial
Z_{\delta}}{\partial\delta}-2\sum_{i=1}^{n}{x_{i}(1-\delta
x_{i})\,\over(1-\delta x_{i})^{2}+1}=0.$

The MLEs $\widehat{\mu},\widehat{\sigma}$ and $\widehat{\delta}$ of
$\mu,\sigma$ and $\delta$ are defined as the values of
$\widehat{\mu},\widehat{\sigma}$ and $\widehat{\delta}$ that maximize the log-
likelihood function in (5.1). There will be, in general, no closed form for
the MLE and their obtention will need, in practice, numerical methods.

Since the second order partial derivatives of $Z_{\delta}$ are

$\displaystyle\begin{array}[]{lllll}&\frac{\partial^{2}Z_{\delta}}{\partial\mu^{2}}=2\delta^{2};&\frac{\partial^{2}Z_{\delta}}{\partial\sigma^{2}}=\frac{\delta\pi^{2}}{3}+2\delta^{2}\gamma^{2};\\\\[5.69046pt]
&\frac{\partial^{2}Z_{\delta}}{\partial\delta^{2}}=\frac{\sigma\pi^{2}}{3}+2(\mu+\sigma\gamma)^{2};&\frac{\partial^{2}Z_{\delta}}{\partial\mu\partial\sigma}=\frac{\partial^{2}Z_{\delta}}{\partial\sigma\partial\mu}=2\delta^{2}\gamma;\\\\[5.69046pt]
&\frac{\partial^{2}Z_{\delta}}{\partial\mu\partial\delta}=\frac{\partial^{2}Z_{\delta}}{\partial\delta\partial\mu}=4\delta(\mu+\sigma\gamma);&\frac{\partial^{2}Z_{\delta}}{\partial\sigma\partial\delta}=\frac{\partial^{2}Z_{\delta}}{\partial\delta\partial\sigma}=\frac{\sigma\pi^{2}}{3}+4\delta\gamma(\mu+\sigma\gamma);\end{array}$

the second-order partial derivatives of the log-likelihood function
$\ell(\Theta;x)$ are given by

$\displaystyle\frac{\partial^{2}\ell(\Theta;x)}{\partial\mu^{2}}$
$\displaystyle=-D_{\mu,\mu}(\Theta;n)-\frac{1}{\sigma^{2}}\sum_{i=1}^{n}\exp\Big{[}-\Big{(}\frac{x_{i}-\mu}{\sigma}\Big{)}\Big{]};$
$\displaystyle\frac{\partial^{2}\ell(\Theta;x)}{\partial\sigma^{2}}$
$\displaystyle=-D_{\sigma,\sigma}(\Theta;n)+\frac{n}{\sigma^{2}}+\frac{1}{\sigma^{3}}\sum_{i=1}^{n}(x_{i}-\mu)\left\\{2-\Big{[}2-\Big{(}\frac{x_{i}-\mu}{\sigma}\Big{)}\Big{]}\exp\Big{[}-\Big{(}\frac{x_{i}-\mu}{\sigma}\Big{)}\Big{]}\right\\};$
$\displaystyle\frac{\partial^{2}\ell(\Theta;x)}{\partial\delta^{2}}$
$\displaystyle=-D_{\delta,\delta}(\Theta;n)-2\sum_{i=1}^{n}{x_{i}^{2}\big{[}(1-\delta
x_{i})^{2}-1\big{]}\over\big{[}(1-\delta x_{i})^{2}+1\big{]}^{2}}.$

Already, the second-order mixed derivatives of $\ell(\Theta;x)$ can be written
as

$\displaystyle\frac{\partial^{2}\ell(\Theta;x)}{\partial\mu\partial\sigma}$
$\displaystyle=\frac{\partial^{2}\ell(\Theta;x)}{\partial\sigma\partial\mu}=-D_{\mu,\sigma}(\Theta;n)-\frac{n}{\sigma^{2}}+\frac{1}{\sigma^{2}}\sum_{i=1}^{n}\Big{[}1-\Big{(}\frac{x_{i}-\mu}{\sigma}\Big{)}\Big{]}\exp\Big{[}-\Big{(}\frac{x_{i}-\mu}{\sigma}\Big{)}\Big{]};$
$\displaystyle\frac{\partial^{2}\ell(\Theta;x)}{\partial\mu\partial\delta}$
$\displaystyle=\frac{\partial^{2}\ell(\Theta;x)}{\partial\delta\partial\mu}=-D_{\mu,\delta}(\Theta;n);$
$\displaystyle\frac{\partial^{2}\ell(\Theta;x)}{\partial\sigma\partial\delta}$
$\displaystyle=\frac{\partial^{2}\ell(\Theta;x)}{\partial\delta\partial\sigma}=-D_{\sigma,\delta}(\Theta;n),$

where

$\displaystyle\textstyle D_{u,v}(\Theta;n)=\frac{\partial}{\partial
u}\big{(}\frac{n}{Z_{\delta}}\frac{\partial Z_{\delta}}{\partial
v}\big{)}=\frac{n}{Z_{\delta}}\big{(}\frac{\partial^{2}Z_{\delta}}{\partial
u\partial v}-{1\over Z_{\delta}}\frac{\partial Z_{\delta}}{\partial
u}\frac{\partial Z_{\delta}}{\partial v}\big{)},\quad
u,v\in\\{\mu,\sigma,\delta\\}.$ (5.2)

Thus, the elements of the Fisher information matrix $I_{X}(\Theta)$ are
defined by

$\displaystyle[I_{X}(\Theta)]_{jk}=\mathbb{E}\biggl{[}{\partial\ln f_{\rm
BG}({X};{\Theta})\over\partial\theta_{j}}\ {\partial\ln f_{\rm
BG}({X};{\Theta})\over\partial\theta_{k}}\biggr{]},$

for $\theta_{j},\theta_{k}\in\\{\mu,\sigma,\delta\\}$ and $j,k=1,2,3$. Under
known regularity conditions, the elements $[I_{X}(\Theta)]_{jk}$ are written
as $\mathbb{E}\big{[}-{\partial^{2}\ln f_{\rm
BG}({X};{\Theta})\over\partial\theta_{j}\partial\theta_{k}}\big{]}$. Hence, by
considering the following notations:

$\displaystyle\begin{array}[]{lllll}F_{1}(X)&=&\exp\big{[}{-\big{(}\frac{X-\mu}{\sigma}\big{)}}\big{]},\\\\[8.5359pt]
F_{2}(X)&=&\big{[}1-\big{(}\frac{X-\mu}{\sigma}\big{)}\big{]}\exp\big{[}{-\big{(}\frac{X-\mu}{\sigma}\big{)}}\big{]},\\\\[8.5359pt]
F_{3}(X)&=&(X-\mu)\big{\\{}2-\big{[}2-\big{(}\frac{X-\mu}{\sigma}\big{)}\big{]}\exp\big{[}{-\big{(}\frac{X-\mu}{\sigma}\big{)}}\big{]}\big{\\}},\\\\[8.5359pt]
F_{4}(X)&=&{X^{2}\big{[}(1-\delta X)^{2}-1\big{]}/\big{[}(1-\delta
X)^{2}+1\big{]}^{2}},\end{array}$

we get

$\displaystyle
I_{X}(\Theta)=\scalebox{0.95}{\mbox{$\displaystyle\begin{bmatrix}D_{\mu,\mu}(\Theta;1)+\frac{1}{\sigma^{2}}\mathbb{E}[F_{1}(X)]&D_{\mu,\sigma}(\Theta;1)+\frac{1}{\sigma^{2}}-\frac{1}{\sigma^{2}}\mathbb{E}[F_{2}(X)]&D_{\mu,\delta}(\Theta;1)\\\\[8.5359pt]
D_{\mu,\sigma}(\Theta;1)+\frac{1}{\sigma^{2}}-\frac{1}{\sigma^{2}}\mathbb{E}[F_{2}(X)]&D_{\sigma,\sigma}(\Theta;1)-\frac{1}{\sigma^{2}}-\frac{1}{\sigma^{3}}\mathbb{E}[F_{3}(X)]&D_{\sigma,\delta}(\Theta;1)\\\\[8.5359pt]
D_{\mu,\delta}(\Theta;1)&D_{\sigma,\delta}(\Theta;1)&D_{\delta,\delta}(\Theta;1)+2\mathbb{E}[F_{4}(X)]\
\end{bmatrix}$}},$

where $D_{u,v}(\Theta;n)$ is as in(5.2). We emphasize that, by using Theorem 2
and Corollary 3.6, closed expressions for the expectations
$\mathbb{E}[F_{1}(X)]$, $\mathbb{E}[F_{2}(X)]$ and $\mathbb{E}[F_{3}(X)]$ can
be found. Furthermore, by using Theorem 3, a simple argument shows that
$F_{4}(X)$ is an integrable random variable.

##  6 Real-world data analysis

In this section, we used the BG model to fit one data set coming from daily
pressure observations of the state of Mato Grosso do Sul-Brazil, during the
period of the year 2015 to 2020. The data series employed were taken from
Corumbá automatic station of INMET (Instituto Nacional de Meteorologia) at
https://portal.inmet.gov.br/dadoshistoricos. The required numerical
evaluations are implemented using the R software [14], [6], and [8]. The
maximum likelihood estimates of the parameters were obtained as in Section 5,
using the bgumbel package [1].

For a random sample of the pressure $\\{X_{i}\\}_{t=1}^{T}$, $T=1774$, we fit
their maximum values by Gumbel and BG models. For this, we used the maximum
block technique. That is, we provide $\tau=29$ non-overlapping subsamples of
length $N=60$; $\\{m_{i}\\}_{i=1}^{\tau}$, where
$m_{i}=\max\\{x_{(i-1)N+1}\dots x_{iN+1}\\}$, $\forall i=1,\ldots,\tau=[T/N]$
($[.]$ denotes the integer part ). We applied Ljung – Box’s test [15] for
verifying the null hypothesis of serial independence of the new sample of
maxima. The test statistics did not reject the null hypothesis at the
significance level of $1,7\%$ for the $28$ blocks of size $N=60$ and the last
of size $34$.

The descriptive statistics of sample $\\{m_{i}\\}_{i=1}^{\tau}$ and the sample
centered by the mean are shown in Table 2.

Table 2:  Descriptive statistics Distribution | Mean | Median | Maximum | Minimum | Std Dev
---|---|---|---|---|---
Extreme values | 1009.2 | 1010.0 | 1017.9 | 1001.1 | 4.9301
Centralized extreme values | 0.0000 | 0.7893 | 8.7393 |  $-$8.0607 | 4.9301

The parameter estimates and their corresponding standard error estimates (SE)
from the BG and Gumbel models are shown in Table 3. The results of goodness-
of-fit tests based on the Kolmogorov–Smirnov (KS) test, Akaike information
criterion (AIC) and Bayesian information criterion (BIC) indicate that the BG
model is a better fit, judging on the basis of $p$-values and the lowest value
of the AIC and BIC, see Table 3. The Figures 4 and 5 also indicates that the
adjustment of the data by the BG model is better than the simple Gumbel model.

Table 3: Parameters estimates of Gumbel and BG models. Model | Parameters | Estimate | SE | ks $p$-value | AIC | BIC
---|---|---|---|---|---|---
| $\mu$ |  $-$3.6972 | 0.4613 |  |  |
BG | $\sigma$ | 2.4661 | 0.2222 | 0.8239 | 169.8228 | 173.8194
| $\delta$ |  $-$0.4128 | 0.1220 |  |  |
Gumbel | $\mu$ |  $-$2.4080 | 0.8667 | 0.6576 | 173.2287 | 175.8932
| $\beta$ | 4.3326 | 0.6405 |  |  |
Figure 4: On the left: centralized extreme values histogram versus fitted BG
and Gumbel models. On the right: empirical cdf versus theorical cdf of BG and
Gumbel models. Figure 5:  QQ-plot empirical versus theoretical of BG model.

##  7 Conclusions

In this paper, we proposed an extension to the Gumbel distribution by using a
quadratic transformation technique used to generate bimodal functions produced
due to using the quadratic expression, with an additional bimodality parameter
which modifies the mode of the distribution, composing as a alternative models
for single maxima events. In this generalization, the Gumbel distribution
appears as a particular case. We provide a mathematical treatment of the new
model including the mode, bimodality, moment generating function and moments.
We performed the modelling under a frequentist approach and the estimation of
the parameters was proposed using the maximum likelihood estimation. Finally,
we have performed a statistical modeling with real data by using the new
proposed model in the article. The application demonstrated the practical
relevance of the new model, which also showed the advantage of Gumbel model.

#### Acknowledgements

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

#### Disclosure statement

There are no conflicts of interest to disclose.

## References

  * [1] P. C. Brom, C. E. G. Otiniano, R. Vila, M. B. Pereira (2021). Bimodal Gumbel Distribution, Package bgumbel. https://CRAN.R-project.org/package=bgumbel.
  * [2] M. N. Çankaya, Y. M. Bulut, F. Z. Doğru, O. Arslan (2015). A bimodal extension of the generalized gamma distribution. Revista Colombiana de Estadística, 38:371–384.
  * [3] D. Elal-Olivero (2010). Alpha-skew-normal distribution. Proyecciones Journal of Mathematics, 29:224–240.
  * [4] R. A. Fisher, L. H. C. Tippett (1928). Limiting forms of the frequency distribution of the largest or smallest member of a sample. Mathematical Proceedings of the Cambridge Philosophical Society, 24:180–190.
  * [5] E. J. Gumbel (1958). Statistics of Extremes. Columbia University Press, New York.
  * [6] W. K. Hastings (1970). Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 57:97–109.
  * [7] S. Klugman, H. Panjer, G. Willmot (1998). Loss models: From data to decisions, Wiley, New York.
  * [8] A. D. Martin, K. M. Quinn, J. H. Park (2011). MCMCpack: Markov Chain Monte Carlo in R. Journal of Statistical Software. 42(9): 1-21. URL https://www.jstatsoft.org/v42/i09/.
  * [9] E. Z. Martinez, J. A. Achcar, A. A. Jacome, J. S. Santos (2013). Mixture and non-mixture cure fraction models based on the generalized modified Weibull distribution with an application to gastric cancer data. Computer Methods and Programs in Biomedicine, 112:343–355.
  * [10] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. (1953). Teller, and E. Teller. Equations of state calculations by fast computing machine. Journal of Chemical Physics, 21:1087–1091.
  * [11] I. E. Okorie, A. C. Akpanta, J. Ohakwe (2016). The Exponentiated Gumbel Type-2 Distribution: Properties and Application. International Journal of Mathematics and Mathematical Sciences, vol. 2016, 10 pages.
  * [12] I. E. Okorie, A. C. Akpanta, J. Ohakwe, D. C. Chikezie, E. O. Obi (2017). The Kumaraswamy G Exponentiated Gumbel type-2 distribution. Afrika Statistika, vol. 12, No. 3.
  * [13] E. C. Pinheiro, S. L. P. Ferrari (2015). A comparative review of generalizations of the Gumbel extreme value distribution with an application to wind speed data. Journal of Statistical Computation and Simulation, 86:2241–2261.
  * [14] R-Team (2020). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.
  * [15] A. Trapletti (2016). R: Box-Pierce and Ljung-Box Tests. stat.ethz.ch. Retrieved 2016-06-05. URL https://stat.ethz.ch/R-manual/R-devel/library/stats/html/box.test.html.
  * [16] R. Vila, J. Leao, H. Saulo, M. N. Shahzad, M. Santos-Neto (2020). On a bimodal Birnbaum-Saunders distribution with applications to lifetime data. Brazilian Journal of Probability and Statistics, 34:495–518
  * [17] R. Vila, L. Ferreira, H. Saulo, F. Prataviera, E. M. M. Ortega (2020). A bimodal gamma distribution: Properties, regression model and applications. Statistics: A Journal of Theoretical and Applied Statistics, 54:469–493.
  * [18] R. Vila, M. N. Çankaya (2021). A Bimodal Weibull Distribution: Properties and Inference. To appear in Journal of Applied Statistics.
  * [19] R. Vila, H. Saulo, J. Roldan (2021). On some properties of the bimodal normal distribution and its bivariate version. Preprint, http://arxiv.org/abs/2106.00097.

## Appendix

Figure 6: TPDF (left) and TCDF (right), n = 1000. Figure 7: TPDF (left) and
TCDF (right), n = 10000. Figure 8: TPDF (left) and TCDF (right), n = 100000.
Figure 9: TPDF (left) and TCDF (right), n = 1000. Figure 10: TPDF (left) and
TCDF (right), n = 10000. Figure 11: TPDF (left) and TCDF (right), n = 100000.
Figure 12: TPDF (left) and TCDF (right), n = 1000. Figure 13: TPDF (left) and
TCDF (right), n = 10000. Figure 14: TPDF (left) and TCDF (right), n = 100000.