# A Bimodal Model for Extremes Data

In extreme values theory, for a sufficiently large block size, the maxima distribution is approximated by the generalized extreme value (GEV) distribution. The GEV distribution is a family of continuous probability distributions, which has wide applicability in several areas including hydrology, engineering, science, ecology and finance. However, the GEV distribution is not suitable to model extreme bimodal data. In this paper, we propose an extension of the GEV distribution that incorporate an additional parameter. The additional parameter introduces bimodality and to vary tail weight, i.e., this proposed extension is more flexible than the GEV distribution. Inference for the proposed distribution were performed under the likelihood paradigm. A Monte Carlo experiment is conducted to evaluate the performances of these estimators in finite samples with a discussion of the results. Finally, the proposed distribution is applied to environmental data sets, illustrating their capabilities in challenging cases in extreme value theory.

## Authors

• 2 publications
• 1 publication
• 18 publications
• 12 publications
08/18/2021

### A Model for Bimodal Rates and Proportions

The beta model is the most important distribution for fitting data with ...
04/22/2022

### Choice of mixture Poisson models based on Extreme value theory

Count data are omnipresent in many applied fields, often with overdisper...
07/23/2020

### The r-largest four parameter kappa distribution

The generalized extreme value distribution (GEVD) has been widely used t...
10/19/2021

### Nonstationary seasonal model for daily mean temperature distribution bridging bulk and tails

In traditional extreme value analysis, the bulk of the data is ignored, ...
08/14/2020

### Uniqueness and global optimality of the maximum likelihood estimator for the generalized extreme value distribution

The three-parameter generalized extreme value distribution arises from c...
04/11/2022

### Empirical Bayes inference for the block maxima method

The block maxima method is one of the most popular approaches for extrem...
08/29/2018

### Extreme Value Theory for Open Set Classification - GPD and GEV Classifiers

Classification tasks usually assume that all possible classes are presen...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The generalized extreme value (GEV) distribution is widely used to model extreme events in several areas such as finance, insurance, hydrology, bioinformatics, among others. In journals of statistics and applied areas, a huge amount of works with applications of the GEV distribution can be found. Theory and applications of extreme value distributions can be found in the books by Kotz and Nadarajah (2000), Coles (2001), Beirlant et al. (2004), Haan and Ferreira (2010) , Embrechts et al. (2013), Longin (2016), Scheirer (2017) , among others.

has a Generalized Extreme Value (GEV) distribution with shape parameter , location parameter , and scale paramer denoted by

, if its probability density function (PDF) is given by

valid for in the case , and in the case .

The cumulated distribution function (CDF) of the GEV distribution is given by

 (1.1)

When the pdf (1.1) is also known as the Gumbel distribution.

Here a GEV random variable with and standar scale we denote by . In this case, the PDF and CDF functios are given, respectively, by

 fG(y;ξ,μ)=[1+ξ(y−μ)](−1/ξ)−1exp{−[1+ξ(y−μ)]−1/ξ} (1.2)

and

 FG(y;ξ,μ)=exp{−[1+ξ(y−μ)]−1/ξ}.

In extreme value modeling, in special, environmental data (see Section 7) data with bimodality has appeared more and more. The statistical modeling of this type of data requires distributions that capture bimodality. In this context, Nascimento et al. (2016) provided some new extended models to the GEV distribution as a baseline function and the transmuted GEV distribution introduced by Aryal and Tsokos (2009) , showing advantages compared with the standard GEV distribution. Eljabri and Nadarajah (2017) studied the Kumaraswamy GEV distribution. However, all the models cited above are not suitable for capturing bimodality. In this sense, the mixture of GEV distributions and the mixture of extreme value distributions has been an alternative ( Otiniano et al. (2016) ). The difficulty with these models is in the estimation process, as six sub-models must be considered.

This paper presents a new model for extremes, based on transformation of the standard GEV, with the hope it yields a “better fit” in certain extremes analysis (see Sect. 5). The inclusion of additional parameter in the GEV role is to introduce bimodality, that is, the new distribution is more flexible than the GEV distribution, and without computational complications for the estimation of its parameters in the case .

In Sect. 2, we define the BGEV model. In Sect. 3, we derive the main properties of the bimodal GEV distribution. In Sect. 4, we provide graphical illustrations. In Sect. 5, inference procedure is carried out under the likelihood paradigm. Already, in Sect. 6 discusses some simulation results for the estimation method. In Sect. 7, two illustrative applications in environmental data sets are investigated. Conclusions are addressed in Sect. 8.

## 2 The Bimodal GEV distribution

In this section, we introduce a bimodal distribution of extreme value (BGEV) as Rathei and Swamme (2007) .

We define a random variable with BGEV distribution and parameters , denoted , if its PDF is given by

 fBGEV(x;ξ,μ,σ,δ) =fG(Tσ,δ(x);ξ,μ)T′σ,δ(x) (2.1)

where

 Tσ,δ(x)=σx|x|δ,x∈R,  δ⩾−1,  σ>0 (2.2)

is an invertible transformation with derivatives

 T′σ,δ(x)=σ(δ+1)|x|δ,T′′σ,δ(x)=sign(x)σ(δ+1)δ|x|δ−1, (2.3)

and

 T(k)σ,δ(x)=[sign(x)]k−1σ[k−2∏i=−1(δ−i)]|x|δ−(k−1),k⩾2,

where denotes the sign function.

Since is non-decreasing for , then, its inverse function, denoted by , exists and is given by

 T−1σ,δ(x)=sign(x)(|x|σ)1/(δ+1). (2.4)

The PDF (2) is valid for in the case , and for in the case .

Note that, for , the function is a PDF, because

 ∫∞T−1σ,δ(μ−1/ξ)fBGEV(x;ξ,μ,σ,δ)dx=∫∞μ−1/ξfG(Tσ,δ(x);ξ,μ)dTσ,δ(x)=1.

Similarly, it is verified that, for , is a PDF.

The expression (2) is equivalent to

 fBGEV(x;ξ,μ,σ,δ)=(FG∘Tσ,δ)′(x),

where

 (FG∘Tσ,δ)(x;ξ,μ)=FBGEV(x;ξ,μ,σ,δ) (2.5)

is the CDF of .

When the function (2) is of type (1.1); .

## 3 Some properties of the BGEV distribution

In this section, some mathematical properties as monotonicity, bimodality property, stochastic representation, moments, quantiles and tail behavior of the BGEV distribution are discussed.

### 3.1 Monotonicity

Let be the unique mode for the GEV distribution (1.2). To state the next result we define the following quantities: and .

###### Proposition 3.1.

The following monotonicity properties it hold:

• If and , then the BGEV PDF is increasing for each ;

• If and , then the BGEV PDF is decreasing for each ;

• If and , then the BGEV PDF is decreasing for each ;

• If and , then the BGEV PDF is increasing for each .

###### Proof.

Since , for , is decreasing (resp. increasing) for each (resp. ). For , is increasing (resp. decreasing) for each (resp. ). On the other hand, since is the mode of the GEV distribution (1.2), is increasing (resp. decreasing) for each (resp. ) and (resp. ). Since is increasing and nonnegative for all , and is increasing for each , and , from definition (2) of BGEV density, we have that is the product of two increasing and nonnegative functions. Thus, the BGEV PDF is increasing for each . This completes the proof of first item.

The proof of the other items follows the same reasoning as the one of Item 1). ∎

### 3.2 Bimodality property

###### Proposition 3.2.

The point is a critical point of BGEV density (2) if it is solution of the following equation:

 T′′σ,δ(x)−1+ξ−[1+ξ(Tσ,δ(x)−μ)]−1/ξ1+ξ(Tσ,δ(x)−μ)[T′σ,δ(x)]2=0,

where , and are given in (2.2) and (2.3).

###### Proof.

The proof is trivial and omitted. ∎

###### Lemma 3.3.

Let be a natural number and . The PDF of the BGEV distribution has at most three real critical points.

###### Proof.

Let . Since is a positive number, . By using definitions of , and in (2.2) and (2.3), a simple algebraic manipulation shows that the equation of Proposition 3.2 is equivalent to

 |x|−2Tσ,δ(x)[δy1+1/ξ−(ξ+1)y1/ξ+1]=0,y>0.

For , is a critical point of . So the number of positive roots of the polynomial , when a natural number, determines the number of remaining roots of the BGEV PDF.

By Descartes’ rule of signs (Griffiths 1947 ; Xue 2012 ), the polynomial has two sign changes (the sequence signs is , , ), meaning that this polynomial has two or zero positive roots. This completes the proof. ∎

###### Theorem 3.4 (Uni- or bimodality).

If is an even natural number and , then the PDF of the BGEV distribution is uni- or bimodal.

###### Proof.

If is even, the polynomial defined in Lemma 3.3

has odd degree. Let

be the extension of polynomial to in . Again, by Descartes’ rule of signs the polynomial has one sign change (the sequence signs is , , ), meaning that this one has one negative root. But as a sub-product of the proof of Lemma 3.3 we get that has two or zero positive roots. Since has odd degree and the complex roots are given in pairs, we only have two possibilities: (a) has one negative root and zero positive roots or (b) has one negative root and two positive roots. From case (a) it follows that , , has no positive roots, and by Lemma 3.3, is the only critical point of BGEV density. Already, from case (b) it follows that , , has exactly two positive roots, and therefore, the BGEV PDF has exactly three real critical points: , and .

Finally, since , the uni- or bimodality of the BGEV density follows. ∎

### 3.3 Stochastic representation of the BGEV distribution

###### Proposition 3.5.

Let .

• If , then .

• If , then .

• If and a constant, then .

###### Proof.

If , by (2.5), follows that

 P(Y⩽y)=P(X⩽T−1σ,δ(y))=FBGEV(T−1σ,δ(y);ξ,μ,σ,δ)=(FG∘Tσ,δ)(T−1σ,δ(y);ξ,μ)=FG(y;ξ,μ).

This prove the Item 1). Analogously, the proof of second item follows.

Already, the proof of Item 3) follows by combining the Item 1), the relation , the following well-known fact: If , then for ; with Item 2). ∎

### 3.4 Moments

Here we give a closed analytical formula for the th moment of a random variable with BGEV distribution.

###### Proposition 3.6.

If and , then

where and are the incomplete gamma functions.

###### Proof.

For , by definition of expectation and the expression (1.2), we have

 E(Xk(δ+1)) = ∫∞T−1σ,δ(μ−1/ξ)xk(δ+1)fBGEV(x;ξ,μ,σ,δ)dx (3.1) = ∫∞T−1σ,δ(μ−1/ξ)xk(δ+1)fG(Tσ,δ(x);ξ,μ)dTσ,δ(x) = ∫∞μ−1/ξ[T−1σ,δ(y)]k(δ+1)fG(y;ξ,μ)dy.

When replacing the expression (2.4) of in (3.1) we obtain two cases:
(i) If ,

 E(Xk(δ+1))=1σk∫∞μ−1/ξykfG(y;ξ,μ)dy=1σkE(Yk),    Y∼FG(⋅;ξ,μ),

where

 E(Yk)=k∑i=1(ki)(μ−1ξ)k−i1ξiΓ(1−ξi),ξ<1/k.

(ii) If ,

 E(Xk(δ+1)) = (−1)k(δ+1)σk∫0μ−1/ξykfG(y;ξ,μ)dy+1σk∫∞0ykfG(y;ξ,μ)dy (3.2) = (−1)k(δ+1)σkI1+1σkI2.

Using the PDF (1.2) and the substitution , we obtain

 I1 = 1ξk+1∫∞(1−ξμ)−1[w−1+(ξμ−1)]k w1ξ−1e−w1/ξdw. (3.3)

To solve the integral of (3.3) we use the Newton’s binomial formula, thus

 I1 = (3.4) = 1ξkk∑i=1(ki)(ξμ−1)k−i∫∞(1−ξμ)−1/ξz−ξie−zdz,

where in the last line we used the new substitution .

Similarly we obtain

 I2 = 1ξkk∑i=1(ki)(ξμ−1)k−i∫(1−ξμ)−1/ξ0z−ξie−zdz. (3.5)

The proof is completed by expressing as integrals of (3.4) and (3.5) in terms of the incomplete gamma functions and then updating the equation (3.2).

For the proof is similar. ∎

###### Corollary 3.7.

Let and , then

 E(X)=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩1(ξσ)⟦1/(δ+1)⟧⟦1/(δ+1)⟧∑i=0(⟦1/(δ+1)⟧i)(ξμ−1)⟦1/(δ+1)⟧−iΓ(1−ξi),μ−1ξ>0,(−1)k(δ+1)(ξσ)⟦1/(δ+1)⟧⟦1/(δ+1)⟧∑i=0(⟦1/(δ+1)⟧i)(ξμ−1)⟦1/(δ+1)⟧−iγ(1−ξi;(1−ξμ)−1/ξ)+1(ξσ)⟦1/(δ+1)⟧⟦1/(δ+1)⟧∑i=0(⟦1/(δ+1)⟧i)(ξμ−1)⟦1/(δ+1)⟧−iΓ(1−ξi;(1−ξμ)−1/ξ),μ−1ξ<0.
###### Proof.

The result follows directly of Proposition 3.6 by considering and . ∎

### 3.5 Quantiles

In the next section, through Monte Carlo simulation studies, we test the behavior of the parameter estimates of the BGEV model. The random samples of the population are simulated using the inverse transformation method based on quantiles. The quantile function of a random variable is defined by

 xQ=inf{x∈R:Q≤FBGEV(x;ξ,μ,σ,δ)}

or equivalent by

 xQ = T−1σ,δ(F−1G(Q)).

From (2.5), we have . Now, by using the expression (2.4) of , we get

 xQ=⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩[μσ+(−lnQ)−ξ−1σξ]1/(δ+1),μ+(−lnQ)−ξ−1ξ>0,(−1)2+δ1+δ[μσ+(−lnQ)−ξ−1σξ]1/(δ+1),μ+(−lnQ)−ξ−1ξ<0. (3.6)

### 3.6 Tail behavior of BGEV

According to Embrechts et al. (2013) , asymptotic estimates of probability models are omnipresent in statistics, insurance mathematics, mathematical finance, and hydrology. In this sense, the theory of regular variation plays a crucial role. Concepts and properties of regularly varying can be found in the book of Bingham et al. (1987) In order to study the asymptotic behavior of the tail of the BGEV distribution, here we show some of the results of regular variation.

###### Definition 3.1.

A positive Lebesgue measurable function on is regularly varying at of index () if

 limt→∞h(tx)h(t)=xα.

Suppose is a distribution function such that for all . Let and for some , then

 ¯¯¯¯F(x)=o(x−αL(x)) (3.7)

for some , a slowly varying function. The expression (3.7) indicates that the tail of the distribution slowly decays, as , that is has heavy (right) tail with index . Thus, is Pareto-like distribution.

###### Proposition 3.8.

Let with and . Then

 ¯¯¯¯FBGEV(x;ξ,μ,σ,δ)=o(x−δ+1ξ),  x>0. (3.8)
###### Proof.

The proof follows of the straightforward calculation of

 limt→∞¯¯¯¯FBGEV(tx;ξ,μ,σ,δ)¯¯¯¯FBGEV(t;ξ,μ,σ,δ)=x−δ+1ξ,

where is given in (2.5). ∎

From (3.8) we have that is Pareto-like distribution with tail index . The new parameter affects the weight of the tail.

## 4 Graphical study

The role that the parameters , , , , and play in the BGEV distribution were investigated graphically by generating different densities with variations in each parameter. Fig. 1 shows that is a scale parameter. According to Fig. 2, the shape of the BGEV PDF changed by varying the parameter. That is, is a shape parameter. Fig. 2 suggests that depending on the combinatation of the parameters values and the BGEV density can be unimodal or bimodal. For example, is unimodal and is bimodal. In the bimodal case, the first mode value is greater than the second mode value for , while for the second mode value is greater than the first mode. Figures 3 and 4 correspond to variations of the new parameter . When the BGEV coincides with the basic GEV distribution. Two scenarios were considered. Scenario 1 for (see Fig. 3), and scenario 2 for (see Fig. 4 ). In scenario 1 , Fig. 3, when the PDF is unimodal and for the PDF is bimodal. In bimodal case, the mode values increases as delta grows. In scenario 2, the Fig. 4 shows that when delta takes positive or negative values close to zero and the PDF is bimodal. In this case, the mode value also increases as delta increases.

## 5 Maximum likelihood estimation

In this section we present the expressions to be calculated to obtain the maximum likelihood estimators for the parameters of a BGEV distribution.

Let be a random variable with PDF defined in (2), where . Let be a random sample of X and the corresponding observed values of the random sample . The log-likelihood function for is given by

 (5.1)

where Note that, for all ,

 ∂Ψi(Θ)∂μ=−ξ,∂Ψi(Θ)∂σ=ξxi|xi|δ,∂Ψi(Θ)∂δ=ξσxi|xi|δln|xi|,∂Ψi(Θ)∂ξ=σxi|xi|δ−μ.

Then, the maximum likelihood estimates of are the solutions of the following system of equations

 ∂ℓn(Θ;x)∂μ=μ∑ni=1Ωi(Θ)=0,∂ℓn(Θ;x)∂σ=nσ−∑ni=1xi|xi|δΩi(Θ)=0,∂ℓn(Θ;x)∂δ=nδ+1+∑ni=1[ln|xi|+σxi|xi|δln|xi|Ωi(Θ)]=0,∂ℓn(Θ;x)∂ξ=ξ−2∑ni=1(σxi|xi|δ−μ){lnΨi(Θ)−ξΩi(Θ)}=0, (5.2)

where Since

 ∂Ωi(Θ)∂θ=−ξ−1Ψ−1i(Θ)[ξΩ−1i(Θ)−Ψ−1/ξi(Θ)]∂Ψi(Θ)∂θ,θ∈{ξ,μ,σ,δ},

the second-order partial derivatives of the log-likelihood function are given by

 ∂2ℓn(Θ;x)∂μ2=∑ni=1{Ωi(Θ)+μΨ−1i(Θ)[ξΩi(Θ)−Ψ−1/ξi(Θ)]},∂2ℓn(Θ;x)∂σ2=−nσ2+∑ni=1x2i|xi|2δΨ−1i(Θ)[ξΩi(Θ)−Ψ−1/ξi(Θ)],∂2ℓn(Θ;x)∂δ2=−n(δ+1)2+σ∑ni=1xi|xi|2δ(ln|xi|)2{Ωi(Θ)−σxiΨ−1i(Θ)[ξΩi(Θ)−Ψ−1/ξi(Θ)]},∂2ℓn(Θ;x)∂ξ2=−2ξ−3∑ni=1(σxi|xi|δ−μ)[lnΨi(Θ)−ξΩi(Θ)] −ξ−2∑ni=1(σxi|xi|δ−μ)2{Ωi(Θ)−Ψ−1i(Θ)[ξΩi(Θ)−Ψ−1/ξi(Θ)+1]}, (5.3)

and, by Schwarz’s theorem or Clairaut’s theorem on equality of mixed partials (James 1966), the second-order mixed derivatives of can be written as

 ∂2ℓn(Θ;x)∂μ∂σ=∂2ℓn(Θ;x)∂σ∂μ=−μ∑ni=1xi|xi|δΨ−1i(Θ)[ξΩi(Θ)−Ψ−1/ξi(Θ)],∂2ℓn(Θ;x)∂μ∂δ=∂2ℓn(Θ;x)∂δ∂μ=−μσ∑ni=1xi|xi|δln|xi|Ψ−1i(Θ)[ξΩi(Θ)−Ψ−1/ξi(Θ)],∂2ℓn(Θ;x)∂μ∂ξ=∂2ℓn(Θ;x)∂ξ∂μ=−μξ−1∑ni=1(σxi|xi|δ−μ)Ψ−1i(Θ)[ξΩi(Θ)−Ψ−1/ξi(Θ)],∂2ℓn(Θ;x)∂σ∂δ=∂2ℓn(Θ;x)∂δ∂σ=−∑ni=1xi|xi|2δln|xi|{Ωi(Θ)−σxiΨ−1i(Θ)[ξΩi(Θ)−Ψ−1/ξi(Θ)]},∂2ℓn(Θ;x)∂σ∂ξ=∂2ℓn(Θ;x)∂ξ∂σ=−ξ−1∑ni=1xi|xi|δ(σxi|xi|δ−μ)Ψ−1i(Θ)[ξΩi(Θ)−Ψ−1/ξi(Θ)],∂2ℓn(Θ;x)∂δ∂ξ=∂2ℓn(Θ;x)∂ξ∂δ=−σξ−1∑ni=1xi|xi|δln|xi|(σxi|xi|δ−μ)Ψ−1i(Θ)[ξΩi(Θ)−Ψ−1/ξi(Θ)]. (5.4)

If , under certain regularity conditions, the Fisher information matrix (FIM) may also be written as

 [I(Θ)]i,j=−E[∂2ℓ1(Θ;X)∂θi∂θj∣∣∣Θ],θi,θj∈{ξ,μ,σ,δ},

where are given in (5.3) and (5.4) by taking . It is clear that the analytical calculations of the FIM above are difficult. But by using the Monte Carlo estimates of the Hessian of the negative log-likelihood function we can form an average of these ones as an estimate of the FIM (Spall 2005).

## 6 Numerical Illustrations

The performance of maximum likelihood estimators calculated was tested by Monte Carlo simulation with 45 combinations of and . We consider the scale parameter always equal to 1; .

We use the following procedure:
(i) Generate random samples of by using the quantile function (3.6). Were generated M=100 random samples (replications) of size N=50, 100, 250, 1000, for each choice of .
(ii) Calculate the maximum likelihood estimates (MLEs) of using the equations (5.1)- (5.2). Here, use the Nelder and Mead’ optimization method implemented in the R statistical program (Nelder and Mead 1965).
(iii) For the initial value, consider the real parameter value added by a value of a auniform variable in .

The estimation results are found in Tables 1 - 5. We have five different configurations. In configuration 1, the shape parameters , and . Its mean estimates, the bias and the mean square error (MSE) of the estimates are present in Table 1. In configuration 2, , and . For this configuration, the mean estimates, the bias and the MSE of the estimates are present in Table 2. In cofigurations 3, 4, and 5, , and , and , respectively. The mean estimates, the bias and the MSE of the estimates for these three configurations are shown in the Tables 3, 4, and 5. In the five configurations, the bias and MSE of the mean estimates are small. That is, the algorithm has obtained good estimates. The bias and the MSE decreased as the sample size was increased from 50 to 1000.