Modelling wind speed with a univariate probability distribution depending on two baseline functions
Fábio V. J. Silveira, Frank Gomes-Silva, Cícero C. R. Brito, Jader S. Jale, Felipe R. S. Gusmão, Sílvio F. A. Xavier-Júnior, João S. Rocha
MModelling wind speed with a univariate probability distributiondepending on two baseline functions
Fábio V. J. Silveira a, * , Frank Gomes-Silva b , Cícero C. R. Brito c , Jader S. Jale b , Felipe R.S. Gusmão b , Sílvio F. A. Xavier-Júnior d , João S. Rocha b a Federal Institute of Education, Science and Technology of Paraíba, João Pessoa, PB, Brazil b Department of Statistics and Informatics, Federal Rural University of Pernambuco, Recife, PE, Brazil c Federal Institute of Education, Science and Technology of Pernambuco, Recife, PE, Brazil d Department of Statistics, Paraíba State University, Campina Grande, PB, Brazil
Abstract
Characterizing the wind speed distribution properly is essential for the satisfactory produc-tion of potential energy in wind farms, being the mixture models usually employed in thedescription of such data. However, some mixture models commonly have the undesirableproperty of non-identifiability. In this work, we present an alternative distribution which isable to fit the wind speed data adequately. The new model, called Normal-Weibull-Weibull,is identifiable and its cumulative distribution function is written as a composition of twobaseline functions. We discuss structural properties of the class that generates the proposedmodel, such as the linear representation of the probability density function, moments andmoment generating function. We perform a Monte Carlo simulation study to investigatethe behavior of the maximum likelihood estimates of the parameters. Finally, we presentapplications of the new distribution for modelling wind speed data measured in five differentcities of the Northeastern Region of Brazil.
Keywords:
Goodness-of-fit; Identifiability; L-BFGS-B algorithm; Maximum likelihood;Monte Carlo simulation
1. Introduction and proposed class
The concern about the emission of greenhouse gases and environmental contaminationfrom conventional energy generation procedures like coal and oil power plants encouragesresearch on alternative resources. A smaller impact on the environment is an advantage of * Corresponding author
Email address: [email protected] (Fábio V. J. Silveira)
Preprint arXiv a r X i v : . [ s t a t . M E ] J a n leaner and sustainable energy production techniques, such as solar, geothermal and wind,over the combustion of fossil fuels.The installed capacity of wind power in Brazil increased from 29 MW in 2005 to roughly16,000 MW (9% of the total capacity of electricity generation) in 2020 [1]. The suitablechoice of the wind turbine must match with the wind behavior at the site of installation.Perkin et al. [2] mention that inadequate turbine selection results in a financially sub-optimalinvestment. Thus, setting the probability distribution appropriately to model the wind speedis essential. Eltamaly [3] used the two-parameter Weibull distribution in a new computerprogram to perform the calculations required to precisely design the wind energy systemand to seek the compatibility between sites and turbines.Ilhan and Kantar [4] declare that despite the wide acceptance of Weibull distribution[5, 6, 7, 8], it may sometimes be poor to model all wind speed data available in nature.Hereupon, they propose using two possible models, namely, the skewed generalized errordistribution [9] and the skewed 𝑡 distribution [10]. Some other distributions used for windspeed and power modelling are Rayleigh [5], gamma [11], normal [6], generalized extremevalue [8] and Birnbaum-Saunders [12]. Additionally, applications of nonparametric methodsto wind speed modelling are also found in the literature [13, 14, 15].Oftentimes, one requires more flexibility from the probability density function (pdf), asin case of bimodality [16] or calm winds regime [17]. In general, finite mixture models aremore flexible than the typical single ones. Akdag et al. [18] compared the usual biparametricWeibull and the two-component mixture Weibull distribution in a study focused on windregimes presenting nearly zero percentage of null speeds; they concluded that the mixture ismore suitable to describe such wind conditions. Carta and Ramírez [19] used three differentmethods to estimate the parameters of the two-component mixture Weibull, namely, themethod of moments, maximum likelihood and least squares; they verified that there is nosignificant difference among them.The mixture density can be written as: 𝑓 ( 𝑥 ; 𝜓 ) = 𝑑 ∑︁ 𝑖 =1 𝑤 𝑖 𝑓 𝑖 ( 𝑥 ; 𝜃 𝑖 ) (1)where the vector 𝜓 = ( 𝑤 , . . . , 𝑤 𝑑 − , 𝜂 ⊤ ) ⊤ contains the unknown parameters of the mixturemodel and the vector 𝜂 contains all the distinct parameters in 𝜃 , . . . , 𝜃 𝑑 . Since 𝑤 , . . . , 𝑤 𝑑 are positive and sum up to one, the presence of 𝑤 𝑑 in 𝜓 is unnecessary.The general definition of identifiability states that a family of densities { 𝑓 ( 𝑥 ; 𝜓 ) : 𝜓 ∈ Ψ }
2s identifiable if: 𝑓 ( 𝑥 ; 𝜓 ) = 𝑓 ( 𝑥 ; 𝜓 ⋆ ) ⇔ 𝜓 = 𝜓 ⋆ . (2)It is not seldom that (2) fails when two or more component densities in (1) belong to thesame parametric family. Such is the case of the mixture of normal densities. Consider 𝑑 = 2, 𝑓 and 𝑓 are normal densities, 𝑤 = 0 . 𝜂 = ( 𝜇 , 𝜎 , 𝜇 , 𝜎 ) ⊤ and 𝜂 = ( 𝜇 , 𝜎 , 𝜇 , 𝜎 ) ⊤ , where 𝜇 ̸ = 𝜇 , 𝜎 ̸ = 𝜎 . We have 𝜓 = ( 𝑤 , 𝜂 ⊤ ) ⊤ ̸ = ( 𝑤 , 𝜂 ⊤ ) ⊤ = 𝜓 ⋆ ⇒ 𝑓 ( 𝑥 ; 𝜓 ) = 𝑓 ( 𝑥 ; 𝜓 ⋆ ). That is, (1) may be invariant under certain permutations ofthe elements in the parametric vector. McLachlan and Peel [20] mention an alternativedefinition of identifiability for mixture models, such that the mixture of 𝑑 normal densitieswould be identifiable under specific conditions. Nonetheless, they remark that it does notovercome the complications due to the interchanging of component labels.Models like the mixture of normal or Weibull densities are quite flexible tools, althoughthe parametric estimation is only credible if the distribution under study is identifiable. Wepresent in this paper a class, whose submodels may be feasible alternatives to mixtures oftwo components from the same parametric family. The class is derived using the method ofgenerating classes of probability distributions of Brito et al. [21]. Its cumulative distributionfunction (cdf) is formulated as a composition of two baselines and under certain conditions,it satisfies (2), even if both baselines belong to the same parametric family.The cdf of the general class is given by: 𝐹 ( 𝑥 ) = 𝜁 ( 𝑥 ) 𝑛 ∑︁ 𝑗 =1 ∫︁ 𝑈 𝑗 ( 𝑥 ) 𝐿 𝑗 ( 𝑥 ) d 𝐻 ( 𝑡 ) − 𝜈 ( 𝑥 ) 𝑛 ∑︁ 𝑗 =1 ∫︁ 𝑉 𝑗 ( 𝑥 ) 𝑀 𝑗 ( 𝑥 ) d 𝐻 ( 𝑡 ) (3)where 𝐻 is a cdf, 𝑛 ∈ N , 𝜁, 𝜈 : R ↦→ R and 𝐿 𝑗 , 𝑈 𝑗 , 𝑀 𝑗 , 𝑉 𝑗 : R ↦→ R ∪ {±∞} are specialfunctions that will be discussed in the next section. ( 𝐺 , 𝐺 ) class and some structural properties The method established by Brito et al. [21] states that if
𝐻, 𝜁, 𝜈 : R ↦→ R and 𝐿 𝑗 , 𝑈 𝑗 , 𝑀 𝑗 , 𝑉 𝑗 : R ↦→ R ∪ {±∞} for 𝑗 = 1 , , , . . . , 𝑛 are monotonic and right continuous functions suchthat:(c1) 𝐻 is a cdf and 𝜁 and 𝜈 are non-negative;(c2) 𝜁 ( 𝑥 ), 𝑈 𝑗 ( 𝑥 ) and 𝑀 𝑗 ( 𝑥 ) are non-decreasing and 𝜈 ( 𝑥 ), 𝑉 𝑗 ( 𝑥 ), 𝐿 𝑗 ( 𝑥 ) are non-increasing ∀ 𝑗 = 1 , , , . . . , 𝑛 ; 3c3) If lim 𝑥 →−∞ 𝜁 ( 𝑥 ) ̸ = lim 𝑥 →−∞ 𝜈 ( 𝑥 ), then lim 𝑥 →−∞ 𝜁 ( 𝑥 ) = 0; or lim 𝑥 →−∞ 𝑈 𝑗 ( 𝑥 ) = lim 𝑥 →−∞ 𝐿 𝑗 ( 𝑥 ) ∀ 𝑗 = 1 , , , . . . , 𝑛 , and lim 𝑥 →−∞ 𝜈 ( 𝑥 ) = 0; or lim 𝑥 →−∞ 𝑀 𝑗 ( 𝑥 ) = lim 𝑥 →−∞ 𝑉 𝑗 ( 𝑥 ) ∀ 𝑗 = 1 , , , . . . , 𝑛 ;(c4) If lim 𝑥 →−∞ 𝜁 ( 𝑥 ) = lim 𝑥 →−∞ 𝜈 ( 𝑥 ) ̸ = 0, then lim 𝑥 →−∞ 𝑈 𝑗 ( 𝑥 ) = lim 𝑥 →−∞ 𝑉 𝑗 ( 𝑥 ) and lim 𝑥 →−∞ 𝑀 𝑗 ( 𝑥 ) =lim 𝑥 →−∞ 𝐿 𝑗 ( 𝑥 ) ∀ 𝑗 = 1 , , , . . . , 𝑛 ;(c5) lim 𝑥 →−∞ 𝐿 𝑗 ( 𝑥 ) ≤ lim 𝑥 →−∞ 𝑈 𝑗 ( 𝑥 ) and if lim 𝑥 →−∞ 𝜈 ( 𝑥 ) ̸ = 0, then lim 𝑥 → + ∞ 𝑀 𝑗 ( 𝑥 ) ≤ lim 𝑥 → + ∞ 𝑉 𝑗 ( 𝑥 ) ∀ 𝑗 =1 , , , . . . , 𝑛 ;(c6) lim 𝑥 → + ∞ 𝑈 𝑛 ( 𝑥 ) ≥ sup { 𝑥 ∈ R : 𝐻 ( 𝑥 ) < } and lim 𝑥 → + ∞ 𝐿 ( 𝑥 ) ≤ inf { 𝑥 ∈ R : 𝐻 ( 𝑥 ) > } ;(c7) lim 𝑥 → + ∞ 𝜁 ( 𝑥 ) = 1;(c8) lim 𝑥 → + ∞ 𝜈 ( 𝑥 ) = 0 or lim 𝑥 → + ∞ 𝑀 𝑗 ( 𝑥 ) = lim 𝑥 → + ∞ 𝑉 𝑗 ( 𝑥 ) ∀ 𝑗 = 1 , , , . . . , 𝑛 and 𝑛 ≥ 𝑥 → + ∞ 𝑈 𝑗 ( 𝑥 ) = lim 𝑥 → + ∞ 𝐿 𝑗 +1 ( 𝑥 ) ∀ 𝑗 = 1 , , , . . . , 𝑛 − 𝑛 ≥ 𝐻 is a cdf without points of discontinuity or all functions 𝐿 𝑗 ( 𝑥 ) and 𝑉 𝑗 ( 𝑥 ) are constantat the right of the vicinity of points whose image are points of discontinuity of 𝐻 , beingalso continuous in that points. Moreover, 𝐻 does not have any point of discontinuity inthe set {︂ lim 𝑥 →±∞ 𝐿 𝑗 ( 𝑥 ) , lim 𝑥 →±∞ 𝑈 𝑗 ( 𝑥 ) , lim 𝑥 →±∞ 𝑀 𝑗 ( 𝑥 ) , lim 𝑥 →±∞ 𝑉 𝑗 ( 𝑥 ) }︂ for some 𝑗 = 1 , , , . . . , 𝑛 ;then Equation (3) is a cdf.Let 𝑛 = 1, 𝐻 ( 𝑡 ) = Φ( 𝑡 ), namely, the standard normal cdf, 𝜁 ( 𝑥 ) = 1, 𝜈 ( 𝑥 ) = 0, 𝑈 ( 𝑥 ) = 𝐺 ( 𝑥 ) / [1 − 𝐺 ( 𝑥 )] and 𝐿 ( 𝑥 ) = log[1 − 𝐺 ( 𝑥 )], where 𝐺 ( 𝑥 ) and 𝐺 ( 𝑥 ) are cdfs. The functionin Equation (3) turns into: 𝐹 𝐺 ,𝐺 ( 𝑥 ) = ∫︁ 𝐺 𝑥 )1 − 𝐺 𝑥 ) log[1 − 𝐺 ( 𝑥 )] dΦ( 𝑡 ) . (4)We took 𝑈 and − 𝐿 from the table of differentiable and monotonically non-decreasingfunctions presented in the well-known paper of Alzaatreh et al. [22], whose method wasused to create generalized distributions of the T-X family. We have intentionally pickedthe two simplest functions from the cited table; alternative (and more complicated) choicesfor 𝑈 and 𝐿 would naturally give rise to different classes. Defining 𝑀 ( 𝑥 ) and 𝑉 ( 𝑥 ) isnot relevant, since 𝜈 ( 𝑥 ) = 0. Also, for obvious reasons, there is no need to verify (c4),(c5) and (c9). The conditions (c1), (c7), (c8) and (c10) are straightforward. As 𝑈 ( 𝑥 )and 𝜁 ( 𝑥 ) are non-decreasing and 𝐿 ( 𝑥 ) is non-increasing, (c2) is true. It is easy to verify4hat lim 𝑥 →−∞ 𝑈 ( 𝑥 ) = 0 = lim 𝑥 →−∞ 𝐿 ( 𝑥 ); and since lim 𝑥 →−∞ 𝜈 ( 𝑥 ) = 0, (c3) is satisfied. Thecondition (c6) is also true because lim 𝑥 → + ∞ 𝑈 ( 𝑥 ) = + ∞ = sup { 𝑥 ∈ R : Φ( 𝑥 ) < } andlim 𝑥 → + ∞ 𝐿 ( 𝑥 ) = −∞ = inf { 𝑥 ∈ R : Φ( 𝑥 ) > } .Thereby, in agreement with the method exposed above, Equation (4) is a cdf. As alreadymentioned, it can be viewed as a composite function of two baselines. Henceforth, let it bedenoted by Normal-( 𝐺 , 𝐺 ) class of probability distributions.Since 𝜑 ( 𝑡 ) = √ 𝜋 𝑒 − 𝑡 / , and Φ( 𝑥 ) = ∫︀ 𝑥 −∞ 𝜑 ( 𝑡 )d 𝑡 , one can write Equation (4) as follows: 𝐹 𝐺 ,𝐺 ( 𝑥 ) = Φ (︃ 𝐺 ( 𝑥 )1 − 𝐺 ( 𝑥 ) )︃ − Φ (log[1 − 𝐺 ( 𝑥 )]) . (5)In case of continuous 𝐺 ( 𝑥 ) and 𝐺 ( 𝑥 ), one can take the derivative of Equation (5) withrespect to 𝑥 to obtain the following pdf: 𝑓 𝐺 ,𝐺 ( 𝑥 ) = 𝜑 (︃ 𝐺 ( 𝑥 )1 − 𝐺 ( 𝑥 ) )︃ 𝑔 ( 𝑥 )[1 − 𝐺 ( 𝑥 )] + 𝜑 (log[1 − 𝐺 ( 𝑥 )]) 𝑔 ( 𝑥 )1 − 𝐺 ( 𝑥 ) , (6)where 𝑔 𝑖 ( 𝑥 ) is the pdf of the random variable whose cdf is 𝐺 𝑖 ( 𝑥 ), for 𝑖 ∈ { , } .At this point, we need to define properly the support of the distributions that emergefrom the new class. Submodels of classes that may be written as a composite function ofone single baseline usually have the same support of the baseline. However, characterizingthe support of a submodel from (4) is not so straightforward, especially if the two baselineshave different supports. As previously mentioned, given that 𝑈 ( 𝐺 ( 𝑥 ) , 𝐺 ( 𝑥 )) = 𝐺 ( 𝑥 ) / [1 − 𝐺 ( 𝑥 )], 𝐿 ( 𝐺 ( 𝑥 ) , 𝐺 ( 𝑥 )) = log[1 − 𝐺 ( 𝑥 )] and 𝑆 𝐻 = R , namely, the support of 𝐻 ( 𝑡 ) is R ,we have that:(a) 𝑆 𝐻 is a convex set;(b) 𝑈 (1 ,
1) = 𝑈 ( 𝐺 (+ ∞ ) , 𝐺 (+ ∞ )) = + ∞ = sup { 𝑥 ∈ R : Φ( 𝑥 ) < } , 𝐿 (1 ,
1) = 𝐿 ( 𝐺 (+ ∞ ) , 𝐺 (+ ∞ )) = −∞ = inf { 𝑥 ∈ R : Φ( 𝑥 ) > } , 𝑈 ( 𝐺 ( 𝑥 ) , 𝐺 ( 𝑥 )) and 𝐿 ( 𝐺 ( 𝑥 ) , 𝐺 ( 𝑥 )) are monotonic functions.According to the Theorem (T4) in [21], (a) and (b) entail that the support of a distributionfrom (4) is the union of the supports of 𝐺 and 𝐺 .In the following lines, we demonstrate that, under specific conditions, the distributionsgenerated by (5) enjoy the attractive property of identifiability. It is important because itassures the uniqueness of the estimates of the parameters.5 heorem 1. Let 𝐺 ( 𝑥 | 𝜃 ) and 𝐺 ( 𝑥 | 𝜃 ) be the baseline cdfs of the normal- ( 𝐺 , 𝐺 ) cdf 𝐹 𝐺 ,𝐺 ( 𝑥 | 𝜃 ) (5), 𝜃 = ( 𝜃 , . . . , 𝜃 𝑟 ) ∈ Θ , 𝜃 = ( 𝜃 𝑟 +1 , . . . , 𝜃 𝑟 + 𝑚 ) ∈ Θ and 𝜃 = ( 𝜃 , . . . , 𝜃 𝑟 , 𝜃 𝑟 +1 , . . . , 𝜃 𝑟 + 𝑚 ) ∈ Θ , where Θ , Θ and Θ are the parametric spaces as-sociated with 𝐺 , 𝐺 and 𝐹 𝐺 ,𝐺 respectively. If 𝐺 and 𝐺 are identifiable, then 𝐹 𝐺 ,𝐺 isidentifiable.Proof. Assume that Φ (︁ 𝐺 ( 𝑥 | 𝜃 )1 − 𝐺 ( 𝑥 | 𝜃 ) )︁ = Φ (︁ 𝐺 ( 𝑥 | 𝜃 ⋆ )1 − 𝐺 ( 𝑥 | 𝜃 ⋆ ) )︁ , where { 𝜃 , 𝜃 ⋆ } ⊂ Θ and 𝜃 ̸ = 𝜃 ⋆ .Since Φ is injective, 𝐺 ( 𝑥 | 𝜃 )1 − 𝐺 ( 𝑥 | 𝜃 ) = 𝐺 ( 𝑥 | 𝜃 ⋆ )1 − 𝐺 ( 𝑥 | 𝜃 ⋆ ) ⇒ 𝐺 ( 𝑥 | 𝜃 ) = 𝐺 ( 𝑥 | 𝜃 ⋆ ); it is a contradiction,because it denies the identifiability of 𝐺 . Therefore, if 𝜃 ̸ = 𝜃 ⋆ then Φ (︁ 𝐺 ( 𝑥 | 𝜃 )1 − 𝐺 ( 𝑥 | 𝜃 ) )︁ ̸ =Φ (︁ 𝐺 ( 𝑥 | 𝜃 ⋆ )1 − 𝐺 ( 𝑥 | 𝜃 ⋆ ) )︁ . Analogously, it is easy to verify that for { 𝜃 , 𝜃 ⋆ } ⊂ Θ , if 𝜃 ̸ = 𝜃 ⋆ thenΦ (log[1 − 𝐺 ( 𝑥 | 𝜃 )]) ̸ = Φ (log[1 − 𝐺 ( 𝑥 | 𝜃 ⋆ )]).Now consider { 𝜃 , 𝜃 ⋆ } ⊂ Θ such that 𝜃 ̸ = 𝜃 ⋆ and assume that 𝐹 𝐺 ,𝐺 ( 𝑥 | 𝜃 ) = 𝐹 𝐺 ,𝐺 ( 𝑥 | 𝜃 ⋆ ).If 𝜃 = 𝜃 ⋆ and 𝜃 ̸ = 𝜃 ⋆ , then we can infer from (5) that 𝐺 ( 𝑥 | 𝜃 ) = 𝐺 ( 𝑥 | 𝜃 ⋆ ), namely, anabsurd. Likewise, if 𝜃 ̸ = 𝜃 ⋆ and 𝜃 = 𝜃 ⋆ , we get to similar contradiction. If 𝜃 ̸ = 𝜃 ⋆ and 𝜃 ̸ = 𝜃 ⋆ , then the assumption fails since 𝐹 𝐺 ,𝐺 ( 𝑥 | 𝜃 ) ̸ = 𝐹 𝐺 ,𝐺 ( 𝑥 | 𝜃 ⋆ ) for almost all values of 𝑥 within the support. Therefore, 𝐹 𝐺 ,𝐺 is identifiable. The normal cdf can be written in terms of the error function erf as follows:Φ( 𝑧 ) = 12 [︃ (︃ 𝑧 √ )︃]︃ , (7)where erf( 𝑧 ) = √ 𝜋 ∫︀ 𝑧 𝑒 − 𝑡 d 𝑡 . Since erf( 𝑧/ √
2) may be linearly represented by:erf (︃ 𝑧 √ )︃ = 2 √ 𝜋 ∞ ∑︁ 𝑛 =0 ( − 𝑛 · ( 𝑧/ √ 𝑛 +1 𝑛 !(2 𝑛 + 1)= √︃ 𝜋 · ∞ ∑︁ 𝑛 =0 (︂ − )︂ 𝑛 𝑧 𝑛 +1 𝑛 !(2 𝑛 + 1) , (8)replacing Equation (8) in Equation (7), we get to:Φ( 𝑧 ) = 12 + 1 √ 𝜋 ∞ ∑︁ 𝑛 =0 (︂ − )︂ 𝑛 𝑧 𝑛 +1 𝑛 !(2 𝑛 + 1) . (9)6ow using the result of Equation (9) in Equation (5), we have: 𝐹 𝐺 ,𝐺 ( 𝑥 ) = ∞ ∑︁ 𝑛 =0 ( − / 𝑛 𝑛 !(2 𝑛 + 1) √ 𝜋 ⎡⎢⎢⎢⎢⎣(︃ 𝐺 ( 𝑥 )1 − 𝐺 ( 𝑥 ) )︃ 𝑛 +1 ⏟ ⏞ A1 − (log[1 − 𝐺 ( 𝑥 )]) 𝑛 +1 ⏟ ⏞ A2 ⎤⎥⎥⎥⎥⎦ . (10)A well-known result on power series raised to a positive integer 𝑁 states that: (︃ ∞ ∑︁ 𝑘 =0 𝑎 𝑘 𝑦 𝑘 )︃ 𝑁 = ∞ ∑︁ 𝑘 =0 𝑐 𝑘 𝑦 𝑘 , (11)where 𝑐 = 𝑎 𝑁 , 𝑐 𝑘 = 𝑘𝑎 ∑︀ 𝑘𝑠 =1 ( 𝑠𝑁 − 𝑘 + 𝑠 ) 𝑎 𝑠 𝑐 𝑘 − 𝑠 for 𝑘 ≥ 𝑁 ∈ N . Setting 𝑁 = 2 𝑛 + 1, 𝑦 = 𝐺 ( 𝑥 ) and 𝑎 𝑘 = 1 for all 𝑘 ≥
0, we can use the result in Equation (11) to rewrite A1 inEquation (10):A1 = 𝐺 ( 𝑥 ) 𝑛 +1 (︃ − 𝐺 ( 𝑥 ) )︃ 𝑛 +1 = 𝐺 ( 𝑥 ) 𝑛 +1 (︃ ∞ ∑︁ 𝑘 =0 𝐺 ( 𝑥 ) 𝑘 )︃ 𝑛 +1 = 𝐺 ( 𝑥 ) 𝑛 +1 ∞ ∑︁ 𝑘 =0 𝑐 ,𝑘 · 𝐺 ( 𝑥 ) 𝑘 = ∞ ∑︁ 𝑘 =0 𝑐 ,𝑘 · 𝐺 ( 𝑥 ) 𝑘 +2 𝑛 +1 , (12)such that 𝑐 , = 1 and 𝑐 ,𝑘 = 𝑘 ∑︀ 𝑘𝑠 =1 (2 𝑠 [ 𝑛 + 1] − 𝑘 ) 𝑐 ,𝑘 − 𝑠 for 𝑘 ≥
1. Equation (11) also allowsus to rewrite A2 in Equation 10 as follows:A2 = (︃ − ∞ ∑︁ 𝑚 =1 𝐺 ( 𝑥 ) 𝑚 𝑚 )︃ 𝑛 +1 = − (︃ ∞ ∑︁ 𝑘 =0 𝐺 ( 𝑥 ) 𝑘 +1 𝑘 + 1 )︃ 𝑛 +1 = − 𝐺 ( 𝑥 ) 𝑛 +1 (︃ ∞ ∑︁ 𝑘 =0 𝐺 ( 𝑥 ) 𝑘 𝑘 + 1 )︃ 𝑛 +1 = − 𝐺 ( 𝑥 ) 𝑛 +1 ∞ ∑︁ 𝑘 =0 𝑐 ,𝑘 · 𝐺 ( 𝑥 ) 𝑘 = − ∞ ∑︁ 𝑘 =0 𝑐 ,𝑘 · 𝐺 ( 𝑥 ) 𝑘 +2 𝑛 +1 (13)where 𝑐 , = 1 and 𝑐 ,𝑘 = 𝑘 ∑︀ 𝑘𝑠 =1 2 𝑠 ( 𝑛 +1) − 𝑘𝑠 +1 𝑐 ,𝑘 − 𝑠 for 𝑘 ≥
1. Now inserting (12) and (13) in(10), we have: 𝐹 𝐺 ,𝐺 ( 𝑥 ) = ∑︁ 𝑖 =1 ∞ ∑︁ 𝑛,𝑘 =0 𝑐 𝑖,𝑛,𝑘 · 𝐺 𝑖 ( 𝑥 ) 𝑘 +2 𝑛 +1 (14)where 𝑐 ,𝑛,𝑘 = ( − / 𝑛 𝑛 !(2 𝑛 +1) √ 𝜋 𝑐 ,𝑘 and 𝑐 ,𝑛,𝑘 = ( − / 𝑛 𝑛 !(2 𝑛 +1) √ 𝜋 𝑐 ,𝑘 . Using Fubini’s theorem on differ-7ntiation we can write the derivative of (14) as follows: 𝑓 𝐺 ,𝐺 ( 𝑥 ) = ∑︁ 𝑖 =1 ∞ ∑︁ 𝑛,𝑘 =0 𝑐 𝑖,𝑛,𝑘 · 𝑔 𝑖,𝑘 +2 𝑛 +1 ( 𝑥 ) (15)where 𝑔 𝑖,𝑘 +2 𝑛 +1 ( 𝑥 ) = ( 𝑘 + 2 𝑛 + 1) 𝑔 𝑖 ( 𝑥 ) 𝐺 𝑖 ( 𝑥 ) 𝑘 +2 𝑛 is the pdf of a random variable from the ex-ponentiated family [23]. Thus, we can say that (15) is the Normal-( 𝐺 , 𝐺 ) pdf (6) expressedas a linear combination of pdfs of exponentiated distributions. Given that 𝑋 is a random variable following a distribution from the normal-( 𝐺 , 𝐺 )class, we can use (15) to write the 𝑟 -th raw moment of 𝑋 as follows: 𝐸 ( 𝑋 𝑟 ) = ∑︁ 𝑖 =1 ∞ ∑︁ 𝑛,𝑘 =0 𝑐 𝑖,𝑛,𝑘 ∫︁ ∞−∞ 𝑥 𝑟 𝑔 𝑖,𝑘 +2 𝑛 +1 ( 𝑥 )d 𝑥 (16)= ∑︁ 𝑖 =1 ∞ ∑︁ 𝑛,𝑘 =0 𝑐 𝑖,𝑛,𝑘 𝐸 ( 𝑌 𝑟𝑖,𝑘 +2 𝑛 +1 ) (17)where 𝑌 𝑖,𝑘 +2 𝑛 +1 follows the exponentiated distribution whose pdf is 𝑔 𝑖,𝑘 +2 𝑛 +1 .Let 𝑄 𝑖 be the quantile function of the baseline 𝐺 𝑖 . Replacing 𝑥 in (16) by 𝑄 𝑖 (︁ 𝑣 /𝑘 +2 𝑛 +1 )︁ we can also represent (17) as: 𝐸 ( 𝑋 𝑟 ) = ∑︁ 𝑖 =1 ∞ ∑︁ 𝑛,𝑘 =0 𝑐 𝑖,𝑛,𝑘 ∫︁ [︁ 𝑄 𝑖 (︁ 𝑣 /𝑘 +2 𝑛 +1 )︁]︁ 𝑟 d 𝑣 . Similarly, one can write the 𝑟 -th incomplete moment of 𝑋 as follows: 𝑚 𝑟 ( 𝑧 ) = ∫︁ 𝑧 −∞ 𝑥 𝑟 𝑓 𝐺 ,𝐺 ( 𝑥 )d 𝑥 = ∑︁ 𝑖 =1 ∞ ∑︁ 𝑛,𝑘 =0 𝑐 𝑖,𝑛,𝑘 𝑚 ⋆𝑟 ( 𝑧 )= ∑︁ 𝑖 =1 ∞ ∑︁ 𝑛,𝑘 =0 𝑐 𝑖,𝑛,𝑘 ∫︁ [ 𝐺 𝑖 ( 𝑧 )] 𝑘 +2 𝑛 +1 [︁ 𝑄 𝑖 (︁ 𝑣 /𝑘 +2 𝑛 +1 )︁]︁ 𝑟 d 𝑣 where 𝑚 ⋆𝑟 ( 𝑧 ) is the 𝑟 -th incomplete moment of 𝑌 𝑖,𝑘 +2 𝑛 +1 mentioned above.The moment generating function (mgf) of 𝑋 is denoted by 𝑀 𝑋 ( 𝑡 ) = 𝐸 (︁ 𝑒 𝑡𝑋 )︁ . It can be8etermined from (15) as: 𝑀 𝑋 ( 𝑡 ) = ∑︁ 𝑖 =1 ∞ ∑︁ 𝑛,𝑘 =0 𝑐 𝑖,𝑛,𝑘 ∫︁ ∞−∞ 𝑒 𝑡𝑥 𝑔 𝑖,𝑘 +2 𝑛 +1 ( 𝑥 )d 𝑥 = ∑︁ 𝑖 =1 ∞ ∑︁ 𝑛,𝑘 =0 𝑐 𝑖,𝑛,𝑘 𝑀 𝑌 𝑘 +2 𝑛 +1 ( 𝑡 ) , where 𝑀 𝑌 𝑘 +2 𝑛 +1 ( 𝑡 ) is the mgf of 𝑌 𝑖,𝑘 +2 𝑛 +1 .Other meaningful quantities, as the characteristic function, the probability-weighted mo-ments, the Rényi entropy and the order statistics can be derived likewise by using (15). Let 𝑋 = ( 𝑥 , . . . , 𝑥 𝑛 ) be a complete random sample of size 𝑛 from the random variable 𝑋 ∼ normal-( 𝐺 , 𝐺 ). Given that 𝜃 = ( 𝜃 , . . . , 𝜃 𝑟 ) ⊤ is the 𝑟 × 𝐺 ( 𝑥 ) = 𝐺 ( 𝑥 | 𝜃 ), 𝜃 = ( 𝜃 𝑟 +1 , . . . , 𝜃 𝑟 + 𝑚 ) ⊤ is the 𝑚 × 𝐺 ( 𝑥 ) = 𝐺 ( 𝑥 | 𝜃 ) and 𝑓 𝐺 ,𝐺 ( 𝑥 ) = 𝑓 𝐺 ,𝐺 ( 𝑥 | 𝜃 ) where 𝜃 = ( 𝜃 , . . . , 𝜃 𝑟 , 𝜃 𝑟 +1 , . . . , 𝜃 𝑟 + 𝑚 ) ⊤ ,we can write the log-likelihood function of 𝑋 as follows: ℓ ( 𝜃 | 𝑋 ) = 𝑛 ∑︁ 𝑖 =1 log {︃ 𝜑 (︃ 𝐺 ( 𝑥 𝑖 )1 − 𝐺 ( 𝑥 𝑖 ) )︃ 𝑔 ( 𝑥 𝑖 )[1 − 𝐺 ( 𝑥 𝑖 )] + 𝜑 (log[1 − 𝐺 ( 𝑥 𝑖 )]) 𝑔 ( 𝑥 𝑖 )1 − 𝐺 ( 𝑥 𝑖 ) }︃ . The solution of the system of equations 𝑈 ( 𝜃 | 𝑋 ) = 𝑟 + 𝑚 provides the maximum likelihoodestimates (MLEs) for 𝜃 , where 𝑟 + 𝑚 is an ( 𝑟 + 𝑚 ) × 𝑈 ( 𝜃 | 𝑋 ) = ∇ 𝜃 ℓ ( 𝜃 | 𝑋 )is the score vector. The elements of 𝑈 ( 𝜃 | 𝑋 ) = ( 𝑢 𝑗 ) ≤ 𝑗 ≤ 𝑟 + 𝑚 are: 𝑢 𝑗 = 𝑛 ∑︁ 𝑖 =1 𝑓 𝐺 ,𝐺 ( 𝑥 𝑖 ) 𝜑 (︃ 𝐺 ( 𝑥 𝑖 )1 − 𝐺 ( 𝑥 𝑖 ) )︃ − 𝐺 ( 𝑥 𝑖 )) [︃ 𝜕𝜕𝜃 𝑗 𝑔 ( 𝑥 𝑖 )+ 𝑔 ( 𝑥 𝑖 )1 − 𝐺 ( 𝑥 𝑖 ) (︃ − 𝐺 ( 𝑥 𝑖 )[1 − 𝐺 ( 𝑥 𝑖 )] )︃ 𝜕𝜕𝜃 𝑗 𝐺 ( 𝑥 𝑖 ) ]︃ , for 1 ≤ 𝑗 ≤ 𝑟 and 𝑢 𝑗 = 𝑛 ∑︁ 𝑖 =1 𝑓 𝐺 ,𝐺 ( 𝑥 𝑖 ) 𝜑 (log[1 − 𝐺 ( 𝑥 𝑖 )])1 − 𝐺 ( 𝑥 𝑖 ) [︃ 𝜕𝜕𝜃 𝑗 𝑔 ( 𝑥 𝑖 )+ (1 + log[1 − 𝐺 ( 𝑥 𝑖 )]) 𝑔 ( 𝑥 𝑖 )1 − 𝐺 ( 𝑥 𝑖 ) 𝜕𝜕𝜃 𝑗 𝐺 ( 𝑥 𝑖 ) ]︃ , for 𝑟 < 𝑗 ≤ 𝑟 + 𝑚 . For testing hypotheses and constructing confidence intervals for 𝜃 , the information ma-9rix 𝐽 ( 𝜃 | 𝑋 ) is needed. The expectation of 𝐽 ( 𝜃 | 𝑋 ), denoted by ℐ 𝜃 , is the expected Fisherinformation matrix. Given that certain conditions of regularity are fulfilled, the quantity √ 𝑛 ( ̂︀ 𝜃 − 𝜃 ) follows approximately a multivariate normal distribution 𝑁 𝑟 + 𝑚 ( 𝑟 + 𝑚 , ℐ − 𝜃 ). Ap-pendix A brings the expression for 𝐽 ( 𝜃 | 𝑋 ).
2. The proposed model
The Weibull cdf is given by 𝐺 𝑊 ( 𝑥 | 𝑘, 𝜆 ) = 1 − 𝑒 − ( 𝑥/𝜆 ) 𝑘 , for 𝑥 ≥ 𝑘 > 𝜆 > 𝐺 and 𝐺 in (5) by 𝐺 𝑊 ( 𝑥 | 𝑘 , 𝜆 ) and 𝐺 𝑊 ( 𝑥 | 𝑘 , 𝜆 ) respectively, we get to thecdf of the Normal-Weibull-Weibull distribution (NWW, for short): 𝐹 𝑁𝑊 𝑊 ( 𝑥 | 𝜃 ) = Φ (︁ 𝑒 ( 𝑥/𝜆 ) 𝑘 − )︁ − Φ (︃ − (︂ 𝑥𝜆 )︂ 𝑘 )︃ , where 𝜃 = ( 𝑘 , 𝜆 , 𝑘 , 𝜆 ) ⊤ . The corresponding pdf can be obtained using (6): 𝑓 𝑁𝑊 𝑊 ( 𝑥 | 𝜃 ) = 𝜑 (︁ 𝑒 ( 𝑥/𝜆 ) 𝑘 − )︁ 𝑘 𝜆 (︂ 𝑥𝜆 )︂ 𝑘 − 𝑒 ( 𝑥/𝜆 ) 𝑘 + 𝜑 (︃ − (︂ 𝑥𝜆 )︂ 𝑘 )︃ 𝑘 𝜆 (︂ 𝑥𝜆 )︂ 𝑘 − . (18)Figure 1 displays some plots of the NWW pdf for different values of the parameters. Thedistribution is able to fit unimodal right-skewed data (top-left chart) and also left-skeweddata (top-right chart). Notice the different shapes of the bimodal curves in the remainingcharts. For instance, in the bottom-left chart, the vertical distance between the modes andthe local minimum in the purple curve is much greater than in the green one. We may alsonotice that 𝜆 and 𝜆 somehow behave like shape parameters, as in the original Weibullbaselines, controlling the shape of the “bells” (compare purple and gray curves to see theeffect of varying 𝜆 ; same for blue and red curves concerning 𝜆 ).10 . . . . x pd f k = 1.3, λ = 2.0, k = 1.5, λ = 1.8 k = 1.4, λ = 2.5, k = 1.3, λ = 0.9 k = 1.5, λ = 3.0, k = 1.5, λ = 1.8 k = 1.3, λ = 3.0, k = 1.2, λ = 3.5 k = 1.2, λ = 1.7, k = 2.5, λ = 2.5 . . . . x pd f k = 1.4, λ = 3.0, k = 6.5, λ = 3.8 k = 2.0, λ = 5.0, k = 8.5, λ = 4.4 k = 2.0, λ = 6.0, k = 4.5, λ = 5.4 k = 1.5, λ = 7.0, k = 4.0, λ = 6.0 k = 1.3, λ = 3.5, k = 6.0, λ = 2.9 k = 2.0, λ = 4.5, k = 7.0, λ = 5.0 . . . . x pd f k = 2.0, λ = 2.2, k = 6.5, λ = 3.8 k = 2.0, λ = 2.2, k = 6.5, λ = 3.5 k = 1.9, λ = 2.5, k = 5.5, λ = 3.1 k = 2.0, λ = 1.5, k = 6.5, λ = 4.1 k = 2.0, λ = 3.0, k = 6.5, λ = 4.1 −2 0 2 4 6 . . . . . . x pd f k = 7.2, λ = 3.1, k = 3.2, λ = 2.5k = 7.0, λ = 3.8, k = 3.2, λ = 2.5k = 6.8, λ = 4.8, k = 3.2, λ = 2.5k = 8.5, λ = 5.2, k = 6.2, λ = 5.9k = 9.4, λ = 7.0, k = 15.9, λ = 6.0 Figure 1: Plots of the Normal-Weibull-Weibull pdf.
3. Simulation
Performing Monte Carlo simulation studies is considerably relevant whenever one wantsto test and confirm assumptions on new statistical methods. In this work, we want toinvestigate the behavior of the estimates of the parameters of the NWW distribution underthe method of maximum likelihood. For this purpose, we used the software R version3.4.4 [24].Initially, we employed the Von Neumann’s acceptance-rejection method [25] to generatepseudo-random samples from the NWW distribution; this simple method requires only thecorresponding pdf 𝑦 = 𝑓 ( 𝑥 ), a minorant and a majorant for 𝑥 and a majorant for 𝑦 . Theprocedure was replicated 10,000 times and at each replication, four different sample sizeswere considered, namely, 𝑛 = 50, 100, 200 and 500. We examined scenarios with fourdifferent values of the parametric vector 𝜃 = ( 𝑘 , 𝜆 , 𝑘 , 𝜆 ) ⊤ , which are presented from thesecond to fifth columns of Tables 1 and 2. 11or each scenario, we calculated the bias and the mean squared error (MSE) as follows:Bias 𝑖 = 110000 ∑︁ 𝑗 =1 (︁ ̂︀ 𝜃 𝑖𝑗 − 𝜃 𝑖 )︁ , MSE 𝑖 = 110000 ∑︁ 𝑗 =1 (︁ ̂︀ 𝜃 𝑖𝑗 − 𝜃 𝑖 )︁ where 𝜃 𝑖 is the 𝑖 -th element of 𝜃 and ̂︀ 𝜃 𝑖𝑗 is the estimate for 𝜃 𝑖 at the 𝑗 -th replication.The global maximum of the log-likelihood function was found by using the L-BFGS-Balgorithm. It is based on the gradient projection and uses a limited memory BFGS matrixto approximate the Hessian of the objective function [26].Besides presenting small values, the desired behavior for both bias and MSE is to decreaseinasmuch as the sample size increases. According to Tables 1 and 2, the values of biasand MSE for all the estimated parameters are small and the greater the sample size, thesmaller the values. Thus, the results presented in this section indicate that the MLEs of theparameters of the NWW distribution are well-behaved when calculated using the L-BFGS-Balgorithm. Table 1: Bias of the estimates under the maximum likelihood method for the NWW model.
Actual value Bias 𝑛 𝑘 𝜆 𝑘 𝜆 ̂︀ 𝑘 ̂︀ 𝜆 ̂︀ 𝑘 ̂︀ 𝜆
50 1.3 2 1.5 1.8 0.54169 0.30227 0.42743 0.291743 1.5 2.8 2.5 0.55544 0.23582 0.63534 0.21462 2.2 6.5 4.1 0.30213 0.1419 1.42262 0.120171.4 1.6 4.8 5.1 0.18208 0.10899 0.86019 0.13781100 1.3 2 1.5 1.8 0.34967 0.24391 0.31488 0.227523 1.5 2.8 2.5 0.43866 0.1199 0.41435 0.138042 2.2 6.5 4.1 0.24521 0.11322 0.89279 0.089231.4 1.6 4.8 5.1 0.13014 0.07167 0.54369 0.09648200 1.3 2 1.5 1.8 0.22345 0.18795 0.24873 0.171873 1.5 2.8 2.5 0.28333 0.04346 0.26341 0.072612 2.2 6.5 4.1 0.161 0.06349 0.57998 0.050341.4 1.6 4.8 5.1 0.09163 0.04958 0.37906 0.069500 1.3 2 1.5 1.8 0.13429 0.11752 0.17774 0.111643 1.5 2.8 2.5 0.16708 0.01806 0.17836 0.049922 2.2 6.5 4.1 0.09747 0.03622 0.34572 0.02771.4 1.6 4.8 5.1 0.06333 0.03001 0.22205 0.0425912 able 2: MSE of the estimates under the maximum likelihood method for the NWW model.
Actual value MSE 𝑛 𝑘 𝜆 𝑘 𝜆 ̂︀ 𝑘 ̂︀ 𝜆 ̂︀ 𝑘 ̂︀ 𝜆
50 1.3 2 1.5 1.8 0.5477 0.13969 0.37049 0.130653 1.5 2.8 2.5 0.43676 0.2146 0.88316 0.114822 2.2 6.5 4.1 0.1976 0.07682 3.60838 0.054331.4 1.6 4.8 5.1 0.05343 0.02546 1.32698 0.03216100 1.3 2 1.5 1.8 0.23088 0.09363 0.17553 0.08093 1.5 2.8 2.5 0.27616 0.09641 0.30237 0.05712 2.2 6.5 4.1 0.16973 0.07868 1.39857 0.064961.4 1.6 4.8 5.1 0.02719 0.00888 0.49768 0.01486200 1.3 2 1.5 1.8 0.0914 0.0588 0.10077 0.045833 1.5 2.8 2.5 0.12419 0.01817 0.11285 0.014722 2.2 6.5 4.1 0.05914 0.02211 0.57197 0.019051.4 1.6 4.8 5.1 0.01352 0.00396 0.23394 0.00751500 1.3 2 1.5 1.8 0.03154 0.02332 0.05049 0.018873 1.5 2.8 2.5 0.04441 0.00052 0.04762 0.003832 2.2 6.5 4.1 0.01723 0.00406 0.19449 0.003151.4 1.6 4.8 5.1 0.00633 0.00143 0.07816 0.00286
4. Results and Discussion
The hourly wind speed data measured at 10 m above ground level were collected by theNational Institute of Meteorology of Brazil (INMET). The anemometers used for measuringthe wind speed (in m/s) are installed in stations located in five cities spread in four statesof the Brazilian Northeastern Region, as illustrated in Figure 2 (blue dots indicate thegeographical position of the stations). Esperantina (denoted by Station 1) is located inthe north part of the State of Piauí. Jaguaruana (denoted by Station 2) is located in themesoregion of Jaguaribe in the State of Ceará. Cabaceiras (denoted by Station 3) andMonteiro (denoted by Station 4) are located in the mesoregion of Borborema, State ofParaíba. Arapiraca (denoted by Station 5) is located in the mesoregion of Agreste, Stateof Alagoas. Table 3 brings further details about the stations and the years of wind dataavailable.As we can see in Table 4, Station 4 has the highest value of the mean among the stationsin the study, whereas the highest value of variance belongs to Station 2. Except for Station3, whose skewness is negative, the remaining stations have different degrees of positive13 igure 2: Northeastern Region of Brazil and geographical position of the stations.Table 3: Details of the regions where the wind speed was measured.
Station Latitude Longitude Altitude (m) Period1 Esperantina 3 ∘ ′ ′′ S 42 ∘ ′ ′′ W 59 2007–20182 Jaguaruana 4 ∘ ′ ′′ S 37 ∘ ′ ′′ W 20 2007–20183 Cabaceiras 7 ∘ ′ ′′ S 36 ∘ ′ ′′ W 382 2008–20184 Monteiro 7 ∘ ′ ′′ S 37 ∘ ′ ′′ W 599 2007–20185 Arapiraca 9 ∘ ′ ′′ S 36 ∘ ′ ′′ W 264 2008–2018skewness. On the other hand, Station 1 has the only positive value of kurtosis and Station 5has the lowest one. Thus, the descriptive statistics indicate that the statistical characteristicsof the wind speed differ from station to station.
Table 4: Descriptive statistics.
St. 𝑛 mean median min max variance skewness kurtosis1 72637 1.56472 1.5 0.1 8.8 0.96082 0.74469 0.615852 71797 3.09401 3 0.1 9.4 2.94352 0.25522 − . − . − . − . − . * ) and Cramér-von Mises (W * ) [28] were also used to com-pare the fitted models. Since these statistics are measures of the difference between theempirical distribution function and the real underlying cdf, it is reasonable to say that thesmaller their values, the better the fit. Table 6 brings the aforementioned goodness-of-fitmeasures for the five cited models fitted to each station.According to Table 6, the distributions NWW and WW present the better fits amongthe competing models for all stations. The four information criteria indicate that NWWpresents a higher performance over WW concerning stations 1, 2, 4 and 5. Regarding station3, the difference between the values of AIC of NWW and WW is not considerable, althoughthe AIC for NWW is slightly smaller in such case; the same behavior states for CAIC, BICand HQIC.Both goodness-of-fit statistics A * and W * (see last ten rows of Table 6) agree with theinformation criteria in pointing NWW and WW as the two better fits. However, they suggestthat NWW outperforms WW in fitting the datasets for all the five stations. To get moreinsight into these results, plots of the fitted densities overlapping the histograms of the windspeed data for stations 1 to 5 and the corresponding cdfs are presented in Figure 3.It is worth pointing out that mixture models are commonly used to fit non-unimodaldatasets, such as those represented by the histograms of stations 2, 3 and 5. Nonetheless,the results attest that the NWW accommodates such data better than the two mixturemodels presented in this study. Furthermore, NWW has one parameter less than NN orWW do. 15 able 5: Estimates and standard errors in parentheses. Distr. Par. St1 St2 St3 St4 St5NWW 𝑘 𝜆 𝑘 𝜆 𝜇 𝜎 𝜇 𝜎 𝑤 𝑘 𝜆 𝑘 𝜆 𝑤 𝜇 𝜎 𝑘 𝜆 able 6: Goodness-of-fit measures. Crit. Distr. St1 St2 St3 St4 St5AIC NWW 189324.5 272976.1 306366.6 286607.1 268213.1NN 197221.2 277168.5 310138.0 289011.9 273287.3WW 189379.2 273075.2 306370.5 286647.3 268244.1N 203235.1 281266.5 315112.0 291251.2 281292.0W 190666.7 278737.7 319347.3 291232.9 277572.1CAIC NWW 189365.2 273016.8 306408.0 286648.1 268253.8NN 197272.2 277219.4 310189.7 289063.1 273338.3WW 189430.2 273126.1 306422.2 286698.5 268295.1N 203255.4 281286.9 315132.7 291271.7 281312.4W 190687.1 278758.0 319368.0 291253.4 277592.5BIC NWW 189361.2 273012.8 306404.0 286644.1 268249.8NN 197267.2 277214.4 310184.7 289058.1 273333.3WW 189425.2 273121.1 306417.2 286693.5 268290.1N 203253.4 281284.9 315130.7 291269.7 281310.4W 190685.1 278756.0 319366.0 291251.4 277590.5HQIC NWW 189335.8 272987.4 306378.0 286618.5 268224.4NN 197235.3 277182.7 310152.3 289026.1 273301.4WW 189393.4 273089.3 306384.8 286661.5 268258.3N 203240.7 281272.2 315117.7 291256.9 281297.7W 190672.4 278743.3 319353.0 291238.6 277577.7A* NWW 71.78 22.43 22.60 19.31 30.98NN 181.13 51.72 54.33 27.79 78.72WW 81.27 31.05 32.44 24.62 50.25N 531.70 230.18 303.82 118.73 494.43W 257.57 471.01 1117.26 337.07 784.70W* NWW 7.73 2.33 2.94 2.58 3.36NN 19.49 5.17 6.60 3.21 7.61WW 9.35 3.59 4.22 3.89 4.82N 72.08 31.04 45.14 16.72 72.83W 33.87 62.25 164.54 46.18 108.58criteria and formal goodness-of-fit statistics, we have good reasons to recommend its useto model similar data in future works. We also encourage practitioners of statistics toinvestigate the modelling benefits of the NWW (and other submodels from the Normal-( 𝐺 , 𝐺 ) class) with respect to data describing different phenomena usually modelled bymixtures. 17 tation 1 wind speed (m/s) pd f . . . . . NWWNNWWNW . . . c d f Station 2 wind speed (m/s) pd f . . . NWWNNWWNW . . . c d f Station 3 wind speed (m/s) pd f . . . NWWNNWWNW . . . c d f Station 4 wind speed (m/s) pd f . . . NWWNNWWNW . . . c d f Station 5 wind speed (m/s) pd f . . . NWWNNWWNW . . . c d f Figure 3: Histograms and fitted densities.
5. Conclusions
An alternative distribution for modelling wind speed data is proposed and some mathe-matical properties of the class that generates it are discussed, like the series representation ofthe pdf, the moments and the moment generating function. The general cdf of the Normal-( 𝐺 , 𝐺 ) class is written as a composition of two baselines and its submodels are identifiableas long as both baseline cdfs are. Such is the case of the NWW distribution.The novel model has four parameters and high flexibility. It is able to fit right-skewedand left-skewed data and its pdf presents unimodal and bimodal shapes.18he Monte Carlo simulation studies indicate that the MLEs of the NWW parametersbehave appropriately when the optimization is performed via the L-BFGS-B algorithm.The modelling gains of the NWW distribution are upheld by the satisfactory resultsconcerning the application to the wind speed data collected in the Northeastern Regionof Brazil. The considered information criteria and the formal goodness-of-fit statistics ofAnderson-Darling and Cramér-von Mises suggest that the proposed model outperforms othercompeting distributions that are commonly employed in wind speed modelling, especiallythe highly competitive mixture model of two Weibull components.We hope that this work may encourage the investigation of the modelling benefits of theNWW (and other identifiable submodels from the Normal-( 𝐺 , 𝐺 ) class) with respect todata describing other natural phenomena usually modelled by mixtures. Appendix A.
The information matrix cited in section 1.4 is given by 𝐽 ( 𝜃 | 𝑋 ) = −∇ 𝜃 ∇ 𝜃 ⊤ ℓ ( 𝜃 | 𝑋 ) = − ( 𝑢 𝑗𝑘 ) ≤ 𝑗 ≤ 𝑟 + 𝑚, ≤ 𝑘 ≤ 𝑟 + 𝑚 where: 𝑢 𝑗𝑘 = 𝑛 ∑︁ 𝑖 =1 𝑓 𝐺 ,𝐺 ( 𝑥 𝑖 ) 𝜑 (︃ 𝐺 ( 𝑥 𝑖 )1 − 𝐺 ( 𝑥 𝑖 ) )︃ − 𝐺 ( 𝑥 𝑖 )) {︃(︃ − 𝐺 ( 𝑥 𝑖 ) − 𝐺 ( 𝑥 𝑖 )[1 − 𝐺 ( 𝑥 𝑖 )] )︃ × (︃ 𝜕𝜕𝜃 𝑘 𝐺 ( 𝑥 𝑖 ) 𝜕𝜕𝜃 𝑗 𝑔 ( 𝑥 𝑖 ) + 𝜕𝜕𝜃 𝑗 𝐺 ( 𝑥 𝑖 ) 𝜕𝜕𝜃 𝑘 𝑔 ( 𝑥 𝑖 ) )︃ + 𝜕 𝜕𝜃 𝑗 𝜕𝜃 𝑘 𝑔 ( 𝑥 𝑖 ) − 𝑔 ( 𝑥 𝑖 )(1 − 𝐺 ( 𝑥 𝑖 )) × [︃ 𝜕𝜕𝜃 𝑗 𝐺 ( 𝑥 𝑖 ) 𝜕𝜕𝜃 𝑘 𝐺 ( 𝑥 𝑖 ) (︃ 𝐺 ( 𝑥 𝑖 )(1 − 𝐺 ( 𝑥 𝑖 )) + 4 𝐺 ( 𝑥 𝑖 )1 − 𝐺 ( 𝑥 𝑖 ) + 1 )︃ + 𝐺 ( 𝑥 𝑖 )1 − 𝐺 ( 𝑥 𝑖 ) × 𝜕 𝜕𝜃 𝑗 𝜕𝜃 𝑘 𝐺 ( 𝑥 𝑖 ) ]︃ + 𝑔 ( 𝑥 𝑖 ) 𝐺 ( 𝑥 𝑖 )(1 − 𝐺 ( 𝑥 𝑖 )) (︃ 𝜕𝜕𝜃 𝑗 𝐺 ( 𝑥 𝑖 ) + 𝜕𝜕𝜃 𝑘 𝐺 ( 𝑥 𝑖 ) )︃ + 2 𝑔 ( 𝑥 𝑖 )(1 − 𝐺 ( 𝑥 𝑖 )) (︃ − 𝐺 ( 𝑥 𝑖 )(1 − 𝐺 ( 𝑥 𝑖 )) )︃ 𝜕𝜕𝜃 𝑘 𝐺 ( 𝑥 𝑖 ) 𝜕𝜕𝜃 𝑗 𝐺 ( 𝑥 𝑖 ) + 𝜕 𝜕𝜃 𝑗 𝜕𝜃 𝑘 𝐺 ( 𝑥 𝑖 ) × 𝑔 ( 𝑥 𝑖 )1 − 𝐺 ( 𝑥 𝑖 ) }︃ − 𝑛 ∑︁ 𝑖 =1 𝑓 𝐺 ,𝐺 ( 𝑥 𝑖 ) 𝜑 (︃ 𝐺 ( 𝑥 𝑖 )1 − 𝐺 ( 𝑥 𝑖 ) )︃ − 𝐺 ( 𝑥 𝑖 )) × [︃ 𝜕𝜕𝜃 𝑗 𝑔 ( 𝑥 𝑖 ) + (︃ 𝑔 ( 𝑥 𝑖 )1 − 𝐺 ( 𝑥 𝑖 ) − 𝑔 ( 𝑥 𝑖 ) 𝐺 ( 𝑥 𝑖 )(1 − 𝐺 ( 𝑥 𝑖 )) )︃ 𝜕𝜕𝜃 𝑗 𝐺 ( 𝑥 𝑖 ) ]︃ × [︃ 𝜕𝜕𝜃 𝑘 𝑔 ( 𝑥 𝑖 ) + (︃ 𝑔 ( 𝑥 𝑖 )1 − 𝐺 ( 𝑥 𝑖 ) − 𝑔 ( 𝑥 𝑖 ) 𝐺 ( 𝑥 𝑖 )(1 − 𝐺 ( 𝑥 𝑖 )) )︃ 𝜕𝜕𝜃 𝑘 𝐺 ( 𝑥 𝑖 ) ]︃ , for 1 ≤ 𝑗 ≤ 𝑟, ≤ 𝑘 ≤ 𝑟 ;19 𝑗𝑘 = 𝑛 ∑︁ 𝑖 =1 − 𝑓 𝐺 ,𝐺 ( 𝑥 𝑖 ) 𝜑 (︃ 𝐺 ( 𝑥 𝑖 )1 − 𝐺 ( 𝑥 𝑖 ) )︃ − 𝐺 ( 𝑥 𝑖 )) 𝜑 (log[1 − 𝐺 ( 𝑥 𝑖 )])1 − 𝐺 ( 𝑥 𝑖 ) (︃ 𝜕𝜕𝜃 𝑘 𝑔 ( 𝑥 𝑖 )+ [︃ 𝑔 ( 𝑥 𝑖 )1 − 𝐺 ( 𝑥 𝑖 ) − 𝑔 ( 𝑥 𝑖 ) 𝐺 ( 𝑥 𝑖 )(1 − 𝐺 ( 𝑥 𝑖 )) ]︃ 𝜕𝜕𝜃 𝑘 𝐺 ( 𝑥 𝑖 ) )︃ (︃ 𝜕𝜕𝜃 𝑗 𝑔 ( 𝑥 𝑖 ) + (1 + log[1 − 𝐺 ( 𝑥 𝑖 )]) × 𝑔 ( 𝑥 𝑖 )1 − 𝐺 ( 𝑥 𝑖 ) 𝜕𝜕𝜃 𝑗 𝐺 ( 𝑥 𝑖 ) )︃ , for 𝑟 < 𝑗 ≤ 𝑟 + 𝑚, ≤ 𝑘 ≤ 𝑟 ; 𝑢 𝑗𝑘 = 𝑛 ∑︁ 𝑖 =1 − 𝑓 𝐺 ,𝐺 ( 𝑥 𝑖 ) 𝜑 (︃ 𝐺 ( 𝑥 𝑖 )1 − 𝐺 ( 𝑥 𝑖 ) )︃ − 𝐺 ( 𝑥 𝑖 )) 𝜑 (log[1 − 𝐺 ( 𝑥 𝑖 )])1 − 𝐺 ( 𝑥 𝑖 ) (︃ 𝜕𝜕𝜃 𝑗 𝑔 ( 𝑥 𝑖 )+ [︃ 𝑔 ( 𝑥 𝑖 )1 − 𝐺 ( 𝑥 𝑖 ) − 𝑔 ( 𝑥 𝑖 ) 𝐺 ( 𝑥 𝑖 )(1 − 𝐺 ( 𝑥 𝑖 )) ]︃ 𝜕𝜕𝜃 𝑗 𝐺 ( 𝑥 𝑖 ) )︃ (︃ 𝜕𝜕𝜃 𝑘 𝑔 ( 𝑥 𝑖 ) + (1 + log[1 − 𝐺 ( 𝑥 𝑖 )]) × 𝑔 ( 𝑥 𝑖 )1 − 𝐺 ( 𝑥 𝑖 ) 𝜕𝜕𝜃 𝑘 𝐺 ( 𝑥 𝑖 ) )︃ , for 1 ≤ 𝑗 ≤ 𝑟, 𝑟 < 𝑘 ≤ 𝑟 + 𝑚 ; 𝑢 𝑗𝑘 = 𝑛 ∑︁ 𝑖 =1 𝑓 𝐺 ,𝐺 ( 𝑥 𝑖 ) 𝜑 (log[1 − 𝐺 ( 𝑥 𝑖 )])1 − 𝐺 ( 𝑥 𝑖 ) [︃ − 𝐺 ( 𝑥 𝑖 )]1 − 𝐺 ( 𝑥 𝑖 ) (︃ 𝜕𝜕𝜃 𝑗 𝐺 ( 𝑥 𝑖 ) 𝜕𝜕𝜃 𝑘 𝑔 ( 𝑥 𝑖 )+ 𝜕𝜕𝜃 𝑘 𝐺 ( 𝑥 𝑖 ) 𝜕𝜕𝜃 𝑗 𝑔 ( 𝑥 𝑖 ) )︃ + 𝜕 𝜕𝜃 𝑗 𝜕𝜃 𝑘 𝑔 ( 𝑥 𝑖 ) + 𝜕𝜕𝜃 𝑗 𝐺 ( 𝑥 𝑖 ) 𝜕𝜕𝜃 𝑘 𝐺 ( 𝑥 𝑖 ) 𝑔 ( 𝑥 𝑖 )(1 − 𝐺 ( 𝑥 )) × (︁ [1 − 𝐺 ( 𝑥 𝑖 )] + 3 log[1 − 𝐺 ( 𝑥 𝑖 )] )︁ + (1 + log[1 − 𝐺 ( 𝑥 𝑖 )]) 𝜕 𝜕𝜃 𝑗 𝜕𝜃 𝑘 𝐺 ( 𝑥 𝑖 ) × 𝑔 ( 𝑥 𝑖 )1 − 𝐺 ( 𝑥 𝑖 ) ]︃ − 𝑛 ∑︁ 𝑖 =1 𝑓 𝐺 ,𝐺 ( 𝑥 𝑖 ) 𝜑 (log[1 − 𝐺 ( 𝑥 𝑖 )])(1 − 𝐺 ( 𝑥 𝑖 )) [︃ (1 + log[1 − 𝐺 ( 𝑥 𝑖 )]) 𝑔 ( 𝑥 𝑖 )1 − 𝐺 ( 𝑥 𝑖 ) × 𝜕𝜕𝜃 𝑗 𝐺 ( 𝑥 𝑖 ) + 𝜕𝜕𝜃 𝑗 𝑔 ( 𝑥 𝑖 ) ]︃ [︃ (1 + log[1 − 𝐺 ( 𝑥 𝑖 )]) 𝑔 ( 𝑥 𝑖 )1 − 𝐺 ( 𝑥 𝑖 ) 𝜕𝜕𝜃 𝑘 𝐺 ( 𝑥 𝑖 ) + 𝜕𝜕𝜃 𝑘 𝑔 ( 𝑥 𝑖 ) ]︃ , for 𝑟 < 𝑗 ≤ 𝑟 + 𝑚, 𝑟 < 𝑘 ≤ 𝑟 + 𝑚. References [1] J. C. H. Araújo, W. F. d. Souza, A. J. d. A. Meireles, C. Brannstrom, Sustainability challenges of windpower deployment in coastal Ceará state, Brazil, Sustainability 12 (2020) 5562.[2] S. Perkin, D. Garrett, P. Jensson, Optimal wind turbine selection methodology: A case-study forBúrfell, Iceland, Renewable Energy 75 (2015) 165–172.[3] A. Eltamaly, Design and implementation of wind energy system in Saudi Arabia, Renewable Energy60 (2013) 42–52.
4] U. Ilhan, Y. M. Kantar, Analysis of some flexible families of distribution for estimation of wind speeddistributions, Applied Energy 89 (2012) 355–367.[5] S. Pishgar-Komleh, A. Keyhani, P. Sefeedpari, Wind speed and power density analysis based on Weibulland Rayleigh distributions (a case study: Firouzkooh county of Iran), Renewable and SustainableEnergy Reviews 42 (2015) 313 – 322.[6] B. Safari, Modeling wind speed and wind power distributions in Rwanda, Renewable and SustainableEnergy Reviews 15 (2011) 925–935.[7] D. Weisser, A wind energy analysis of Grenada: an estimation using the ‘Weibull’ density function,Renewable Energy 28 (2003) 1803 – 1812.[8] R. Kollu, S. Rayapudi, S. Narasimham, K. Pakkurthi, Mixture probability distribution functions tomodel wind speed distributions, International Journal of Energy and Environmental Engineering 3(2012).[9] T. G. Bali, P. Theodossiou, Risk measurement performance of alternative distribution functions, Jour-nal of Risk and Insurance 75 (2008) 411–437.[10] B. E. Hansen, Autoregressive conditional density estimation, International Economic Review 35 (1994)705–730.[11] E. C. Morgan, M. Lackner, R. M. Vogel, L. G. Baise, Probability distributions for offshore wind speeds,Energy Conversion and Management 52 (2011) 15 – 26.[12] K. Mohammadi, O. Alavi, J. Mcgowan, Use of Birnbaum-Saunders distribution for estimating windspeed and wind power probability distributions: A review, Energy Conversion and Management 143(2017) 109–122.[13] Z. Qin, W. Li, X. Xiong, Estimating wind speed probability distribution using kernel density method,Electric Power Systems Research 81 (2011) 2139–2146.[14] B. Hu, Y. Li, H. Yang, H. Wang, Wind speed model based on kernel density estimation and itsapplication in reliability assessment of generating systems, Journal of Modern Power Systems andClean Energy 5 (2016).[15] Q. Han, S. Ma, T. Wang, F. Chu, Kernel density estimation model for wind speed probability dis-tribution with applicability to wind energy assessment in China, Renewable and Sustainable EnergyReviews 115 (2019) 109387.[16] O. Jaramillo, M. Borja, Wind speed analysis in La Ventosa, Mexico: A bimodal probability distributioncase, Renewable Energy 29 (2004) 1613–1630.[17] T. Chang, Estimation of wind energy potential using different probability density functions, AppliedEnergy 88 (2011) 1848–1856.[18] S. Akdag, H. Bagiorgas, G. Mihalakakou, Use of two-component Weibull mixtures in the analysis ofwind speed in the Eastern Mediterranean, Applied Energy 87 (2010) 2566–2573.[19] J. Carta, P. Ramírez, Analysis of two-component mixture Weibull statistics for estimation of windspeed distributions, Renewable Energy 32 (2007) 518–531.[20] G. McLachlan, D. Peel, Finite Mixture Models, Wiley Interscience, 2000.[21] C. R. Brito, L. C. Rego, W. R. Oliveira, F. Gomes-Silva, Method for generating distributions and classesof probability distributions: the univariate case, Hacettepe Journal of Mathematics and Statistics 48(2019) 897–930.
22] A. Alzaatreh, C. Lee, F. Famoye, A new method for generating families of continuous distributions,Metron 71 (2013) 63–79.[23] G. S. Mudholkar, D. K. Srivastava, Exponentiated Weibull family for analyzing bathtub failure-ratedata, IEEE transactions on reliability 42 (1993) 299–302.[24] R Core Team, R: A Language and Environment for Statistical Computing, R Foundation for StatisticalComputing, Vienna, Austria, 2018. URL: .[25] J. Von Neumann, Various techniques used in connection with random digits, Applied MathematicsSeries 12, National Bureau of Standards, Washington, DC, USA, 1951.[26] R. H. Byrd, P. Lu, J. Nocedal, C. Zhu, A limited memory algorithm for bound constrained optimization,SIAM Journal on Scientific Computing 16 (1995) 1190–1208.[27] H. D. Nguyen, D. Wang, G. J. McLachlan, Randomized mixture models for probability density ap-proximation and estimation, Information Sciences 467 (2018) 135 – 148.[28] G. Chen, N. Balakrishnan, A general purpose approximate goodness-of-fit test, Journal of QualityTechnology 27 (1995) 154–161..[25] J. Von Neumann, Various techniques used in connection with random digits, Applied MathematicsSeries 12, National Bureau of Standards, Washington, DC, USA, 1951.[26] R. H. Byrd, P. Lu, J. Nocedal, C. Zhu, A limited memory algorithm for bound constrained optimization,SIAM Journal on Scientific Computing 16 (1995) 1190–1208.[27] H. D. Nguyen, D. Wang, G. J. McLachlan, Randomized mixture models for probability density ap-proximation and estimation, Information Sciences 467 (2018) 135 – 148.[28] G. Chen, N. Balakrishnan, A general purpose approximate goodness-of-fit test, Journal of QualityTechnology 27 (1995) 154–161.