A new versatile modification of the Rayleigh distribution for modeling COVID-19 mortality rates.

The aim of this paper is to specify a new flexible statistical model to analyze COVID-19 mortality rates in Italy and Canada. A new versatile lifetime distribution with four parameters is proposed by using the exponentiated generalized class of distributions and the gull alpha power Rayleigh distribution to form the exponentiated generalized gull alpha power Rayleigh (EGGAPR) distribution. This new distribution is characterized by a tractable cumulative distribution function. To estimate the unknown parameters of the proposed distribution the maximum likelihood estimation method is used. In evaluating the effectiveness of the MLE method graphical displays of the Monte Carlo simulation are presented. The EGGAPR distribution is compared to its sub-models which include the exponentiated gull alpha Rayleigh distribution, the gull alpha Rayleigh distribution, exponentiated generalized Rayleigh distribution, exponentiated Rayleigh distribution and the Rayleigh distribution. Different measures of goodness-of-fit are used to investigate whether the EGGAPR distribution is more flexible and fit than its sub-models in modeling COVID-19 mortality rates.

Introduction

The main objective of statistics is to be able to develop models that can be used to model life phenomenon. Probability distributions are at the heart of modeling the phenomenon that are characterized by uncertainties especially in biomedical sciences. Much effort has been devoted in developing new statistical models. Recently, the approach of extending a family of distribution by adding extra parameters has gained the attention of researchers. This technique has benefits in that it leads to flexible models that can capture all aspects of real life data.

The Rayleigh distribution is best suited choice for analyzing data characterized by increasing hazard rates. However in real life situations data is characterized by monotonic hazard rates such as decreasing and non-monotonic failure rates such as bathtub or modified bathtub hazard rates. Due to this shortcoming of the Rayleigh distribution, many extensions of the Rayleigh distribution has been developed in literature. The generalized Rayleigh distribution was proposed by [1] where they studied different estimation procedures for the unknown parameters and compared the performances by the use of Monte Carlo simulation. The type I half logistic Rayleigh distribution by [2] and studied statistical properties. The Rayleigh–Rayleigh distribution by [3], the slashed Rayleigh distribution by [4], the slashed generalized distribution by [5], the Rayleigh–weibull and the Rayleigh–inverted-weibull distribution by [6], the kumaraswamy generalized distribution by [7], the Gompertz-exponentiated Rayleigh distribution by [8], the inverse lomax-Rayleigh distribution by [9], the extended odd weibull Rayleigh distribution by [10].

Many distributions have been proposed to model COVID-19 data. For further reading on this the reader is referred to [10], [11], [12], [13], [14], [15], [16], [17], [18] among others.

A new model with four parameters called the EGGAPR is introduced in this paper. The EGGAPR distribution is obtained by following the following methodology:

Let and be the PDF and CDF of a distribution with parameter vector respectively. Recently, [19] introduced the Gull Alpha Power Family of distributions with cumulative distribution function given by:

denotes the shape parameter and is . The PDF of the CDF in Eq. (1) is: To be able to develop the new generator of continuous distributions, the study used the methodology that was developed by [20]. Let be a CDF of a given model and , [20] proposed the exponentiated generalized class of distributions which has two extra shape parameters and with the CDF defined as: and are two extra shape parameters. Of importance, the function in Eq. (3) is not complicated and has tractable properties especially in simulations. In addition, it is easy to obtain the inverse of the CDF of Eq. (3). The purpose of the two extra shape parameters is that they can control both tail weights which allows more flexible distributions to be generated which have heavier tails. The probability density function (PDF) associated with the equation in (3) is: Eq. (4) has many applications because it offers greater flexibility on its tails. The Cumulative Density Function and the probability density function of the Exponentiated-Generalized gull alpha power family of distribution will be obtained by using and given in (1), (2) as the baseline CDF and PDF in Eqs. (3), (4).

The resulting new distribution is called the Exponentiated-Generalized Gull Alpha Power Family(EGGAPF) of distributions and its CDF is given by:

The PDF associated with Eq. (5) is given by: is the baseline survival function, is the parameters of the baseline distribution and is the baseline CDF.

The primary objective of this paper is to propose and study a new modification of the Rayleigh distribution using it as the baseline distribution in the EGGAPF family. The proposed four parameter distribution may be more flexible compared to its sub-models and it would accommodate data sets that have a monotonic and a non-monotonic failure rates.

The remainder of the paper is presented as described, the proposed distribution is presented in Section “EGGAPR distribution”. Statistical properties of the proposed distribution are given in Section “Mathematical Properties”, the estimation of the parameters and results for Monte Carlo simulation is explained in Section “Estimation of Parameters”, Section “Application to COVID-19 data” presents the application to real life data and comparison to competing models. The concluding remarks are presented in Section “Conclusion”.

EGGAPR distribution

The four parameter EGGAPR distribution is a special case of the Exponentiated Generalized Gull Alpha Power Family of distribution with the Rayleigh distribution serving as the baseline distribution. The CDF of the Rayleigh distribution with scale parameter is given by: and the corresponding PDF is

By using the CDF and the PDF of the Rayleigh distribution as the baseline in Eqs. (5), (6) respectively we obtain the EGGAP-Rayleigh distribution which has CDF given as: where is the scale parameter, are the shape parameters. The corresponding PDF is given by: The PDF of EGGAPR distribution can be decreasing, unimodal, inverted bathtub, left-skewed, almost symmetrical or right-skewed. Fig. 1 gives the different shapes of the EGGAPR density function.

Fig. 1

EGGAPR density function for selected parameter values.

The survival distribution of this distribution: and the hazard function is given by: where

As observed in Fig. 2, the hazard plot for the Rayleigh distribution shows the flexibility of the distribution. The various shapes exhibited for the various parameter values include the J-shape, bathtub shape, decreasing, decreasing–increasing–decreasing shape. This demonstrates the flexibility of the distribution in modeling data with monotone and non-monotone hazard shape.

Fig. 2

EGGAPR hazard function shapes.

The EGGAPR distribution consists of sub-models that can be widely used in modeling real life data. The include:

The Rayleigh distribution When the EGGAPR reduces to the Rayleigh distribution with CDF given by

The Gull Alpha Power Rayleigh distribution When the EGGAPR distribution is reduced to Gull Alpha Power Rayleigh distribution with CDF:

Exponentiated Gull Alpha Power Rayleigh distribution When , EGGAPR reduces to exponential Gull Alpha Power Rayleigh whose CDF is:

Exponentiated Rayleigh distribution When the EGGAPR will reduce to the exponentiated Rayleigh distribution with CDF:

Exponentiated generalized Rayleigh distribution When the EGGAPR will reduce to the exponentiated generalized Rayleigh distribution with CDF:

A summary of the various sub-models of the EGGAPW distribution is displayed in Table 1

Table 1

Summary of submodels from the EGGAPR distribution.

Distribution	σ	α	a	b
Rayleigh	–	1	1	1
GAPR	–	–	1	1
EGAPR	–	–	1	–
ER	–	1	1	–
EGR	–	1	–	–

Mathematical properties

Here some mathematical properties of the EGGAPR distribution such as quantiles, moments, moment generating function, renyi entropy and order statistics are derived.

Linear representation

The EGGAPR can be expressed in form of the mixture representation which is useful in the derivation of the mathematical properties of the distribution. This is important since it is usually difficult or even impossible to obtain explicit analytical solutions for certain mathematical quantities.

The EGGAPR distribution PDF ca be written in mixture form as:

For any real integer , a series representation of for —z— ¡ 1 is given by

By applying Eq. (14) in Eq. (10) twice we have:

Substituting Eq. (15) into Eq. (9) yields

Which can be rewritten as:

Upon further simplification Eq. (17) can be written in terms of the exponentiated gull alpha power Rayleigh distribution as: where is the PDF of the Exponentiated Gull Alpha Power Rayleigh distribution with parameters

Quantile function

The quantile function is useful in describing the distribution of the random variable. It is important in the simulation of random numbers and also in the computation of the median, kurtosis and the skewness.

EGGAPR distribution quantile function for is given by:

By definition the quantile function always returns the value such that:

Thus

Letting in Eq. (20) and solving for gives:

By substituting , gives the first quartile, median, third quartile of the EGGAPR distribution respectively. In order to simulate random variables from the EGGAPR distribution the following relation is used:

Moments

For a newly proposed distribution, it is always vital to derive the moments. They are important in deriving measures of shape, dispersion and central tendency measures.

The th non-central moment of the EGAPR distribution is given by: where ,

By definition

Replacing the EGAPR PDF into Eq. (22) we get:

Expanding Eq. (23) and utilizing the following representations and

Thus the th moment can be formed as follows:

Let

Using the binomial expansion thrice we have: which can be written as: where

Let

Utilizing the expansions and substituting we have:

Which can be rewritten as: where

Combining Eqs. (25), (27) completes the proof.

It can be observed that Eq. (21) has no closed form hence the moments are obtained through numerical integration. The following parameters were used for the computation. , , and . As observed in Table 2 it can be observed that the EGGAPR distribution is numerically versatile in means and variance. In addition, the values of the CS reveal that the distribution can be right skewed, left skewed and almost symmetrical. For the values of CK, the EGGAPR distribution can be almost mesokurtic, platykurtic and leptokurtic.

Table 2

Summary statistics of moments of the EGGAPR distribution.

r	I	II	III	IV	V
μ1′	0.9320	1.9493	1.8913	2.1964	0.4102
μ2′	1.1308	5.0975	4.7793	5.4815	0.4071
μ3′	1.1617	15.5039	14.6953	14.9212	0.5504
μ4′	2.6046	52.3216	52.4291	43.4588	0.8926
μ5′	4.6114	191.4388	210.7946	133.9179	1.6440
SD	0.5120	1.3919	1.0965	0.8108	0.4887
CV	0.5494	0.5844	0.5798	0.3692	1.1914
CS	0.5525	0.3435	0.8406	−0.0117	1.6066
CK	3.0175	2.5748	3.7671	2.8066	5.5296

The MGF of a random variable that follows the EGGAPR distribution is given by the proposition:

The MGF of the random variableis given by:whereandare defined in Eq. (21)

Using the identity given by: then the MGF of the EGGAPR distribution is given by: replacing with Eq. (21) completes the proof.

Skewness and Kurtosis

To demonstrate the effects of the additional shape parameters on the measures of shape skewness and kurtosis, the Bowley skewness coefficient based on the quartiles and the Moors Kurtosis based on the octiles are calculated. The measures are robust and less sensitive to outliers and even exists for distribution with no moments.

The Bowley Skewness is given by: and the Moors Kurtosis Coefficient is defined by:

The plot for Bowley’s skewness is given in Fig. 3 for fixed baseline parameter values and . As observed the additional parameters have an effect on the shapes of the skewness coefficient. This enhances the flexibility of the EGGAPR distribution and reinforces the importance of the additional parameters.

Fig. 3

Plots of Bowley skewness for EGGAPR distribution for fixed baseline parameter values of and .

The plot for Moors kurtosis is given in Fig. 4 for fixed baseline parameter values and . As observed the additional parameters have an effect on the shapes of the kurtosis coefficient. This enhances the flexibility of the EGGAPR distribution and reinforces the importance of the additional parameters in controlling for kurtosis.

Fig. 4

Plots of Moors Kurtosis for EGGAPR distribution for fixed baseline parameter values of and .

Entropy

The Renyi Entropy of the EGGAPR distribution is derived in this section.

For a random variable having the EGGAPR distribution the Renyi Entropy is given by where ,

The Renyi entropy is defined by:

Applying the technique for expanding the moments we have:

Hence, which completes the proof. The values of the entropy of the EGGAPR are obtained through numerical integration and displayed in Table 3 gives the Renyi Entropy for different parameter values .

Table 3

Values of Renyi Entropy for different parameter values.

δ	I	II	III
0.1	1.4976	0.7346	1.6265
0.5	0.9227	0.3069	1.0020
0.9	0.7551	0.2039	0.8284
1.5	0.6322	0.1321	0.7047
2.5	0.5307	0.0721	0.6042

Estimation of parameters

The method of maximum likelihood estimation is used to derive the maximum likelihood estimators for the unknown parameters of the EGGAPR distribution. Let be a random sample from the EGGAPR distribution with unknown parameter vector . Then the likelihood function is given by:

Replacing Eq. (10) into Eq. (33) and taking the logarithm we get the log-likelihood function given by:

The maximum likelihood estimates for the parameters are the values which maximize the log-likelihood function given in Eq. (34). Therefore, by setting the first partial derivatives of the log likelihood function in Eq. (34) with respect to to be zero are as follows:

Let and such that The partial derivatives of with respect to each parameter and equating each to zero is given by:

Explicit forms of the Eqs. (35), (36), (37), (38) are not available and numerical solutions will be provided in practice.

In order to estimate the parameters of the EGGAPR distribution using the BFGS algorithm both the gradient vector of the log likelihood function and the Hessian matrix are required. The Hessian matrix is defined as a matrix of which its elements involve the second derivatives of the log likelihood function with respect to the parameters.

The observed information matrix is given by: where . The expressions for the observed information can be availed upon further request.

Monte Carlo Simulation

Simulation study was conducted to evaluate the average bias, the root mean square errors and the mean of the parameter estimates. Different sample sizes and different sets of the parameter values were used in the simulation study. The quantile function given in Eq. (19) was used in simulation. The algorithm used is as follows:

A random sample of sizes are generated from the EGGAPR distribution using selected initial values for parameters and .

For each sample size the estimates are obtained

For each sample size Steps(1–2) are repeated 2000 times and the MLE’s of the parameters, their average biases and the root mean square errors are recorded.

The Average Biases (AB) and the Root Mean Squared Error (RMSE) were calculated using the following Equations: and

Figs. 5 and 7 and 6 displays the Means of the parameter estimates, the RMSE and AB for the maximum likelihood estimators of for increasing sample sizes. As observed, the parameter values converge to the true parameter value as the sample size increases. For the RMSE and AB they decrease as the sample size increases. Parameter is sensitive to changes in the sample size.

Fig. 5

Plots of the estimated MLEs against sample size.

Fig. 7

RMSE for estimators.

Fig. 6

Average bias for estimators.

Figs. 8 and 9 and 10 and displays the Means of the parameter estimates, the RMSE and AB for the maximum likelihood estimators of for increasing sample sizes. As observed, the parameter values converge to the true parameter value as the sample size increases. For the RMSE and AB they decrease as the sample size increases. Parameter is sensitive to changes in the sample size.

Fig. 8

Estimated MLEs against sample size.

Fig. 9

RMSE for estimators.

Fig. 10

Average bias for estimators.

Concluding remarks on simulation results

The values of the parameter estimates approach the true value as sample size increase.

The RMSE of the parameters decrease with increase in sample size.

The AB of the parameter estimates decrease with increase in sample size.

All the estimators perform well and they provide small values of RMSE and the mean value of the estimators tend to the initial value of the parameters.

Application to COVID-19 data

Data set I: Italy COVID-19 mortality rates data

The first data set represents COVID-19 mortality rates Data belonging to Italy for a period of 59 days recorded from 27 February to 27 April 2020. The data was obtained from https://covid19.who.int/. The same data set also has been used by [10]. The data is given in Table 4

Table 4

Italy COVID-19 mortality rates data.

4.571	7.201	3.606	8.479	11.410	8.961	10.919	10.908	6.503	18.474
11.010	17.337	16.561	13.226	15.137	8.697	15.787	13.333	11.822	14.242
11.273	14.330	16.046	8.646	8.905	8.906	7.407	7.445	7.214	6.194
4.640	5.542	5.073	4.416	4.859	4.408	4.639	3.148	4.040	4.253
4.011	3.564	3.827	3.134	2.780	2.881	3.341	2.686	2.814	2.508
2.450	1.518	11.950	10.282	11.775	10.644	10.138	9.037	12.396

The descriptive statistics for the data are displayed in Table 5. The data is right skewed because of the positive sign of the skewness coefficient. The kurtosis value is less than three implying that the data is platykurtic.

Table 5

Summary statistics for Italy COVID-19 mortality rates data.

Statistic	Min	Max	Mean	Std.dev	Median	Kurtosis	Skewness
Value	1.518	18.474	8.156	4.527	7.445	−0.943	0.441

The data set have an increasing failure rate as displayed by the TTT transform plot which has a concave shape as depicted in Fig. 11. This reinforces the fact that the shape of the hazard function displayed in Fig. 2 is fit for modeling this type of data.

Fig. 11

TTT-transform plot for Italy COVID-19 mortality rates data.

The Maximum Likelihood Estimates of the parameters of the proposed distribution EGGAPR and its sub models are presented in Table 6

Table 6

MLEs and SE (in parenthesis) for Italy COVID-19 mortality rates data.

Model	aˆ	bˆ	αˆ	σˆ
EGGAPR	0.1397(0.0220)	1.3639(0.2415)	2.497 (0.052)	2.480(0.0022)
EGAPR	–	1.0073(0.2287)	1.218(0.586)	6.910(0.8792)
GAPR	–	–	1.206(0.445)	6.908(0.8728)
EGR	6.388(28.571)	0.9478(0.1598)	–	16.929(70.39)
ER	–	1.030(48.250)	–	6.681(33.868)
R	–	–	–	6.583(0.4285)

The EGGAPR distribution provides a better fit than its sub models. As indicated in Table 7, the EGGAPR distribution has the highest log-likelihood and the smallest values of K–S and W* compared to the other models. Considering the formal tests of goodness of fit tests, in order to verify which distributions better fits the jet airplane data, since the EGGAPR distribution has the lowest values for the K–S, Anderson–Darling and W* we then conclude that the distribution provides a better fit than the sub models. The Exponentiated Generalized Rayleigh Distribution, Exponentiated Rayleigh and the Rayleigh distributions do not fit the data as shown by the small p-values.

Table 7

Log-likelihood, information criteria and goodness-of-fit statistics for Italy COVID-19 mortality rates data.

Model	l	AIC	A*	K–S	p-value	W*
EGGAPR(Proposed)	−165.480	338.961	0.560	0.0892	0.7019	0.0952
EGAPR(New)	−167.636	341.708	0.1346	0.1222	0.3156	0.1346
GAPR(New)	−167.637	339.279	0.8039	0.1215	0.3225	0.1344
EGR	−167.710	341.421	98.787	0.9999	<0.05	12.530
ER	−167.766	339.533	0.7993	0.2931	<0.05	0.1328
R	−167.767	337.533	0.8044	0.2931	<0.05	0.1328

The plots of the densities of the fitted distributions are shown in Fig. 12.

Fig. 12

Fitted densities plot for Italy COVID-19 mortality rates.

Data set II: Canada COVID-19 mortality rate data

The second data set consists of Mortality rates data belonging to Canada for a period of 56 days ranging from 1 November to 26 December 2020. The data sets formed a drought mortality rates. The data are 0.1622 0.1159 0.1897 0.1260 0.3025 0.2190 0.2075 0.2241 0.2163 0.1262 0.1627 0.2591 0.1989 0.3053 0.2170 0.2241 0.2174 0.2541 0.1997 0.3333 0.2594 0.2230 0.2290 0.1536 0.2024 0.2931 0.2739 0.2607 0.2736 0.2323 0.1563 0.2677 0.2181 0.3019 0.2136 0.2281 0.2346 0.1888 0.2729 0.2162 0.2746 0.2936 0.3259 0.2242 0.1810 0.2679 0.2296 0.2992 0.2464 0.2576 0.2338 0.1499 0.2075 0.1834 0.3347 0.2362.

The descriptive statistics for the data are displayed in Table 8. The data is right skewed because of the positive sign of the skewness coefficient. The kurtosis value approximately three implying that the data is mesokurtic.

Table 8

Summary statistics for Canada COVID-19 mortality rates data.

Statistic	Min	Max	Mean	Std.dev	Median	Kurtosis	Skewness
Value	1.515	6.868	3.281	0.998	3.177	2.815	1.163

The data has an increasing failure rate as depicted in the TTT plot in Fig. 13

Fig. 13

TTT-transform plot for Canada COVID-19 mortality rates data.

The Maximum Likelihood Estimates of the parameters of the proposed distribution EGGAPR and its sub models are presented in Table 9

Table 9

Maximum likelihood estimates of parameters and standard errors (in parenthesis) for Canada COVID-19 mortality rates data.

Model	aˆ	bˆ	αˆ	σˆ
EGGAPR	3.003(3.166)	1.894(4.838)	5.175 (2.624)	2.173(2.1026)
EGAPR	–	2.2833(0.336)	5.600(2.047)	2.560(0.256)
GAPR	–	–	1.577(0.111)	0.0095(0.0143)
EGR	3.683(20.23)	4.853(53.30)	–	4.015(1.185)
ER	–	1.155(14.759)	–	0.227(5.807)
R	–	–	–	2.422(0.202)

The EGGAPR distribution provides a better fit than its sub models. As indicated in Table 10, the EGGAPR distribution has the highest log-likelihood and the smallest values of K–S and W* compared to the other models. Considering the formal tests of goodness of fit tests, in order to verify which distributions better fits the jet airplane data, since the EGGAPR distribution has the lowest values for the K–S, Anderson–Darling and W* we then conclude that the distribution provides a better fit than the sub models. The Exponentiated Generalized Rayleigh Distribution, Exponentiated Rayleigh and the Rayleigh distributions do not fit the data as shown by the small p-values.

Table 10

Log-likelihood and goodness-of-fit statistics for Italy COVID-19 mortality rates data.

Model	l	A*	K–S	p-value	W*
EGGAPR(Proposed)	−47.926	0.5279	0.1048	0.8234	0.0919
EGAPR(New)	−47.845	0.5231	0.1092	0.7837	0.0892
GAPR(New)	−48.609	0.5533	0.1300	0.5769	0.0975
EGR	−48.380	70.916	0.999	<0.05	10.825
ER	−58.466	0.9330	0.8827	<0.05	0.1632
R	−58.466	0.8044	0.2802	0.007	0.1186

The plots of the densities of the fitted distributions are shown in Fig. 14.

Fig. 14

Fitted densities plot for Canada COVID-19 mortality rates.

Concluding remarks on real data

Concerning data set I we can conclude that EGGAPR provides the highest p value and the lowest W*, A and K–S compared to its sub-models.

From Fig. 12 we can observe that EGGAPR was the best fitting model for the data set.

Concerning data set II we can conclude that EGGAPR provides the highest p value and the lowest W*, A and K–S compared to its sub-models.

From Fig. 14 we can conclude that EGGAPR was the best fitting model for data set II.

From Table 7, EGR, ER and R provides a poor fitting for the data set I.

From Table 10, EGR, ER and R provides a poor fitting for the data set II.

From both applications we can conclude that EGGAPR distribution provides the best fitting among all its competitive distributions which gives it a superiority in fitting this kind of COVID-19 mortality rate data.

Conclusion

The Rayleigh distribution has been widely used in the practice of statistical sciences especially in reliability, engineering and financial sciences. In this study we suggested a new four parameter model called the EGGAPR distribution which can be denoted as EGGAPR distribution. The survival function, hazard function, linear representation, quantile function, moments, moment generating function and entropy are given. The maximum likelihood estimators of the EGGAPR parameters are obtained and a Monte Carlo simulation conducted. The Average Bias, the Root Mean Square Error and the Average of the parameter estimates were calculated. Generally, as the sample size increases the maximum likelihood estimates converge to the true value, the Average Bias and the Root Mean Square Error decrease with increase in sample size. In the sense of statistics, we have implemented the EGGAPR distribution for COVID-19 data. The data sets showed that EGGAPR provides the best fit for the data sets compared with the competitive sub-models.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.