A Bivariate Weibull-Epsilon Distribution with Application to Wind Speed Data

Gongsin Isaac Esbond; Funmilayo Westnand Oshogboye Saporu

doi:doi:10.11648/j.ijsda.20261203.11

Review Article |

| Peer-Reviewed

A Bivariate Weibull-Epsilon Distribution with Application to Wind Speed Data

Gongsin Isaac Esbond^*

, Funmilayo Westnand Oshogboye Saporu

Published in International Journal of Statistical Distributions and Applications (Volume 12, Issue 3)

Received: 25 July 2026 Accepted: 8 August 2026 Published: 18 September 2026

Views: Downloads:

Download PDF

Share This Article

Twitter
Linked In
Facebook

Abstract

Modelling dependence between non-normal environmental variables is important because conventional multivariate models may not adequately represent asymmetric distributions, nonlinear associations, and complex dependence structures commonly observed in environmental data. Wind speed, in particular, exhibits substantial temporal variability, making appropriate modelling of its dependence across different times of the day important for forecasting and renewable energy planning. This study develops a flexible Bivariate Weibull-Epsilon Distribution (BWED) for modelling dependence between two non-normal random variables. The study derives the joint probability density, cumulative distribution, conditional distributions, and copula representation of the proposed model, together with a copula-based measure of dependence. Model parameters are estimated using the Exact Maximum Likelihood (EML) method, while Gibbs sampling is employed to generate observations and evaluate the finite-sample performance of the estimators. Simulation experiments are conducted for sample sizes ranging from 20 to 5,000, the results show that parameter estimates become increasingly accurate and precise as sample size increases, with bias and standard errors decreasing substantially. A comparison with two existing bivariate Weibull models showed the BWED provides a superior fit, with an AIC of 3468.56 compared with 4257.98 and 3886.72 for the competing models. The proposed model is subsequently applied to 517 paired observations of wind speed recorded at 9:00 a.m. and 3:00 p.m. on the same day at Christmas Island, Australia, covering October 2018 to March 2020. The estimated Kendall’s tau of 0.583 indicates moderate positive monotonic dependence between morning and afternoon wind speeds. The study concludes that the BWED provides a flexible, robust, and tractable framework for modelling dependence in wind-speed and other environmental data, with potential applications in renewable energy forecasting and risk management.

Published in	International Journal of Statistical Distributions and Applications (Volume 12, Issue 3)
DOI	10.11648/j.ijsda.20261203.11
Page(s)	46-60
Creative Commons	This is an Open Access article, distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution and reproduction in any medium or format, provided the original work is properly cited.
Copyright	Copyright © The Author(s), 2026. Published by Science Publishing Group

Keywords

Copula Dependence, Exact Maximum Likelihood, Gibb’S Sampling, Multivariate Distribution, Renewable Energy Modelling, Wind Speed

1. Introduction

Modelling dependence among non-normal random variables has become increasingly important in modern statistical and applied research, largely because real-world datasets often violate the normality assumption. Multivariate Gaussian models are often inadequate when variables exhibit heavy tails, asymmetries, or nonlinear dependencies. Copula methods, which separate marginal behavior from joint dependence, offer a flexible framework for constructing valid joint distributions even when marginal distributions are non-Gaussian

[1, 2]

The modern appeal to copulas stems from the work of Embrechts et al.

[3]

, whose seminal exposition showed how to model non-normal dependence via copula theory. As a result, copula models have rapidly gained traction across many fields. In finance, directional dependence in foreign exchange markets has been studied using copulas

[4]

, and copulas have been used in pricing collateralized debt obligations

[5]

. In hydrology and climatology, copulas have enabled multivariate drought assessment

[6]

. In genomics, they have been used to capture dependence structures between genes

[7]

. And in reliability and engineering, copulas have been applied to model dependence in wind speeds or other system components

[8]

. Recently, dependence among anthropometric variables were also modeled using copulas

[9]

In recent literature, many further developments in copula theory and applications have been witnessed. For example, large skew-t copula models have been proposed to capture asymmetric and heavy-tail dependence in financial time series

[10]

. New reviews highlight evolving trends such as vine copula simplification, factor copulas, and applications in censored survival settings

[11, 12]

. Copula-based joint regression frameworks and survival models are receiving more attention, especially for non-normal correlated outcomes

[13]

Beyond copulas, the literature also contains direct methods for constructing multivariate non-normal distributions. Early expositions include Ref.s

[14-16]

. These approaches generally specify joint densities or cumulative functions in closed form, often by extending univariate distributions or using conditioning constructs. In the conditional approach, one builds a bivariate (or multivariate) distribution by specifying a marginal and then a conditional distribution. For example, Gongsin & Saporu

[17]

developed a bivariate conditional Weibull model, and Sens et al.

[18]

proposed a bivariate distribution combining conditional gamma constructs. Another path uses the probability integral transform, for example, Ref.

[19]

constructed multivariate Weibull distribution by mapping the Weibull variables into standard normal variables for wind speed and wind power modelling. In Ref.

[20]

the power law transformations of vector components with Gaussian, isotropic, mean-zero variability, was used to construct a bivariate distribution for modelling wind speed from different altitudes.

Nevertheless, many of the existing constructions are limited in flexibility, that is, they may only produce simple hazard shapes, or be difficult to generalize. In environmental and reliability contexts, where variables like wind speed often exhibit complex dependence and nonlinear tails, richer multivariate models are needed. Indeed, recent joint modelling of wind and waves, or wind and other environmental covariates, relies heavily on vine copula frameworks, see for example,

[21]

. The Clayton copula was used to model the joint distribution of extreme wind speed and air density for structural wind load assessment

[22]

. However, which bivariate model to use for a given situation is not generally known a priori and is usually determined empirically through statistical model selection criteria.

In light of these developments, we propose the bivariate Weibull-epsilon distribution, defined via a joint survival distribution method. This approach combines features of the epsilon (baseline) distribution and Weibull margins, and seeks to integrate dependence directly through the joint survival function rather than relying solely on a copula plug-in. The resulting model builds on and extends the univariate Weibull-epsilon architecture to the bivariate setting, offering more flexible joint behavior and tail dependence features. Specifically, our objectives are to derive explicit expressions for the joint density and cumulative distribution of the bivariate Weibull-epsilon model; analyze its dependence properties; compare the flexibility and fit of the model against other bivariate distributions in the Weibull family; and, apply the proposed model to actual wind speed datasets, demonstrating its usefulness in modelling dependence in environmental applications. The motivation for deriving the new distribution is the complexity in the expression of the bivariate Weibull distributions in Villanueva et al.

[19]

and Monahan

[20]

. The expressions of the dependence relationship of the new distribution is simpler, which makes for easier parameter estimation.

The rest of the paper is organized as follows: Section 2 briefly looks at the bivariate survival function; in Section 3 the new distribution is derived along with the methods for estimating the co-dependence between the bivariate datasets, and also estimating the parameters of the distribution. In section 4, the robustness of the new distribution is examined through comparison with other bivariate distributions in the Weibull family. Modelling wind speed data using the new distribution is also presented in this section. Section 5 concludes the study.

2. Theoretical Framework

The joint survival function of any two random variables,

X_{1}

and

X_{2}

, can be obtained by letting

h (x_{1})

and

h (x_{2})

be their respective failure rate functions on

[0, \infty)

, with corresponding cumulative failure rate functions

H (x_{1})

and

H (x_{2})

. Then according to Ref.

[23]

, their joint survival function is given by

\overset{̅}{F} (x_{1}, x_{2}) = \exp \{- {[{(H (x_{1}))}^{\frac{1}{φ}} + {(H (x_{2}))}^{\frac{1}{φ}}]}^{φ}\},

(1)

where

H (x_{i}) = - \log (1 - F (x_{i})) = - logS (x_{i}), i = 1, 2

and

φ \in (0, 1]

is the dependence parameter. When

φ = 1

, it implies

X_{1}

and

X_{2}

are independent.

The derivative of Eq. (1) with respect to the variables

X_{1}

and

X_{2}

yields a bivariate probability density function. Eq. (1) has been used in the literature to construct the bivariate Weibull probability density function and generalized to its multivariate form

[24]

3. Methodology

3.1. The Bivariate Weibull-epsilon Distribution

The derivation of the BWED from Eq. (1), for random variables

X_{1}

and

X_{2}

, is presented in Appendix Ⅰ. However, the new distribution is given by

f (x_{1}, x_{2}) = \frac{α_{1} λ_{1} {δ_{1}}^{2}}{{δ_{1}}^{2} - {x_{1}}^{2}} \frac{α_{2} λ_{2} {δ_{2}}^{2}}{{δ_{2}}^{2} - {x_{2}}^{2}} z (x_{1}) z (x_{2}) {[z (x_{1}) - 1]}^{\frac{α_{1}}{φ} - 1} {[z (x_{2}) - 1]}^{\frac{α_{2}}{φ} - 1} {({[z (x_{1}) - 1]}^{\frac{α_{1}}{φ}} + {[z (x_{2}) - 1]}^{\frac{α_{2}}{φ}})}^{φ - 2}

\{{({[z (x_{1}) - 1]}^{\frac{α_{1}}{φ}} + {[z (x_{2}) - 1]}^{\frac{α_{2}}{φ}})}^{φ} + \frac{1}{φ} - 1\} \exp \{- {({[z (x_{1}) - 1]}^{\frac{α_{1}}{φ}} + {[z (x_{2}) - 1]}^{\frac{α_{2}}{φ}})}^{φ}\},

(2)

where

z (x_{i}) = {(\frac{x_{i} + δ_{i}}{δ_{i} - x_{i}})}^{\frac{λ_{i} δ_{i}}{2}}

α_{i} > 0

λ_{i} > 0

δ_{i} > 0

0 < φ \leq 1

, and

0 < x_{i} < δ_{i}

i = 1, 2

The corresponding bivariate Weibull-epsilon cumulative distribution function is given by

F (x_{1}, x_{2}) = \exp \{- {({[z (x_{1}) - 1]}^{\frac{α_{1}}{φ}} + {[z (x_{2}) - 1]}^{\frac{α_{2}}{φ}})}^{φ}\} - \exp (- {[z (x_{1}) - 1]}^{α_{1}}) - \exp (- {[z (x_{2}) - 1]}^{α_{2}}) + 1 .

(3)

The parameter

φ

governs the strength of dependence, where values close to 1 indicates weak dependence and smaller values indicate strong dependence (or association).

The validity of the BWED is presented in Appendix Ⅱ. Plots of the BWE density and the corresponding contour for varying parameter values are presented in Figures 1-3.

Download: Download full-size image

Figure 1. Joint probability density and corresponding contour plots of the BWE density for varying values of

λ_{1} = λ_{2}

with fixed parameters

α_{1} = α_{2} = 3

δ_{1} = δ_{2} = 15.5

φ = 0.5

Download: Download full-size image

Figure 2. Joint probability density and corresponding contour plots of the BWE density for varying values of

α_{1} = α_{2}

;, with fixed parameters

λ_{1} = λ_{2} = 0.1

δ_{1} = δ_{2} = 15.5

φ = 0.5

Download: Download full-size image

Figure 3. Joint probability density and corresponding contour plots of the BWE density for varying values of

φ

, with fixed parameters

α_{1} = α_{2} = 3

λ_{1} = λ_{2} = 0.1

δ_{1} = δ_{2} = 25

From Figure 1, the plots illustrate how changes in the scale parameters influence the concentration, spread, and dependence structure of the joint distribution. In Figure 2, the plots demonstrate the effect of the shape parameters on the skewness, tail behavior, and peak concentration of the joint density. As the shape parameters vary, the contour lines change in curvature and orientation, reflecting alterations in marginal behavior and dependence structure. In Figure 3, the plots illustrate the role of the dependence parameter

φ

in regulating the strength and structure of association between the two variables. As

φ

varies, the density surface and contour configurations exhibit noticeable changes in concentration, symmetry, and dependence intensity, ranging from stronger (top-left panel) to weaker (bottom-right panel) joint interactions.

The BWE density plots together with their corresponding contour representations for the varying parameter configurations show diverse shapes observed under different parameter settings, illustrating the potential for adaptability of the BWE model in representing a wide range of bivariate data-generating processes. The contour structures further emphasize the model’s flexibility in capturing varying degrees of dependence, association patterns, and tail behavior. Collectively, these graphical features demonstrate the suitability of the BWED for modelling complex bivariate relationships encountered in practical applications, while maintaining the prescribed marginal characteristics. The generalization of the BWED to its multivariate form is given in Appendix Ⅲ.

3.2. Conditional Distributions

The conditional distribution of

X_{1}

given

X_{2} = x_{2}

can be expressed as

f (X_{1} | X_{2} = x_{2}) = \frac{α_{1} λ_{1} {δ_{1}}^{2}}{{δ_{1}}^{2} - {x_{1}}^{2}} z (x_{1}) {[z (x_{1}) - 1]}^{\frac{α_{1}}{φ} - 1} \frac{{[z (x_{2}) - 1]}^{α_{2} (\frac{1}{φ} - 1)}}{\exp (- {[z_{2} (x_{2}) - 1]}^{α_{2}})} {({[z (x_{1}) - 1]}^{\frac{α_{1}}{φ}} + {[z (x_{2}) - 1]}^{\frac{α_{2}}{φ}})}^{φ - 2}

∙ \{{({[z (x_{1}) - 1]}^{\frac{α_{1}}{φ}} + {[z (x_{2}) - 1]}^{\frac{α_{2}}{φ}})}^{φ} + \frac{1}{φ} - 1\} \exp \{- {({[z (x_{1}) - 1]}^{\frac{α_{1}}{φ}} + {[z (x_{2}) - 1]}^{\frac{α_{2}}{φ}})}^{φ}\} .

(4)

Conditional distributions are used in simulating the BWED via Gibb’s sampling technique in this study. They can also be used in regression and risks modelling.

3.3. Copula-based Dependence Measure

Sklar’s theorem

[25]

states that any multivariate joint distribution can be written in terms of univariate marginal distribution functions and a copula (link function) which describes the dependence structure between the variables. That is, for marginal cumulative distribution functions denoted by

F (x_{1}), F (x_{2}), \dots, F (x_{n})

and a joint cumulative distribution function

F (∙)

of random variables

x_{1}, x_{2}, \dots, x_{n}

, there exist a copula function,

C (∙)

, such that

F (x_{1}, x_{2}, \dots, x_{n}) = C (F (x_{1}), F (x_{2}), \dots, F (x_{n})) .

(5)

Hence, given random variables

X_{1}

and

X_{2}

from the BWED with respective marginal distributions

F (x_{1})

and

F (x_{2})

, their copula (link) distribution and density functions are given, respectively, by

C_{φ} (u_{1}, u_{2}) = \exp \{- {[{(- ζ (u_{1}))}^{\frac{1}{φ}} + {(- ζ (u_{2}))}^{\frac{1}{φ}}]}^{φ}\} + u_{1} + u_{2} - 1,

(6)

and

c_{φ} (u_{1}, u_{2}) = \frac{{(- ζ (u_{1}))}^{\frac{1}{φ} - 1}}{1 - u_{1}} \frac{{(- ζ (u_{2}))}^{\frac{1}{φ} - 1}}{1 - u_{2}} {[{(- ζ (u_{1}))}^{\frac{1}{φ}} + {(- ζ (u_{2}))}^{\frac{1}{φ}}]}^{φ - 2}

∙ \{{[{(- ζ (u_{1}))}^{\frac{1}{φ}} + {(- ζ (u_{2}))}^{\frac{1}{φ}}]}^{φ} + \frac{1}{φ} - 1\} \exp \{- {[{(- ζ (u_{1}))}^{\frac{1}{φ}} + {(- ζ (u_{2}))}^{\frac{1}{φ}}]}^{φ}\},

(7)

where

u_{i} = F (x_{i})

ζ (u_{i}) = \log (1 - u_{i}), i = 1, 2

and

0 < φ \leq 1

One of the most desired outcomes in the study of multivariate distributions is the understanding of co-dependence structure among the variables. Unlike Pearson’s correlation, which measures straight-line relationships, monotonic dependence (often measured by the Spearman’s rho and Kendall’s tau) captures nonlinear but consistent trends. The strength of monotonic dependence is measured by copula-based Kendall’s tau, denoted by

τ_{φ}

. This can be obtained

[26]

from Eq. (5), using

τ_{φ} = 4 \int_{0}^{1} \int_{0}^{1} C_{φ} (u_{1}, u_{2}) c_{φ} (u_{1}, u_{2}) d u_{1} d u_{2} - 1,

(8)

τ_{φ}

is the measure of nonlinear monotonic dependence relationship that will be used in this study. The strength of this measure as classified in www.statisticshowto.com and used in Saporu & Gongsin (2020)

[9]

is employed to aid the understanding of our interpretation.

3.4. Parameter Estimation

The EML method is used for estimation purpose. The EML method enables simultaneous estimation of all the seven parameters in the bivariate Weibull-epsilon distribution at once, using the log-likelihood function of the BWE density in Eq. (2). This function is coded in R and executed via the optim package. The log-likelihood function for the bivariate Weibull-epsilon distribution is given by

L (θ | x_{1}, x_{2}) = \sum_{j = 1}^{n} \log (\frac{α_{1} λ_{1} {δ_{1}}^{2}}{{δ_{1}}^{2} - {x_{1 j}}^{2}} \frac{α_{2} λ_{2} {δ_{2}}^{2}}{{δ_{2}}^{2} - {x_{2 j}}^{2}} z (x_{1 j}) z (x_{2 j})) + \sum_{j = 1}^{n} \log ({[z (x_{1 j}) - 1]}^{\frac{α_{1}}{φ} - 1} {[z (x_{2 j}) - 1]}^{\frac{α_{2}}{φ} - 1})

+ \sum_{j = 1}^{n} \log {({[z (x_{1 j}) - 1]}^{\frac{α_{1}}{φ}} + {[z (x_{2 j}) - 1]}^{\frac{α_{2}}{φ}})}^{φ - 2} + \sum_{j = 1}^{n} \log \{{({[z (x_{1 j}) - 1]}^{\frac{α_{1}}{φ}} + {[z (x_{2 j}) - 1]}^{\frac{α_{2}}{φ}})}^{φ} + \frac{1}{φ} - 1\}

+ \sum_{j = 1}^{n} {({[z (x_{1 j}) - 1]}^{\frac{α_{1}}{φ}} + {[z (x_{2 j}) - 1]}^{\frac{α_{2}}{φ}})}^{φ},

(9)

where,

θ

is the parameter vector space and

x_{i} = (x_{i 1}, x_{i 2}, ∙ ∙ ∙, x_{in}), i = 1, 2

, are random samples of size

n

The execution of the log-likelihood function (Eq. (9)), simulation and other computations are carried out in R, and the script will be made available upon request.

τ_{φ}

will be computed from Eq. (8) through numerical integration in R package, curvature.

4. Results and Discussion

4.1. Simulation

Here, we simulate data from the bivariate Weibull-epsilon distribution for given parameter values using Gibb’s sampling method. Letting

X = X_{1}

and

Y = X_{2}

, we define the conditional distributions of

y / X = x

and that of

x / Y = y

from Eq. (4). Using initial value for

x

we simulate

y

from

f (y / X = x)

, employ the simulated value of

y

into

f (x / Y = y)

to simulate for

x

. The process is repeated using the previous simulated value of one variable to simulate the next other value until the required sample size is obtained. The performance of the parameter estimates for 20 - 5,000 realizations of sampled values are presented in Table 1.

From Table 1, the simulation results indicate that parameter estimates improve steadily with increasing sample size, reflecting the maximum likelihood estimates’ asymptotic consistency. For small samples (

n < 100

), estimates show noticeable bias and higher variability, while for larger samples (

n \geq 1000

), they converge closely to the true parameter values with markedly reduced standard errors.

Table 1. Parameter Estimates Performance of the BWED from the Sampled Bivariate Datasets.

$n$	Parameters
$n$	$α_{1} (se)$	$α_{2} (se)$	$λ_{1} (se)$	$λ_{2} (se)$	$δ_{1} (se)$	$δ_{2} (se)$	$φ (se)$
Actual values	1.5	3.0	0.05	0.03	17.5	20.0	0.25
20	0.961 (0.173)	3.830 (1.167)	0.156 (0.025)	0.048 (0.006)	67.024 (105.64)	15.261 (2.215)	0.646 (0.143)
50	1.452 (0.216)	2.097 (0.355)	0.111 (0.010)	0.031 (0.002)	31.274 (41.043)	15.159 (0.226)	0.564 (0.077)
70	1.768 (0.209)	2.456 (0.373)	0.108 (0.007)	0.030 (0.002)	44.989 (88.176)	15.193 (0.219)	0.584 (0.067)
100	2.031 (0.196)	4.414 (0.973)	0.100 (0.005)	0.037 (0.004)	27.012 (12.952)	16.844 (1.446)	0.628 (0.058)
150	2.201 (0.205)	3.660 (0.598)	0.088 (0.004)	0.032 (0.002)	22.019 (7.077)	16.702 (0.622)	0.543 (0.044)
200	1.938 (0.166)	4.321 (0.556)	0.076 (0.003)	0.036 (0.002)	28.545 (10.902)	20.122 (1.307)	0.386 (0.029)
250	1.989 (0.112)	5.004 (0.585)	0.076 (0.002)	0.041 (0.002)	48.695 (23.661)	26.819 (4.793)	0.388 (0.025)
500	1.845 (0.105)	4.094 (0.362)	0.068 (0.002)	0.038 (0.002)	28.591 (6.791)	24.151 (2.117)	0.304 (0.015)
750	1.754 (0.063)	3.589 (0.187)	0.076 (0.002)	0.039 (0.001)	63.484 (45.204)	23.869 (1.227)	0.254 (0.010)
1000	1.672 (0.055)	3.481 (0.182)	0.078 (0.002)	0.039 (0.001)	39.607 (10.938)	23.408 (1.217)	0.252 (0.009)
1500	1.471 (0.042)	2.941 (0.098)	0.071 (0.001)	0.036 (0.001)	21.672 (1.308)	20.723 (0.341)	0.232 (0.007)
2000	1.452 (0.035)	2.968 (0.083)	0.066 (0.001)	0.035 (0.0004)	19.216 (0.602)	20.755 (0.265)	0.236 (0.006)
3000	1.420 (0.028)	2.847 (0.066)	0.063 (0.0008)	0.033 (0.0003)	18.347 (0.342)	20.139 (0.171)	0.233 (0.005)
5000	1.501 (0.023)	3.006 (0.054)	0.059 (0.0006)	0.032 (0.0003)	18.456 (0.254)	20.181 (0.123)	0.248 (0.004)

Specifically, the shape parameters (

α_{1}, α_{2}

) show gradual reduction in bias;

α_{1}

improves from 0.961 (bias = − 0.359) at

n = 20

to 1.501 (bias = 0.0007) at

n = 5000

, while

α_{2}

decreases from 3.83 (bias = 0.277) to 3.006 (bias = 0.002). Standard errors shrink by nearly an order of magnitude, from nearly between 0.17 and 1.17 to nearly between 0.02 and 0.05, indicating increasing precision. Secondly, the scale parameters (

λ_{1}

λ_{2}

) exhibit rapid stability;

λ_{1}

converges from 0.156 (bias = 2.12) at

n = 20

to 0.059 (bias = 0.18) at

n = 5000

, and

λ_{2}

from 0.048 (bias = 0.6) to 0.032 (bias = 0.067). Variability in the scale parameters decline sharply after

n = 200

, with standard errors falling below 0.001 for large samples. Thirdly, the location parameters (

δ_{1}, δ_{2}

) initially fluctuate widely, for example,

δ_{1} = 67.02

versus the true value 17.5 at

n = 20

which progressively stabilize around the true values for

n \geq 1500

, where at

n = 5000

δ_{1} = 18.46

and

δ_{2} = 20.18

with biases of 0.055 and 0.009, respectively. Finally, the dependence parameter (

φ

) displays consistent downward adjustment with increasing

n

; from 0.646 (bias = 1.58) at

n = 20

to 0.248 (bias = − 0.008) at

n = 5000

, confirming convergence toward the true value 0.25. Generally, bias and variance decline monotonically across parameters, and by

n = 1000

, all estimates closely approximate their true values. These confirm the consistency, efficiency, and robustness of the EML estimators of the BWE model under moderate to large sample sizes.

4.2. Robustness of the Bivariate Weibull-epsilon Distribution

To show the robustness of the bivariate Weibull-epsilon distribution, we compare with other bivariate distributions in the Weibull family applied in wind speed modelling, namely the bivariate Weibull distribution (BWD₁) in Villanueva et al.

[19]

and the bivariate Weibull distribution (BWD₂) in Monahan

[20]

. We simulated data from these distributions and fitted each to the parent distribution and the BWED. The Akaike Information Criterion (AIC) is used to determine the more compatible model.

The bivariate Weibull distributions in

[19]

(after modification) and

[20]

are given in Eq.s (10) and (11), respectively.

f (v_{1}, v_{2}) = \frac{1}{\sqrt{1 - ρ}} \frac{k_{1}}{c_{1}} \frac{k_{2}}{c_{2}} {(\frac{v_{1}}{c_{1}})}^{k_{1} - 1} {(\frac{v_{2}}{c_{2}})}^{k_{2} - 1} \exp \{- [{(\frac{v_{1}}{c_{1}})}^{k_{1}} + {(\frac{v_{2}}{c_{2}})}^{k_{2}}]\} ∙ \exp \{- \frac{ρ}{2 (1 - ρ^{2})} [ρ {x_{1}}^{2} + ρ {x_{2}}^{2} - 2 x_{1} x_{2}]\}

(10)

where

x_{i} = \sqrt{2} \erf^{- 1} (1 - 2 \exp [- {(\frac{v_{i}}{c_{i}})}^{k_{i}}]), i = 1, 2

and

\erf^{- 1} (∙)

is the inverse of the error function, and

f (w_{1}, w_{2}) = \frac{1}{1 - ρ^{2}} \frac{b_{1}}{a_{1}} \frac{b_{2}}{a_{2}} {(\frac{w_{1}}{a_{1}})}^{b_{1} - 1} {(\frac{w_{2}}{a_{2}})}^{b_{2} - 1} \exp \{- \frac{1}{1 - ρ^{2}} [{(\frac{w_{1}}{a_{1}})}^{b_{1}} + {(\frac{w_{2}}{a_{2}})}^{b_{2}}]\} ∙ I_{0} [\frac{2 ρ}{1 - ρ^{2}} {(\frac{w_{1}}{a_{1}})}^{\frac{b_{1}}{2}} {(\frac{w_{2}}{a_{2}})}^{\frac{b_{2}}{2}}]

(11)

where

I_{0} (∙)

is the Bessel function of the order zero expressed as

I_{0} (z) = \sum_{κ = 0}^{\infty} \frac{z^{2 κ}}{2^{2 κ} κ! Γ (κ + 1)}

[27]

. In both Eq.s (10) and (11),

ρ

represents dependence, with

ρ \to 1

indicating strong to perfect correlation. Estimated results from 500 realizations of the simulated values from the two distributions are given in Table 2.

Table 1. Comparison of model performance for datasets simulated from two bivariate Weibull models in Eq.s (10) and (11).

Density	Parameters							$- LL$	$AIC$
${BWED}^{֎}$	$α_{1}$	$α_{2}$	$λ_{1}$	$λ_{2}$	$δ_{1}$	$δ_{2}$	$φ$
${BWD}_{1}^{*}$	$c_{1} = 5$	$c_{2} = 9$	$k_{1} = 3$	$k_{2} = 2$	$ρ = 0.8$
${BWD}_{1}^{* *}$	4.866 (0.046)	8.761 (0.131)	4.878 (0.157)	3.103 (0.108)	0.206 (0.045)			1907.37	3824.74
${BWED}^{* *}$	3.401 (0.112)	0.142 (0.004)	110.12 (119.80)	2.184 (0.012)	0.078 (0.001)	102.44 (116.77)	0.700 (0.029)	1892.12	3798.24
${BWD}_{2}^{*}$	$a_{1} = 7.5$	$a_{2} = 8.5$	$b_{1} = 4$	$b_{2} = 3$	$ρ = 0.8$
${BWD}_{2}^{* *}$	9.178 (0.100)	11.005 (0.135)	4.122 (0.033)	3.505 (0.022)	0.698 (0.003)			2192.92	4395.84
${BWED}^{* *}$	4.127 (0.164)	0.087 (0.001)	53.138 (33.952)	3.066 (0.279)	0.073 (0.004)	21.839 (7.039)	0.993 (0.024)	1942.04	3898.08

^*actual parameter values, ^**estimated parameter values (standard errors), ^֎parameters of the BWED

Table 2 compares the proposed BWED with two existing Weibull-based models developed by Villanueva et al. (

{BWD}_{1}

)

[19]

and Monahan (

{BWD}_{2}

)

[20]

. The comparison uses simulated datasets generated from each of these two models and evaluates performance using the AIC, where smaller AIC values indicate a better model fit. When data were generated from

{BWD}_{1}

, the BWED produced a lower AIC (3798.24) than the fitted

{BWD}_{1}

model (3824.74), showing that BWED describes the joint behavior better. The difference in AIC (

∆ AIC \approx 26.5

) is large enough to indicate strong support for BWED over the original model. When data were generated from BWD₂, the advantage of BWED became even more pronounced, with an AIC of 3898.08 compared to 4395.84 for

{BWD}_{2}

(

∆ AIC \approx 498

). The consistent superiority of BWED across both scenarios suggests that it remains reliable even when the true dependence structure differs from its own formulation. Overall, the results demonstrate that BWED is a flexible and robust model for representing joint wind-speed behavior and is therefore suited for practical environmental and renewable energy applications.

4.3. Application to Wind Speed Data

4.3.1. Motivation

Wind energy is one of the most promising alternatives to fossil fuels, offering a path to decarbonization and sustainable power generation. The growth of wind power globally has been extraordinary: for example, installed wind capacity increased from 158.5 GW in 2009 to over 542 GW by 2018

[28, 29]

. In 2020, wind power additions surged by approximately 93 GW, representing a 53% increase in new installations alone

[30]

, with onshore wind capacity hitting about 699 GW

[29]

. More recent data show that global capacity continued to grow to a total (onshore + offshore) exceeding 1000 GW by 2023

[31]

. This accelerating deployment is driven by urgent climate targets, falling renewable technology costs, and increased policy support.

Despite this promising growth, one of the persistent challenges in wind energy planning, forecasting, and management is the intermittent and variable nature of wind speed - both in time and across spatial locations. Wind speeds fluctuate due to meteorological variability, terrain effects, and atmospheric dynamics; and capturing their dependence structure is critical for improved prediction, site layout design, and reliability assessment of wind farms. In particular, understanding temporal dependence of wind speeds (e.g. at two different times of the same day at the same location) can reduce forecasting error and enhance wind farm modelling.

A substantial body of literature have investigated spatial dependence in wind speeds. For example, Shen et al.

[32]

examined correlations among turbines under varying meteorological conditions; Liu et al.

[33]

explored correlation structures, both within and between wind farms; and Kadhem et al.

[34]

associated wind variability with seasonal weather patterns. These works, however, focus predominantly on dependence across locations at the same time instant. In contrast, our interest is temporal dependence, that is modelling wind speeds at two distinct times (e.g. morning vs afternoon) recorded at the same location. This type of dependence is especially relevant in forecasting, diurnal cycle analysis, and operational decision-making.

To this end, we adopt the BWED (Eq. (2)) to model the joint behavior of wind speeds at two time points of the same day. The key goals are: (1) fit the BWED to empirical wind speed data and quantify the degree of co-dependence between the two time-point variables; (2) compare model performance with two alternative bivariate Weibull constructions (Eqs. (10) and (11)); and (3) assess whether the enhanced flexibility of the BWED improves fit and/or dependence capture, in the context of wind speed data application.

4.3.2. The Data and Parameter Estimation

We analyzed daily wind speed records from Christmas Island, Australia, collected at two time points: 9:00 a.m. and 3:00 p.m of the same day. The sampling period spans from October 1, 2018 to March 25, 2020. Days with missing observations at either time point were removed, yielding 517 complete bivariate observations out of an initial 542 days. The original measurements reported in km/h were transformed to m/s using a conversion factor of 0.278. The data source is 错误!超链接引用无效。

The descriptive statistics of the wind speed data are given in Table 3.

Table 3. Descriptive Statistics of the Wind Speed.

Time	Min	$Q_{1}$	$M$	$μ$	$Q_{3}$	Max
9 a. m.	0.56	3.61	5.28	5.01	6.12	9.17
3 p. m.	1.11	3.61	4.73	4.84	6.12	8.34

cor = 0.807

The descriptive summaries in Table 3 indicates positive correlation between the wind speeds at 9 a. m. and 3 p. m.

We fitted three models to the data, namely: the proposed BWED model (Eq. (2)); the bivariate Weibull model in Villanueva et al.

[19]

(Eq. (10)); and the bivariate Weibull model in Monahan

[20]

(Eq. (11)). For all the models, parameters were estimated using maximum likelihood method. In the BWED and Villanueva

[19]

model cases, we maximized the log-likelihood (Eq. 9) using the optim package in R. For Monahan’s model

[20]

, estimation encountered identifiability or likelihood surface irregularities, thus, the three parameters (

a_{1}, a_{2}, b_{1}

) were estimated via direct maximization, while the remaining two parameters (

b_{2}, ρ

) were located by examining the profile of the log-likelihood surfaces, holding other parameters fixed. The models’ fit and parameter estimates are summarized in Tables 4 and 5.

Table 4. Maximum likelihood estimates and Wald significance tests for the parameters of the BWED fitted to Christmas Island wind speed data.

Parameter	${\hat{α}}_{1}$	${\hat{λ}}_{1}$	${\hat{δ}}_{1}$	${\hat{α}}_{2}$	${\hat{λ}}_{2}$	${\hat{δ}}_{2}$	$\hat{φ}$
Estimate	1.890	0.107	10.988	2.199	0.118	12.874	0.417
Std error	0.096	0.003	0.655	0.128	0.004	1.872	0.019
Wald stat	19.618	36.267	16.780	17.201	32.065	6.877	22.064
p-value	0.000	0.000	0.000	0.000	0.000	0.000	0.000
Remark	Sig	Sig	Sig	Sig	Sig	Sig	Sig

Sig = significant,

LL = - 1727.28

AIC = 3468.56

Table 5. Comparative fit results for competing bivariate Weibull-type models fitted to the Christmas Island wind speed data.

Density	Parameters					$- LL$	$AIC$
${BWD}_{1}$	$a_{1}$	$a_{2}$	$b_{1}$	$b_{2}$	$ρ$
	5.930 (0.081)	5.726 (0.082)	3.243 (0.114)	3.431 (0.132)	0.927 (0.003)	2123.99	4257.98
${BWD}_{2}$	$c_{1}$	$c_{2}$	$k_{1}$	$k_{2}$	$ρ$
	5.606 (0.076)	5.387 (0.067)	3.197 (0.112)	3.477 (0.120)	0.390 (0.036)	1938.36	3886.72

All estimated parameter values are statistically significant. The marginal shape parameters (

{\hat{α}}_{1} = 1.89

{\hat{α}}_{2 .} = 2.199

) are highly significant and confirm that the wind-speed distributions are not memoryless and exhibit pronounced peaks, a characteristic commonly observed in Weibull-type wind speed behavior. The scale parameters (

{\hat{λ}}_{1} = 0.107

{\hat{λ}}_{2} = 0.118

) are also statistically significant, indicating their importance in shaping the tail behavior and variability of wind speeds. The location parameters (

{\hat{δ}}_{1} = 10.99

{\hat{δ}}_{2} = 12.87

) reflect the physical limits of observed wind speeds and provide realistic bounds for extreme values. The estimated dependence parameter (

\hat{φ} = 0.42

) is precise and significant, indicating moderate nonlinear dependence between the paired wind speed measurements of the same day. This suggests that wind speeds observed at different times of the day exhibit some noticeable degree of persistence, while still retaining variability.

Compared to the dependence measure formulations proposed by Villanueva et al.

[19]

and Monahan

[20]

, with estimates presented in Table 5, the BWED’s presentation captures the nonlinear dependence characteristics of the two times of the day much better. Their dependence parameters (

\hat{ρ} = 0.93

and

\hat{ρ} = 0.39

, respectively) imply strong and relatively weak linear dependence.

Model comparison using the AIC clearly favours the BWED (AIC = 3468.56) over

{BWD}_{1}

(4257.98) and

{BWD}_{2}

(3886.72). The large AIC differences (

∆ AIC = 789.42

and

∆ AIC = 418.16

, respectively) provide strong evidence that the BWED offers a substantially better representation of the joint wind-speed behavior. This improved fit likely arises from its flexible survival-based dependence structure and the ability to capture nonlinear association.

Given its superior fit and realistic dependence characterization, the BWED is adopted for subsequent dependence visualization, including the copula-based plots shown in Figure 4 and drawing of inference from application of model to data.

Download: Download full-size image

Figure 4. Copula-based visualization of wind speed dependence at Christmas Island, Australia. Panel.

From Figure 4, panel (a) presents the copula-based joint density contours for wind speeds measured at 9 a.m. (x-axis) and 3 p.m. (y-axis) at Christmas Island. The contour orientation along the 45° diagonal axis signifies some degree of positive dependence between the two time periods. This indicates that higher morning wind speeds are likely associated with higher afternoon wind speeds, reflecting to some extent temporal persistence of wind regimes over the location. The dense yellow-to-orange region (4 - 6 m/s) marks the area of highest joint probability, implying that moderate wind speeds dominate the daily wind profile. The smooth gradient and absence of outlying contour distortion suggest that the BWE copula adequately captures the underlying dependence structure without overestimating tail co-movement. Panel (b) shows the empirical copula representation of the same dataset, with axes

u_{1} = F (x_{1})

and

u_{2} = F (x_{2})

denoting the pseudo-observations transformed to the unit square

{[0, 1]}^{2}

. The color intensity represents the frequency of paired observations, and the overlaid contour lines correspond to the BWE copula density in the copula domain. The clustering of points along the main diagonal confirms monotonic co-movement, and the contour overlay demonstrates a close alignment between the empirical dependence pattern and the theoretical BWE copula density; an indicator of good model calibration. The few sparse points near the corners (high-high and low-low regions) suggest mild tail dependence.

Finally, the good overlap (Panel a) between the empirical and theoretical contours provides a supporting visual evidence of the goodness-of-fit of the BWED. Also, Figure 4b is indicative of a consistent graphical evidence that the BWED offers an adequate representation of the co-dependence structure in the Australian wind speed data.

4.3.3. Dependence Measure

The copula-based measure of monotonic dependence used here is Kendall’s tau

({\hat{τ}}_{φ})

. Monotonic dependence refers to how two variables consistently move in the same relative direction (positive) or opposite direction (negative), but not necessarily at a constant rate. Unlike linear relationships, monotonic trends can curve, so long as they never reverse direction. This estimate is presented in Table 6 together with the classification of the strength of this measure as used in

[9]

Table 6. Estimate of co-dependence measure of BWED for Wind speed data.

Method	${\hat{τ}}_{φ}$	Strength
EML	0.583	Moderate

From Table 6, the value obtained indicates that there is a moderate positive monotonic dependence between the wind speeds at the two time points of the day. This implies that the tendency of one variable to increase with the other is visible but there is a significant variation in the data which is clearly evident in Figure 4b.

4.4. Discussion of Results

The new bivariate non-normal probability density function generated in this study, referred to as the BWED, describes the simultaneous occurrences of two random phenomena, each assumed to be characterized by the univariate Weibull-epsilon distribution

[41]

. The shapes of the bivariate density function, shown in Figures 1-3, at varying parameter values, depict different random simultaneous occurrences in nature, thus providing visual evidence of its flexibility to describe the intermittence in natural processes.

The simulation results confirm that the EML estimators yields consistent and efficient estimates for the parameters of the new bivariate distribution. Across sample sizes from 20 to 5,000, the estimates approach the true values monotonically, with both bias and standard error diminishing as

n

increases. This convergence behavior validates the asymptotic properties of the EML estimators and demonstrates its reliability even under moderate sample conditions typical of environmental data. The accurate recovery of the dependence parameter

φ

, which converged from 0.646 at

n = 20

to approximately 0.25 at

n = 5000

, underscores the method’s effectiveness in capturing dependence structures crucial for joint risk and reliability modelling. Similar improvements in convergence and variance reduction have been reported in recent MCMC-based estimation frameworks for copula and compound distributions applied to wind and hydrological data

[35, 36]

Comparative analysis with existing bivariate Weibull models further emphasizes the robustness of the BWED. The log-likelihood and AIC results show that the proposed BWED is consistently more amenable to data compared to the bivariate Weibull models in Villanueva et al.

[19]

and Monahan

[20]

, confirming a better flexibility in modelling both central and mild tail dependence. This advantage arises from the epsilon component, which adjusts the joint tail thickness without distorting marginal behavior, thereby offering a more adaptable dependence structure. Such flexibility has also been highlighted in recent studies where hybrid or copula-based Weibull models yielded improved fits in meteorological and renewable energy datasets

[37, 38]

. These findings suggest that the BWED can serve as a general framework for multivariate modelling where marginal Weibull behavior is evident but dependence deviates from conventional Gaussian or Gumbel copula assumptions.

The empirical application to wind speed data from Christmas Island further demonstrates the practical usefulness of the proposed BWED model. Graphical validation through contour density and empirical copula plots shows strong agreement between the theoretical BWED copula and the empirical pseudo-observations, indicating a well-calibrated model with little evidence of misfit. Similar copula-based diagnostics have shown that such approaches effectively capture nonlinear dependence structures in energy-related variables

[39]

. The marginal shape parameter estimates (

{\hat{α}}_{1} = 1.89, {\hat{α}}_{2} = 2.20

) indicate typical Weibull-type behavior with clear nonzero modes, consistent with observed wind speed distributions.

The dependence parameter (

\hat{φ} \approx 0.42

) is significant and reflects moderate nonlinear co-movement driven by shared atmospheric conditions. The estimate of the copula-based Kendall’s tau,

{\hat{τ}}_{φ} = 0.583

, suggests that this co-movement is a moderate positive monotonic dependence. Consequently, other meteorological and atmospheric factors (humidity, temperature, etc.) are influencing morning wind speed alongside that of the afternoon of the same day causing more noise in the data than a strong relationship.

There are practical insights provided by the nature of this dependence. These are:

1) Morning wind speed contain to some extent predictive information for corresponding afternoon conditions. This can improve the skill of short-term wind power forecasts and enable more reliable scheduling of maintenance or storage operations.

2) Wind power variability in the morning is likely to be mirrored later in the day, a feature critical for grid stability assessment.

3) The dependence relationship captured by the BWED copula allows realistic simulation of the joint wind profiles, aiding probabilistic forecasting and uncertainty quantification in renewable energy integration.

These insights align with recent developments emphasizing dependence-focused stochastic modelling in power systems and atmospheric sciences

[35, 36]

From a methodological standpoint, the results confirm that the BWED combines parsimony with flexibility, providing a tractable yet expressive model for wind-speed dependence. Its superior AIC performance over traditional bivariate Weibull formulations suggests that accounting for mild tail dependence significantly improves model fit. Furthermore, the copula representation inherent in the BWED framework allows its direct extension to higher dimensions, making it suitable for spatiotemporal modelling of wind data across multiple sites or time intervals. Given the increasing complexity of renewable energy systems, such scalable dependence models are becoming essential for integrated forecasting and risk management.

These findings are also timely given the rapid global expansion of wind power. According to recent reports by the Global Wind Energy Council

[40]

and the International Energy Agency

[31]

, global installed wind capacity surpassed 1,000 GW by 2023, with continued double-digit annual growth. This accelerating deployment heightens the need for accurate probabilistic models of wind dependence to support grid balancing, energy market forecasting, and investment risk analysis. As renewable penetration increases, the cost of misestimating joint variability also rises, making advanced dependence models like the BWE particularly valuable for operational planning and reliability studies.

Notwithstanding its strengths, several limitations suggest directions for further research. The present analysis assumes temporal stationarity and constant dependence; extending the BWED to time-varying or regime-switching frameworks could better capture evolving dependence structures. Although AIC-based selection establishes relative fit, out-of-sample validation using cross-validated forecast accuracy would provide stronger evidence of predictive performance. Future studies may incorporate meteorological covariates - such as atmospheric pressure, humidity, or stability indices - into the marginal or copula components and evaluate the BWED alongside alternative copula families to assess robustness. For extreme-event analysis, tail-adaptive or vine copula extensions could improve characterization of joint wind extremes across spatially correlated sites. Further work could also explore spatio-temporal formulations and applications in hydrology and reliability engineering, where modeling dependent extremes and component lifetimes is essential.

5. Conclusion

A new bivariate distribution, termed the BWED, has been developed and explored. Theoretical and graphical analyses, including density and contour plots under varying parameter settings, demonstrate its flexibility and wide applicability to bivariate lifetime systems. The model was also extended to a multivariate form, reinforcing its generality. Robustness studies indicate that the BWED is more amenable to data than the existing members of the bivariate Weibull family in benchmark applications

[19, 20]

. This provides strong motivation for its use. Empirical evaluation using paired wind-speed measurements at 9 a.m. and 3 p.m. of the same day from Christmas Island, Australia, further confirms its suitability. The BWED achieved the best fit relative to competing bivariate Weibull models, capturing the co-dependence between wind speeds at different times of the day. The copula-based Kendall’s τ estimate of 0.58 reveals a moderate monotonic dependence between the two time periods. The fact that a measure of co-dependence is derivable from the model gives credence to a comparative advantage in data application.

Generally, both simulation and empirical evidence affirm that the BWED is a robust and flexible framework for modelling dependence in environmental and reliability data. It provides efficient parameter estimation and realistic joint behavior, making it a valuable addition to stochastic wind-energy modelling and broader environmental applications.

Abbreviations

BWED	Bivariate Weibull Epsilon Distribution
BWE	Bivariate Weibull Epsilon
$BWD$	Bivariate Weibull distribution
EML	Exact Maximum Likelihood
AIC	Akaike Information Criterion
MCMC	Markov Chain Monte Carlo
Eq.	Equation

Author Contributions

Gongsin Isaac Esbond: Conceptualization, Formal Analysis, Methodology, Resources, Writing – original draft

Funmilayo Westnand Oshogboye Saporu: Supervision, Resources, Writing – review & editing

Conflicts of Interest

The authors declare no conflicts of interest.

Appendix: The Bivariate Weibull-epsilon Distribution and Its Multivariate Form

Appendix Ⅰ: Derivation of the Bivariate Weibull-epsilon Distribution

The Weibull-epsilon distribution was introduced

[41]

as a life time distribution with flexible shapes within its parameter range that can be used as density and hazard rate functions. For a random variable

X

, its density and distribution functions are, respectively, given by

f (x; α, δ, λ) = αλ \frac{δ^{2}}{δ^{2} - x^{2}} z (x) {[z (x) - 1]}^{α - 1} \exp \{- {[z (x) - 1]}^{α}\},

(A-1)

and

F (x; α, δ, λ) = 1 - \exp \{- {[z (x) - 1]}^{α}\},

(A-2)

where

z (x) = {(\frac{x + δ}{δ - x})}^{\frac{λδ}{2}}

0 < x < δ

α, λ, δ > 0

The bivariate Weibull-epsilon distribution is constructed using the joint survival function method defined in Eq. (1). The Weibull-epsilon cumulative failure rate function,

H (x_{i}), i = 1, 2

, is obtained as

H (x_{i}) = - \log \exp \{- {[{(\frac{x_{i} + δ_{i}}{δ_{i} - x_{i}})}^{\frac{λ_{i} δ_{i}}{2}} - 1]}^{α_{i}}\}

= {[{(\frac{x_{i} + δ_{i}}{δ_{i} - x_{i}})}^{\frac{λ_{i} δ_{i}}{2}} - 1]}^{α_{i}} .

(A-3)

For simplification of notation, let

z (x_{i}) = {(\frac{x_{i} + δ_{i}}{δ_{i} - x_{i}})}^{\frac{λ_{i} δ_{i}}{2}}, i = 1, 2

. Then the Weibull-epsilon joint survival rate function is given by

\overset{̅}{F} (x_{1}, x_{2}) = \exp \{- {({[z (x_{1}) - 1]}^{\frac{α_{1}}{φ}} + {[z (x_{2}) - 1]}^{\frac{α_{2}}{φ}})}^{φ}\} .

(A-4)

Differentiating

Eq . (A . 4)

with respect to

x_{1}

and

x_{2}

gives the bivariate Weibull-epsilon probability density function in Eq. (2). The corresponding cumulative distribution function is derived as

The bivariate Weibull-epsilon distribution function can be determined from the relationship between bivariate cumulative distribution and survival functions. This is given

[42]

\overset{̅}{F} (x_{1}, x_{2}) = 1 - F (x_{1}) - F (x_{2}) + F (x_{1}, x_{2})

\Rightarrow F (x_{1}, x_{2}) = \overset{̅}{F} (x_{1}, x_{2}) + F (x_{1}) + F (x_{2}) - 1 .

(A-5)

Substituting

(A . 2)

and

(A . 4)

into

(A . 5)

, the bivariate Weibull-epsilon distribution function is obtained, and given in Eq.(3).

Appendix Ⅱ: Validity of the BWED

Theorem

The bivariate Weibull-epsilon density function,

f (x_{1}, x_{2})

, given in Eq. (2) is a true probability density function.

Proof

It suffices to show that

0 \leq F (x_{1}, x_{2}) \leq 1

F (x_{1}, x_{2}) = P (X_{1} \leq x_{1}, X_{2} \leq x_{2})

= \int_{0}^{x_{2}} \int_{0}^{x_{1}} f (t_{1}, t_{2}) d t_{1} d t_{2}

= \exp \{- {({[z (x_{1}) - 1]}^{\frac{α_{1}}{φ}} + {[z (x_{2}) - 1]}^{\frac{α_{2}}{φ}})}^{φ}\} - \exp (- {[z (x_{1}) - 1]}^{α_{1}}) - \exp (- {[z (x_{2}) - 1]}^{α_{2}}) + 1 .

Since

z (0) = {(\frac{0 + δ_{i}}{δ_{i} - 0})}^{\frac{λ_{i} δ_{i}}{2}} = 1

and

z (δ_{i}) = {(\frac{δ_{i} + δ_{i}}{δ_{i} - δ_{i}})}^{\frac{λ_{i} δ_{i}}{2}} = \infty, i = 1, 2

, then

F (0, 0) = \exp \{- {({[1 - 1]}^{\frac{α_{1}}{φ}} + {[1 - 1]}^{\frac{α_{2}}{φ}})}^{φ}\} - \exp (- {[1 - 1]}^{α_{1}}) - \exp (- {[1 - 1]}^{α_{2}}) + 1

= \exp \{- 0\} - \exp (- 0) - \exp (- 0) + 1

= 1 - 1 - 1 + 1 = 0

Also

F (δ_{1}, δ_{2}) = \exp \{- {({[\infty - 1]}^{\frac{α_{1}}{φ}} + {[\infty - 1]}^{\frac{α_{2}}{φ}})}^{φ}\} - \exp (- {[\infty - 1]}^{α_{1}}) - \exp (- {[\infty - 1]}^{α_{2}}) + 1

= \exp \{- \infty\} - \exp (- \infty) - \exp (- \infty) + 1

= 0 - 0 - 0 + 1

= 1

Since

F (x_{1}, x_{2})

is a strictly monotonic increasing function, the proof is complete.

Appendix Ⅲ: Multivariate Weibull-epsilon distribution

The bivariate Weibull-epsilon distribution can be generalized to its multivariate form by the expansion of the expression of the bivariate survival function

[24]

. That is, given an

m -

dimemsional random vector

X = {[X_{1}, X_{2}, \dots, X_{m}]}^{T}

, the joint survival function is given by

{\overset{̅}{F}}_{X} (x) = \exp \{- {({[H_{1} (x_{1})]}^{\frac{1}{φ}} + {[H_{2} (x_{2})]}^{\frac{1}{φ}} + \dots + {[H_{m} (x_{m})]}^{\frac{1}{φ}})}^{φ}\} .

(A-6)

And for a Weibull-epsilon random vector,

X

, the joint survival function is given by

{\overset{̅}{F}}_{X} (x) = \exp \{- {({[z (x_{1}) - 1]}^{\frac{α_{1}}{φ}} + {[z (x_{2}) - 1]}^{\frac{α_{2}}{φ}} + \dots + {[z (x_{m}) - 1]}^{\frac{α_{m}}{φ}})}^{φ}\}

= \exp \{- {(\sum_{i = 1}^{m} {[z (x_{i}) - 1]}^{\frac{α_{i}}{φ}})}^{φ}\},

(A-7)

where

0 < φ \leq 1

0 < x_{i} < δ_{i}

x = {[x_{1}, x_{2}, \dots, x_{m}]}^{T}

and

α_{i}, λ_{i}, δ_{i} > 0, i = 1, 2, \dots, m

The multivariate Weibull-epsilon probability density function is derived from the joint survival function (

A . 7

) by differentiation

[24]

. That is,

f_{X} (x_{1}, x_{2}, \dots, x_{m}) = {(- 1)}^{m} \frac{\partial^{m}}{\partial x_{1} \partial x_{2} \dots \partial x_{m}} {\overset{̅}{F}}_{X} (x_{1}, x_{2}, \dots, x_{m})

= {(- 1)}^{m} \frac{\partial^{m}}{\partial x_{1} \partial x_{2} \dots \partial x_{m}} \exp \{- {(\sum_{i = 1}^{m} {[z (x_{i}) - 1]}^{\frac{α_{i}}{φ}})}^{φ}\}

= {(\frac{- 1}{φ})}^{m} \prod_{i = 1}^{m} \{α_{i} λ_{i} \frac{δ_{i}^{2}}{δ_{i}^{2} - x_{i}^{2}} z (x_{i}) {[z (x_{i}) - 1]}^{\frac{α_{i}}{φ} - 1}\}

\sum_{j = 1}^{N (m)} \{{(- 1)}^{n_{j}} N_{s} (m, j) (\prod_{k = 1}^{n_{j}} φ^{\underset{̲}{m_{k}}}) {(\sum_{i = 1}^{m} {[z (x_{i}) - 1]}^{\frac{α_{i}}{φ}})}^{n_{j} φ - m}\} \exp \{- {(\sum_{i = 1}^{m} {[z (x_{i}) - 1]}^{\frac{α_{i}}{φ}})}^{φ}\},

(A-8)

where

n_{j} =

number of summands of the

j^{th}

partition of

m

such that

m_{1} + m_{2} + \dots + m_{n_{j}} = m

m_{1} \geq m_{2} \geq \dots \geq m_{n_{j}} > 0

0 \leq n_{j} \leq m

φ^{\underset{̲}{m_{k}}} = φ (φ - 1) \times \dots \times (φ - m_{k} + 1)

is the falling factorial of

φ

N (m) =

total number of partitions of

m

, and

N_{s} (m, j) =

total number of set partitions of the set

S_{m} = \{1, 2, \dots, m\}

corresponding to the

j^{th}

partition of

m

References

[1]	C. Genest, O. Okhrin & T. Bodnard, “Copula modelling from Abe Sklar to the present day”, Journal of Multivariate Analysis 201 (2024) 105278, https://doi.org/10.1016/j.jmva.2023.105278
[2]	J. Größer & O. Okhrin (2021) “Copulae: An overview and recent developments”, WIREs Comput Stat. 2022; 14: e1557. https://doi.org/10.1002/wics.1557
[3]	P. Embrechts, A. McNeil & D. Straumann, Correlation and dependence in risk management: Properties and pitfalls. In M. Dempster (Ed.), Risk management: Value at risk and beyond, 1999), pp. 176-223. Cambridge University Press. https://doi.org/10.1017/CBO9780511615337.008
[4]	R. C. Jung, S. Kim & J. V. K. Rombouts, “Directional dependence in foreign exchange markets: Evidence from copula models”, Journal of Empirical Finance, (2008) 15(2), 287-298. https://doi.org/10.1016/j.jempfin.2007.05.002
[5]	M. Hofert & M. Scherer, “CDO pricing with nested Archimedean copulas”, Quantitative Finance, (2011) 11(5), 775-787. https://doi.org/10.1080/14697680903410073
[6]	C. de Michele, G. Salvadori, R. Vezzoli, & S. Pecora, Multivariate assessment of droughts: Copula-based analysis of severity and duration. Water Resources Research, 2013, 49(10), 6985 - 6994, https://doi.org/10.1002/wrcr.20548
[7]	S. Kim, R. C. Jung & J. Chung, “Dependence modelling in gene expression data with copulas”, Computational Statistics & Data Analysis (2009), 53(9), 3610-3620. https://doi.org/10.1016/j.csda.2009.02.020
[8]	X. Xie, Y. Ding, & J. Chen, “Copula-based wind speed dependence modelling for system reliability”, Renewable Energy, 2012, 40(1), 50-56. https://doi.org/10.1016/j.renene.2011.09.007
[9]	F. W. O. Saporu & I. E. Gongsin, “Modelling Dependence Relationships of Anthropometric Variables Using Copula Approach”, American Journal of Theoretical and Applied Statistics, 2020, Vol. 9, No. 5, pp. 245 - 255, https://doi.org:10.11648/j.ajtas.20200905.18
[10]	L. Deng, M. S. Smith & W. Maneesoonthon, “Large skew-t copula models and Asymmetric Dependence in intraday Equity Returns”, Journal of Business & Economic Statistics, 43(2) 269 - 285, (2025) https://doi.org/10.1080/07350015.2024.2360592
[11]	O. Crommen, P. Janssen & N. Veraverbeke, “Vine copula simplification: A comprehensive review”, SORT, 2025, 49(1), 67 - 102, https://raco.cat/index.php/SORT/article/view/433164
[12]	A. I. N. Miguel, “A Comprehensive Review of Copulas in Multivariate Dependency Modelling”, (September 09, 2024). Available at SSRN: https://dx.doi.org/10.2139/ssrn.4950425
[13]	Y. Wei, S. Yang, & Z. Sun, “Copula-based joint regression models for multivariate non-normal outcomes”, BMC Medical Research Methodology, 23(1) (2023), 185, https://doi.org/10.1186/s12874-023-01921-3
[14]	S. Kotz, N. Balakrishnan & N. L. Johnson, Continuous multivariate distributions, (2000). Volume 1: Models and applications (2nd ed.). Wiley. https://doi.org/10.1002/0471722065
[15]	J. M. Sarabia & G. Emilio, “Construction of multivariate distributions with given marginal and correlation matrix”, Journal of Multivariate Analysis, 99(6) (2008), 1262 - 1275. https://doi.org/10.1016/j.jmva.2007.07.005
[16]	N. Balakrishnan, & C. D. Lai, Continuous bivariate distributions (2nd ed.), 2009. Springer. https://doi.org/10.1007/b101765
[17]	I. E. Gongsin & F. W. O. Saporu, A bivariate conditional Weibull distribution with application, Afrika Matematika, 31:565-583, 2020. https://doi.org/10.1007/s13370-019-00742-8
[18]	S. Sens, R Lamichhane & N. Diawara, “A Bivariate Distribution with Conditional Gamma and its Multivariate Form”, Journal of Modern Applied Statistical Methods, vol. 13, No 2 Article 9, 169 - 184. 2014, https://doi.org/10.22237/jmasm/1414814880
[19]	D. Villanueva, A. Feijoo & J. L.Pazos, “Multivariate Weibull Distribution for Wind Speed and Wind Power Behavior Assessment”, Resources, 2013, 2, 370 - 384, https://doi.org:10.3390/resources2030370
[20]	A. H. Monahan, “Idealized models of the joint probability distribution of wind speeds”, Nonlin. Processes Geophys., 2018, 25, 335-353. https://doi.org/10.5194/npg-25-335-2018
[21]	J. Wu, X. Li & T. Zhou, “Copula-based joint modelling of wind and wave parameters for marine engineering” Journal of Marine Science and Engineering, 2025, 13(3), 396. https://doi.org/10.3390/jmse13030396
[22]	H. Zhang, Y. Liu, & K. Xu, “Modelling dependence of extreme wind speed and air density using Clayton copula”, Atmosphere, 2024, 15(12), 1437. https://doi.org/10.3390/atmos15121437
[23]	J. C. Lu & G. Bhattacharyya, “Some New Constructions of Bivariate Weibull Models” Annals of the Institute of Statistical Mathematics, 1990, 42, 543-559. https://www.ism.ac.jp/editsec/aism/pdf/042_3_0543.pdf
[24]	C. Lee & M. J. Wen, “A Multivariate Weibull Distribution”, Pakistan Journal of Statistics and Operations Research, 2009, 5(2), pp 55 - 66. https://doi.org/10.18187/pjsor.v5i2.120
[25]	A. Sklar, “Fonctions de Répartition à n Dimensions et Leurs Marges”, Publications de l'Institut de Statistique de L'Université de Paris, 1959, 8, 229-231.
[26]	R. B. Nelsen, An Introduction to Copulas. Springer Verlag, New York, 2006. http://dx.doi.org/10.1007/0387-28678-0
[27]	M. Abramowitz & I. A. Stegun, Handbook of Mathematical Functions, Dover, New York (1965)
[28]	GWEC. Global wind power boom continues despite economic woes. In: Global Wind Report 2019, pp. 8, (GWEC). www.gwec.net
[29]	IRENA. Future of wind: Deployment, investment, technology, grid integration and socio-economic aspects (A Global Energy Transformation paper). International Renewable Energy Agency (2019), Abu Dhabi. www.irena.org/publications
[30]	F. Zhao, “A record year for the wind industry, 2020”, In: Global Wind Report 2021, pp.6, GWEC.
[31]	IEA. Wind Power in 2023: Global Overview and Trends. International Energy Agency (2023), Paris.
[32]	X. Shen, C. Zhou & X. Fu, “Study of Time and Meteorological Characteristics of Wind Speed Correlation in Flat Terrains Based on Operation Data”, Energies, 2018, 11, 219. https://doi.org/10.3390/en11010219
[33]	S. Liu, G. Li, H. Xie & X. Wang, “Correlation Characteristic Analysis for Wind Speed in Different Geographical Hierarchies”, Energies, 2017, 10, 237. https://doi.org/10.3390/en10020237
[34]	A. A. Kadhem, N. I. Abdul-Wahab, I. Aris, J. Jasni & A. N. Abdalla, “Advanced Wind Speed Prediction Model Based on a Combination of Weibull Distribution and an Artificial Neural Network”, Energies (2017), 10, 1744. https://doi.org:10.3390/en10111744
[35]	S. Huang, Q. S. Li, Z. Shu & P. W. Chan, “Copula-based estimation of directional extreme wind speeds: Application for wind-resistant structural design”, Structures 60(4):105845, 2024. https://doi.org/10.1016/j.istruc.2023.105845
[36]	H. Qasem, N-E. Joergensen, A. Rahman, H. A. Samman, S. Al-Malki & A. S. Al-Ansari, “Vine Copula-Based Multivariate Distribution of Rainfall Intensity, Wind Speed, and Wind Direction for Optimizing Qatari Meteorological Stations”, Water, 16(9):1257, 2024, https://doi.org/10.3390/w16091257
[37]	W. Wang, F. Chen, Y. Li & L. Weng, “A hybrid WOA-KDE and mixed copula framework for directional wind assessment in complex terrain”, Sustainable Energy Technologies and Assessments, 83, 2025, 104590. https://doi.org/10.1016/j.seta.2025.104590
[38]	S. Huang, Q. Li, Z. Shu, & P. W. Chan, “Copula-based joint distribution analysis of wind speed and wind direction: Wind energy development for Hong Kong”, Wind Energy, 26(9), 2023, 900 - 922. https://doi.org/10.1002/we.2847
[39]	T. E. B. de Carvalho Silva, P. Maçaira, F. L. Cyrino Oliveira, & G. A. A. Pereira, “Enhancing energy planning with copula-based dependency models”, Energy Reports, Vol. 14, (2025), Pages 43-55. https://doi.org/10.1016/j.egyr.2025.05.028
[40]	GWEC. Global Wind Report 2024. Global Wind Energy Council (2024), Brussels.
[41]	I. E. Gongsin & F. W. O. Saporu, “The Weibull-epsilon distribution - its properties and applications”, International Journal of Mathematical and Computational Methods, volume 5, 2020, pp 80 - 88. http://www.iaras.org/iaras/journals/ijmcm
[42]	D. Kundu & R. D. Gupta, “Bivariate generalized exponential distribution”, Journal of Multivariate Analysis (2009), 100, 581 - 593. https://doi.org/10.1016/j.jmva.2008.06.012

Cite This Article

Plain Text BibTeX RIS

APA Style

Esbond, G. I., Saporu, F. W. O. (2026). A Bivariate Weibull-Epsilon Distribution with Application to Wind Speed Data. International Journal of Statistical Distributions and Applications, 12(3), 46-60. https://doi.org/10.11648/j.ijsda.20261203.11

Copy | Download

ACS Style

Esbond, G. I.; Saporu, F. W. O. A Bivariate Weibull-Epsilon Distribution with Application to Wind Speed Data. Int. J. Stat. Distrib. Appl. 2026, 12(3), 46-60. doi: 10.11648/j.ijsda.20261203.11

Copy | Download

AMA Style

Esbond GI, Saporu FWO. A Bivariate Weibull-Epsilon Distribution with Application to Wind Speed Data. Int J Stat Distrib Appl. 2026;12(3):46-60. doi: 10.11648/j.ijsda.20261203.11

Copy | Download

@article{10.11648/j.ijsda.20261203.11,
  author = {Gongsin Isaac Esbond and Funmilayo Westnand Oshogboye Saporu},
  title = {A Bivariate Weibull-Epsilon Distribution with Application to Wind Speed Data},
  journal = {International Journal of Statistical Distributions and Applications},
  volume = {12},
  number = {3},
  pages = {46-60},
  doi = {10.11648/j.ijsda.20261203.11},
  url = {https://doi.org/10.11648/j.ijsda.20261203.11},
  eprint = {https://article.sciencepublishinggroup.com/pdf/10.11648.j.ijsda.20261203.11},
  abstract = {Modelling dependence between non-normal environmental variables is important because conventional multivariate models may not adequately represent asymmetric distributions, nonlinear associations, and complex dependence structures commonly observed in environmental data. Wind speed, in particular, exhibits substantial temporal variability, making appropriate modelling of its dependence across different times of the day important for forecasting and renewable energy planning. This study develops a flexible Bivariate Weibull-Epsilon Distribution (BWED) for modelling dependence between two non-normal random variables. The study derives the joint probability density, cumulative distribution, conditional distributions, and copula representation of the proposed model, together with a copula-based measure of dependence. Model parameters are estimated using the Exact Maximum Likelihood (EML) method, while Gibbs sampling is employed to generate observations and evaluate the finite-sample performance of the estimators. Simulation experiments are conducted for sample sizes ranging from 20 to 5,000, the results show that parameter estimates become increasingly accurate and precise as sample size increases, with bias and standard errors decreasing substantially. A comparison with two existing bivariate Weibull models showed the BWED provides a superior fit, with an AIC of 3468.56 compared with 4257.98 and 3886.72 for the competing models. The proposed model is subsequently applied to 517 paired observations of wind speed recorded at 9:00 a.m. and 3:00 p.m. on the same day at Christmas Island, Australia, covering October 2018 to March 2020. The estimated Kendall’s tau of 0.583 indicates moderate positive monotonic dependence between morning and afternoon wind speeds. The study concludes that the BWED provides a flexible, robust, and tractable framework for modelling dependence in wind-speed and other environmental data, with potential applications in renewable energy forecasting and risk management.},
 year = {2026}
}

Copy | Download

TY  - JOUR
T1  - A Bivariate Weibull-Epsilon Distribution with Application to Wind Speed Data
AU  - Gongsin Isaac Esbond
AU  - Funmilayo Westnand Oshogboye Saporu
Y1  - 2026/09/18
PY  - 2026
N1  - https://doi.org/10.11648/j.ijsda.20261203.11
DO  - 10.11648/j.ijsda.20261203.11
T2  - International Journal of Statistical Distributions and Applications
JF  - International Journal of Statistical Distributions and Applications
JO  - International Journal of Statistical Distributions and Applications
SP  - 46
EP  - 60
PB  - Science Publishing Group
SN  - 2472-3509
UR  - https://doi.org/10.11648/j.ijsda.20261203.11
AB  - Modelling dependence between non-normal environmental variables is important because conventional multivariate models may not adequately represent asymmetric distributions, nonlinear associations, and complex dependence structures commonly observed in environmental data. Wind speed, in particular, exhibits substantial temporal variability, making appropriate modelling of its dependence across different times of the day important for forecasting and renewable energy planning. This study develops a flexible Bivariate Weibull-Epsilon Distribution (BWED) for modelling dependence between two non-normal random variables. The study derives the joint probability density, cumulative distribution, conditional distributions, and copula representation of the proposed model, together with a copula-based measure of dependence. Model parameters are estimated using the Exact Maximum Likelihood (EML) method, while Gibbs sampling is employed to generate observations and evaluate the finite-sample performance of the estimators. Simulation experiments are conducted for sample sizes ranging from 20 to 5,000, the results show that parameter estimates become increasingly accurate and precise as sample size increases, with bias and standard errors decreasing substantially. A comparison with two existing bivariate Weibull models showed the BWED provides a superior fit, with an AIC of 3468.56 compared with 4257.98 and 3886.72 for the competing models. The proposed model is subsequently applied to 517 paired observations of wind speed recorded at 9:00 a.m. and 3:00 p.m. on the same day at Christmas Island, Australia, covering October 2018 to March 2020. The estimated Kendall’s tau of 0.583 indicates moderate positive monotonic dependence between morning and afternoon wind speeds. The study concludes that the BWED provides a flexible, robust, and tractable framework for modelling dependence in wind-speed and other environmental data, with potential applications in renewable energy forecasting and risk management.
VL  - 12
IS  - 3
ER  -

Copy | Download

Author Information

Gongsin Isaac Esbond

Department of Statistics, University of Maiduguri, Maiduguri, Nigeria

Contact Email

http://orcid.org/0000-0003-2448-5543
Funmilayo Westnand Oshogboye Saporu

Department of Mathematics and Statistics, Confluence University of Science and Technology, Osara, Nigeria

Contact Email

Download PDF

Submit an Article

Abstract
Keywords
Document Sections
Show Full Outline
Abbreviations
Author Contributions
Conflicts of Interest
Appendix: The Bivariate Weibull-epsilon Distribution and Its Multivariate Form
References
Cite This Article
Author Information

Plain Text BibTeX RIS

APA Style

Esbond, G. I., Saporu, F. W. O. (2026). A Bivariate Weibull-Epsilon Distribution with Application to Wind Speed Data. International Journal of Statistical Distributions and Applications, 12(3), 46-60. https://doi.org/10.11648/j.ijsda.20261203.11

Copy | Download

ACS Style

Esbond, G. I.; Saporu, F. W. O. A Bivariate Weibull-Epsilon Distribution with Application to Wind Speed Data. Int. J. Stat. Distrib. Appl. 2026, 12(3), 46-60. doi: 10.11648/j.ijsda.20261203.11

Copy | Download

AMA Style

Esbond GI, Saporu FWO. A Bivariate Weibull-Epsilon Distribution with Application to Wind Speed Data. Int J Stat Distrib Appl. 2026;12(3):46-60. doi: 10.11648/j.ijsda.20261203.11

Copy | Download

@article{10.11648/j.ijsda.20261203.11,
  author = {Gongsin Isaac Esbond and Funmilayo Westnand Oshogboye Saporu},
  title = {A Bivariate Weibull-Epsilon Distribution with Application to Wind Speed Data},
  journal = {International Journal of Statistical Distributions and Applications},
  volume = {12},
  number = {3},
  pages = {46-60},
  doi = {10.11648/j.ijsda.20261203.11},
  url = {https://doi.org/10.11648/j.ijsda.20261203.11},
  eprint = {https://article.sciencepublishinggroup.com/pdf/10.11648.j.ijsda.20261203.11},
  abstract = {Modelling dependence between non-normal environmental variables is important because conventional multivariate models may not adequately represent asymmetric distributions, nonlinear associations, and complex dependence structures commonly observed in environmental data. Wind speed, in particular, exhibits substantial temporal variability, making appropriate modelling of its dependence across different times of the day important for forecasting and renewable energy planning. This study develops a flexible Bivariate Weibull-Epsilon Distribution (BWED) for modelling dependence between two non-normal random variables. The study derives the joint probability density, cumulative distribution, conditional distributions, and copula representation of the proposed model, together with a copula-based measure of dependence. Model parameters are estimated using the Exact Maximum Likelihood (EML) method, while Gibbs sampling is employed to generate observations and evaluate the finite-sample performance of the estimators. Simulation experiments are conducted for sample sizes ranging from 20 to 5,000, the results show that parameter estimates become increasingly accurate and precise as sample size increases, with bias and standard errors decreasing substantially. A comparison with two existing bivariate Weibull models showed the BWED provides a superior fit, with an AIC of 3468.56 compared with 4257.98 and 3886.72 for the competing models. The proposed model is subsequently applied to 517 paired observations of wind speed recorded at 9:00 a.m. and 3:00 p.m. on the same day at Christmas Island, Australia, covering October 2018 to March 2020. The estimated Kendall’s tau of 0.583 indicates moderate positive monotonic dependence between morning and afternoon wind speeds. The study concludes that the BWED provides a flexible, robust, and tractable framework for modelling dependence in wind-speed and other environmental data, with potential applications in renewable energy forecasting and risk management.},
 year = {2026}
}

Copy | Download

TY  - JOUR
T1  - A Bivariate Weibull-Epsilon Distribution with Application to Wind Speed Data
AU  - Gongsin Isaac Esbond
AU  - Funmilayo Westnand Oshogboye Saporu
Y1  - 2026/09/18
PY  - 2026
N1  - https://doi.org/10.11648/j.ijsda.20261203.11
DO  - 10.11648/j.ijsda.20261203.11
T2  - International Journal of Statistical Distributions and Applications
JF  - International Journal of Statistical Distributions and Applications
JO  - International Journal of Statistical Distributions and Applications
SP  - 46
EP  - 60
PB  - Science Publishing Group
SN  - 2472-3509
UR  - https://doi.org/10.11648/j.ijsda.20261203.11
AB  - Modelling dependence between non-normal environmental variables is important because conventional multivariate models may not adequately represent asymmetric distributions, nonlinear associations, and complex dependence structures commonly observed in environmental data. Wind speed, in particular, exhibits substantial temporal variability, making appropriate modelling of its dependence across different times of the day important for forecasting and renewable energy planning. This study develops a flexible Bivariate Weibull-Epsilon Distribution (BWED) for modelling dependence between two non-normal random variables. The study derives the joint probability density, cumulative distribution, conditional distributions, and copula representation of the proposed model, together with a copula-based measure of dependence. Model parameters are estimated using the Exact Maximum Likelihood (EML) method, while Gibbs sampling is employed to generate observations and evaluate the finite-sample performance of the estimators. Simulation experiments are conducted for sample sizes ranging from 20 to 5,000, the results show that parameter estimates become increasingly accurate and precise as sample size increases, with bias and standard errors decreasing substantially. A comparison with two existing bivariate Weibull models showed the BWED provides a superior fit, with an AIC of 3468.56 compared with 4257.98 and 3886.72 for the competing models. The proposed model is subsequently applied to 517 paired observations of wind speed recorded at 9:00 a.m. and 3:00 p.m. on the same day at Christmas Island, Australia, covering October 2018 to March 2020. The estimated Kendall’s tau of 0.583 indicates moderate positive monotonic dependence between morning and afternoon wind speeds. The study concludes that the BWED provides a flexible, robust, and tractable framework for modelling dependence in wind-speed and other environmental data, with potential applications in renewable energy forecasting and risk management.
VL  - 12
IS  - 3
ER  -

Copy | Download