Qualitative results and numerical approximations of the $(k, \psi)$-Caputo proportional fractional differential equations and applications to blood alcohol levels model

1.   Introduction

Modeling industrial data is now considered a substantial area of interest for researchers across different disciplines, including biometrics, engineering, survival, lifetime, reliability sciences, and numerous other areas. Different probability distributions are available, but they have some limitations in fitting these types of data sets. For example, see Verevka et al. ^[1] and Barskov et al. ^[2]. There is a constant increase in the range of count data, and the constraints of existing models make modeling these data challenging. Therefore, in the last few decades, many researchers have sought to propose adaptable models for modeling these types of data sets using different generalized approaches. For more details, see Alzaatreh and Famoye ^[3], Ieren et al. ^[4], Riad et al. ^[5], EL-Helbawy et al. ^[6], Altun et al. ^[7], Alotaibi et al. ^[8], Maya et al. ^[9], Meraou et al. ^{[10,11,12,13]}, Alrweili et al. ^[14], Alrweili ^[15,16].

On the other side, several situations exist where the suggested extension distributions are unsuitable for analyzing different data sets. Additionally, the role and importance of generating a new family of distribution using various generators becomes a very important and compulsory task for researchers to accommodate the variety of data patterns being generated in every field of life. Data coming from different fields of study, specifically in the industrial field, all need a better model to fit their diverse data patterns and be motivated by the urgency of highly flexible statistical models. Consequently, Nofal et al. ^[17] introduced a new methodology to create a new distribution called the generalized Kumaraswamy (GK) family. It is a new concept of generalizing a given distribution, which introduces three additional parameters in a baseline distribution, and it has wider applications in industrial, engineering, survival, and other fields. The cumulative distribution function (cdf) and the probability density function (pdf) of the GK class of distributions can be defined as

and

where $f(x)$ and $F(x)$ are, respectively, the pdf and cdf of the basic distribution.

It is well documented that the Rayleigh distribution (RD) is frequently used to model diverse data sets drawn from different areas, especially for analyzing industrial data. Consequently, the RD is suitable for modeling industrial data. This is vital in the reliability analysis of industrial devices, such as 3D printing, drones, and robots. Let us consider the random variable $X$ has RD with positive parameter $\theta$, so its cdf and pdf, respectively are

and

The RD is widely used in reliability and survival analysis for mortality rates, especially when studying extreme events. Since it captures tail behavior effectively, it is particularly also useful in understanding the upper quantiles of life expectancy or survival time. Also, it has undoubtedly established itself as a crucial tool for data modeling across nearly all sectors, including survival, hydrology, insurance, and energy theory. However, despite its widespread use and advantages, the RD is constrained by its inherent limitations. One of the primary constraints of the RD is its capacity to represent solely monotonically increasing forms of hazard functions, as it can only model data where the hazard rate increases or decreases consistently over time. Also, the RD is regarded as a limiting model for residual lifetimes. For this, we have seen increased interest in studying the RD and its applications in various fields such as medicine, engineering, insurance, industry, and risk management. Some of these efforts are listed below by Chukwudi et al. ^[18], Muzamil et al. ^[19], Aijaz et al. ^[20], Abdulsalam et al. ^[21], Javed et al. ^[22].

The basic motivations for the recommended GKRD in practice are:

$(1)$ The GKRD provides a crucial important role in analyzing numerous kinds of data sets. Its parameters provide a flexible way to manipulate the shape and characteristics of a probability distribution. This adaptability allows researchers and analysts to tailor the distribution to better fit real-world data, making it a valuable tool in diverse fields such as statistics, engineering, biology, etc. Further, the four parameters of the proposed GKRD make the underlying patterns more interpretable. This enhanced interpretation ability can lead to deeper insights and a better understanding of the factors influencing the data.

$(2)$ With adding three additional parameters, the GKRD has the ability to represent the unimodal or bimodal probability distribution.

$(3)$ Another motivation is the ability to induce skewness in symmetrical and asymmetrical distributions. This capability is particularly valuable in fields where skewed distributions are prevalent, such as finance, economics, and insurance.

The following are the key objectives:

● The primary objective that must be fulfilled is introducing a novel model, and the new distribution is named the generalized Kumaraswamy Rayleigh distribution (GKRD). We also determined its various statistical properties, including moments, the moment-generating function, and order statistics.

● Derive and discuss its reliability characteristics.

● Estimate the model parameters using the maximum likelihood and Bayesian approaches under the square error loss function (SELF) and illustrate the pattern of these derived estimators using a comprehensive simulation study.

● Check the validity and flexibility of the GKRD using industrial and financial data sets.

The rest of the study is organized as follows: The recommended GKRD is defined in Section 2 with some distributional properties. In Section 3, we present several mathematical properties including moments, the quantile function, and the moment generating function. Section 4 demonstrates the estimation of the model parameters based on two different proposed methods. The effectiveness of the proposed estimation tools are studied using some simulation studies in Section 5, and three distinct real data sets are applied to show the results of the application of the GKRD. The concluding report is given in Section 7.

2.   Derivation of the generalized Kumaraswamy Rayleigh distribution

Here, in this part of the study, several distributional properties such as the pdf, cdf, survival, and hazard rate function of the GKRD are derived.

2.1.   Model and assumption

Based on Eqs (1.1)–(1.2) and by replacing the classical distribution with the RD, the cdf and pdf of the proposed GKRD are

The pdf curves of the GKRD are explored using several parametric values of parameters and displayed in Figure 1. As shown in Figure 1, the proposed GKRD is unimodal and has decreasing behavior.

The survival function (sf) and hazard rate function (hrf) can be obtained from the following equations:

and

Figure 2 reports the hrf curves of the GKRD using several parameter values. It is upside down and increasing depending on the parameter values.

2.2.   Special cases

The proposed GKRD has several specialized sub-models, which confirm its importance in modeling various types of data sets. The specialized sub-models are displayed in Table 1.

3.   Theoretical characteristics of the GKRD

In this section, some key mathematical characteristics of the GKRD are investigated.

3.1.   Quantile function

The quantile function $\mathcal{Q}_p$ of the proposed GKRD is given by

where $p\in (0, 1)$ represents the probability level. Further, with $p = \dfrac{1}{2}$, the value of the median is obtained, and it is

The skewness ($\mathcal{S}$) and kurtosis ($\mathcal{K}$) coefficients can be obtained using the formula:

and

3.2.   Useful expansion

Let us define the following series as:

Thus, by applying Eq (3.4) in (2.2), the pdf of the GKRD becomes

where $\eta_{l}(\alpha, \beta, \gamma) = (-1)^{2l}\dbinom{\beta-1}{l}\dbinom{\gamma-1}{l} \ \alpha^{l+1}$.

3.3.   Moments and associated measures

The $k^{th}$ ordinary moment of $T$ that follows the GKRD is defined as follows:

with $\omega_{k, l}(t, \beta, \theta) = \int_{0}^{\infty} t^{k+1} e^{-\theta(l+1)t^2} \ \left(1-e^{-\theta t^2}\right)^{\beta l} \ dt$.

Taking $k = 1$ and $2$ in Eq (3.6), the first and second moments of origin of the GKRD can be obtained as

and

Next, the variance and coefficient of variance (CV) for the GKRD can be found as follows:

and

3.4.   Moment-generating function

The moment-generating function (MGF) of the GKRD can be derived as

We have numerically assessed several statistical summary measures such as Mean, Variance, CV($T$), $\mathcal{S}(T)$, and $\mathcal{K}(T)$ for different parametric values and posted them in Tables 2 and 3. The same can easily be observed for these quantities from the plots presented in Figures 3 and 4. Henceforth, the GKRD is an option to model the positively skewed and leptokurtic data sets.

3.5.   Order statistic

Let a random sample $t_1, t_2, \ldots, t_m$ represent a continuous GKRD. Based on Arnold et al. ^[23] and David et al. ^[24], the density function of the $k^{th}$ order statistic is as follows:

Next,

The $\phi_{(1:m)}(t)$ minimum-order statistics is obtained by substituting $m = 1$ in equation:

Similarly, we will get the expression of the $m^{th}$ order-statistic by replacing $k = m$,

4.   Statistical inference

In this estimation section, we estimate the parameters of the GKRD using two estimation methods. These methods, which include maximum likelihood (MLE) and Bayesian estimators under SELF, are crucial in enhancing ecological studies.

4.1.   Maximum likelihood method (MLE)

Consider $\{t_1, t_2, \ldots, t_m\}$ a random sample of size $m$ is taken from the GKRD and its associated log likelihood function $\mathcal{LL}$ is

Now differentiate the above equation for $\rho = (\alpha, \beta, \gamma, \theta)$

and

The parameter estimates of $\rho = (\alpha, \beta, \gamma, \theta)$ are the solution of the above non-linear equations (4.2)–(4.5). Because these normal equations lack closed-form solutions, we use numerical methods to effectively solve them and derive ML estimates such as the Newton-Raphson, fixed point, or secant methods. To achieve this goal, we used the optim function in R software for the estimation process.

4.2.   Bayesian estimation

Compared to the maximum likelihood estimation approach, Bayesian estimation is a more current and efficient approximation. Considering past data and samples, we can make the Bayesian estimation.

We consider the independent informative type of priors for the parameters $\rho = (\alpha, \beta, \gamma, \theta)$ as

The posterior density of $\rho$ has the below form:

Hence, the Bayes estimation based on SELF

is obtained to be:

By obtaining the joint prior, the posterior function can be determined, and it can be applied to the Metropolis-Hasting method.

5.   Simulation experiments

Here, we discuss the performance of the two proposed estimators, MLE and Bayes, considering a finite number of samples. We do a simulation study with various samples $(\alpha, \beta, \gamma, \theta)$ (Scenario 1: $\rho = (0.75, 1, 2, 1.5)$, Scenario 2: $\rho = (0.8, 1.1, 2.3, 1.8)$, Scenario 3: $\rho = (0.9, 1.2, 2.5, 2)$) from the GKRD. We calculated the mean estimates (Mean), average bias (Bias), root mean square error (RMSE), and the efficiency (Eff) to weigh the MLEs and Bayes accuracy. Additionally, we present the mean of the number of iterations (NIT) required for convergence in each method (the Metropolis-Hasting technique for Bayes estimator and Newton Raphson for the MLE technique), showing that convergence occurs within $number = 1000$ steps.

The computations were obtained employing the R program with the function optim for Newton Raphson technique and optim for the Metropolis-Hasting procedure by taking the values of $\rho$ as Scenario 1, Scenario 2, and Scenario 3 respectively. Recall that, for the Bayesian estimation, we choose the gamma informative prior to obtaining the final estimate $\hat{\rho}$. The algorithm for computing the unknown parameters for the GKRD is presented in details in Appendix. The results of the simulation are presented in Tables 4–6. The following expression is utilized to generate random samples from the suggested model:

Based on the findings presented in Tables 4–6, we can concluded the following points:

● As $m$ increases the parameter estimates become closer to the true parameter. It appears obvious that the estimates of $\rho$ are generally unbiased for the two methods of estimates.

● The RMSEs also show a decreasing pattern with an increase in $m$ for the two methods of estimates.

● The results show that the Bayes estimator under SELF achieves an excellent performance among ML estimators. This is evident from the consistently low values of RMSE observed across all cases.

● Based on the NIT for the two proposed estimation procedures, it indicates the excellent performance of the Bayes estimator under SELF compared to the MLE method.

● The performance Bayes under SELF estimator is better than the MLE technique in all procedure scenarios because all efficiency values are greater than 1.

6.   Real application data sets

The most important part of statistical inference is the application of actual data and mathematical modeling. Data is frequently modeled under established probability distributions in the manufacturing and quality control industries in order to discover faults, guarantee product quality, and keep standards consistent. So we always need a new distribution.

Here in this study, the proposed GKRD is used to model two industrial data sets taken from Saudi Arabia (KSA) and one from financial data, and the resulting fits are compared to the competitive continuous distributions including the Kumaraswamu Gul alpha power transformed Rayleigh distribution (KGAPRD), Poisson generalized Rayleigh distribution (PGRD), new generalized Rayleigh distribution (NGRD), generalized Rayleigh distribution (GRD), Rayleigh distribution (RD), gamma distribution (GD), Weibull distribution (WD), and log-normal distribution (LND). Using the conventional criteria of the lowest values of the Akaike information criterion (AIC) and Bayesian information criterion (BIC), the fits given by the GKRD and the other examined models were compared. Furthermore, the comparison of the fitted distributions was assessed using the Kolmogorov-Smirnov ($\mathcal{KS}$) test with its associated $\mathcal{P}$-values and Anderson-Darling test (AD). For more information about recently financial applications see Atchadé et al. ^[25], and Kamal et al. ^[26].

6.1.   First data

The first real life data set used in this study was studied by Yu et al. ^[27], and it is based on the efficiency of the construction industry and their pure technical between 2013 and 2022 in the KSA. The values of the proposed data set are shown in Table 7.

6.2.   Second data

Here are the values provided about the construction industry and their scale efficiency in the KSA from 2013 to 2022. The suggested data set was considered by Yu et al. ^[27] and is shown in Table 8.

6.3.   Third data

In this subsection, the data set is drawn from group medical insurance. It is defined as the total loss for all the claim amounts exceeding 25,000 USD during 1991. The values of the data set are given in euros (EUR), and it is available at http://www.soa.org as well as used by Meraou et al. ^[10].

Based on the three proposed data sets, various non-parametric plots including the kernel density, fitted histogram, scaled total time on the test (TTT), probability-probability (PP), QQ normal, and box plots are plotted in Figures 5–7.

Table 9 summarized the obtained results of goodness-of-fit test with the ML estimates of all fitted model using the industrial and financial data sets. It is well documented that the GKRD is recognized as the optimal choice for the three data sets, with the following:

$(1)$ Data 1: $\alpha = 0.2, \ \ \beta = 2.4, \ \ \gamma = 10, \ \ \theta = 0.05$.

$(2)$ Data 2: $\alpha = 0.9, \ \ \beta = 3.75, \ \ \gamma = 20, \ \ \theta = 0.01$.

$(3)$ Data 3: $\alpha = 0.1, \ \ \beta = 1.8, \ \ \gamma = 30, \ \ \theta = 0.02$.

6.4.   Concluding remarks on data analysis

Clearly, from the obtained results in Table 9, for the three data sets, the GKRD is an efficient, superior model among competing models based on all AIC and BIC measures. This ensures that the GKRD is the most appropriate model among the choices. Specifically, in terms of $P$-value and AD, the GKRD outperforms all other models considered, confirming its status as the most optimal distribution for the three data sets when compared to alternative models. Additionally, the plots for the estimated pdf versus fitted histogram and estimated cdf versus empirical cdf are plotted in Figures 8–10. The various plots presented also confirm a good fit for the GKRD to the considered data sets.

7.   Concluding remarks

For modeling industrial and financial data sets that many models lack, to model, we defined a novel four-parameter probabilistic model. The new model was created using the generalized Kumaraswamy technique, resulting in the generalized Kumaraswamy Rayleigh model. Several important distributional and statistical characteristics have been determined and analyzed. By using a wide range of methods, including the classical MLE and Bayesian techniques, we are able to handle statistical analysis of the GKRD distribution and its unknown parameters. Therefore, we came to the conclusion that the Bayes approach is superior to the conventional estimating method since it consistently produces lower values for the MSE. Additionally, using three real-life data sets taken from industrial and financial domains, the results indicate that the proposed GKRD distribution effectively analyzes both data sets compared to competing distributions.

8.   Future work

In future work, this study may attract the bivariate case of the GKRD. In addition, this study may contribute to the estimation of model parameters in censored samples based on several cases. In addition, the proposed model can attract a wider set of applications, such as in engineering and environmental fields.

Author contributions

Alanazi Talal Abdulrahman, Tariq S. Alshammari and Ramlah H Albayyat worked on mathematics; Eslam Hussam, Amirah Saeed Alharthi and Khudhayr A. Rashedi worked on english and programming. All the authors have read and approved the final version of the manuscript for publication.

Acknowledgments

This research has been funded by the Scientific Research Deanship at the University of Ha'il, Saudi Arabia, through project number (RG-24 068).

Availability of data and materials

All data exists in the paper with their related references.

Funding

This research has been funded by the Scientific Research Deanship at the University of Ha'il, Saudi Arabia, through project number RG-24 068.

Use of Generative-AI tools declaration

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

Conflict of interest

All authors declare no conflicts of interest in this paper.

Appendix

pdf.GKRD=function(star,x){

alpha=star[1]

beta=star[2]

delta=star[3]

theta=star[4]

2*beta*delta*alpha*theta*x*exp(-theta*x^2)/(1-(1-alpha)^delta)*

(1-exp(-theta*x^2))^(beta-1)*(1-alpha*(1-exp(-theta*x^2))^beta)^(delta-1)

}

t=seq(0,2,len=1000)

plot(t,pdfGKRD(c(0.95,0.45,3,1.2),t),col="red",type="l",lwd=2,lty=4)

cdf.GKRD=function(star,x){

alpha=star[1]

beta=star[2]

delta=star[3]

theta=star[4]

(1-(1-alpha*(1-exp(-theta*x^2))^beta)^delta)/(1-(1-alpha)^delta)

}

t=seq(0,2,len=1000)

plot(t,cdf.GKRD(c(0.2,1.3,1.5,4.5),t),col="red",type="l",lwd=2,lty=4)

hrf.GKRD <- function(star,x){

alpha=star[1]

beta=star[2]

delta=star[3]

theta=star[4]

2*beta*delta*alpha*theta*x*exp(-theta*x^2)/(1-(1-alpha)^delta)*

(1-exp(-theta*x^2))^(beta-1)*(1-alpha*(1-exp(-theta*x^2))^beta)^(delta-1)/

(1-(1-(1-alpha*(1-exp(-theta*x^2))^beta)^delta)/(1-(1-alpha)^delta))

}

t=seq(0,2,len=1000)

plot(t,hrf.GKRD(c(0.95,0.3,0.6,3.5),t),col="red",type="l",lwd=2,lty=4)

## Estimation

fMLE<-function(star,x){

alpha=star[1]

beta=star[2]

delta=star[3]

theta=star[4]

-sum(log(beta*delta*alpha*theta*exp(-theta*x^2)/(1-(1-alpha)^delta)*

(1-exp(-theta*x^2))^(beta-1)*(1-alpha*(1-exp(-theta*x^2))^beta)^(delta-1)

))

}

NB=100; nb=100; res.alpha=numeric(NB);res.beta=numeric(NB)

res.delta=numeric(NB);res.theta=numeric(NB)

for(i in 1:NB){

alpha=0.75;delta=2; beta=1.0;theta=1.5

u=runif(nb,0,1)

X=sqrt(-1/theta*log(1-(1-(1-u*(1-(1-alpha)^delta))^(1/delta))^(1/beta)))

res.alpha[i]=optim(c(alpha,beta,delta,theta),fMLE,method="N",x=X)＄par[1]

res.beta[i]=optim(c(alpha,beta,delta,theta),fMLE,method="N",x=X)＄par[2]

res.delta[i]=optim(c(alpha,beta,delta,theta),fMLE,method="N",x=X)＄par[3]

res.theta[i]=optim(c(alpha,beta,delta,theta),fMLE,method="N",x=X)＄par[4]

}

AEMLE.alpha=mean(res.alpha); AEMLE.beta=mean(res.beta)

AEMLE.delta=mean(res.delta); AEMLE.theta=mean(res.theta)

AB.alpha=abs(mean(res.alpha-alpha)); AB.beta=abs(mean(res.beta-beta))

AB.delta=abs(mean(res.delta-delta)); AB.theta=abs(mean(res.theta-theta))

MSEMLE.alpha=mean((alpha-res.alpha)**2); MSEMLE.beta=mean((beta-res.beta)**2)

MSEMLE.delta=mean((delta-res.delta)**2); MSEMLE.theta=mean((theta-res.theta)**2)

## Application

result=goodness.fit(pdf = pdf.GKRD, cdf =cdf.GKRD, method = "BFGS",

starts = c(alpha,beta,delta,theta), data = data, domain = c(0,Inf),mle = NULL)