Approximating reproduction numbers: a general numerical method for age-structured models

1.   Introduction

Classical probability distribution models play a pivotal role in statistics, providing essential theoretical tools for data modeling, parameter estimation, hypothesis testing, and statistical inference. They facilitate a comprehensive understanding and analysis of random variables and uncertain phenomena. However, when dealing with bounded random variables within the range (0, 1), such as proportions, probabilities, and percentages, classical probability distribution models encounter limitations in accurately describing and predicting data. To address this issue, we introduce the concept of unit distribution to generate novel models that enhance the flexibility of existing approaches for better adaptation to real-world data. The unit distribution is typically derived by redefining traditional analytical methods or transforming classical continuous distributions. By avoiding the introduction of new parameters while effectively translating the flexibility inherent to classical distributions into unit intervals, it offers improved efficiency compared to baseline methods. In recent years, research on unit distribution has made significant progress across various fields. For instance, Alvarez et al. ^[1] investigated the correlation properties and statistical characteristics of a novel distribution in the unit interval derived from transforming random variables with a semi-normal distribution. They employed maximum likelihood (ML) estimation to simulate correlation statistics and validate the superiority of this new distribution over other distributions defined within the unit interval. Mazucheli et al. ^[2] proposed a unit-Gompertz distribution model based on the unit interval, which can replace existing models such as unit-Birnbaum-Saunders, unit-Weibull, L-Logistic, Kumaraswamy, and Johnson SB distributions. The authors reparametrized the unit-Gompertz model to incorporate covariate effects across the response distribution. Okorie et al. ^[3] examined the regression model for the upper truncated Weibull distribution in the unit interval and extended it to include 0-1 inflation. They estimated parameters using ML estimation and demonstrated improved fitting performance of their proposed distribution using real data analysis. Shakhatreh et al. ^[4] explored Bayesian estimation (BE) for the logarithmic logical distribution in the unit interval by considering non-informative priors for parameter estimation and employing Monte Carlo simulations to evaluate Bayesian estimates' performance.

The compound Rayleigh distribution (CRD), being a significant probability distribution model in statistics, finds extensive applications in reliability testing, survival analysis, communication systems, and various other fields. It represents a novel type of distribution derived by fixing one parameter of the three-parameter Burr-XII distributions. Hence, studying the CRD holds immense significance. Currently, considerable progress has been made in statistical inference research on the CRD. For instance, Shao et al. ^[5] investigated the BE of CRD parameters using progressively type II censored data. Wang et al. ^[6] obtained ML estimates of CRD parameters based on complete samples and employed pivot parameter construction to derive inverse moment estimates for shape and scale parameters while utilizing Monte Carlo simulation for obtaining corresponding statistics. Badr ^[7] utilized the modified Kolmogorov-Smirnov test, Cramer-Von Mises test, and Anderson-Darling test to assess the goodness of fit of the CRD model under complete samples as well as type II censored samples. Barot and Patel ^[8] analyzed the BE of CRD based on progressively type II censored data.

The probability density function (PDF) and cumulative distribution function (CDF) of the CRD are respectively represented as follows:

Here, $\beta$ is the shape parameter, and $\theta$ is the scale parameter.

In order to address the limitations of the conventional distribution model in describing bounded random variables within the interval (0, 1), this study proposes a novel probability distribution model known as the unit compound Rayleigh distribution (UCRD). The objective of this research is to analyze the mathematical and statistical properties of UCRD and explore its potential applications. This paper initially derives fundamental mathematical properties of UCRD, including PDF, CDF, survival function (SF), and hazard function (HF). Furthermore, it analyzes quantile function, moments, and order statistics of UCRD. Additionally, detailed expressions for five entropy measures under UCRD are derived. These entropy measures effectively quantify uncertainty in random variables and hold significant importance for data analysis and model selection. To accurately estimate parameters in the UCRD model, ML estimation and BE are employed due to their wide adaptability range, solid theoretical foundation, excellent properties, and practical value in statistics. Through Monte Carlo simulation experiments, this study evaluates the performance of different entropy measures based on ML estimation to verify their validity and reliability. Although some progress has been made regarding unit distribution models in existing literature, the research on CRD within a unit interval remains relatively scarce. The UCRD model proposed in this study possesses the advantages of flexibility, excellent fitting performance, easy implementation, and a solid theoretical foundation, thereby rendering it highly promising for extensive applications across various fields. In comparison with other models, UCRD exhibits superior capability in describing and predicting random variables within the interval (0, 1), thus not only enriching the theory of related models in probability theory and statistics but also providing a potent analytical tool for researchers and practitioners in relevant domains. Moreover, the UCRD model holds great potential for application in reliability analysis, survival analysis, and communication systems, among others. This study aims to introduce novel research ideas and methodologies to these fields to foster their further advancement.

The rest of this article is described below. In Section 2, we introduce the PDF, CDF, survival function, and HF of the UCRD. In Section 3, we derive the quantile function, moment-generating function, K-order incomplete moment, and other mathematical properties of UCRD. In Section 4, five entropy expressions of UCRD are derived. In Section 5, two methods of ML estimation and BE are used to obtain the parameter estimates of UCRD. In Section 6, the parameters of the UCRD model are estimated by Monte Carlo simulation combined with the estimation method used, and the estimated values, average bias (AB), corresponding mean square error (MSE), and mean relative estimate (MRE) of the parameters are obtained. Based on ML estimation, the average entropy estimate (AEE), average entropy bias (AEB), MSE, and MRE of each entropy are obtained to evaluate the performance of each entropy. In Section 7, two sets of real data are presented to test the validity of the UCRD model and various entropy measures. Section 8 provides relevant conclusions and future prospects.

2.   Related models of UCRD

Assuming T follows a CRD distribution with parameters $\beta$ and $\theta$, the new UCRD can be obtained by transforming T to $T = X/(1 - X)$ using Eqs (1) and (2), resulting in the following PDF and CDF, respectively:

When one parameter is held constant, the PDF plot of UCRD exhibits distinct variations as the value of another parameter undergoes transformation. Specifically, when the $\beta$ value remains fixed, altering the $\theta$ value leads to an increasing trend within the range of (0, 0.5), followed by a decreasing trend within the range of (0.5, 1). Notably, this PDF plot demonstrates unimodality and overall symmetry throughout, as depicted in Figure 1.

The SF and HF are crucial analytical tools in the field of reliability analysis and survival analysis, utilized for examining the system's reliability, lifetime distribution, and event occurrence probability. Analyzing SF and HF enables us to assess system reliability and conduct risk assessment, among other applications. The expressions for SF and HF of UCRD are as follows:

As the $\beta$ and $\theta$ values change, it is evident from Figure 2 that the HF curve of UCRD always shows a clear upward trend, thus demonstrating the remarkable flexibility of UCRD.

The cumulative HF and the reverse HF are two important concepts in the fields of reliability engineering and survival analysis. They provide assessments and predictions of the probability of system failure and failure times, which aid in optimizing system design, maintenance, and decision-making processes. Below, we present the CDF and the reverse HF for the UCRD, respectively:

3.   Mathematical properties of UCRD

In this section, we present various statistical properties of UCRD, including the quantile function, K-order moment, mean value, variance, etc. Through analyzing these statistics, we can gain a deeper understanding of the characteristics and properties of the dataset and evaluate the performance of the UCRD model. These statistical properties have significant applications in data analysis, modeling, and inference fields; they provide us with robust tools and guidance for interpreting and utilizing data more effectively.

3.1.   Quantile function of UCRD

In the field of statistics and probability theory, the quantile function is extensively utilized to assess the distribution pattern, central tendency, and dispersion of data. It serves as a crucial tool in characterizing positional information within a probability distribution, enabling us to determine the values taken by a random variable at a given probability. The quantile function plays an integral role in data analysis by facilitating our understanding of dataset distribution patterns, identifying outliers, and executing statistical inference and prediction operations. By analyzing the quantile function, we can extract vital insights regarding data location and distribution. The derivation process for obtaining the quantile function pertaining to UCRD involves employing Eq (4) as outlined below:

where w represents the probability value, which is within the range of (0, 1). The first quantile, median, and third quantile represent the quantiles corresponding to $w = \frac{1}{4}, \frac{1}{2}, \frac{3}{4}$, respectively.

Skewness and kurtosis play a crucial role in data analysis, providing valuable insights into the shape of the data distribution. They are utilized to assess the degree of asymmetry, peakedness, and deviation from theoretical distributions. Moreover, they serve as indicators for testing the normality assumption of the data distribution and selecting appropriate statistical models to gain a comprehensive understanding of its distributional characteristics. Let CSk represent skewness and CKu represent kurtosis in the UCRD.

3.2.   Moments and moment-generating function

Assuming the random variable X is distributed with UCRD, the k-th moment of X can be expressed as:

Let $z = \frac{x}{{1 - x}}$, then $x = \frac{z}{{1 + z}} = 1 - \frac{1}{{1 + z}}.$ Hence ${x^k} = {(1 - \frac{1}{{1 + z}})^k} = \sum\limits_{i = 0}^k {{{(- 1)}^i}\left({\begin{array}{*{20}{c}} k \\ i \end{array}} \right)} {(1 + z)^{ - i}}.$

Thus

Here, ${\lambda _i}(z; \beta, \theta) = \int_0^\infty {z{{(1 + z)}^{ - i}}(1 + \frac{{{z^2}}}{\theta }} {)^{ - (\beta + 1)}}dz.$ As can be seen from Eq (11), the mean and variance of X are:

Therefore, the coefficient of variation (CV) for Z is $CV = \frac{\sigma }{{{\mu _1}}}$. The expression for the moment-generating function of X is as indicated below:

Table 1 presents the measured values of different statistical properties of the UCRD for various parameter values, and Figure 3 provides a 3D image of the statistical properties of the UCRD.

3.3.   Incomplete moments

The computation of moments for a random variable in the field of statistics and probability theory considers only a subset of possible values, rather than encompassing all potential values. This is particularly relevant when dealing with random variables that have an extensive or infinite range. However, practical applications often face limitations in data collection or computational complexity, making it impractical to calculate or estimate moments of all orders. Consequently, incomplete moments are sometimes considered for random variables to simplify the problem and reduce computational complexity. The k-th incomplete moment of the random variable X can be defined using the following formula:

Therefore, the k-th incomplete moment of X is:

Here, ${\eta _i}(z; \beta, \theta) = \int_0^{\frac{x}{{1 - x}}} {z{{(1 + z)}^{ - i}}(1 + \frac{{{z^2}}}{\theta }} {)^{ - (\beta + 1)}}dz.$

3.4.   Lorentz curve

The Lorenz curve is a graphical tool commonly used in economics and statistics to describe the degree of inequality in the distribution of income or wealth. Proposed by the American economist Lorenz in 1905, it is widely used in the fields of economics, sociology, and policy analysis to study and compare the income or wealth distribution of different groups, countries, or regions, and to assess the severity of inequality. The Lorenz curve can provide quantitative information and intuitive graphical presentations to help researchers and policymakers better understand and deal with inequality. The Lorenz curve for X with UCRD is given below:

3.5.   Bonferroni curve

The Bonferroni curve is also a relative index reflecting income inequality, and it is a slight modification of the Lorenz curve. Its definition is the ratio of the Lorenz curve to the CDF, that is:

3.6.   Order statistics

Let $({X_1}, {X_2}, \cdots, {X_n})$ be a random sample drawn from the population X and denote $({x_1}, {x_2}, \cdots, {x_n})$ as the observed values of the sample. When the sample is arranged in ascending order from smallest to largest, we obtain $({X_{(1)}}, {X_{(2)}}, \cdots, {X_{(n)}})$, and $({x_{(1)}}, {x_{(2)}}, \cdots, {x_{(n)}})$ are the observed values of $({X_{(1)}}, {X_{(2)}}, \cdots, {X_{(n)}})$. These are referred to as the order statistics of the sample $({X_1}, {X_2}, \cdots, {X_n})$, where ${X_{(k)}}$ is the k-th order statistic of the random sample. Consequently, the PDF for the k-th order statistic ${X_{(k)}}$ of the UCRD is given by:

The CDF for the k-th order statistic ${X_{(k)}}$ of the UCRD is given by:

However, ${X_{(n)}} = \max \{ {x_1}, {x_2}, \cdots, {x_n}\}$, thus the PDF of ${X_{(n)}}$ is:

Due to ${X_{(1)}} = \min \{ {x_1}, {x_2}, \cdots, {x_n}\}$, the PDF of ${X_{(1)}}$ is:

4.   Expressions for entropy measures

Entropy measure, as a fundamental concept in information theory, is utilized to quantify the uncertainty or informational content of random variables. It serves as a metric that characterizes the probability distribution of a random variable and is employed to depict the level of uncertainty associated with said distribution. A higher entropy measure indicates greater uncertainty and thus more information contained within the random variable. Conversely, a lower entropy measure signifies reduced uncertainty and consequently less information present. This notion finds extensive applications in diverse fields such as information theory, data analysis, pattern recognition, signal processing, and machine learning. Currently, numerous researchers have made notable advancements in exploring entropy measures. For instance, Kashyap et al. ^[9] proposed an innovative entropy measure specifically tailored for Pythagorean fuzzy sets by axiomatically defining it and introducing key properties related to it. Sayyari and Barsam ^[10] investigated novel practical inequalities for extended entropy and relative entropy across various parameters. Abd El-Raouf and AbaOud ^[11] analyzed different entropy indices of unit generalized Rayleigh distributions while employing Monte Carlo simulation techniques to estimate their performance under this particular model.

4.1.   Shannon entropy

Shannon entropy, also known as information entropy, is extensively utilized in various fields including communication, data compression, channel capacity, and cryptography to quantify the average amount of information, data compression rate, and security of cryptographic systems. Moreover, Shannon entropy plays a crucial role in statistics, machine learning, and data analysis where it finds applications in tasks such as feature selection, cluster analysis, and pattern recognition. Wang et al. ^[12] comprehensively investigated the impact of temperature and quantum confidence on the Shannon entropy of hydrogen impurity states in gallium arsenide quantum wells. Nascimento et al. ^[13] employed Shannon entropy to explore electrons confined in a double quantum dot system by mapping changes in spatial entropy parameters as an indicator of decoupling/coupling extent between twin quantum dots. Piga et al. ^[14] proposed a fast semi-analytic estimator for sparse sampling distributions using a hierarchical Bayesian approach to estimate Shannon entropy. Formentin et al. ^[15] discussed the concept of entropy along with degree-based indices for characterizing the chemical structure of iron chloride while utilizing these indices to calculate Shannon's entropy. Below, we present the expression for Shannon entropy under the UCRD:

4.2.   Rényi entropy

Rényi entropy, proposed by Rényi in 1960 as a generalization of Shannon entropy, measures the uncertainty and diversity of random variables ^[16]. When the parameter $\alpha = 1$, Rényi entropy transforms into the form of Shannon entropy. Recent progress has been made on Rényi entropy. Zhong ^[17] investigated the replication technique and homogenization mapping for calculating the heat of Rényi entropy in individual intervals on cylinder calibration. Baharith ^[18] discussed the basic statistical properties of the Chweeble inverse Gompertz distribution, including the reliability function, moments, Rényi entropy, and order statistics. Below is the expression for Rényi entropy under UCRD:

4.3.   Havrda-Charvat entropy

Havrda-Charvat entropy, also known as Tsallis entropy in the context of non-extensive thermodynamics, is a generalized concept that expands on Shannon entropy. It was introduced by Havrda and Charvat in 1967 and has found applications in various fields such as pattern recognition, image processing, and data clustering. Compared to Shannon entropy, Havrda-Charvat entropy offers a more flexible and comprehensive set of measures that effectively describe complex systems with non-standard probability distributions. Shi et al. ^[19] proposed a novel complexity measurement method called weighted Havrda-Charvat entropy based on permutation patterns and the Havrda-Charvat entropy itself to differentiate uncertainty among time series possessing identical pattern orders. Brochet et al. ^[20] conducted a quantitative comparison of loss functions by implementing parametric Tsallis-Havrda-Charvat entropy alongside conventional Shannon entropy for training deep learning networks in medical imaging tasks constrained by limited data volumes. Wang and Shang ^[21] utilized the Havrda-Charvat entropy plane to analyze complexity features within time series data while Brochet et al. ^[22] designed a classifier using convolutional neural networks with a novel loss function based on Havrda-Charvat entropy for analyzing various types of data. Below, we present the expression for Havrda-Charvat entropy under the UCRD:

4.4.   Arimoto entropy

The Arimoto entropy is a measure in information theory that quantifies the correlation or mutual information between random variables. It was introduced by the Japanese scientist Arimoto in 1972 and is based on relative entropy, which measures the discrepancy between probability distributions. The Arimoto entropy has wide applications in fields such as information theory, pattern recognition, machine learning, and communications. It can be utilized for tasks including feature selection, cluster analysis, information encoding, and channel capacity to comprehend and quantify relationships and information transfer between random variables. Li et al. ^[23] proposed a novel objective non-reference metric for evaluating image fusion by leveraging the properties of Arimoto entropy in their calculations. Additionally, Li et al. ^[24] introduced a new method for non-rigid registration of medical images that employs the Arimoto entropy of gradient distribution as an information-theoretic metric. Almarashi ^[25] and Abd El-Raouf and AbaOud ^[26] conducted research on the statistical properties of Rényi entropy, Shannon entropy, Havrda-Charvat entropy, and Arimoto entropy under different distributional settings. Below, we present the expression for Arimoto entropy under the UCRD:

4.5.   Mathai-Haubold entropy

Almanjahie et al. ^[27] studied the fundamental properties of Mathai-Haubold entropy through order statistics and evaluated the performance of magnitude-based Mathai-Haubold entropy via simulation. Additionally, Asgharzadeh et al. ^[28] analyzed various properties of a generalization of the Lindley distribution and derived estimates of the Rényi entropy and Mathai-Haubold entropy for this distribution. Below, we present the expression for Mathai-Haubold entropy under the UCRD:

5.   Estimation of UCRD model parameters

5.1.   ML estimation

Let $({x_1}, {x_2}, \cdots, {x_n})$ be an observation of a random sample $({X_1}, {X_2}, \cdots, {X_n})$ drawn from the UCRD. The corresponding likelihood function is:

The log-likelihood function is:

By taking the partial derivatives of the log-likelihood function with respect to the parameters $\beta$ and $\theta$, we obtain the likelihood Eqs (29) and (30):

By solving the likelihood Eqs (29) and (30), we can obtain the ML estimates for $\beta$ and $\theta$. However, these equations cannot be directly solved analytically, so we consider using numerical methods such as the Newton-Raphson iteration or dichotomy method to obtain the ML estimates for $\beta$ and $\theta$. Due to the invariance of ML estimation, the parameter estimates $\hat \beta$, $\hat \theta$ obtained through numerical methods can be substituted back into Eqs (22)−(26) to obtain the ML estimates for the corresponding entropies ${\hat H_S}, \, {\hat H_R}, \, {\hat H_H}, \, {\hat H_A}, \, {\hat H_M}$.

5.2.   BE

In this section, we use BE to estimate the parameters of the URED model under the squared error loss function. BE, as a traditional estimation method, is based on Bayes theorem to estimate the probability of unknown parameters in the case of given observation data. It has a better property than ML estimation, using prior knowledge and sample size to get more accurate parameter estimation. Assuming that $\beta$ and $\theta$ are independent random variables and follow $\Gamma ({\sigma _1}, {\omega _1})$ and $\Gamma ({\sigma _2}, {\omega _2})$, respectively, the density functions of $\beta$ and $\theta$ are:

Then, the joint prior distribution of $\beta$ and $\theta$ is:

According to Bayes' theorem, the joint posterior density of $\beta$ and $\theta$ is:

In order to quantify the cost loss resulting from errors in the estimation process, it is customary to introduce a loss function aimed at minimizing the expected loss. In this study, we explore the utilization of both square error loss functions (SELF) for parameter estimation in UCRD models. The SELF is formally defined as ^[29]:

Where $\hat \phi$ represents the ML estimate for the parameter $\phi$.

Thus, the Bayesian estimator under the SELF is:

The non-explicit form of Eq (35) makes direct calculation more intricate, leading us to consider employing the Lindley approximation algorithm for parameter estimation. The definition of Lindley's approximation is obtained from ^[30]:

Where $\phi (\beta, \theta)$ represents the parameter vector, $h(\beta, \theta) = \ln \pi (\beta, \theta)$, if the sample size is sufficiently large, Eq (36) can be further simplified to the following expression:

Where $\hat \beta$ and $\hat \theta$ are ML estimates of $\beta$ and $\theta$, and the subscripts represent partial derivatives of variables, such as ${\phi _\beta }$ representing the first derivative of $\beta$ in $\phi (\beta, \theta)$; similarly, the others are similar representations, ${\lambda _{i, j}}$ representing the (i, j) element of ${[- \frac{{{\partial ^2}l(\left. {\beta, \theta } \right|x)}}{{\partial \beta \partial \theta }}]^{ - 1}}$, (i, j = 1, 2). Then, there is:

When we calculate the estimate of $\beta$ under SELF, where $\phi (\beta, \theta) = \beta$, we have

Bring the Eqs (38)−(44) into the Eq (37) to get the following formula:

When we calculate the estimate of $\theta$ under SELF, where $\phi (\beta, \theta) = \theta$, we have

${\phi _\theta } = 1, {\phi _\beta } = {\phi _{\beta \beta }} = {\phi _{\beta \theta }} = {\phi _{\theta \beta }} = {\phi _{\theta \theta }} = 0.$ (45)

Bring the Eqs (38)−(43) and Eq (45) into the Eq (37) to get the following formula:

6.   Simulation study

In this section, we utilize the Monte Carlo simulation method in conjunction with the estimation technique employed in this study and MATLAB software to compute the mean square error (MSE), mean bias (AB), and mean relative estimate (MRE) of the parameters of the UCRD model. Initially, we set $\beta = 0.3, \theta = 0.1$ as the initial parameter values, while ${\sigma _1} = {\omega _1} = {\sigma _2} = {\omega _2} = 1$ serve as hyperparameters. We consider sample sizes of n = 30, 50, 80, and 100 for which we conduct 1000 repeated tests per sample size; the main code can be found in Appendix A. The specific results are presented in Table 2.

Based on ML estimation, we estimate different entropy measures under UCRD to evaluate and compare their performance. First, we set the initial values of the parameters as $\beta = 0.35, \, 0.75, \, 1.25, 1.65$, $\theta = 0.5, \, 1, \, 1.5, \, 2$, and the entropy parameter $\alpha = 0.5, 1.5$, and obtain the estimates of different entropy measures under UCRD (see Tables 3 and 4). Then, we further set the initial values of the parameters to $\beta = 0.5, \theta = 0.3$ and $\beta = 0.75, \theta = 0.4$, and set the entropy parameters $\alpha = 0.8, 1.2, 1.5, 1.7$. For the sample size n = 30, 50, 80, 100, we conducted 1000 repetitions. Based on these experiments, we used MATLAB software to calculate the AEEs, AEBs, and corresponding MSEs and MREs of the parameters and five entropy measures under UCRD (see Tables 5−8). The main code is in Appendix B.

Based on the above tables, we can draw the following conclusions:

(1) The ML estimation exhibits superior parameter accuracy compared to the BE when considering the SELF. Moreover, as the sample size increases, the MSEs, AEBs, and MREs of parameter estimation all demonstrate a decreasing trend, indicating an improvement in estimation accuracy with larger data volumes.

(2) When $\theta$ is fixed, different entropy measures exhibit a decrease with increasing $\beta$. Similarly, when $\beta$ is fixed, various entropy measures tend to decrease as $\theta$ increases.

(3) With fixed parameters $\beta$ and $\theta$, Rényi entropy, Havrda-Charvat entropy, and Arimoto entropy measures tend to decrease as the entropy parameters $\alpha$ increase. Conversely, Mathai-Haubold entropy measures tend to increase with higher values of the entropy parameters $\alpha$.

(4) As the sample size increases, the AEs of parameters and different entropies gradually approach the initial values, and the MSEs, AEBs, and MREs of parameters and entropies decrease.

(5) Increasing the entropy parameter $\alpha$ generally leads to higher MSEs and AEBs for entropies.

7.   Real data analysis

In this section, we employ two sets of actual data to examine the versatility and applicability of the UCRD across the unit interval in comparison with other classical probability distribution models. Additionally, we investigate the feasibility of utilizing ML estimation for estimating various entropy measures. The initial dataset was provided by Bjerkedal ^[31] and pertains to the survival time of 72 guinea pigs infected with virulent tuberculosis bacillus. The second dataset consists of infection times (in months) among dialysis patients as proposed by Klein and Moeschberger ^[32], with specific details presented in Table 8. It is important to note that we need to normalize dataset 2 within the range (0, 1), resulting in converted values as follows: 0.08333, 0.08333, 0.11667, 0.11667, 0.11667, 0.15000, 0.18333, 0.21667, 0.21667, 0.25000, 0.25000, 0.25000, 0.25000, 0.28333, 0.31667, 0.35000, 0.38333, 0.41667, 0.41667, 0.45000, 0.48333, 0.48333, 0.71667, 0.71667, 0.75000, 0.75000, 0.85000, 0.91667.

In this study, to comprehensively evaluate the goodness of fit of the UCRD model to both datasets, we selected various probability distribution models for analysis. These included the Kumaraswamy (Kw) distribution, CRD, generalized Rayleigh distribution (GRD), inverse exponential Rayleigh distribution (IERD), and two-parameter Rayleigh distribution (RD), which were compared with the UCRD model. Using MATLAB software, we applied Akaiike information criteria (AIC), Kolmogorov-Smirnov (KS) statistics, and Anderson-Darling (AD) statistics to these models. The P-value of KS statistic was adopted as the criterion for selecting the best-fitting model; the main code can be found in Appendix C and D, and corresponding numerical results are presented in Table 10. Additionally, empirical distributions of both datasets were plotted and visually compared with CDF diagrams corresponding to each model. These comparison results are illustrated in Figure 4. It can be observed from Table 10 and Figure 4 that the UCRD exhibits superior fitting performance for both datasets. Based on this finding, we further analyze the proposed statistics using these two datasets; relevant analysis results are shown in Table 11. In Table 11, it is evident that metric values of Rényi entropy, Havrda-Charvat entropy, and Arimoto entropy decrease with increasing entropy parameters, while metric values of Mathai-Haubold entropy increase with increasing entropy parameters. This conclusion aligns with the findings obtained in the simulation section and further validates the effectiveness and reliability of our proposed method for estimating entropy measures.

8.   Conclusions

In this paper, we construct a novel UCRD model, which is defined on the bounded interval (0, 1). First, the basic characteristics of the UCRD model are briefly summarized. Then, the core statistical properties such as quantile function, K-moment, expectation, and variance of the distribution are derived in detail, and the application potential and adaptability of the model in practice are deeply discussed. In addition, five entropy measures under the framework of the model are systematically analyzed, and ML estimation and BE are used to estimate the model parameters. Through Monte Carlo simulation, the AEEs, AEBs, MSEs, and MREs of these entropy measures are further calculated. Finally, two sets of real data are used to verify the performance of the UCRD model, and the results show that compared with other traditional distributions, the UCRD model presents a better fitting effect. In future work, we plan to further expand the UCRD model to enhance its adaptability to different types of data. We are considering combining the UCRD model with other models to improve the predictive power of the model. At the same time, we will explore the development of the UCRD model in other subject areas, promote knowledge exchange between disciplines, and solve interdisciplinary problems. In addition, we will deepen theoretical research and model validation to enhance the theoretical support and practical application value of the UCRD model.

Author contributions

Qin Gong: Writing-original draft, Method, Software, Writing-review & editing; Lijun Luo: Methods, Writing-original draft, Writing-review & editing; Haiping Ren: Methods, Writing-review & editing. All authors have read and approved the final version of the manuscript for publication.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This research was funded by National Natural Science Foundation of China, grant number 71661012 and the Foundation of Jxust (No.XJG-2020-27).

Conflict of interest

The authors declare no conflict of interest.

Appendix

A.Estimating the parameters of the UCRD model using maximum likelihood estimation and Bayesian estimation

clear, clc;

n = 30;%nIndicates the total number of samples

beta0 = 0.3;theta0 = 0.1;%Set the initial value of parameters

a1 = 1;b1 = 1;a2 = 1;b2 = 1;%Set the value of hyperparameterss

for i = 1:1000

G = rand(1, m);

for j = 1:m

H(j) = G(j).^(1/(j+sum(R(m-j+1:m))));

end

for j = 1:m

Z(j) = 1-prod(H(m-j+1:m));

x(j) = 1-(1+(theta0.*(1-Z(j)).^(-1./beta0)-theta0).^(1/2)).^(-1);%Generate random numbers that follow the UCRD distribution

end

thetaML(i) = ER(n, x);

betaML(i) = n./sum(log(1+(1./thetaML(i)).*(x./(1-x)).^2)); %Obtaining maximum likelihood estimates of parameters from the dichotomy method

[beta_BE(i), theta_BE(i)] = Lindley(n, x, a1, b1, a2, b2, betaML(i), thetaML(i));

end

theta = mean(thetaML); beta = mean(betaML); %Maximum likelihood estimation mean of parameters

theta_B = mean(theta_BE); beta_B = mean(beta_BE); %Bayesian estimation mean of parameters

Attached call function:

(1) Dichotomy method

function thetaML = ER(n, x)

f = @(n, x, theta)-n./theta+((n./sum(log(1+(1./theta).*(x./(1-x)).^2)))+1).*sum(((x./(1-x)).^2)./(theta.^2+theta.*(x./(1-x)).^2)); %Merge of likelihood functions

theta_lower = 0;%The lower bound of theta

theta_upper = 5;%The upper bound of theta

tolerance = 1e-5;%Error tolerance

while theta_upper-theta_lower > tolerance

theta_mid = (theta_lower+theta_upper)/2;%The midpoint of theta

f_mid = f(n, x, theta_mid);

if f_mid < 0

theta_upper = theta_mid;

theta_mid = (theta_lower+theta_upper)/2;

else

theta_lower = theta_mid;

theta_mid = (theta_lower+theta_upper)/2;

end

end

thetaML = theta_mid;

end

(2) Lindly approximation

function [beta_BE, theta_BE] = Lindley(n, x, a1, b1, a2, b2, betaML, thetaML)%Lindley approximation for Bayesian estimation

L11 = @(a, b)-n./(a.^2);

L22 = @(a, b)n./(b.^2)-(a+1).*sum(((x./1-x).^2.*(2.*b+(x./1-x).^2))./(b.^2+b.*(x./1-x).^2).^2);

L12 = @(a, b)sum(((x./1-x).^2)./(b.^2+b.*(x./1-x).^2));

L21 = L12;

L111 = @(a, b)2.*n./(a.^3);

L121 = 0;

L211 = L121;

L221 = @(a, b)-sum(((x./1-x).^2.*(2.*b+(x./1-x).^2))./(b.^2+b.*(x./1-x).^2).^2);

L122 = L221;

L212 = L221;

L222 = @(a, b)(-2.*n)./(b.^3)-2.*(a+1).*sum(((x./1-x).^2)./(b.^2+b.*(x./1-x).^2).^2)+2.*(a+1).*sum(((x./1-x).^2.*(2.*b+(x./1-x).^2).^2)./(b.^2+b.*(x./1-x).^2).^3);%On the differentiation of parameters of logarithmic likelihood function

P1 = @(a, b)((a1-1)./a)-b1;

P2 = @(a, b)((a2-1)./b)-b2;%The logarithmic derivative of the logarithmic joint prior distribution of parameters

%Next, we will apply maximum likelihood estimates to each function

L11 = L11(betaML, thetaML); L12 = L12(betaML, thetaML); L21 = L21(betaML, thetaML); L22 = L22(betaML, thetaML);

L111 = L111(betaML, thetaML); L122 = L122(betaML, thetaML); L212 = L122;L221 = L122;L222 = L222(betaML, thetaML);

P1 = P1(betaML, thetaML); P2 = P2(betaML, thetaML);

O = inv([-L11, -L12;-L21, -L22]); O11 = O(1, 1);O12 = O(1, 2);O21 = O(2, 1);O22 = O(2, 2);%Using O to represent Fisher inverse matrix

%Next, calculate the Bayesian estimation of beta under the squared error loss function

U1 = betaML;

U1_1 = 1;

U1_2 = 0;U1_11 = 0;U1_22 = 0;U1_12 = 0;U1_21 = 0;

%Substitute maximum likelihood estimates into U1

beta_BE = U1+0.5*((U1_11+2*U1_1*P1)*O11+(U1_21+2*U1_2*P1)*O21+(U1_12+2*U1_1*P2)*O12+(U1_22+2*U1_2*P2)*O22)+0.5*((U1_1*O11+U1_2*O12)...

*(L111*O11+L221*O22)+(U1_1*O21+U1_2*O22)*(L122*O12+L212*O21+L222*O22));

%Next, calculate the Bayesian estimation of theta under the squared error loss function

U2 = thetaML;

U2_1 = 1;

U2_2 = 0;U2_11 = 0;U2_22 = 0;U2_12 = 0;U2_21 = 0;

%Substitute maximum likelihood estimation into U2

theta_BE = U2+0.5*((U2_11+2*U2_1*P1)*O11+(U2_21+2*U2_2*P1)*O21+(U2_12+2*U2_1*P2)*O12+(U2_22+2*U2_2*P2)*O22)+0.5*((U2_1*O11+U2_2*O12)*...

(L111*O11+L221*O22)+(U2_1*O21+U2_2*O22)*(L122*O12+L212*O21+L222*O22));

end

B.Calculate the average estimate of parameters and entropy, the average deviation, as well as the corresponding mean square error and average relative estimate

clear, clc;

n = 30;%Indicates the total number of samples

beta0 = 0.5;theta0 = 0.3;%Set the initial value of parameters

alpha = 0.8;%Set entropy parameter values

HR = @(a, b)(1./(1-alpha)).*log(integral(@(t)((2*a)./b).^alpha.*(t./(1-t).^3).^alpha.*(1+(1./b).*(t./(1-t)).^2).^(-alpha.*(a+1)), 0, 1));

%Expressions for Renyi entropy

HR0 = HR(beta0, theta0);

%The initial values of Renyi entropy

for i = 1:1000

G = rand(1, m);

for j = 1:m

H(j) = G(j).^(1/(j+sum(R(m-j+1:m))));

end

for j = 1:m

Z(j) = 1-prod(H(m-j+1:m));

x(j) = 1-(1+(theta0.*(1-Z(j)).^(-1./beta0)-theta0).^(1/2)).^(-1);

end

thetaML(i) = ER(n, x);

betaML(i) = n./sum(log(1+(1./thetaML(i)).*(x./(1-x)).^2)); %Estimating parameters through dichotomy method

HR_ML(i) = HR(betaML(i), thetaML(i));

%Obtain estimates of Renyi entropy

end

theta_AEE = mean(thetaML); beta_AEE = mean(betaML); %Mean of parameters

HR_AEE = mean(HR_ML);

%The mean of Renyi entropy

theta_MSE = sum((thetaML-theta0).^2)./1000;

beta_MSE = sum((betaML-beta0).^2)./1000;%Mean squared error of parameters

HR_MSE = sum((HR_ML-HR0).^2)./1000;

%Mean squared error of Renyi entropy

theta_AEB = sum(abs(thetaML-theta0))./1000;

beta_AEB = sum(abs(betaML-beta0))./1000;%The average bias of parameters

HR_AEB = sum(abs(HR_ML-HR0))./1000;

%The average bias of Renyi entropy

theta_MRE = sum((thetaML-theta0)./theta0)./1000;

beta_MRE = sum((betaML-beta0)./beta0)./1000;%The mean relative estimate of parameters

HR_MRE = sum((HR_ML-HR0)./HR0)./1000;

%Mean relative estimates of Renyi entropy

C. Calculate the KS and AD statistics of UCRD on real data:

clear, clc;

n = 72;%Indicates the total number of samples

w = [0.010 0.033 0.044 0.056 0.059 0.072 0.074 0.077 0.092 0.093 0.096 0.100 0.100 0.102 0.105 0.107 0.107 0.108 0.108 0.108...

0.109 0.112 0.113 0.115 0.116 0.120 0.121 0.122 0.122 0.124...

0.130 0.134 0.136 0.139 0.144 0.146 0.153 0.159 0.160 0.163...

0.163 0.168 0.171 0.172 0.176 0.183 0.195 0.196 0.197 0.202...

0.213 0.215 0.216 0.222 0.230 0.231 0.240 0.245 0.251 0.253...

0.254 0.254 0.278 0.293 0.327 0.342 0.347 0.361 0.402 0.432…

0.458 0.555];%Real dataset

thetaML = ER_2(n, w);

betaML = n./sum(log(1+(1./thetaML).*(w./(1-w)).^2)); %Estimated values of UCRD parameters obtained through dichotomy method

F_u = @(x)1-(1+(1./thetaML).*(x./(1-x)).^2).^(-betaML); %The CDF of UCRD

Fu_values = F_u(w); %Bring this set of real data into the CDF

z = length(w); % Number of data points

ecdf_values = (1:z) / z; % The value of the empirical distribution function

KS_statistic_U = max(abs(Fu_values - ecdf_values)); %Formula for calculating KS test (cumulative distribution function - empirical distribution function)

p_value_U = exp(-2 * (KS_statistic_U)^2 * n); %Formula for calculating p-value

i = 1:z;

AD_u = -g-sum(((2*i-1)./z).*(log(Fu_values)+log(1-Fu_values(end:-1:1)))); %Formula for calculating AD test

D.Calculate the AIC of UCRD on real data:

clear, clc;

n = 72;%Indicates the total number of samples

w = [0.010 0.033 0.044 0.056 0.059 0.072 0.074 0.077 0.092 0.093 0.096 0.100 0.100 0.102 0.105 0.107 0.107 0.108 0.108 0.108...

0.109 0.112 0.113 0.115 0.116 0.120 0.121 0.122 0.122 0.124...

0.130 0.134 0.136 0.139 0.144 0.146 0.153 0.159 0.160 0.163...

0.163 0.168 0.171 0.172 0.176 0.183 0.195 0.196 0.197 0.202...

0.213 0.215 0.216 0.222 0.230 0.231 0.240 0.245 0.251 0.253...

0.254 0.254 0.278 0.293 0.327 0.342 0.347 0.361 0.402 0.432…

0.458 0.555];%Real dataset

thetaML = ER_2(n, w);

betaML = n./sum(log(1+(1./thetaML).*(w./(1-w)).^2)); %Estimated values of UCRD parameters obtained through dichotomy method

f_u = ((2.*betaML.*w)./(thetaML.*(1-w).^3)).*(1+(1./thetaML).*(w./1-w).^2).^(-betaML-1);%The PDF of UCRD

L_u = prod(f_u); %Likelihood function of UCRD

% Calculate AIC

k = 2; % The number of model parameters

AIC_u = -2.* log(L_u) + 2.* k.