Stochastic Extreme Wind Speed Modeling and Bayes Estimation under the Inverse Rayleigh Distribution

: Inverse Rayleigh probability distribution is shown in this paper to constitute a valid model for characterization and estimation of extreme values of wind speed, thus constituting a useful tool of wind power production evaluation and mechanical safety of installations. The first part of this paper illustrates such a model and its validity to interpret real wind speed field data. The inverse Rayleigh model is then adopted as the parent distribution for assessment of a dynamical “risk index” defined in terms of a stochastic Poisson process, based upon crossing a given value with part of the maximum value of wind speed on a certain time horizon. Then, a novel Bayes approach for the estimation of such an index under the above model is proposed. The method is based upon assessment of prior information in a novel way which should be easily feasible for a system engineer, being based upon a model quantile (e.g., the median value) or, alternatively, on the probability that the wind speed is greater than a given value. The results of a large set of numerical simulation—based upon typical values of wind-speed parameters—are reported to illustrate the efficiency and the precision of the proposed method, also with hints to its robustness. The validity of the approach is also verified with respect to the two different methods of assessing the prior information.

Before introducing the above model, its generalization, and their application to EWS analyses, it is recalled that studies on wind speed statistical distribution are more and more largely widespread, since wind energy production strongly depends on such distribution. For instance, the annual electricity production (AEP) indicator is frequently used to quantify the electricity output expected in a particular system. The following expression can be considered for the evaluation of the AEP indicator [15]: , 0 x F x e x        (2) to which the following probability density function (PDF) corresponds: In particular, in this paper, a practical method for the Bayes estimation [7][8][9][10][11][12] of the IRD model for EWS is proposed, in the framework of the assessment of a dynamical "risk index", to be outlined in the following. The method is also extended to a generalization of the IRD obtained by means of a continuous mixture of its parameter. The mixture of the scale parameter has been recently adopted also for other models [19,20] in order to allow the randomness on environmental conditions of the windfarm under study.
Such a model is adopted as the parent distribution for assessment by means of a novel Bayes approach of a risk index model accounting for the most severe EWS amplitudes as disturbances for wind tower safety [18]. For the definition of such a risk index, this paper follows a probabilistic approach established in some basic literature on stochastic processes but probably not so widespread as it would be opportune in wind energy applications. Under such a stochastic approach-originated by classical studies on risk analysis as the authoritative book by Thompson [25]-safety is evaluated probabilistically in terms of a suitable stochastic counting process, measuring risk in terms of the occurrence of disturbance levels (i.e., EWS amplitude, in the case under study) higher than a prefixed threshold value. The risk index (RI) here adopted is defined with reference to a given time interval (0,t) as the probability that a given EWS amplitude is With reference to such an RI, this paper also discusses the problem of its statistical estimation in the presence of a few field data with the purpose of efficient assessment of wind tower reliability and safety, not forgetting that such kinds of analyses also allow to obtain useful information for wind power generation statistics, as above outlined.
In brief, the novelty of the proposed procedure discussed in this manuscript is twofold: 1. The use of a new wind-speed model, the IRD one, to characterize extreme wind-speed distribution; such a model is proposed for the first time in the literature on the object, and its validity is illustrated with reference to real wind-speed data. 2. The use of a new suitable and practical Bayes approach is introduced for the estimation of the above RI with reference to the above IRD wind-speed model. The numerical part of the manuscript shows that such a Bayesian statistical inference approach for the estimation of the above RI from available data possesses a very good level of efficiency and precision, as proven by the very low Root Mean Square errors (RMSE).
It is apparent that such results constitute novelties with respect to the current state-of-art on the subject.
The application of the proposed wind-speed model to real field data is reported in Appendix B. Moreover, in order to better highlight the application of the proposed wind-speed model to real field data, in Section 4, a brief summary of the steps which are necessary to perform the application of the proposed model to real field data is reported.
The rest of this paper is organized as follows. In Section 2, the basic steps leading to the RI of (4) are illustrated, together with the key features of the IRD model, characterizing the CDF F(z) of (4). In Section 3, a summary of the methods for IRD classical estimation is outlined; then, in Section 4, a suitable and practical Bayes approach is introduced for the estimation of the RI, with reference to the IRD model. The method constitutes a "practical" Bayesian statistical inference approach for the estimation of the EWS distribution from available data by using the model quantile or the model CDF as the proper input for the prior assessment of the model. The feasibility, efficiency, and precision of such a Bayes estimation method is finally illustrated in the successive and final Sections 5-7 by numerical examples and a broad amount of Monte Carlo simulations performed on the basis of typical values of real EWS data. Such results, illustrated by a large number of tables reported in Sections 5 and 6, assess the good performance of the proposed approach, especially when compared to "classical" estimation methods, such as the maximum likelihood estimation [26][27][28][29][30][31][32][33]. Appendix A reports some details on the derivation of the RI of (6), while Appendix B reports some details on the IRD model, also with examples of their applications to real EWS data. Other Appendixes are also reported for recalling some basic properties of the statistical distributions mentioned in this paper (Appendixes C-E), while Appendix F recalls useful relations concerning the gamma Probability density function (PDF) as a prior model for the scale parameter of the IRD distribution.

A Stochastic Risk Analysis for Extreme Wind Speed Data under the IRD Model
For the purpose of developing a stochastic risk analysis for extreme wind speed data, a stochastic approach established in some literature on extremal processes for safety applications [25] is adopted, also referring to specific application of such extremal processes to wind data [28][29][30]. Let a threshold value z be assigned to the EWS values for safety reasons. If z is large enough, the occurrences of wind speed values which exceed such thresholds (or "exceedances" in the following), are rare events which, under general hypotheses, can be described by the well-known Poisson process, i.e., by the probability law given, as well known, by the following [25,31]: where -N(t) is a process which counts the number of EWS values of exceeding z in the time interval (0,t); -p(k,t) is the probability that the number of exceedances in (0,t), N(t), and attains the nonnegative integer value k; - > 0 is the mean frequency of the process N(t), i.e., the expected number of exceedances, over the threshold value z for a unit time or Indeed, mean (or expected) value and variance of N(t) are expressed, as well known [25,31], by the following: Moreover, a stronger result holds, almost surely (a.s.) in terms of the so called "strong law" of the renewal processes [25,31]: The following developments, leading to the analytical expression of the RI, are discussed with more details in Appendix A. Let u(t,z) denote the probability that z is never exceeded over (0,t), i.e., a "safety" index.
Let the random variable "EWS amplitude occurring at time Tk of the kth exceedance" (k = 1,2, …) be denoted as Zk. If the RV Zk is (as typically assumed in the literature [18]) assumed to be statistically independent and identically distributed with the common, time-independent CDF, then, the following compact expression can be obtained for the above function u(t,z), as shown in Appendix A: Therefore, with the RI defined-as above reported-as the probability of the event that the EWS amplitude z is exceeded at least once, since such event is complementary of the event "z is never exceeded over (0,t)", the above defined RI is expressed in terms of the function u(t,z) in (9) by the following: By its very definition and clearly shown in its expression, such RI is intrinsically dynamic, being dependent on the chosen value of time horizon. If the above limit value z is chosen as the maximum Appl. Sci. 2020, 10, 5643 5 of 35 value of the EWS that the structure can resist, then the RI is an efficient measure of the failure probability the structure. As discussed in [30], this kind of measure appears to be very significant for predicting the risk values associated with gust occurrence and their effects on real wind power systems.
In this framework, this paper discusses the problem of the statistical estimation of wind probabilistic distribution with the primary purpose of an efficient assessment of wind tower reliability and safety and of wind power generation statistics. This is performed by assuming an IRD model as the reference model of the CDF F(z) in (10) and then by performing an estimation of the RI based upon the Poisson exceedances process and the above IRD model.
Let us first briefly illustrate the use of the IRD model for the EWS statistics and then the expression of the RI in (10) in terms of such IRD model. In Figure 1, three typical IRD PDF and CDF curves corresponding to the analytical expressions of (1) and (2)   The expression of the RI assuming the IRD as the reference model becomes the following: A set of possible curves of the RI in time are depicted in Figure 2, assuming a value α = 90 (corresponding to a median value of 11 m/s for the EWS); a mean frequency of the process N(t):  = 10/year; and three different threshold values z of (22, 33, and 44) m/s, i.e., two, three, and four times the median value. From the three curves, it is apparent-as simply deducible from the meaning of the RI as well as from the general expression of (10)-that the RI always increases with time (given the scale parameter value of α and the threshold value z), while it decreases for a given time interval if z increases (allowing for a greater maximum value z implies decreasing the probability that such a value is exceeded). Finally, as some figures will confirm in a following section-holding fixed the values of time t and threshold value z-the RI always increases with increasing values of scale parameter α, since both median and mean values (as also apparent from Figure 1) increase with α. Therefore, the PDF curve shifts to the right for increasing values of the parameter α, thus yielding a higher probability that a given threshold value z is exceeded, which means a higher value of the risk index in any given time interval.

Hints to the Classical Inference Methods of Estimation and Identification of the IRD Model
While the present paper is focused upon new Bayes estimation methods, the topic of classical estimation has been already discussed in the literature [1][2][3][4][5][6]31], so in this brief section, only hints to the classical ML (Maximum Likelihood) estimation method of the IRD distribution are illustrated with some further discussion on the model identification. The extension to the RI estimation is straightforward and is briefly illustrated here for the ML method: by virtue of the well-known property of invariance of the ML estimation, if  * and F * are respectively the Maximum Likelihood estimates (MLE) of the quantities  and F in (10), then the MLE of the RI in (10) is expressed for any fixed value of time t and threshold z by the following: The MLE of the Poisson frequency is also well known in terms of the number of exceedances in a given interval (0,t) denoted by the random variable M (which counts the number of WS values in the given sample which exceed the value z) and is given by the following [31]: Then, let us consider the MLE of the CDF F(z) appearing in (12) for an IRD model. First, we consider the MLE of parameter α of the IRD model in Equations (1) (14) Equating to zero the partial derivative of the above function with respect to parameter α *, the MLE of such parameter is as follows: (15) It is remarked that each RV 1/Xk follows a Raleigh distribution, so its square obeys an exponential law [1][2][3][4][5][6]31].
It is known that a sum of independent exponentially distributed random variates gives a gamma distributed variable, giving the probability density function of the sum in the denominator of (15) as a gamma PDF with parameters: N (shape) and α (scale).
Finally, since-as already recalled-the MLE is invariant, the MLE of the CDF F(z) appearing in (12) is provided, looking at (2), by the following: Let us give a brief hint to the IRD model identification, which is very simple given the basic properties of the inverse Rayleigh model, i.e., the fact that it describes the distribution of X = 1/Y, with Y being a Rayleigh RV, for which identification is a well-established topic. Therefore, the goodness of fit tests for assessing the validity of an IRD model for a given set of EWS data where yk = 1/xk.

A Premise. Bayes Estimation of the Frequency of Exceedances
Bayes method of inference or estimation-as illustrated in its general theory in [26,27] and for some EWS applications in [12,18,20]-is proposed here for the RI under the IRD model. This kind of inference implies the assessment of prior information on the model parameters. In the cases under study, there are two parameters to be estimated: 1) The mean frequency  of the process N(t) in (5), i.e., the expected number of exceedances of the EWS per unit time according (6). Briefly,  is denoted as frequency of exceedances.
2) The parameter a of the IRD model of (2) and (3) or any other parameter related to these ones.
The above phrase "any other parameter related to the parameter a of the IRD model" means, for instance, the CDF of the above model, denoted as F(x|α) in Section 1.
Indeed, the mentioned relation represents a one-to-one correspondence between the parameter α vs. the IRD CDF F(x) for any given EWS value x. In the proposed approach, the prior information can be assessed equivalently both in terms of the "basic" parameters (α) or of the corresponding CDF of their respective model.
It is remarked that assessing the prior information in terms of the basic parameters (α) is also equivalent to assessing the prior information in terms of the of the quantiles of the model, which are related to α.
Whatever the choice, depending on the available information or past experience, for assessing the above information, the ultimate purpose of the inference is to evaluate the Bayes estimates of the above defined RI, in which time t and threshold z are to be considered known constants: It may be convenient in order to simplify some of the following formulations and considering that the time horizon t and threshold z are to be held fixed to express more concisely the RI as follows: Appl. Sci. 2020, 10, 5643 are respectively the "expected number of exceedances" (ENE) in (0,t), and the complementary CDF (CCDF) (or "survival function") of the EWS amplitudes evaluated at z, i.e., S represents the probability that the threshold value z is exceeded by any single EWS values. Since S(z) = 1 − F(z|α) for the IRD model, the RI to be estimated can be considered as functions of two RVs, R = p(, αMoreover, the two input RVs ( and α) for the RI estimation will be considered as statistically independent of each other.
In general terms, denoting by T = (X,Y) the 2-dimensional vector of the two input RV, the problem of the Bayes estimation of the function R = p(X,Y) = p(T) is briefly outlined in the following, leaving the details to the well-established literature on the subject [26,27].
The Bayes estimator (BE) of R requires a joint prior PDF of the RV (X,Y), denoted by p(x,y), expressing prior belief or information about the above parameters. Such a prior PDF on (X,Y) implies-as will be shown for the cases under study-a prior PDF on the RI R and is based upon the well-known Bayes formula used to derive the posterior PDF q(r|D), conditional to data D (the observed EWS data in this application), which has been collected: where L is the Likelihood function [26,27], given the observed data w, and C has the following expression: Finally, as well known [26,27], the best Bayes estimate-in the mean square error sense-of R is given by its "posterior mean", evaluable through the posterior PDF of R, q(r|D): In practice, due to the nature of R, i.e., the fact that it is a probability, the above integration is confined to the interval (0,1).
Analogous formulations applied to the estimation of the single parameters ( and α will be illustrated in the following two subsections. Therefore, before considering the relevant features of the above general formulations to the practical cases of the IRD model, the Bayes estimation of  is illustrated here in a brief summary since, differently from the estimation of α, it is well established in the literature [25,31]. As above recalled, f is the mean frequency of the process N(t) counting the number of exceedances: these are random events which occur at random time,: T1, T2, …,Tn, by well-known properties of the Poisson process; such random times instants are gamma random variables [25,31], and the times between two subsequent events, k = Tk+1 − Tk (assuming and T0 = 0) are Exponential random variables with mean value  = 1/. It is well-known [28] that a gamma PDF, also recalled in Appendix C, apart from being very flexible for describing positive random variables is a valid choice for the prior PDF of the RV , since it is a "conjugate" PDF [26,27], i.e., the posterior PDF is still a gamma PDF. Thus, in this paper, a gamma prior PDF-denoted as "Gampdf (;p0,0)"-for  is assumed: Appl. Sci. 2020, 10, 5643 9 of 35 Then, let a sample of EWS data be observed in a given time interval (0,t), and the number of exceedances of z be denoted by M as in the previous section. The likelihood function of the sample, i.e., of data D, is the probability that M exceedances are observed in time interval (0,t), computed by mean of the above Poisson PDF, It is then straightforward to deduce that the product of the gamma prior PDF of (23) and the likelihood function of (24) yields a posterior PDF of , q(|D), which is still a gamma PDF, i.e., a Gampdf (;p0,0) with updated parameters [26,27]: Therefore, the BE of  being expressed, the mean of its posterior PDF, is computed as follows (see also Appendix C for the mean of a gamma PDF): Of course, the same arguments apply to the Expected number of exceedances (ENE) RV, W = t. It is easily proved [28] that, if  is a gamma RV with parameters  and , then the RV W = t is still a gamma RV with parameters p' and 'simply related to p and as follows: It is also straightforward to verify that the gamma conjugate PDF for  also implies a gamma conjugate PDF for W. Thus, in the paper, a gamma prior PDF Gampdf (; p0, 0 t) for W is assumed with parameters deduced from the gamma prior PDF of As above said, the posterior PDF of W is still a gamma prior PDF Gampdf(; p1, 1), and the BE of W is as follows: A final remark seems opportune before concluding this subsection devoted to the estimation of  The assessment of the prior PDF on S is dealt with in the following subsections. Once the posterior PDF of S, qS(s), is also obtained, it can be used together with the posterior PDF of W, qW(w), for determining the joint posterior PDF of W and S, which-assuming the statistical independence of W and S-will be given by the following: Then, an alternative (but equivalent) formulation of the BE of R, instead of (22), can be made in terms of joint posterior PDF of W and S as follows: It is remarked that, in some cases, the above double integration can be avoided as well as the explicit evaluation of the single integration of (22). First, let us consider the prior PDF assessment.
An opportune method which will be pursued in the following for a feasible estimation of the RI is the one assessing the prior PDF on the CCDF S of (18) in order for-at least in principle-easy evaluation of the implied prior PDF of the product, V = WS, and consequently of the RI: Therefore, an alternative way to evaluate the BE of R may be obtained by means of a single integrations, as follows: if the PDF of V is computed, then the BE of R can be computed in terms of such a posterior PDF of V, q(v|D), with It is noticed that the above formulation establishes the BE of the RI in terms of a "Laplace Transform" of the posterior PDF of V.
A useful analytical tool, given that a gamma prior PDF is assumed for W, is the assessment of the prior PDF for the probability S(z) = 1 − F(z) by mean of a model that possesses hopefully the two following properties: 1) First, the obvious property that such a model is very manageable analytically and rather flexible for a RV with support in the interval (0,1). 2) Moreover, in view of (32), it is opportune that it can render the above product V = WS manageable analytically, i.e., that the PDF of the RV V = WS and R = 1−exp(−V) could be easy computable, at least approximatively, by analytical formulations or simple numerical tools.
This seems to be the case if a beta prior PDF [28] is assumed for S. It is indeed by far the most widely adopted model in Bayes estimation for a RV which is a probability [26,27] since it is well-known that it constitutes the most practical and flexible model to express a prior PDF for a RV confined in the (0,1) interval. It is recalled (see Appendix D for details) that the beta PDF has the following form, in which s is a generic value of the RV S, and a and b are positive parameters: The above formulation is denoted by the symbol "p(s)", as this is adopted in the paper for all prior PDFs. It is possible-as it will be shown in the following-that also the posterior PDF q(s|D) is a beta PDF (or it may be approximated by a beta PDF) which simplifies all the subsequent computations. However, let us continue by examining the prior assessment.
Under the above assumption, the product V = WS has a computable PDF here denoted as a "beta-gamma product" (BGP) PDF, as reported in Appendix E. In some cases, such a product still has a gamma PDF, and in this case also, the RI R possesses a computable prior PDF, i.e., a Negative Log-Gamma distribution (NLGD) PDF (in Appendix C). In some cases, discussed in the following, this also holds for the posterior PDF, at least approximately, as illustrated in the course of the numerical application. This implies that the BE of R, evaluated as the posterior mean of R, is obtained by well-known formulations (reported in the Appendix C), thus avoiding the above integrations.

Bayes Estimation of the Parameters of the IRD Model
The Bayes estimation of the IRD model implies to estimate the parameter α in the expression of the IRD CDF: For this estimation, two methods for are considered in the paper in order to establish appropriate prior distributions: 1) A gamma prior distribution on α, which is a conjugate distribution (see Appendix F and [26,27]). 2) A negative exponential beta distribution (NEBD) on α, (see Appendix D) which is not a conjugate distribution, but has the advantage of implying a beta distribution on S.
The choice of a gamma prior distribution is already reported in the literature and will not be further discussed here but will directly be used in the course of the numerical applications. Before discussing such method, it is recalled that, in the Bayesian setting, also the CDF and other parameter (e.g., moments and quantiles), being functions of α, are themselves RV, described by appropriate distributions. The same also holds for the CCDF: for which estimation, as outlined above, can directly allow for the estimation of the RI. The assessment of the prior information in terms of the parameter α, a quantile, the CDF, or any other parameter is a matter of convenience, depending on which information is available. Generally, it is advisable to base the estimation upon a parameter which possesses a clear physical meaning, which is not the case for α, which appears to be difficult to interpret. The information on a given quantile should be more easily available and also interpretable, so it is chosen here first how to use this kind of information. In particular, for sake of illustration, the median is chosen as the input RV for assessing the estimation, with no loss of generality, since any quantile Xp is proportional to the median through the following relations, which are simply deducible (see Appendix B): Therefore, an assessment of the prior information in terms of any quantile is essentially the same as an assessment of the prior information in terms of the median, which in turn implies a prior information in terms of α, since The above relations also imply prior information in terms of IRD CCDF, expressible in terms of the median m as well as in terms of the parameter α: It is remarked that the function S(z) reported here, being relevant to a fixed threshold value, z, of EWS has to be regarded as a constant, i.e., a parameter, in view of the estimation process. In the above Equation (34), instead, the CDF F(x) reported there denotes a function of the generic EWS value as a variable argument.
Far from what is above discussed, the peculiar method suggested in this paper is that S, a RV with support in (0,1), is assumed as a beta RV, which implies that the parameter α has an NEBD PDF, as can be deduced from Appendix D, since In practice, assuming an NEBD prior PDF for α implies that the negative exponential function of α in (38) has a beta distribution. This assumption has the advantages already discussed in the previous section and may be put in practice, starting from prior information in terms of the median as follows: Let a prior guess on the median be given in terms of a "credible interval" (m1, m2), i.e., a Bayesian (prior) confidence interval on the median, i.e., an "interval" (m1, m2) such that, for a given probability value p, 1) Since α =  (m) 2 , the corresponding quantiles (α, α on α are readily be evaluated as follows: 2) On the basis of such quantiles of α, the corresponding quantiles of S are readily computed by (38) with S being a one-to-one increasing function of α, as follows: This information is sufficient (see [27]) for assessing the beta PDF on S. In any case, the posterior PDF on S must be evaluated numerically. In principle, this also holds for the posterior PDF on R, which is the basic tool for assessing the BE of R, i.e., its posterior mean. Some simplifications are however possible, which will be shown to be numerically acceptable. It was already discussed that, in cases in which S has a beta PDF, V has a computable PDF, the BGP PDF. Apart from such particular cases, the independence of W and S can be used as a basis for evaluating the posterior mean of V: where p and 1 are the posterior gamma parameters computed in the previous subsection. This shows that only the posterior mean of S is needed in order to evaluate the posterior mean, i.e., the BE, of V°. This does not imply a direct estimation of the RI, which theoretically requires the integration of (30) or (32).
However, at least as a first-order approximation, with R = e −V , the following approximation of the BE of R can be evaluated: Finally, if V can be described by a gamma distribution, then R has a "negative log-gamma" distribution (summarized in Appendix D and recalled in the previous subsection) and its posterior mean is evaluable analytically.

A Framework for the Numerical Applications of the Bayes Method of Inference for the Risk Index with Reference to Practical Applications
Before illustrating (in the next sections) extensive numerical applications of the Bayes method aimed at evaluating its performances, a framework for such numerical applications is established here in order to better highlight the rationale behind the computations which will be shown afterwards. For the sake of brevity, only the case of a gamma prior PDF on the IRD parameter a will be illustrated here, leaving the case of a of a beta prior PDF on the IRD CCDF S(α) to be illustrated directly in the course of the numerical applications. Section 5.1 here is devoted to the prior assessment analysis, and the following Section 5.2 is devoted to the posterior analysis, still in the case of a gamma prior PDF on the IRD parameter α. the assumption of a gamma prior PDF on the Poisson parameter  is always maintained instead. It is remarked that the two subsections both referred to practical applications, being based upon typical wind-speed data as those reported in the Appendix B.

The Case of a Gamma Prior PDF on the IRD Parameter α: Prior Analysis
First, for sake of numerical illustration in Table 1 are shown with reference to a fixed time interval of 1 year, a threshold value z = 22 m/s., and nine couples of values (, α), shown respectively in the headings of the rows and columns of Table 1. Also, the following numerical applications will refer to these typical values of the RI and parameters  and α. As already discussed, the values of the RI increase as both  and α increase, so they increase from left to right, looking through rows, and from top to bottom, looking through columns. In particular, for the sake of brevity, the three cases corresponding to boldface characters in the principal diagonal of Table 1 will be those taken as references in the course of the following numerical applications. Considering that the parameters  and α are RVs in the Bayes method, their values are to be considered expected values of their prior distributions. Therefore, the following three cases are taken as references:  Table 1, are merely indicative; indeed they are approximate expected values since the RI is a nonlinear function of (, α): Indeed, for a given prior PDF p(, α) of (, α) which-as long as t = 1 year and W =  t-coincides with the prior PDF p(w,α) of (w,α), the mean of R, as already discussed before since However, hints on a simpler way to obtain such an expected value of R have already been hinted at before and will be shown here in the course of computations, together with numerical or simulation methods, which are always possible to give a better approximation if the above ones are generally satisfactory, but it is remarked that the primary meaning of the above approximations (such as E[R] ≈ 0.0735, and more) is not to serve as approximate expected values of the true means but only as typical values that one should expect from computations, so that the experimenter can choose the typical values of the parameters to be assigned a prior PDF with a better knowledge of its implications.
As for the choice of the prior PDF p(, α), i.e., the gamma PDF (on both  and α) and/or beta PDF on S(, α), their parameters are chosen in order to obtain the expected values of the cases A, B, and C above; moreover, in each of the three cases, different values of the prior coefficient of variation (CV)-i.e., the ratio between standard deviation and mean values of  and a-are chosen in order to show different degrees of prior uncertainty (a larger prior CV value implies a higher degree of prior uncertainty about the random parameter values). Therefore, denoting for a given RV Q, its CV is as follows in terms of its mean value E[Q] and standard deviation STD[Q]: The following values of CV for both  and α have been chosen: 0.05 and 0.10, giving rise to four different subcases of the three cases A, B, and C above, each labeled with (i,j) indexes, "i" standing for 1 or 2 according to CV[] = 0.05 or CV[] = 0.10 and "j" having the same meaning for CV[α].
In detail, 12 = 3 × 4 sets of combinations of prior values (corresponding to different simulation cases) among the many more performed are shown in the following sections with respect to the mean value and CV of the input RV(, α), as listed with more details in the next section, while here, they are only exemplified in Table 2: Details on the results of numerous simulations performed for each one of the above cases will be shown in the following sections, while here, some highlights are reported which give some insight on the process of the prior and posterior PDF assessments and on the BE of the RI. For the sake of exemplification, case C11 is examined. Therefore, it is assumed that, concerning the prior PDF assessment, The corresponding prior gamma PDF with such parameters (0 and 0), i.e., with mean E[α] = 90 and CV: CV[α] = 0.05, is depicted in Figure 3. Due to the high value of the shape parameter n0 = 400, such a PDF closely resembles a Gaussian one [28].  Of course, the primary interest here is not the estimation of the above parameters per se, but they are the input RVs which allow the BE of the RI. Let us first evaluate the prior mean of RI implied by the above assessments on  and α.
By simulation, numerical integration, or following the approximated method illustrated in the previous section, it is not difficult to evaluate such a prior mean and afterwards also the posterior mean based upon the following relations: First, notice that the RV has a prior PDF which can be easily assessed in terms of the NLGD PDF (Appendix C). Indeed, since the RV α has a gamma PDF with parameters 0 = 400 and 0 = 0.2250, it is easy to verify that also the RV, has a prior gamma PDF, with parameters  and  expressed by   0 = 400 and   0/z 2 = 0.2250/z 2 = 4.6488e-04. Therefore, by definition of the NLGD model (Appendix C), U = 1 − S = e −β has an NLGD prior PDF, here denoted by p(u): This implies that S = 1 − U has a prior PDF which can be easily computed from the above NGL PDF, evaluated in u = 1 − s, as follows: In Figure 4, two prior gamma PDFs are superimposed: 1) The one corresponding to the theoretical approximating prior gamma PDF of V, and 2) The one estimated by means of the "smoothed kernel PDF estimator" [31], obtained after generating by Monte Carlo simulation (The Monte Carlo simulation of samples from RV W and S was obtained through the function "gamrnd" of Matlab® used for generating samples of the gamma RV (, a), and then, the corresponding sample of S was obtained through the function S = S(α) = 1 − exp(−α /z 2 )) 5000 values of W and S and so of V = WS, according to their above prior PDF.
The prior mean of V, as estimated from the simulated sample data, is equal to 1.6955, while the one obtained from the approximating prior gamma is equal to 1.6965, which confirms the adequacy of the approximation.
Finally (as long as the priori assessment is concerned) the prior mean of R can be evaluated by the simulated sample, i.e., by the sample mean (average value), 0.8155.
If, instead, the above approximating prior gamma is adopted for V, then 1 − R = e −V has an NLGD prior PDF. Indeed, since the prior gamma PDF of V, as already computed, has parameters 0V = 218.6171 and 0V = 0.0078, then analogously to what was before illustrated, S − R has a prior PDF which can be approximated by a NGL PDF, as follows: Also, in this case, the theoretical approximating prior NLGD PDF of R and the smoothed kernel PDF estimator of R are perfectly superimposed. Assuming the approximating NLGD prior PDF of R, its mean is equal to with   0V and   0V, which yields 0.8155, the same value computed by simulation.
Also, as hinted at in a previous section, it can be roughly computed by means of the "first order" approximation: E Also, this value is close to the previous one and very close to what reported is in Table 1.
It is well known, however, that the Bayes inference method, differently from the classical ones, produces not only point estimates of the parameters but also the probability distributions describing our uncertainty about their values. In this case, we can assess a whole prior PDF, which, for what concerns its numerical estimation, is shown in Figure 5 in the form of the smoothed kernel PDF estimator of R. Such a numerical prior PDF, i.e., the quantile estimates associated with the simulated sample, yields, for instance, the following prior 90% credibility interval: (0.7804 0.8483). P (0.7804 < R < 0.8483) = 0.90 Using the NLGD approximation, taking into account that U = (1−R) is a NLGD RV with parameters (  0V;   0V), it is easy to see that the p-quantile of R is equal to 1, the (1−p)-quantile of U: where nlginv (q, , and ) is the q-quantile of an NLGD RV with parameters  and . (For such a computation, the Matlab® function computing the quantile of a gamma RV, gaminv, was used, since Y = e (−X) is an NLGD RV when X is a gamma RV, so the p-quantile of Y is given by Yp = e (−X(1−p)) ) This yields the following prior 90% credibility interval (0.7796, 0.8489), which is practically the same interval above found numerically (0.7804, 0.8483), with precision up to the third decimal number. Of course, since the numerical computations are very reliable, the above gamma or NLGD approximation, as those dealt with in the following, are useful only because they allow easier and faster computations, a sensitivity or robustness analysis with respect to the assumed values of the parameters or the choice of the basic distributions (e.g., one could be interested in evaluating which are the changes in the range of results if the prior gamma PDF, here assumed for the parameters W and S, is changed into, say, a lognormal or a uniform PDF. This is an issue discussed under the topic of "robustness analysis", briefly dealt with in the final part of the paper). Apart from this remark, in the following, more emphasis will be given to the numerical computations.

The Case of a Gamma Prior PDF on the IRD Parameter α-Posterior Analysis
Obviously, the basic purpose of Bayes inference is to obtain, once the sample data are observed, the posterior estimates and, also, the probability distributions describing our uncertainty about their values. This can be performed in a relatively straightforward way, which can take advantage of the numerical tools shown in the course of the prior analysis illustrated above. Here, a simple example of how the posterior analysis can be accomplished is illustrated, based upon a small data set produced by simulation from an IRD model, according to the relation (see Appendix F): which yields an IRD RV if U is a uniform RV in (0,1). Of course, the posterior analysis requires first that the relevant parameters are generated according their prior distributions; then, that their PDFs are updated in view of the above EWS sample; and then, that the BE of the RI is obtained. This is briefly illustrated here, while the performances of such a BE and their comparison with the one obtained by classical methods are examined in the following sections.
Reference is still made to the above case C11, so that In order to illustrate a simple and easily reproducible numerical example and to highlight how the Bayes inference allows a significant update of information also with very few data, a small sample of only 3 EWS values is produced by a single simulation exemplified here as follows: --a value of α is generated by its gamma prior PDF; its value is 88.00. --a value of  is generated by its gamma prior PDF; its value is 10.00.
The corresponding "true" value of the RI to be estimated is as follows: A sample of three EWS values obtained from three uniform variates according to Equation (59), and the above generated value of α is obtained through simulation, resulting in the following sample: x = (12.07, 45.01, 49.62).
To such a sample, the following observed number of exceedances M is M = 2 (the number of EWS values exceeding the threshold 22), and the Statistics n (sample size) and  of Equation (15) and (A31) in Appendix F for evaluating the MLE and the BE for n equals to 3 is 0.0078 and for M equals to 2 is which is of course (given also the very little sample size) much closer to the true value (0.8103) of R than the MLE (0.6672). It is also worth mentioning that, as in the previous subsection, the "rough" computation of such BE by means of the "first-order" approximation E[R|D] ≈ 1 − e −E[V|D] = 1 − e −1.6955 = 0.8165.
Finally, the numerical computations, both when performed by simulation or numerical computations, give essentially the same result: E[R|D] ≈ 1-e −E[V|D] = 1-e −1.6955 = 0.8120. Of course, the merits of the BE are not appreciable by a single simulation, here shown only to highlight the mathematics, but require an overall assessment, which is the object of the next sections.

Numerical Application for Performances of the Bayes Approach in the Case of a Gamma Prior PDF on the IRD Parameter α
This and the next section follow the framework illustrated in the previous section, in which only the features of an estimation procedure for a given sample have been reported. Since in practice the Bayes estimate and their relevant features (e.g., their mean square errors) are computable only numerically, in order to assess the overall performance of the proposed Bayes estimation, a large series of Monte Carlo simulations have been carried out, relevant to different sample sizes and extending the same approach. In detail, the Bayesian methodology has been evaluated numerically, as illustrated here, by means of a large set of simulated experiments aiming at assessing the estimation efficiency of the Bayes approach. This is illustrated here in the case of a gamma prior PDF on the IRD parameter α, while the alternative assumption of a beta prior PDF on the IRD parameter S(α) will be briefly illustrated in the next section. Both in this and in the next section, the assumption of a gamma prior PDF on the Poisson parameter (α, ) is maintained, as customary [26,27].
After a sound investigation, different sets of typical values have been selected for the mean and the standard deviation (or the CV) of the main random parameters (α and ).
Different CV prior values were selected with the same values of prior means, since-as intuitive-the estimation performances are rather sensitive to variations in CV prior values while variations in the prior means are less effective in terms of estimation efficiency. As above said, 12 sets of simulation cases are illustrated here, among the many more performed, and shown in Tables 3-5. For each of the above cases, N random samples of size n (with n varying as specified below) of EWS values Xj have been generated according to the assumed IRD PDF of stress X with parameter α and with a Poisson parameter  for the process of exceedances. The parameters ( and α) are in turn RVs generated according to their prior gamma distributions, as shown in a previous section (different priors were also investigated, with similar results to those shown in the following).
For each sample size n, a number N = 10 4 of replications has been performed; let us refer to each considered couple of values of n and N as a single simulation case study; different simulation case studies were considered by testing different sample sizes. Here, in particular, the results for small (n = 3, n = 10) or moderate (n = 30, n = 50) sample sizes are reported in the tables that follow, among the many more performed. The choice to show also the results for very small sample sizes as n = 3 takes its origin in view of the previously discussed features of Bayesian information updated also in the case of a limited amount of data has mainly a theoretical interest; nonetheless, it can also have a practical aspect when dealing with very high values of WS, so that their occurrence is very rare and only very few data are available.
With reference to the generic parameter  to be estimated (which in this case is basically the RI R, but one could be interested also in the BE of and α), it will be denoted as j°, the estimate of the "true" value j of  relevant to the particular n-sized sample generated at the jth simulation cycle. The basic statistics-estimated at the end of each simulation case study-which describe the efficiency of the proposed estimates are the following: The REFF index is the most widely adopted measure of the "relative efficiency" [26,27] of the Bayes estimator with respect to the ML estimator: the more it exceeds unity, the more efficient is the Bayes estimate with respect to the ML estimate.
The above-defined quantities are based on the concept of "RMSE". Given an estimator, °, of the parameter , its RMSE is-as well known-defined as with both ° and  as random variables in the Bayesian approach, while the latter is an unknown, deterministic constant in the classical statistical approach. The RMSE (which is not analytically computable) is evaluated here at the end of each simulation case study by means of the ordinary large-sample estimator [31] as follows: While the above measures concern the efficiency of the Bayes estimates, other significant "performance indexes" were also evaluated, as measure of their precision: With the same notations, the Maximum Relative Error, MRE, of an estimator over N trials is given by the following: The values of RMSE, ARE, and MRE can be calculated for each simulation case study for both the ML estimator and the proposed Bayesian approach, thus obtaining the values of quantities RMSEB, RMSEL, AREB, ARELB, MREB, and MREL (defined above). In other words, RMSEB denotes the RMSE value for a BE, RMSEL denotes that for an MLE; AREB denotes the ARE value for a BE; AREL denotes that for an MLE; and so on. Note that, for an easier comparison, both ARE and MRE are evaluated in terms of absolute values.
In addition, the following (dimensionless) ratios are computed as useful synthetic measures in evaluating the precision of the Bayes estimates with respect to the MLE: As for REFF, RARE and RMREL are also greater than 1 if the Bayes estimate is "better"-in the specified sense-than the ML one.
A summary of the simulation results, in Tables 3-5, reported the average value of R, α and f; the average value of the MLE and BE of R; and the measures RMSEL, RMSEB, RARE, RMRE, and REFF.
Any average value is always computed as for indexes ARE and MRE above as the sample mean computed over the N simulations performed for each generated sample. The average values of a random parameter, such as α, should be ideally equal to its mean value (they would be equal if an infinite number of simulations would be possible), e.g., 15 in case A11; in practice, they should be close to the theoretical mean values.
All the above values are reported as a function of the sample size considered in the simulations (see row headings). It is noticed that the values of the RI average value, AV(R), and the ratio REFF (relative efficiency) are boldface, being the most significant measures.
Since they are evaluated through simulation, all the reported quantities must be regarded as sample estimates. The number of digits used for the data representation is purely indicative.
All the results show that the Bayes estimates errors are always acceptably small and much smaller-in particular for very small sample sizes-than those of the traditional ML method. In all cases (as in the many more cases omitted here), the reported results show the good performances of the Bayes estimates, assessing the Bayes estimation robustness, with very high REFF values, especially for small sample sizes. As could be anticipated, cases implying a higher confidence in the prior values (i.e., with minor CV values) are more favorable in terms of efficiency of the Bayes estimates: it is remarkable, however, that also cases corresponding to a relatively high degree of uncertainty in the prior estimates (CV = 10%) are characterized by high values of efficiency. The same conclusions can be drawn about the precision of the Bayes estimators, as measured by the RARE and RMRE index.
In summary, the relative efficiency with respect to the ML estimate-as measured by REFF-is always larger than 1; in particular, it is much larger for small sample sizes, for which the ML estimates are definitively outperformed by the Bayes ones. For example, it can be noticed in Table 5, case C11, that the REFF assumes a value around 7 for n = 50 but increases until 33 about for n = 3. This latter feature is reinforced by the examination of ratios RARE and RMRE, which behave, as could be expected, less regularly from sample to sample.
It is again remarked that the above Bayes estimate must be computed by numerical methods, as above discussed. In particular, both numerical integration of the posterior PDF and simulation were performed, giving essentially the same results.
Finally, it is remarked that the numerical simulations are not (only) theoretical but are based upon typical values of wind-speed parameters which have been found in actual practice. Moreover, a larger series of Monte Carlo simulations carried out for a much wider range of parameter values (prior means and CV) were performed, all with similar results, which are often highlighted in the applications of Bayes estimation [32][33][34][35][36].

Numerical Application for Performances of the Bayes Approach in the Case of a Beta Prior PDF on the IR Parameter S(α)
In this section, the topic of the previous section is continued, assuming a beta prior PDF on the IR parameter S(α) while maintaining the customary assumption of a gamma prior PDF for the random parameter . Also, in this case, interest is focused on the evaluation of the performances of the Bayes approach by numerical methods since, as highlighted in Section 4, also in this case, generally an analytical solution is not available or is very cumbersome [37,38].
In practice, as illustrated in Section 4, a prior guess on the parameter S (α can, as in the case of a prior guess on α, be originated by a prior guess on the median or may be directly expressed on the basis of "credible interval" (S1, S2), i.e., a Bayesian (prior) confidence interval on S).
As already discussed (see the discussion after eqn. (38)), , it is convenient that S, a RV with support in (0,1), is assumed as a beta RV (already introduced in Section 4 and with more details in Appendix D), since this is a very flexible distribution in the interval (0,1), which implies that the parameter α has an NEBD PDF (Appendix D), since: For brevity, here, only three cases are examined. In detail, we start from the same cases examined in previous section, here recalled only in the "most unfavorable" case of both CV = 0. For easier comparison with the cases examined above, we assume that the beta(r,s) prior PDF for S has the mean value corresponding to the above mean values of α, with the value CV = 0.10 assigned to S, so to obtain the following three cases, denoted with apexes A'22, B'22, and A'22 corresponding respectively to A22, B22, and C22. For example, the curve of the prior beta(r,s) PDF of case C'22 (S has mean value 0.1697 and CV = 0.10, implying: r = 82.86, s = 405.4) is depicted in Figure 6. The results of these 3 cases are reported in Tables 6-8, illustrated as in the previous section, except that, of course, in this section in the 2nd row of the tables, the average value of S, AV(S) is shown instead of that of α, AV(α). It is apparent that the results are quite comparable and indeed very similar to those of the previous section, confirming the validity of the approach with respect to different methods of assessing the researcher's prior information.
Finally, it is remarked that an analysis of the robustness of the Bayesian Estimation with respect to the assumed prior hypotheses has been successfully carried out following the same methodology illustrated in [27]. In particular, many different prior PDFs were chosen instead of the assumed ones, e.g., lognormal or uniform prior PDF with the same mean value and CV of the assumed gamma and beta PDF. For reasons of length, these results are not reported here but are available on request from the authors.

Conclusions
The correct assessment of the probabilistic distribution for wind and extreme wind speeds is very important for a technical and economical point of view and gives the possibility to optimize the electricity production of a wind farm. In particular, the occurrence of very high values of wind speed is an undesirable event, since wind power systems provide zero power output for values greater than their cutoff thresholds. Also, the mechanical safety of the installations can be seriously compromised by extreme values of wind speed. The IRD is proposed in the paper as a valid model for extreme wind speeds; such a model is proposed for the first time in the literature on the object, and its validity is illustrated with reference to real wind-speed data. Then, this paper investigates the problem of a proper and efficient Bayes estimation of a risk index described in terms of a stochastic process of extreme values of wind speed based upon a practical parameter estimation of the IRD using the gamma and/or the prior beta distribution according to which kind of prior information is available. The methodology proposed is based upon the use of real data of EWS and using different estimation procedures, also comparing such estimation procedures to the ML procedure. The proposed model is verified by means of many numerical simulations. The obtained results show the efficiency and precision of the analyzed model. The considered Bayes estimator appears also to be robust with respect to incorrect prior assumptions.
It is highlighted that, beyond the theoretical aspects which have been developed in the paper, its basic purpose is to address a real engineering problem, since it deals both with wind power production evaluation and the mechanical safety of the installations. This is a clearly basic topic to be addressed both in the planning and the operation of a wind farm and is afforded in the manuscript by a novel risk index which the authors deem to be an important tool in helping a decision making procedure for deciding if a wind farm is worth being installed and whether its operation will be safe or reliable enough in time.
Author Contributions: conceptualization, E.C. and L.P.D.N.; methodology E.C.; software, E.C.; validation, E.C. and L.P.D.N.; writing-original draft preparation, E.C. and L.P.D.N. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Acknowledgments:
The authors express their most sincere thanks to the reviewers for their insightful and very helpful comments, which considerably improved the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature
while the model variance does not exist. The application to EWS data is exemplified in this study by data analysis related to EWS data in [19], summarized in Table A1, where the 52 weekly maxima (z1,… zn, with n = 52) of wind speed data in (m/s) are reported (These data represent the 52 maximum WS values registered in each of the 52 weeks of year 2014 in Boulder, CO, USA at latitude 39°54′ north and longitude 105°14′ west; they are part of large datasets reported in [19]). Some basic statistical parameters of the sample in Table A1 are summarized in Table A2. Table A2. Some basic statistical parameters (quartiles, mean, and SD) of the EWS data in Table 1. 25 The goodness of fit tests for assessing the validity if an IRD model for the above data were all successfully satisfied. A simple way to get evidence of such a fitting is drawing a "Weibull probability plot" with shape parameter  = 2 for the data on the RV Y = 1/Z as shown in Figure A1 (The Weibull probability plot in Figure A1 was obtained by the Matlab® function "wblplot (Z)", with Z being the data sample vector of Table A1). For the (classical) estimation of parameter α with the purpose of finding the best IRD fitting PDF, it is readily seen that the quantile estimation methods [24], the method of moments [31], and the maximum likelihood one yield essentially the same result. For brevity, here, the first one is reported. Since the median is equal to 33.9 m/s, the parameter α can be estimated as follows: a * = 33.9 [−log(0.5)] = 28.22

x0.25 (m/s) x0.50 (m/s) x0.75 (m/s) Mean Value (m/s) Standard Deviation (m/s)
The graph of the corresponding IRD fitting PDF is shown in Figure A2 superimposed to the histogram. Figure A2. IRD fitting PDF for the data histogram of EWS data in Table A1.