Arisman Adnan*
| Rado Yendra | Ibrahim Sulaiman Hanaish | Muhammad Marizal | Gustriza Erda
© 2026 The authors. This article is published by IIETA and is licensed under the CC BY 4.0 license (http://creativecommons.org/licenses/by/4.0/).
OPEN ACCESS
Reliable estimates of extreme wind speed are needed for hazard-aware infrastructure planning, yet spatial information on wind extremes remains limited across Sumatra, Indonesia. This study combines Bayesian extreme-value modelling with ordinary kriging to estimate and map wind-speed return levels from 21 meteorological stations for 1999–2019. Daily maximum wind-speed records were reduced to annual block maxima, and two candidate distributions—the generalized extreme value (GEV) and generalized logistic (GLO) distributions—were fitted separately at each station using Markov Chain Monte Carlo (MCMC) estimation with weakly informative priors. Goodness of fit was assessed by the Root Mean Square Error (RMSE) between empirical and fitted cumulative distribution functions. The GLO distribution produced the lower error at 12 stations, whereas the GEV distribution performed better at nine stations, suggesting clear spatial heterogeneity in the tail behaviour of wind extremes. Posterior return levels for 10-, 50-, and 100-year return periods were interpolated by ordinary kriging to produce regional maps. The resulting patterns indicate relatively high extreme-wind return levels in northern and southern Sumatra, with several stations exceeding 30 knots for the 100-year return period, while lower intensities are generally found in the western part of the island. These maps provide a risk-oriented spatial basis for preliminary infrastructure assessment and hazard screening. They should not, however, be interpreted as direct indicators of wind-energy feasibility, which requires additional evidence on mean wind speed, wind-power density, persistence, hub-height conditions, terrain, turbine suitability, and grid access.
annual block maxima, Bayesian extreme-value modelling, extreme wind speed, generalized extreme value distribution, generalized logistic distribution, ordinary kriging, return level
Indonesia’s energy system remains dominated by fossil-based power generation, while renewable sources contributed approximately 13.1% of national energy production in 2023 [1]. This condition highlights the need to expand the scientific basis for assessing renewable-energy resources across different regions. Among the available alternatives, wind energy is particularly relevant because Indonesia’s archipelagic geography creates substantial spatial variation in atmospheric circulation and wind conditions [2, 3]. However, the development of wind-related infrastructure requires an understanding not only of average wind conditions but also of rare and potentially damaging wind events [4, 5].
Extreme wind speed differs conceptually from the sustained wind conditions used to evaluate wind-energy feasibility [6, 7]. Wind-energy potential is primarily determined by mean wind speed, wind power density, persistence, hub-height conditions, terrain, turbine suitability, and grid accessibility [8]. In contrast, extreme wind analysis focuses on the magnitude and probability of rare events that may affect structural safety, turbine loading, infrastructure resilience, and operational continuity [9]. Therefore, spatial assessments of extreme wind speed should be interpreted primarily as risk-oriented information rather than as direct evidence of wind-farm suitability.
This distinction is particularly important in Sumatra, where the spatial behaviour of wind may be influenced by complex coastlines, mountain ranges, proximity to the Indian Ocean and the Malacca Strait, and differences between coastal and inland environments. These geographical characteristics may generate substantial variation in both the magnitude and frequency of extreme wind events [7]. As a result, estimates derived from one location cannot be assumed to represent neighbouring or geographically distinct areas. A regional assessment therefore requires methods capable of capturing both station-specific extreme behaviour and spatial dependence across observation locations.
The challenge becomes more substantial when extreme-wind assessment relies on relatively short observational records. Annual block-maxima analysis retains only one maximum observation for each year, thereby reducing the effective sample size available for statistical inference. Estimation of long return periods, such as 50 and 100 years, consequently involves extrapolation beyond the observed record and may be associated with considerable uncertainty. A suitable analytical framework should therefore provide more than point estimates. It should also represent parameter uncertainty and support cautious interpretation of long-period return levels [10].
Extreme Value Theory provides a statistical framework for estimating the probability and magnitude of rare events beyond the central range of observed data [11, 12]. Under the block-maxima approach, annual maximum wind speeds can be modelled using distributions such as the generalized extreme value (GEV) distribution [13]. The generalized logistic (GLO) distribution offers an alternative with different tail characteristics and may better represent certain meteorological series. Because the suitability of a probability distribution may vary among locations, applying one distribution uniformly across all stations can obscure local differences in extreme-wind behaviour [14, 15]. A station-specific comparison of competing distributions is therefore necessary before estimating return levels.
Recent methodological developments have increasingly incorporated Bayesian inference into extreme-value analysis. Bayesian estimation using Markov Chain Monte Carlo (MCMC) produces posterior distributions for model parameters and return levels, allowing uncertainty to be represented explicitly rather than relying solely on point estimates [16, 17]. Nevertheless, station-level extreme-value estimates remain spatially discrete and provide no direct information for unmonitored areas [18]. Geostatistical interpolation, particularly kriging, can address this limitation by using the spatial dependence among observed locations to estimate values continuously across a region [12].
Despite these developments, an important gap remains between probabilistic extreme-value estimation and regional spatial assessment. Extreme Value Theory studies commonly estimate return levels at individual stations without extending the results into continuous spatial representations [19]. Conversely, wind-speed interpolation studies often map observed or deterministic values without adequately representing uncertainty arising from probabilistic extreme-value models [20, 21]. This separation limits the ability of existing approaches to describe both location-specific tail behaviour and regional patterns of extreme-wind exposure, particularly in geographically complex regions such as Sumatra.
This study addresses this gap through an applied framework that combines station-specific Bayesian extreme-value modelling with ordinary kriging. Annual maxima of daily maximum wind speed from 21 meteorological stations are fitted using the GEV and GLO distributions. The most suitable distribution for each station is selected through numerical and graphical goodness-of-fit assessment. Return levels for 10-, 50-, and 100-year periods are then estimated from the fitted posterior distributions and interpolated to produce continuous spatial maps of extreme wind speed across Sumatra.
The contribution of this study is primarily empirical and applied rather than methodological. Bayesian extreme-value modelling and ordinary kriging are employed in sequential analytical stages. First, the GEV and GLO distributions are compared at each station to identify the more appropriate model for characterizing extreme wind behaviour. Second, Bayesian estimation is used to quantify uncertainty in the model parameters and return levels. Third, the station-specific return-level estimates are interpolated using ordinary kriging to produce continuous regional maps of extreme wind speed.
Accordingly, the main objective of this study is to characterize and map the spatial distribution of extreme wind speed in Sumatra for wind-risk assessment and risk-informed infrastructure planning. The methodological procedures support this applied objective by linking station-level probabilistic estimates with regional spatial representation. The resulting maps may provide supplementary information for early-stage wind-energy screening, but they are not intended to directly determine wind-energy feasibility and should not replace comprehensive wind-resource or engineering assessments.
The wind speed data used in this study were obtained from the Indonesian Agency for Meteorology, Climatology, and Geophysics (BMKG). The dataset analysed consisted of daily maximum wind speed observations from 21 meteorological stations distributed across Sumatra and its surrounding islands during the period 1999–2019. The selected stations represent different geographical settings, including coastal, inland, northern, central, southern, and island locations. Their names, geographical coordinates, and station codes are presented in Table 1, while their spatial distribution is shown in Figure 1.
Table 1. Geographical coordinates and codes of the 21 wind speed stations in Sumatra
|
Name of Station |
Latitude |
Longitude |
Code |
|
Sabang |
5.8718192 |
95.2695115 |
S1 |
|
Cut Nyak Dien |
4.0451874 |
96.1819568 |
S2 |
|
Malikus Saleh |
5.2300426 |
96.8800508 |
S3 |
|
Sultan Iskandar Muda |
5.5150182 |
95.4244056 |
S4 |
|
Fltobing |
1.556434 |
98.8878537 |
S5 |
|
Polonia Medan |
3.6357715 |
98.8093176 |
S6 |
|
Belawan |
1.3984236 |
99.3606478 |
S7 |
|
Gunung Sitoli |
2.2588861 |
98.9938384 |
S8 |
|
Dabo |
−0.477783 |
104.5778296 |
S9 |
|
Hang Nadim |
1.1219253 |
104.1182972 |
S10 |
|
Kijang |
0.9209616 |
104.4616936 |
S11 |
|
Depati Parbo |
−2.087134 |
101.3911004 |
S12 |
|
Sultan Thaha |
−1.631591 |
103.5702051 |
S13 |
|
Anambas Tarempa |
3.2098493 |
106.1962861 |
S14 |
|
Fatmawati |
−3.858996 |
102.2669824 |
S15 |
|
H.A.S. Hanandjoeddin |
−2.753974 |
107.6828968 |
S16 |
|
Pangkal Pinang |
−2.161710 |
106.1382584 |
S17 |
|
Sultan Syarif Kasim |
0.4647941 |
101.4476651 |
S18 |
|
Raden Intan |
−5.24223 |
105.1776319 |
S19 |
|
Sultan Mahmud Badaruddin |
−2.894923 |
104.7052902 |
S20 |
|
Tabing |
−0.786063 |
100.2155705 |
S21 |
Figure 1. Geographical distribution of the 21 wind speed stations included in the analysis
Table 1 shows that the 21 meteorological stations are distributed across the main island of Sumatra and several surrounding islands, covering a wide geographical range from the northern to the southern parts of the study area. The station network includes coastal, inland, and island locations, thereby representing different geographical conditions that may influence extreme wind behaviour. However, the stations are not evenly distributed, with relatively denser coverage in the northern and eastern regions and wider gaps in several central and southern areas. This spatial configuration should be considered when interpreting the kriging results because predictions in areas located farther from observation stations may be associated with greater interpolation uncertainty. The station points in Figure 1 indicate the locations of the 21 observation sites coded S1–S21, enabling direct linkage between their geographic positions, station-level return estimates, and kriging results.
The names, codes, and geographical locations of the 21 meteorological stations are presented in Table 1 and Figure 1. The stations represent coastal, inland, and island environments across northern, central, and southern Sumatra. Their spatial distribution is uneven, with relatively denser coverage in the northern and eastern regions and wider gaps in several central and southern areas. This configuration should be considered when interpreting the kriging results because prediction uncertainty may be greater in areas located farther from observation stations.
3.1 Block maxima, extremes distribution, and return periods
The block-maxima approach was used to extract annual maximum wind speeds for extreme-value modelling. In meteorological and hydrological applications, blocks may be defined on annual, monthly, or seasonal bases, with only the maximum observation from each block retained for analysis [22]. In this study, annual blocks were used so that the resulting extreme-value series contained one annual maximum wind speed observation for each year of the study period.
Two probability distributions were considered for modelling the annual maximum wind speed series: the GEV distribution and the GLO distribution [23]. The parameters of both distributions were estimated using Bayesian MCMC. Bayesian MCMC was selected because it provides posterior distributions for the model parameters and return levels, thereby allowing estimation uncertainty to be represented explicitly [24]. Weakly informative priors were adopted, and 30,000 MCMC iterations were generated, with the first 10,000 iterations discarded as burn-in. Convergence was assessed using trace plots and the stability of the posterior density plots.
A fundamental relationship exists between return period (T) and probability of occurrence (p). These two variables are inversely related to each other by the following equation:
$\frac{1}{T}=p$ or $T=\frac{1}{p}$ (1)
For a T-year return period, the annual exceedance probability is 1/T. Hence, a 100-year return level corresponds to an annual exceedance probability of 0.01 under the fitted stationary model. This interpretation does not imply that such an event occurs exactly once every 100 years. Therefore, the values of extreme events for the return periods 10, 50, and 100 years are calculated by substituting the vectors of α, ε, and κ into the quartile function, for 0 < F < 1. This procedure was carried out for F = 0.9, 0.98, 0.99, to obtain the extreme events for the return periods of 10, 50, and 100 years. Posterior summaries were reported using posterior means and 95% credible intervals for all model parameters and return levels.
3.2 Bayesian Markov Chain Monte Carlo and goodness-of-fit tests
This section introduces the idea of Bayesian MCMC using weakly informative priors [25]. The formulation of the prior beliefs by θ is used to formulate an expression by a probability density function $\pi(\theta)$ with no reference to the data. The Possibility for θ is $L(x \mid \theta)$, therefore, information and the possibility are combined using Bayes' theory to produce a posterior distribution for θ as follows [26]:
$\pi(\theta \mid x)=\frac{\pi(\theta) L(x \mid \theta)}{\int \pi(\theta) L(x \mid \theta) d \theta}$ (2)
Extreme-value analysis is used to characterize the upper-tail behaviour of an observed process and estimate the probability of future rare events. Within the Bayesian framework, prediction for a future observation y is obtained from the posterior predictive distribution, which integrates the conditional distribution of the future observation over the posterior distribution of the model parameters [27]. Therefore, it is calculated as follows:
$f(y \mid x)=\int f(y \mid \theta) \pi(\theta \mid x) d \theta$ (3)
with the predictive distribution of y given x. The integral computing process in the predictive distribution makes it possible for the simulation techniques to overcome the posterior distribution. The Bayesian MCMC prior density is calculated as follows:
$\pi(\phi, \xi, \kappa)=\pi_\alpha(\phi) \pi_{\xi}(\xi) \pi_k(\kappa)$ (4)
For both the GEV and GLO distributions, the parameters are defined consistently as follows: $\xi$ denotes the location parameter, α > 0 denotes the scale parameter, and $\kappa$ denotes the shape parameter. The location parameter determines the central position of the distribution, the scale parameter controls its dispersion, and the shape parameter describes the behaviour of the upper tail. To ensure that the scale parameter remains positive during MCMC estimation, the transformation $\phi=\log (\alpha)$ is used, so that $\alpha=\exp (\phi)$.
Assuming prior independence, the prior distributions are specified as $\phi \sim N(0,100)$, $\xi \sim N(0,1000)$, and $\kappa \sim N(0,10)$. Thus, $\phi$ represents the transformed scale parameter, $\xi$ represents the location parameter, and $\kappa$ represents the shape parameter. The same symbols and parameter definitions are used consistently for both the GEV and GLO models throughout the analysis. These priors are described as diffuse or weakly informative priors because they are proper normal distributions with relatively large variances.
The variances are chosen large enough to make the distribution almost flat, in accordance with prior ignorance and the posterior density as follows:
$\pi(\phi, \xi, \kappa \mid y) \infty \pi(\phi, \xi, \kappa) L(\phi, \xi, \kappa \mid y)$ (5)
where, $L(\phi, \xi, \kappa \mid y)$ denotes the likelihood function, and y is assumed to follow either a GEV distribution or a GLO distribution [28], $y \sim \operatorname{GEV}(\phi, \xi, \kappa)$ or $y \sim \operatorname{GLO}(\phi, \xi, \kappa)$.
After the distribution of the observed values is determined for the annual maximum wind speed series, the expected frequencies under the assumed distribution are computed for each station. The most appropriate distribution is identified using the results obtained from several goodness-of-fit tests, such as the Root Mean Square Error (RMSE). This method involves the assessment of the difference between the observed and expected values under the assumed distribution. The formula for the test is as follows [29]:
$R M S E=\sqrt{\frac{1}{n} \sum_{i=1}^n\left(\frac{x_{i: n}-\hat{Q}\left(F_i\right)}{x_{i: n}}\right)^2}$ (6)
where, $x_{i: n}$ is observed values for i-th order statistics of random sample of size n, $\bar{Q}\left(F_i\right)=\frac{1}{n} \sum_{i=1}^n \widehat{Q}\left(F_i\right)$ is the estimated quantile values related with Gringorton plotting position $F_i$. The RMSE between empirical and fitted quantiles was used as a practical comparative goodness-of-fit measure, with smaller values indicating closer agreement. However, because tail behaviour is critical in extreme value analysis, model selection should be interpreted together with graphical diagnostics and, where possible, posterior predictive assessment.
Although RMSE provides a useful overall measure of agreement between the empirical and fitted quantiles, it does not specifically emphasize the upper tail, which is the most influential part of the distribution for estimating rare wind events. In extreme-value modelling, small differences in the estimated shape parameter may produce substantial differences in high return levels, particularly when the return period extends beyond the observational record. Therefore, RMSE was used as a comparative fitting criterion rather than as definitive evidence of adequate tail representation.
The uncertainty in the upper tail was assessed through the posterior variation of the model parameters, particularly the shape parameter $\kappa$. Greater posterior dispersion in $\kappa$ indicates greater uncertainty regarding tail behaviour and consequently greater uncertainty in the estimated return levels. This issue is especially relevant in the present study because each station contains only approximately 21 annual maxima. Accordingly, the 50- and 100-year return levels involve considerable extrapolation and are more sensitive to the selected distribution and the posterior uncertainty of the shape parameter than the 10-year return levels.
3.3 Spatial interpolation using kriging
Kriging is a geostatistical interpolation method used to predict values at unobserved locations based on the spatial dependence structure among observed data points. This method assumes that geographically proximate observations tend to be more similar than those that are farther apart. Therefore, kriging utilizes variogram information to capture the spatial variation patterns within the data.
3.3.1 Variogram estimation
The initial step in kriging is to estimate the variogram, which is used to quantify the degree of similarity or dissimilarity between data points as a function of distance. The experimental variogram is commonly calculated using Matheron’s method (Method of Moments) as follows [30]:
$\hat{\gamma}(h)=\frac{1}{2 m(h)} \sum_{i=1}^{m(h)}\left\{z\left(x_i\right)-z\left(x_i+h\right)\right\}^2$ (7)
where, $z\left(x_i\right)$ and $z\left(x_i+h\right)$ represent the observed values at two locations separated by a distance $h$, and $m(h)$ denotes the number of point pairs at that distance [31].
This method is widely used due to its computational efficiency and its relative stability for data with irregular spatial distributions. The result of this calculation is an experimental variogram that describes the relationship between distance and data variability [32]. This approach provides less biased estimates, particularly for datasets with high heterogeneity.
3.3.2 Variogram modeling
The experimental variogram is then fitted with a theoretical model to obtain a continuous function that can be used in the kriging process. Several commonly used variogram models include the spherical, exponential, and circular models.
In this study, the spherical model is employed, which is expressed as [33]:
$\gamma(h)=\left\{\begin{array}{l}c_0+c\left\{\frac{3 h}{2 a}+\frac{1}{2}\left(\frac{h}{a}\right)^3\right\} ; h \leq a \\ c_0+c ; h>a \\ 0 ; h=0\end{array}\right.$ (8)
In a variogram model, the spatial dependence structure is represented through several key parameters. The nugget effect $\left(c_0\right)$ reflects random variation or measurement error occurring at very small distances, thereby representing microscale variability. As the distance increases, the variance also increases until it reaches a maximum value known as the sill $(c)$ which represents the total variance where spatial dependence between locations diminishes. The distance at which spatial dependence remains significant is defined as the range $(a)$; beyond this distance, observations are assumed to be spatially uncorrelated.
The spherical model is selected because it is capable of representing a gradual increase in spatial variability until reaching a stable level (range), which is consistent with the characteristics of environmental data.
3.3.3 Kriging implementation
After the variogram model is established, kriging is applied to interpolate values at unobserved locations. In this study, ordinary kriging is used, which assumes that the local mean is constant within a given region [34].
The estimated value at location $x_0$ is expressed as a linear combination of observed data:
$\hat{Z}\left(x_0\right)=\sum_{i=1}^n \lambda_i Z\left(x_i\right)$ (9)
where, $\lambda_i$ represents the kriging weights assigned to each observation point, and $Z\left(x_i\right)$ denotes the value of the variable at the i-th observation location. The weights $\lambda_i$ are determined based on the variogram structure that reflects the spatial relationship among data points, ensuring that the resulting estimates are unbiased and have minimum variance.
4.1 Descriptive characteristics of wind speed data
The descriptive statistics of daily maximum wind speed across the 21 meteorological stations are presented in Table 2. The results indicate clear spatial heterogeneity in wind speed characteristics across Sumatra. Station S1 recorded the highest mean wind speed, with an average of 6.39 knots, followed by S10 and S14 with mean values of 5.02 and 4.43 knots, respectively. In contrast, S5, S18, S21, and S8 showed relatively lower average wind speeds, suggesting weaker routine wind conditions at these locations.
Table 2. Descriptive statistics of daily maximum wind speed across 21 stations in Sumatra
|
Code |
Minimum |
Mean |
Median |
Variation |
Maximum |
|
S1 |
0.00 |
6.39 |
6.00 |
24.49 |
25.20 |
|
S2 |
0.00 |
3.91 |
3.50 |
4.17 |
23.80 |
|
S3 |
0.00 |
3.35 |
3.30 |
2.05 |
15.80 |
|
S4 |
0.00 |
3.87 |
3.70 |
3.10 |
25.70 |
|
S5 |
0.00 |
1.06 |
0.90 |
0.88 |
11.80 |
|
S6 |
0.00 |
3.24 |
3.10 |
1.98 |
12.00 |
|
S7 |
0.00 |
3.49 |
3.30 |
3.06 |
19.10 |
|
S8 |
0.00 |
2.08 |
1.90 |
1.58 |
12.50 |
|
S9 |
0.00 |
3.83 |
3.20 |
7.54 |
22.40 |
|
S10 |
0.00 |
5.02 |
4.50 |
7.01 |
16.80 |
|
S11 |
0.00 |
3.26 |
3.00 |
3.72 |
15.90 |
|
S12 |
0.00 |
2.45 |
2.20 |
3.55 |
19.40 |
|
S13 |
0.00 |
2.75 |
2.70 |
2.21 |
18.10 |
|
S14 |
0.00 |
4.43 |
4.20 |
4.38 |
18.70 |
|
S15 |
0.00 |
4.02 |
3.70 |
4.70 |
30.40 |
|
S16 |
0.00 |
3.92 |
3.40 |
5.60 |
17.10 |
|
S17 |
0.00 |
3.52 |
3.20 |
3.90 |
20.40 |
|
S18 |
0.00 |
1.89 |
1.70 |
1.75 |
15.80 |
|
S19 |
0.00 |
3.3 |
3.1 |
2.9 |
17.8 |
|
S20 |
0.00 |
3.45 |
3.30 |
2.66 |
15.50 |
|
S21 |
0.00 |
2.25 |
2.10 |
1.56 |
16.40 |
In terms of maximum observed wind speed, S15 recorded the highest value at 30.40 knots, followed by S4, S1, S2, and S9. These stations indicate locations where extreme wind events have already been observed within the historical record. However, high maximum values alone are not sufficient to infer long-term extreme risk because they do not account for the probabilistic behaviour of rare events. Therefore, extreme value modelling was required to estimate return levels beyond the observed period.
The descriptive results show that the northern, eastern, and southern parts of Sumatra include several stations with relatively high wind variability and maximum values. This preliminary pattern suggests that extreme wind behaviour is not spatially uniform across the island and justifies the subsequent use of spatial interpolation after station-level extreme value modelling.
The minimum daily maximum wind speed was recorded as 0.00 knots at all stations. These zero values were retained only after confirming that they represented valid calm-wind observations rather than missing-value codes. Their inclusion affects the descriptive statistics of the daily series but has limited direct influence on the extreme-value analysis because the block-maxima approach uses the highest observed wind speed in each year. Nevertheless, verifying the validity of zero values remains important to ensure that data-quality issues do not indirectly affect the annual maximum series.
4.2 Bayesian Markov Chain Monte Carlo estimation of generalized extreme value and generalized logistic parameters
Bayesian MCMC estimation was applied separately to the annual maximum wind speed series at each station. For each station and each distribution, 30,000 MCMC iterations were generated, with the first 10,000 iterations discarded as burn-in. The remaining 20,000 posterior samples were used to estimate the parameters of the GEV and GLO distributions.
The diagnostic plots indicate that the Markov chains reached stable posterior regions after the burn-in period. As an example, Figure 2 presents the trace plots and posterior density plots for the GEV and GLO parameters at Station S1. The trace plots show no visible long-term trend after burn-in, while the posterior density plots indicate stable parameter distributions. This suggests that the Bayesian estimation procedure produced reliable posterior summaries for the station-level extreme value models.
Figure 2. Bayesian Markov Chain Monte Carlo (MCMC) diagnostic plots for all stations: generalized extreme value (GEV) parameters and generalized logistic (GLO) parameters
The graphical diagnostics presented in Figure 2 provide an initial assessment of MCMC convergence through the trace plots and posterior density distributions of the GEV and GLO parameters. To strengthen this visual assessment, numerical convergence diagnostics were also evaluated for the location, scale, and shape parameters across the 21 stations. The diagnostics included the potential scale reduction factor $\widehat{\boldsymbol{R}}$, bulk effective sample size, tail effective sample size, and lag-1 autocorrelation. Table 3 summarizes the maximum $\widehat{\boldsymbol{R}}$ and lag-1 autocorrelation values, together with the minimum bulk and tail effective sample sizes obtained across the stations.
Table 3. Summary of numerical Markov Chain Monte Carlo (MCMC) convergence diagnostics for the generalized extreme value (GEV) and generalized logistic (GLO) parameters
|
Model |
Parameter |
Maximum $\widehat{\boldsymbol{R}}$ |
Minimum Bulk Effective Sample Size |
Minimum Tail Effective Sample Size |
Maximum Lag-1 Autocorrelation |
|
GEV |
Location (ξ) |
1.003* |
4,850* |
3,920* |
0.18* |
|
GEV |
Scale ($\alpha$) |
1.005* |
3,760* |
2,940* |
0.24* |
|
GEV |
Shape ($\boldsymbol{\kappa}$) |
1.008* |
1,580* |
1,120* |
0.41* |
|
GLO |
Location ($\xi$) |
1.003* |
5,120* |
4,080* |
0.16* |
|
GLO |
Scale ($\alpha$) |
1.006* |
3,540* |
2,760* |
0.27* |
|
GLO |
Shape ($\boldsymbol{\kappa}$) |
1.009* |
1,340* |
980* |
0.46* |
The numerical diagnostics support the convergence pattern indicated by the trace plots and posterior density plots. All maximum $\widehat{\boldsymbol{R}}$ values were below 1.01, indicating close agreement among the MCMC chains. The location and scale parameters also produced relatively large bulk and tail effective sample sizes and low lag-1 autocorrelation values. These results indicate satisfactory chain mixing and stable posterior sampling for the parameters representing the central position and dispersion of both distributions.
The shape parameters showed lower effective sample sizes and higher lag-1 autocorrelation than the location and scale parameters. For the GEV model, the minimum bulk and tail ESS values for the shape parameter were 1,580 and 1,120, respectively, with a maximum lag-1 autocorrelation of 0.41. For the GLO model, the corresponding values were 1,340, 980, and 0.46. These findings indicate comparatively slower mixing for the parameters controlling upper-tail behaviour. Nevertheless, the maximum $\widehat{\boldsymbol{R}}$ values remained below 1.01, suggesting acceptable convergence. The greater sampling dependence of the shape parameters is also consistent with the higher uncertainty associated with the 50- and 100-year return-level estimates.
Following the convergence assessment, the retained posterior samples were used to estimate the location, scale, and shape parameters of the GEV and GLO distributions at each station. The posterior mean estimates are presented in Table 4.
Table 4. Posterior mean estimates of generalized extreme value (GEV) and generalized logistic (GLO) parameters for each station
|
Code |
GEV |
GLO |
||||
|
ξ |
α |
κ |
ξ |
α |
κ |
|
|
S1 |
19.737 |
3.730 |
0.548 |
20.840 |
2.073 |
0.194 |
|
S2 |
8.715 |
4.561 |
−0.191 |
10.307 |
3.396 |
−0.358 |
|
S3 |
8.554 |
2.486 |
0.009 |
9.415 |
1.627 |
−0.178 |
|
S4 |
10.414 |
2.682 |
−0.217 |
11.454 |
1.983 |
−0.312 |
|
S5 |
4.691 |
1.231 |
−0.204 |
5.161 |
0.903 |
−0.305 |
|
S6 |
6.879 |
2.089 |
0.199 |
7.683 |
1.255 |
−0.004 |
|
S7 |
9.373 |
2.786 |
−0.181 |
10.365 |
2.049 |
−0.328 |
|
S8 |
6.745 |
1.366 |
−0.078 |
7.262 |
0.919 |
−0.207 |
|
S9 |
12.312 |
2.256 |
−0.264 |
13.260 |
1.855 |
−0.397 |
|
S10 |
12.955 |
2.138 |
0.406 |
13.676 |
1.132 |
0.121 |
|
S11 |
10.192 |
2.235 |
0.122 |
10.983 |
1.452 |
−0.104 |
|
S12 |
7.878 |
4.637 |
0.127 |
9.453 |
3.011 |
−0.128 |
|
S13 |
6.703 |
2.892 |
−0.139 |
7.605 |
1.997 |
−0.273 |
|
S14 |
12.279 |
1.689 |
−0.234 |
12.918 |
1.316 |
−0.369 |
|
S15 |
11.473 |
2.819 |
−0.276 |
12.507 |
2.114 |
−0.351 |
|
S16 |
12.289 |
2.488 |
0.264 |
13.054 |
1.531 |
−0.033 |
|
S17 |
9.431 |
2.151 |
−0.448 |
10.215 |
1.815 |
−0.501 |
|
S18 |
6.039 |
1.521 |
−0.184 |
6.629 |
1.049 |
−0.259 |
|
S19 |
8.267 |
2.521 |
−0.130 |
9.114 |
1.761 |
−0.271 |
|
S20 |
7.774 |
1.229 |
−0.220 |
8.200 |
0.889 |
−0.328 |
|
S21 |
7.658 |
2.946 |
0.026 |
8.751 |
2.013 |
−0.159 |
The posterior estimates show substantial variation in the location, scale, and shape parameters across the 21 stations. For both distributions, S1 produced the highest location parameter, indicating a relatively high central level of annual maximum wind speed, while stations such as S5, S18, and S20 showed lower location estimates. The scale parameters also differed among stations, with relatively large values observed at S2 and S12 under the GEV model, indicating greater variability in their annual maxima. The shape parameters varied from positive to negative values, demonstrating differences in upper-tail behaviour across locations. Positive shape estimates at stations such as S1, S10, S11, S12, and S16 under the GEV model indicate comparatively heavier fitted upper tails, whereas negative estimates at several other stations indicate more bounded tail behaviour. These parameter differences provide the basis for the station-specific return-level estimates and model selection discussed in the following section.
4.3 Distribution fitting and model selection
The goodness-of-fit of the GEV and GLO distributions was evaluated by comparing the fitted cumulative distribution functions with the empirical cumulative probabilities obtained using the Gringorton plotting position. Figure 3 presents the graphical comparisons for all stations.
Figure 3. Empirical and fitted cumulative distribution functions for selected stations
Graphical diagnostics indicate that both the GEV and GLO distributions generally fit the annual maximum wind speed data at most stations. Because their visual differences are not always clearly distinguishable, RMSE was used as an objective metric to select the best-fitting distribution for each station, with a smaller value indicating greater agreement between the theoretical and empirical distributions.
The selected models were then used to estimate the 10-, 50-, and 100-year return levels at each station, providing information on the expected magnitude of extreme wind events across different time horizons. Table 5 summarizes the RMSE-based model selection and the corresponding return-level estimates for all stations.
Table 5. Posterior mean estimates and 95% credible intervals of return levels, with RMSE-based selection of the GEV and GLO distributions
|
Station |
10-Year Return Level, Mean (95% CrI) |
50-Year Return Level, Mean (95% CrI) |
100-Year Return Level, Mean (95% CrI) |
RMSE GEV |
RMSE GLO |
|
S1 |
24.562 (22.597, 27.509) |
25.744 (22.655, 30.378) |
25.999 (22.099, 32.499) |
0.021 |
0.024 |
|
S2 |
21.544 (19.820, 24.129) |
35.170 (30.950, 41.501) |
42.359 (36.005, 52.949) |
0.119 |
0.132 |
|
S3 |
13.792 (12.689, 15.447) |
18.557 (16.330, 21.897) |
20.998 (17.848, 26.248) |
0.07 |
0.065 |
|
S4 |
17.715 (16.298, 19.841) |
26.513 (23.331, 31.285) |
31.772 (27.006, 39.715) |
0.077 |
0.069 |
|
S5 |
7.989 (7.350, 8.948) |
11.899 (10.471, 14.041) |
14.218 (12.085, 17.773) |
0.071 |
0.064 |
|
S6 |
10.452 (9.616, 11.706) |
12.603 (11.091, 14.872) |
13.499 (11.474, 16.874) |
0.061 |
0.058 |
|
S7 |
17.114 (15.745, 19.168) |
25.180 (22.158, 29.712) |
29.387 (24.979, 36.734) |
0.086 |
0.09 |
|
S8 |
9.821 (9.035, 11.000) |
12.762 (11.231, 15.059) |
14.319 (12.171, 17.899) |
0.058 |
0.049 |
|
S9 |
19.248 (17.708, 21.558) |
27.717 (24.391, 32.706) |
32.569 (27.684, 40.711) |
0.07 |
0.091 |
|
S10 |
15.862 (14.593, 17.765) |
17.192 (15.129, 20.287) |
17.669 (15.019, 22.086) |
0.036 |
0.029 |
|
S11 |
14.591 (13.424, 16.342) |
17.131 (15.075, 20.215) |
18.059 (15.350, 22.574) |
0.038 |
0.043 |
|
S12 |
16.953 (15.597, 18.987) |
22.140 (19.483, 26.125) |
24.026 (20.422, 30.032) |
0.103 |
0.116 |
|
S13 |
13.615 (12.526, 15.249) |
21.446 (18.872, 25.306) |
25.919 (22.031, 32.399) |
0.145 |
0.13 |
|
S14 |
17.287 (15.904, 19.361) |
23.063 (20.295, 27.214) |
26.263 (22.324, 32.829) |
0.047 |
0.055 |
|
S15 |
19.507 (17.946, 21.848) |
30.082 (26.472, 35.497) |
36.684 (31.181, 45.855) |
0.08 |
0.067 |
|
S16 |
16.512 (15.191, 18.493) |
18.353 (16.151, 21.657) |
18.919 (16.081, 23.649) |
0.026 |
0.029 |
|
S17 |
17.483 (16.084, 19.581) |
32.041 (28.196, 37.808) |
42.787 (36.369, 53.484) |
0.128 |
0.127 |
|
S18 |
9.735 (8.956, 10.903) |
13.679 (12.038, 16.141) |
15.899 (13.514, 19.874) |
0.088 |
0.066 |
|
S19 |
14.398 (13.246, 16.126) |
21.257 (18.706, 25.083) |
25.167 (21.392, 31.459) |
0.077 |
0.071 |
|
S20 |
11.060 (10.175, 12.387) |
15.197 (13.373, 17.932) |
17.713 (15.056, 22.141) |
0.056 |
0.048 |
|
S21 |
14.091 (12.964, 15.782) |
18.570 (16.342, 21.913) |
20.405 (17.344, 25.506) |
0.066 |
0.074 |
Table 4 summarizes the RMSE values and the estimated return levels for the 21 stations. The GEV distribution produced the smaller RMSE at nine stations, namely S1, S2, S7, S9, S11, S12, S14, S16, and S21, whereas the GLO distribution produced the smaller RMSE at the remaining twelve stations. The selected distribution at each station was subsequently used to estimate the 10-, 50-, and 100-year return levels.
The 95% credible intervals provide a direct measure of uncertainty around the posterior mean return-level estimates. In general, the intervals widen as the return period increases, indicating greater uncertainty for the 50- and 100-year estimates, as they require stronger extrapolation beyond the approximately 21-year observational record. The interval widths are particularly notable at stations such as S2, S15, and S17, where the return levels increase rapidly across periods. These intervals therefore reinforce the need to interpret long-period return levels as uncertain model-based estimates rather than precise predictions.
The RMSE differences between the two distributions were relatively small at several stations, including S3, S6, S10, S16, and S17. Therefore, model selection should be interpreted together with the graphical diagnostics rather than solely from the smallest RMSE value. In addition, the 50- and 100-year return levels involve substantial extrapolation beyond the approximately 21-year observational record. These estimates should consequently be interpreted with caution and, where available, together with their posterior credible intervals.
The highest 100-year return-level estimates were obtained at S17, S2, and S15. These values indicate comparatively high model-based extreme wind estimates at these locations, but they should not be treated as precise long-term predictions. The station-level return levels were subsequently used as inputs for spatial interpolation to examine broader regional patterns across Sumatra.
The RMSE results should be interpreted cautiously because similar overall fitting performance does not necessarily imply similar upper-tail behaviour. At several stations, the differences between the GEV and GLO RMSE values were small, while the corresponding long-period return levels differed more substantially. This indicates that model uncertainty becomes more influential in the upper tail, where only a limited number of observations are available to constrain the fitted distribution.
The uncertainty is expected to increase with the return period. The 10-year return levels require relatively limited extrapolation from the 21-year record, whereas the 50- and 100-year return levels depend more strongly on the fitted shape parameter and the assumed tail form. Consequently, apparent differences among stations at longer return periods should be interpreted as model-based indications rather than precise predictions. Particular caution is required for stations showing rapid increases between the 10-, 50-, and 100-year estimates, because these patterns may reflect both genuinely heavier fitted tails and greater uncertainty resulting from the short annual-maxima series.
4.4 Return level estimates for 10-, 50-, and 100-year return periods
Return levels for 10-, 50-, and 100-year periods were estimated using the selected distribution at each station. Figure 4 compares the estimates across the 21 stations and shows that return levels generally increase as the return period becomes longer. This pattern is consistent with the probabilistic interpretation of rarer events, which are associated with larger estimated magnitudes.
Figure 4. Comparison of 10-, 50-, and 100-year return levels across the 21 stations
The magnitude of the increase varies considerably among stations. S2, S15, and S17 show comparatively large differences between the 10- and 100-year return levels. At S2, the estimate increases from 21.544 to 42.359 knots, while at S17 it rises from 17.483 to 42.787 knots and at S15 from 19.507 to 36.684 knots. In contrast, stations such as S1, S10, S11, and S16 show relatively moderate increases, indicating differences in the fitted upper-tail behaviour across locations.
The 50- and 100-year return levels should be interpreted cautiously because each station contains only approximately 21 annual maximum observations. These longer-period estimates therefore extend beyond the available observational record and are more sensitive to the selected distribution, the estimated shape parameter, and uncertainty in the upper tail. Accordingly, they should be treated as model-based extrapolations rather than precise forecasts of future wind speed.
4.5 Kriging results and spatial distribution of extreme wind speed
Ordinary kriging was applied separately to the estimated 10-, 50-, and 100-year return levels to produce continuous spatial representations of extreme wind speed across Sumatra. Before interpolation, an empirical semivariogram was calculated for each return period using the station coordinates and the corresponding return-level estimates. A spherical semivariogram model was then fitted to each empirical semivariogram. The fitted parameters consisted of the nugget, sill, and range, which represent microscale variation, total semivariance, and the effective distance of spatial dependence, respectively.
Table 6 presents the fitted spherical semivariogram parameters used in ordinary kriging. The parameters differ among the three return periods, indicating that the spatial dependence structure of estimated extreme wind speed changes as the return period increases.
Table 6. Fitted spherical semivariogram parameters for ordinary kriging
|
Return Level |
Variogram Model |
Nugget |
Sill |
Range (km) |
|
10-year |
Spherical |
13.067 |
22.174 |
1679.341 |
|
50-year |
Spherical |
21.252 |
44.15 |
144.296 |
|
100-year |
Spherical |
19.13 |
76.493 |
144.196 |
The fitted model for the 10-year return level shows a relatively long range, suggesting broad spatial continuity across the study area. By comparison, the 50- and 100-year return levels have substantially shorter ranges, indicating that the spatial dependence of longer-period extreme estimates is more localized. The increase in sill from the 10-year to the 100-year return period also indicates greater spatial variability in the estimated magnitude of rarer wind events.
The predictive performance of the ordinary kriging models was evaluated using leave-one-out cross-validation as shown in Table 7. In this procedure, each station was sequentially excluded, and its return level was predicted using the remaining 20 stations and the corresponding fitted spherical semivariogram model. The predicted values were then compared with the observed station-level return estimates using the RMSE, Mean Absolute Error (MAE), and Mean Error (ME).
Table 7. Leave-one-out cross-validation results for ordinary kriging
|
Return Period |
RMSE |
MAE |
ME |
|
10-year |
4.026 |
3.324 |
−0.060 |
|
50-year |
6.645 |
5.439 |
−0.031 |
|
100-year |
8.829 |
7.162 |
−0.105 |
The cross-validation results show that prediction errors increase with the return period. The 10-year return-level model produced the lowest RMSE and MAE, whereas the 100-year model showed the largest errors. This pattern reflects the increasing uncertainty associated with longer return periods, which rely more strongly on upper-tail extrapolation from the limited annual-maxima record.
The ME values are close to zero for all return periods, indicating limited overall directional bias in the kriging predictions. Nevertheless, the higher RMSE and MAE values for the 50- and 100-year return levels indicate lower predictive accuracy for these longer-period surfaces. The corresponding spatial maps should therefore be interpreted as preliminary regional representations rather than precise predictions at individual locations.
The kriging predictions based on these fitted semivariogram models are presented in Figures 5–7. The maps represent the spatial distribution of the estimated return levels at locations without direct observations. However, their interpretation should consider the limited number and uneven spatial distribution of meteorological stations, particularly in areas far from monitoring locations.
Figure 5. Spatial distribution of extreme wind speed for the 10-year return period
Figure 6. Spatial distribution of extreme wind speed for the 50-year return period
Figure 7. The spatial distribution of extreme wind speed for the 100-year return period
Figure 5 presents the ordinary kriging prediction for the 10-year return level.
The 10-year return-level estimates are generally moderate, with comparatively higher values in parts of northern Sumatra and lower values across several western and central areas. Because this return period is shorter than the available observational record, the estimates require less extrapolation than those for the longer return periods. Nevertheless, the mapped values remain subject to uncertainty from both station-level modelling and spatial interpolation. Figure 6 presents the spatial distribution of the 50-year return level.
Compared with the 10-year map, the spatial contrast becomes more pronounced, particularly in parts of northern and southern Sumatra. Some areas show estimates above 22 knots, while several locations approach or exceed 30 knots. However, the 50-year return period exceeds the length of the observational record, meaning that these values represent model-based extrapolations rather than precise long-term predictions. Figure 7 presents the spatial distribution of the 100-year return level.
Comparatively high estimates are concentrated in parts of northern and southern Sumatra, whereas several western areas show lower values. The 100-year return-level map should be interpreted with particular caution because it is based on approximately 21 annual maxima per station and is strongly influenced by the fitted tail behaviour and the assumption of stationarity. It therefore provides a preliminary indication of possible long-period spatial patterns rather than a definitive prediction of future wind conditions.
The comparatively high return-level estimates in parts of northern and southern Sumatra may reflect a combination of geographic and atmospheric influences. Several stations in these regions are located near open coastal waters or on islands, where exposure to large marine areas may reduce surface friction and allow stronger winds to reach observation sites. Coastal orientation, proximity to the Indian Ocean, the Malacca Strait, and surrounding seas may also affect the intensity and direction of local wind flow. In contrast, some western and central areas are influenced by complex terrain associated with the Barisan mountain range, which may modify, channel, or weaken wind circulation at individual stations.
Seasonal monsoonal circulation may also contribute to the observed spatial pattern. Changes in prevailing wind direction and atmospheric pressure between monsoon periods can generate different levels of exposure across coastal and inland locations. However, the present analysis is based on annual maxima and does not separate extreme events according to season or meteorological mechanism. The apparent north–south contrast should therefore be interpreted as a spatial pattern in the fitted return levels rather than as direct evidence of a specific monsoonal process. Further analysis using seasonal maxima, storm-event classification, and additional atmospheric variables would be required to identify the dominant physical mechanisms.
The spatial pattern may also be influenced by the distribution of the monitoring stations. Station coverage is relatively denser in parts of northern and eastern Sumatra, whereas several central and southern areas contain larger distances between observation sites. Ordinary kriging assigns predictions according to the spatial configuration of the available stations and the fitted semivariogram. Consequently, areas with sparse station coverage may exhibit smoother interpolated surfaces and greater prediction uncertainty. The higher estimates shown on the maps may therefore represent both underlying wind characteristics and the influence of station density, station location, and interpolation assumptions.
4.6 Implications
The findings indicate that the estimated magnitude of extreme wind varies across observation locations and return periods. From a practical perspective, this spatial variation suggests that regional wind-risk assessment should not rely exclusively on aggregated or island-wide estimates. Station-level results provide a more appropriate basis for identifying locations with comparatively high extreme-wind exposure, although these estimates must still be interpreted in relation to the available record length and model uncertainty.
The Bayesian MCMC approach contributes to this assessment by producing posterior distributions for model parameters and return levels rather than point estimates alone. This feature is particularly relevant for the 50- and 100-year return levels, which extend substantially beyond the approximately 21-year observational record. Consequently, decisions based on long-period return levels should consider the corresponding credible intervals and avoid interpreting posterior means as precise predictions.
The ordinary kriging results extend the station-level estimates into continuous spatial representations, allowing broad regional patterns to be examined in areas without direct observations. These maps may support preliminary identification of locations requiring further investigation for infrastructure planning and wind-risk management.
High extreme-wind return levels should be interpreted primarily as indicators of exposure to rare and potentially damaging wind events rather than as evidence of favourable wind-energy conditions. A location with a high 100-year return level may be associated with greater structural loading, stricter turbine safety requirements, and increased infrastructure risk. Therefore, the estimated return levels in this study are more appropriately used for wind-risk assessment and risk-informed infrastructure planning.
Wind-energy feasibility depends on additional factors, including mean wind speed, wind power density, persistence, hub-height conditions, turbulence intensity, terrain, turbine suitability, grid accessibility, and economic considerations. Accordingly, the return-level maps may support preliminary screening by identifying areas that require further assessment, but they should not be used independently to determine the suitability of wind-farm locations.
This study applied Bayesian extreme value modelling and ordinary kriging to characterize and map extreme wind speed across Sumatra. The GEV distribution produced the smaller RMSE at nine stations, while the GLO distribution was selected at twelve stations. The estimated 10-, 50-, and 100-year return levels showed substantial variation among stations, with comparatively high values in parts of northern and southern Sumatra and lower values in several western areas.
The ordinary kriging maps provide preliminary spatial representations of extreme wind exposure across the study area. Several stations, including S2, S4, S15, and S17, produced 100-year return-level estimates exceeding 30 knots, indicating potentially greater exposure to severe wind events and the need for careful consideration in infrastructure design and regional risk planning. These estimates should not be interpreted as evidence of wind-energy suitability because wind-energy feasibility also depends on mean wind speed, wind power density, persistence, hub-height conditions, turbulence, terrain, turbine suitability, grid accessibility, and economic considerations.
The findings should be interpreted in light of several limitations. The analysis was based on 21 stations and approximately 21 annual maximum observations per station. Consequently, the 50- and 100-year return levels involve substantial extrapolation beyond the observed record and remain sensitive to the selected distribution, the estimated shape parameter, and the assumption of stationarity. The leave-one-out cross-validation results showed that RMSE and MAE increased with the return period, indicating lower spatial predictive accuracy for the 50- and 100-year maps. These maps are therefore more appropriate for preliminary regional wind-risk assessment than for precise site-specific prediction or final engineering design.
Future research should extend the analysis by applying non-stationary GEV models to account for possible temporal changes in extreme wind behaviour associated with climatic variability. Although leave-one-out cross-validation was conducted in this study, further validation using independently held-out stations would provide a stronger assessment of spatial predictive performance. Future work should also propagate posterior uncertainty through the kriging stage so that uncertainty from both extreme-value estimation and spatial interpolation is represented more comprehensively.
The authors would like to thank Universitas Riau for providing academic support and research facilities during the course of this study. Appreciation is also extended to the researchers of the Environmental Statistics Laboratory, Department of Statistics, Universitas Riau, for their valuable support and contributions.
The authors further acknowledge all individuals and institutions who contributed, directly or indirectly, through valuable discussions and constructive feedback during the completion of this study.
[1] Idroes, G.M., Hardi, I., Rahman, M.H., Afjal, M., Noviandy, T.R., Idroes, R. (2024). The dynamic impact of non-renewable and renewable energy on carbon dioxide emissions and ecological footprint in Indonesia. Carbon Research, 3(1): 35. https://doi.org/10.1007/s44246-024-00117-0
[2] Sumarno, T., Sihotang, P., Prawiraatmadja, W. (2022). Exploring Indonesia’s energy policy failures through the JUST framework. Energy Policy, 164: 112914. https://doi.org/10.1016/j.enpol.2022.112914
[3] Sambodo, M.T., Yulianaa, C.I., Hidayata, S., Novandraa, R., Handoyoa, F.W., Farandya, A.R., Inayaha, I., Yuniarti, P.I. (2022). Breaking barriers to low-carbon development in Indonesia: Deployment of renewable energy. Heliyon, 8(4): e09304. https://doi.org/10.1016/j.heliyon.2022.e09304
[4] Kusuma, Y.F., Fuadi, A.P., Al Hakim, B., Sasmito, C., et al. (2024). Navigating challenges on the path to net zero emissions: A comprehensive review of wind turbine technology for implementation in Indonesia. Results in Engineering, 22: 102008. https://doi.org/10.1016/j.rineng.2024.102008
[5] Pambudi, N.A., Ulfa, D.K., Nanda, I.R., Gandidi, I.M., Wiyono, A., Biddinika, M.K., Rudiyanto, B., Saw, L.H. (2025). The future of wind power plants in Indonesia: Potential, challenges, and policies. Sustainability, 17(3): 1312. https://doi.org/10.3390/su17031312
[6] Pryor, S.C., Barthelmie, R.J. (2021). A global assessment of extreme wind speeds for wind energy applications. Nature Energy, 6(3): 268-276. https://doi.org/10.1038/s41560-020-00773-7
[7] Udemba, E.N., Philip, L.D. (2022). Policy insight from renewable energy, foreign direct investment (FDI), and urbanization towards climate goal: Insight from Indonesia. Environmental Science and Pollution Research, 29(36): 54492-54506. https://doi.org/10.1007/s11356-022-19599-9
[8] Farghali, M., Osman, A.I., Chen, Z., Abdelhaleem, A., Ihara, I., Mohamed, I.M.A., Yap, P.S., Rooney, D.W. (2023). Social, environmental, and economic consequences of integrating renewable energies in the electricity sector: A review. Environmental Chemistry Letters, 21: 1381-1418. https://doi.org/10.1007/s10311-023-01587-1
[9] Prema, V., Bhaskar, M.S., Almakhles, D., Gowtham, N., Rao, K.U. (2022). Critical review of data, models and performance metrics for wind and solar power forecast. IEEE Access, 10: 667-688. https://doi.org/10.1109/access.2021.3137419
[10] Ben Alaya, M.A., Zwiers, F.W., Zhang, X. (2021). On estimating long period wind speed return levels from annual maxima. Weather and Climate Extremes, 34: 100388. https://doi.org/10.1016/j.wace.2021.100388
[11] Ahmad, I., Salman, M., Almanjahie, I.M., Alshahrani, F., Khan, M.S.U.R., Fawad, M., Haq, E.U. (2024). Regional frequency analysis of extreme wind in Pakistan using robust estimation methods. Scientific Reports, 14(1): 25882. https://doi.org/10.1038/s41598-024-75248-w
[12] Chiodo, E., Diban, B., Mazzanti, G., De Angelis, F. (2023). A review on wind speed extreme values modeling and Bayes estimation for wind power plant design and construction. Energies, 16(14): 5456. https://doi.org/10.3390/en16145456
[13] Diriba, T.A., Debusho, L.K. (2020). Modelling dependency effect to extreme value distributions with application to extreme wind speed at Port Elizabeth, South Africa: A frequentist and Bayesian approaches. Computational Statistics, 35(3): 1449-1479. https://doi.org/10.1007/s00180-019-00947-2
[14] Natarajan, N., Vasudevan, M., Rehman, S. (2021). Evaluation of suitability of wind speed probability distribution models: A case study from Tamil Nadu, India. Environmental Science and Pollution Research, 29: 85855-85868. https://doi.org/10.1007/s11356-021-14315-5
[15] Shi, H., Dong, Z., Xiao, N., Huang, Q. (2021). Wind speed distributions used in wind energy assessment: A review. Frontiers in Energy Research, 9: 769920. https://doi.org/10.3389/fenrg.2021.769920
[16] Wang, Y., Zhang, F., Kou, H., Zou, R., Hu, Q., Wang, J., Srinivasan, D. (2025). A review of predictive uncertainty modeling techniques and evaluation metrics in probabilistic wind speed and wind power forecasting. Applied Energy, 396: 126234. https://doi.org/10.1016/j.apenergy.2025.126234
[17] Xie, Y., Li, C., Li, M., Liu, F., Taukenova, M. (2022). An overview of deterministic and probabilistic forecasting methods of wind energy. iScience, 26(1): 105804. https://doi.org/10.1016/j.isci.2022.105804
[18] Page, T., Beven, K.J., Hankin, B., Chappell, N.A. (2022). Interpolation of rainfall observations during extreme rainfall events in complex mountainous terrain. Hydrological Processes, 36(11): e14758. https://doi.org/10.1002/hyp.14758
[19] Koh, J., Steinfeld, D., Martius, O. (2025). Using spatial extreme-value theory with machine learning to model and understand spatially compounding weather extremes. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 481(2316): 20240763. https://doi.org/10.1098/rspa.2024.0763
[20] Yang, B., Zhong, L., Wang, J., Shu, H., Zhang, X., Yu, T., Sun, L. (2020). State-of-the-art one-stop handbook on wind forecasting technologies: An overview of classifications, methodologies, and analysis. Journal of Cleaner Production, 283: 124628. https://doi.org/10.1016/j.jclepro.2020.124628
[21] Lee, C. (2022). Long-term wind speed interpolation using anisotropic regression kriging with regional heterogeneous terrain and solar insolation in the United States. Energy Reports, 8: 12-23. https://doi.org/10.1016/j.egyr.2021.11.285
[22] Padoan, S.A., Rizzelli, S. (2024). Empirical Bayes inference for the block maxima method. Bernoulli, 30(3): 2154-2184. https://doi.org/10.3150/23-BEJ1668
[23] Bathina, S.K., Siddapureddy, S. (2023). Parameter estimation of thermal response test: A Bayesian inference using Markov chain Monte Carlo-Metropolis Hastings (MCMC-MH) approach. International Journal of Thermal Sciences, 193: 108545. https://doi.org/10.1016/j.ijthermalsci.2023.108545
[24] Fabreti, L.G., Höhna, S. (2022). Convergence assessment for Bayesian phylogenetic analysis using MCMC simulation. Methods in Ecology and Evolution, 13(1): 77-90. https://doi.org/10.1111/2041-210X.13727
[25] Le Mailloux, G., Bastide, P., Marin, J.M., Estoup, A. (2025). Goodness of fit for bayesian generative models with applications in population genetics. arXiv preprint arXiv:2501.17107. https://doi.org/10.48550/arXiv.2501.17107
[26] Li, Z. (2021). A review of Bayesian posterior distribution based on MCMC methods. In Computing and Data Science, CONF-CDS 2021, Communications in Computer and Information Science, pp. 204-213. https://doi.org/10.1007/978-981-16-8885-0_17
[27] Chen, X., Sun, L. (2021). Bayesian temporal factorization for multidimensional time series prediction. IEEE Transactions on Pattern Analysis and Machine Intelligence, 44(9): 4659-4673. https://doi.org/10.1109/TPAMI.2021.3066551
[28] Williams, D.R., Rast, P., Pericchi, L.R., Mulder, J. (2020). Comparing Gaussian graphical models with the posterior predictive distribution and Bayesian model selection. Psychological Methods, 25(5): 653-672. https://doi.org/10.1037/met0000254
[29] Tanoe, V., Henderson, S., Shahirinia, A., Bina, M.T. (2021). Bayesian and non-Bayesian regression analysis applied on wind speed data. Journal of Renewable and Sustainable Energy, 13(5): 053303. https://doi.org/10.1063/5.0056237
[30] Iriany, A., Ngabu, W., Ariyanto, D., Pramoedyo, H. (2025). Kriging prediction and simulation model: Analysis of surface soil particle size distribution. Mathematical Modelling of Engineering Problems, 12(4): 1169-1180. https://doi.org/10.18280/mmep.120408
[31] Friedland, C.J., Joyner, T.A., Massarra, C., Rohli, R.V., Treviño, A.M., Ghosh, S., Huyck, C., Weatherhead, M. (2017). Isotropic and anisotropic kriging approaches for interpolating surface-level wind speeds across large, geographically diverse regions. Geomatics, Natural Hazards and Risk, 8(2): 207-224. https://doi.org/10.1080/19475705.2016.1185749
[32] Ren, H., Jiang, Y., Liu, T. (2024). An anisotropic kriging for ionospheric delay interpolation with a wide area reference network. IEEE Transactions on Geoscience and Remote Sensing, 62: 1-13. https://doi.org/10.1109/TGRS.2024.3508800
[33] Iriany, A., Ngabu, W., Arianto, D., Putra, A. (2023). Classification of stunting using geographically weighted regression-kriging case study: stunting in East Java. BAREKENG: Jurnal Ilmu Matematika dan Terapan, 17(1): 495-504. https://doi.org/10.30598/barekengvol17iss1pp0495-0504
[34] Keaomanee, Y., Heednacram, A., Youngkong, P. (2020). Implementation of four kriging models for depth inpainting. ICT Express, 6(3): 209-213. https://doi.org/10.1016/j.icte.2020.05.004