A difference scheme for a triangular system of conservation laws with discontinuous flux modeling three-phase flows

1.   Introduction

Record statistics are interest and importance in many areas of real-life applications including data relating to meteorology and weather, economics, sport, athletic events, oil, mining surveys and life testing studies and so on. In the midst of practical life there are many situations that depend entirely on the record statistics (record values) involving: Guinness World Records where only record values are observed (lower or upper record values). Shortest ever tennis matches both in terms of number of games (number of runs) and duration in terms of time (lower record values). News items like fastest time taken to recite the periodic table of the elements (lower record values). Weightlifting competitions where the highest weight is lifted in the fewest number of attempts (upper record values). Swimming competitions where the swimmer travels a certain distance in the least time of swimming (lower record values). Mountaineering, fastest marathon, longest time to hop, etc are of great importance to sports federations and athletes.

From the above, we note that to create a record, several attempts are made, and when the attempt is successful, the record is created, and accordingly we do not get the data for all the attempts that were made to break records. The data that we have are records. Many scientists specially the statisticians have become interested in record values and discussed the statistical inferences for different models based on lower and upper record values including Chandler ^[1] studied the distribution and frequency of lower record values, Resnick ^[2] compared the behavior of the sequence of record values, Shorrock ^[3] proposed asymptotic results on record values and inter-record times, Glick ^[4] discussed the breaking records and breaking boards, Nevzorov ^[5] and Nagaraja ^[6] listed some theoretical characteristics for record values for various univariate model, Kumar ^[7] developed some recurrence relations for generalized logistic distribution based on lower record values, Mahmoud et al. ^[8] studied the estimation problem for the behavior of the coefficient of variation of Lomax distribution based on upper record values, Kumar and Saran ^[9] investigated the ratio and inverse moments of record values from Marshall-Olkin log-logistic distribution, and EL-Sagheer ^[10] studied the estimation of the parameters for generalized logistic distribution based on the lower record values. Arnold et al. ^[11] summarized the study of upper and lower record values and presented many of the basic properties of them as follows: Assume $X_{1}, X_{2}, X_{3}, ...$ a sequence of independent and identically distributed random variables with probability density function (PDF) $f(x)$ and cumulative distribution function (CDF)$\ F\left(x\right)$. Setting $Y_{j} = \max\left(X_{1}, X_{2}, X_{3}, ..., X_{j}\right),$ for $j\geq1$, we say $X_{j}$ an upper record and denoted by $X_{U\left(j\right) }$ if $Y_{j}$ $>$ $Y_{j-1}$ $,$ $j > 1$. Also, let $Y_{j} = \min\left(X_{1}, X_{2}, X_{3}, ..., X_{j}\right),$ for $j\geq1$ we say that $X_{j}$ is a lower record and denoted by $X_{L\left(j\right) }$ if $Y_{j}$ $<$ $Y_{j-1}$ $,$ $j > 1$. Thus the joint density function of the first $n$ upper record values is given by

and the joint density function of the first $n$ lower record values is given by

Dagum ^[12] introduced a statistical distribution closely fitting income distributions which accommodating the heavy tails, called the Dagum distribution (DD). In the actuarial literature, is also called the inverse Burr. Furthermore, DD admits a mixture representation in terms of generalized Gamma and inverse Weibull distributions. The DD can also be obtained as a compound generalized Gamma distribution whose scale parameter follows an inverse Weibull distribution with identical shape parameters. Since Dagum introduced his model as income distribution, many authors have discussed its properties in economics and financial fields, reliability, and survival data, etc. See for instance, Kleiber and Kotz ^[13], Kleiber ^[14], Domma et al. ^[15], Domma and Condino ^[16], among others. Recently, Alotaib et al. ^[17] discussed parameter estimation for the DD based on progressive type-Ⅰ interval censored data. A random variable $X$ has the Dagum distribution and denoted by $X\sim D(\lambda, \theta, \beta)$, if it has PDF, CDF, reliability function $r\left(x\right)$ and hazard rate function $h\left(x\right)$, given by

and

where $\beta,$ $\theta > 0$ are the shape parameters and $\lambda > 0$ is the scale parameter.

In this article, we investigate the Bayes and maximum likelihood estimates for the unknown quantities of the DD using lower record values. The Lindley's approximation and the Metropolis-Hasting (M-H) algorithm within Gibbs sampler are proposed to construct the Bayes estimates. To this end, the results of the Bayes estimates are obtained under both squared error (SE) and linear exponential (LINEX) loss functions. Also, the corresponding credible intervals (CRIs) are constructed under MCMC technique.

The rest of this article is organized as follows. Section $2$ deals with the maximum likelihood estimates and asymptotic confidence intervals. Bayesian estimates using Lindley's approximation and MCMC technique are provided in Section $3$. In Section $4$ a simulation study is provided to examine the performance of the proposed estimation methods. Section $5$ discusses an application to real life data. Finally, a brief conclusion is given in Section $6$.

2.   Maximum likelihood estimation

In this section, we discuss the maximum likelihood estimates (MLEs) and approximate confidence intervals (ACIs) of the parameters as well as reliability and hazard rate functions of DD given in (3) when the available data are lower record values. Suppose that $\underline{x} = x_{L\left(1\right) }, x_{L\left(2\right) }, x_{L\left(3\right) }, ..., x_{L\left(n\right) }$ be the lower record values of size $n$ from $D(\lambda, \theta, \beta)$, then the likelihood function for observed record $\underline{x}$ is given by

The logarithm of the likelihood function may then be written as

The estimators $\hat{\lambda},$ $\hat{\theta}$ and $\hat{\beta}$ of the parameters $\lambda$, $\theta$, and $\beta$, respectively, can be then obtained as the solution likelihood equations

and

Since (9)–(11) cannot be solved analytically for $\hat{\lambda},$ $\hat {\theta}$ and $\hat{\beta}$, some numerical methods such as Newton-Raphson method must be employed to obtain the desired MLEs in such situations. Once we get desired MLEs of $\lambda,$ $\theta$ and $\beta$ using the invariant property of MLEs, for given mission time $t,$ the MLEs of $r\left(t\right)$ and $h\left(t\right)$ can be obtained after replacing $\lambda,$ $\theta$ and $\beta$ by $\hat{\lambda},$ $\hat{\theta}$ and $\hat{\beta}$ as

2.1.   Asymptotic confidence intervals

From the log-likelihood function $\ell(\lambda, \theta, \beta{\Large |} \underline{x})$ in (8), we have

The Fisher information matrix $I\left(\lambda, \theta, \beta\right)$ is then obtained by taking expectation of negative Eqs (13)–(17). Based on some mild regularity conditions, $\left(\hat{\lambda}, \hat{\theta}, \hat{\beta }\right)$ is approximately bivariate normal with mean $\left(\lambda, \theta, \beta\right)$ and covariance matrix $I^{-1}\left(\lambda, \theta, \beta\right)$. Practically, we usually estimate $I^{-1}\left(\lambda, \theta, \beta\right)$ by $I^{-1}\left(\hat{\lambda }, \hat{\theta}, \hat{\beta}\right)$. Moreover, a simpler and equally valid procedure is to use the approximation

where

ACIs for $\lambda,$ $\theta$ and $\beta$ can be obtained by to be bivariate normal distributed with mean $\left(\lambda, \theta, \beta\right)$ and variance-covariance matrix $I_{\circ}^{-1}\left(\lambda, \theta, \beta\right)$. Therefore, the $100(1-\gamma)\%$ ACIs for $\lambda,$ $\theta$ and $\beta,$ respectively, are

where $\Lambda_{11},$ $\Lambda_{22}$ and $\Lambda_{33}$ are the elements on the main diagonal of the variance-covariance matrix $I_{\circ}^{-1}\left(\lambda, \theta, \beta\right)$, and $z_{\gamma/2}$ is the percentile of the standard normal distribution with right-tail probability $\gamma/2$.

In order to find the approximate estimates of the variance of $r\left(t\right)$ and $h\left(t\right)$, we use the delta method, see Agresti ^[18]. Assume that $G_{1} = (\frac{\partial r\left(t\right) } {\partial\lambda}, \frac{\partial r\left(t\right) }{\partial\theta}, \frac{\partial r\left(t\right) }{\partial\beta})$ and $G_{2} = (\frac{\partial h\left(t\right) }{\partial\lambda}, \frac{\partial h\left(t\right) }{\partial\theta}, \frac{\partial h\left(t\right) }{\partial\beta }),$ where $\frac{\partial r\left(t\right) }{\partial\lambda},$ $\frac{\partial r\left(t\right) }{\partial\theta},$ $\frac{\partial r\left(t\right) }{\partial\beta},$ $\frac{\partial h\left(t\right) }{\partial\lambda},$ $\frac{\partial h\left(t\right) }{\partial\theta},$ and $\frac{\partial h\left(t\right) }{\partial\beta}$ are the first partial derivatives of $r(t)$ and $h(t)$ with respect to the parameters $\lambda,$ $\theta$ and $\beta$, respectively. The approximate asymptotic variances of $\hat{r}(t)$ and $\hat{h}(t)$ are given by

where and $G_{i}^{T}$ is the transpose of $G_{i}$, $i = 1, 2$. Thus, the $(1-\gamma)100\%$ ACIs for $\hat{r}(t)$ and $\hat{h}(t)$ can be written as

3.   Prior information and Bayesian estimation

Besides being an accurate analysis and not easy to handle mathematically, the use of the Bayesian method allows the incorporation of previous knowledge of the parameters through informative prior densities. In the Bayesian approach, the information referring to the model parameters is obtained through posterior marginal distributions. In this section, two Bayesian inference procedures (Lindley's approximation and MCMC technique) are proposed to estimate the parameters $\lambda,$ $\theta$ and $\beta$ as well as the reliability $r(t)$ and hazard rate $h(t)$ functions. For estimating these quantities, we assume mainly both SE and LINEX loss functions. However, any other loss function can be easily implemented. From the perspective of frequentist, a natural choice for the prior distributions of $\lambda,$ $\theta$ and $\beta$ would be to assume that the three quantities are independent gamma distributions $G(a_{1}, b_{1}),$ $G(a_{2}, b_{2})$ and $G(a_{3}, b_{3})$ respectively, where $a_{i}$ and $b_{i},$ $i = 1, 2, 3$ reflect the knowledge of prior about $\left(\lambda, \theta, \beta\right)$ and they are assumed to be known and nonnegative hyper-parameters. Further, the joint prior distribution can be written as

Subsequently, via Bayes' theorem across combining the likelihood function given (7) with the joint prior given in (20), the joint posterior density function of $\lambda,$ $\theta$ and $\beta$ can be written as follows

Thus, the conditional posterior densities of $\lambda,$ $\theta$ and $\beta$ can be computed and written, respectively, as

and

Therefore, for any function, say $u\left(\lambda, \theta, \beta\right),$ the Bayes estimates under SE and LINEX loss functions can be written as

There is a difficulty to compute (25) analytically. Therefore, we propose to approximate it by using both Lindley's approximation and MCMC technique with Metropolis-Hasting algorithm within Gibbs sampler, to obtain the Bayes estimates of $\lambda,$ $\theta$, $\beta$, $r(t)$ and $h(t)$.

3.1.   Lindley's approximation

There are various methods to approximate the ratio of integrals of the above form, may be the simplest one is Lindley approximation method, which was introduced by Lindley ^[19] to approximate the Bayes estimates into a form containing no integrals. This method has been used with various lifetime distributions by several authors for instance, Singh et al. ^[20], Ying and Gui ^[21], Ramin and Kohansal ^[22]. Briefly, it can be described as follows: If the sample size is sufficiently large, according to Lindley ^[19], any ratio of the integral of the form

where $u\left(\lambda, \theta, \beta\right)$ is the function of $\lambda,$ $\theta$ and $\beta$ only, $\ell\left(\lambda, \theta, \beta\right)$ is the log-likelihood function and $\rho\left(\lambda, \theta, \beta\right) = \log g(\lambda, \theta, \beta)$. Hence, $I(x)$ can be estimated as

where $\hat{\lambda},$ $\hat{\theta},$and $\hat{\beta}$ are the MLEs of $\lambda,$ $\theta$ and $\beta$ respectively, and

where $\phi_{1} = \lambda$, $\phi_{2} = \theta$, $\phi_{3} = \beta$, $i, j, l = 1, 2, 3$ and $\sigma_{ij}$ is $\left(i, j\right) th$ elements of the matrix $\hat {I}_{\circ}^{-1}\left(\hat{\lambda}, \hat{\theta}, \hat{\beta}\right)$ in Subsection (2.1).

From (20), $\rho(\lambda, \theta, \beta)$ can be written as

Hence, we get

Now, under SE and LINEX loss function, the approximate Bayes estimators $\lambda$ are given by

and

Similarly, we can obtain the approximate Bayes estimators of $\theta$, $\beta$, $r(t)$ and $h(t)$.

3.2.   MCMC technique

It is known that there are several procedures of MCMC technique available in which samples are generated from the conditional posterior densities. One of the simplest MCMC procedures is the Gibbs sampling procedure. Another procedure is considered the Metropolis-Hastings (M-H) algorithm and later extended by Hastings ^[23]. A more general procedure of MCMC procedures which we will use is considered the M-H within Gibbs sampling. Gibbs sampler is required to decompose the joint posterior distribution into full conditional distributions for each parameter and then sample from them. Since, $g_{2}^{\ast}(\theta|\lambda, \beta, \underline{x})$ follows gamma $\left[n+a_{2}, b_{2}+\log\left(1+\lambda x_{L\left(n\right) }^{-\beta}\right) \right]$, it is quite simple to generate samples of $\theta$ from $g_{2}^{\ast}(\theta|\lambda, \beta, \underline{x})$ by implementing any gamma generating routine. In addition, the conditional posterior densities $g_{1}^{\ast}(\lambda|\theta, \beta, \underline{x})$ and $g_{3}^{\ast} (\beta|\theta, \lambda, \underline{x})$ can't be reduced analytically to well-known distributions. So, according to Tierney ^[24] M-H algorithm within Gibbs sampling with normal proposal distribution is used to conduct the MCMC methodology. The hybrid M-H algorithm and Gibbs sampler works as follows:

(1) Use the MLEs as the initial value, denoted by $\hat{\lambda}^{\left(0\right) }, \hat{\theta}^{\left(0\right) }$and $\hat{\beta}^{\left(0\right) }$.

(2) Set $j = 1.$

(3) Generate $\theta^{(j)}$ from Gamma$\left(n+a_{2}, b_{2}+\log\left(1+\lambda x_{L\left(n\right) }^{-\beta}\right) \right)$.

(4) Using M-H algorithm, generate $\lambda^{\left(j\right) }$ and $\beta^{\left(j\right) }$ from $g_{1}^{\ast}(\lambda^{\left(j-1\right) }|\theta^{\left(j\right) }, \beta^{\left(j-1\right) }, \underline{x}),$ and $g_{3}^{\ast}(\beta^{\left(j-1\right) }|\theta^{\left(j\right) }, \lambda^{\left(j-1\right) }, \underline{x})$, respectively, with the normal distributions $N\left(\lambda^{\left(j-1\right) }, \Lambda_{11}\right)$ and $N\left(\beta^{\left(j-1\right) }, \Lambda_{33}\right)$.

(a) Generate $\lambda^{\ast}$ from $N\left(\lambda^{\left(j-1\right) }, \Lambda_{11}\right)$ and $\beta^{\ast}$ from $N\left(\beta^{\left(j-1\right) }, \Lambda_{33}\right)$.

(b) Evaluate the acceptance probabilities

(c) Generate a $\rho_{1}$ and $\rho_{2}$ from a Uniform $\left[0, 1\right]$ distribution.

(d) If $\rho_{1} < \Omega_{\lambda}$ accept the proposal and set $\lambda^{\ast} = \lambda^{\left(j\right) }$, else set $\lambda^{\left(j\right) } = \lambda^{\left(j-1\right) }$.

(e) If $\rho_{2} < \Omega_{\beta}$ accept the proposal and set $\beta^{\ast } = \beta^{\left(j\right) }$, else set $\beta^{\left(j\right) } = \beta^{\left(j-1\right) }$.

(5) Compute $r\left(t\right)$ and $h\left(t\right)$ as

(6) Set $j = j+1$.

(7) Repeat Steps $(3)-(6)$ $M$ times. The first $M_{0}$ simulated varieties are ignored to remove the affection of the selection of initial value and to guarantee the convergence. Then the selected sample are $\lambda^{(j)}, \theta^{(j)}, \beta^{\left(j\right) }, r^{^{(j)}}(t)$ and $h^{^{(j)}}(t),$ $j = M_{0}+1, ..., M,$ for sufficiently large $M$, forms an approximate posterior sample which can be used to develop the Bayesian inference.

(8) Under SE and LINEX Loss functions the approximate Bayes estimate of $\eta$ (where $\eta =$ $\lambda,$ $\theta,$ $\beta,$ $r\left(t\right)$ and $h\left(t\right)$) can be obtained by

where $M_{0}$ is the burn-in period and $\eta^{\left(j\right) } =$ $\lambda^{\left(j\right) },$ $\theta^{\left(j\right) },$ $\beta^{\left(j\right) },$ $r^{\left(j\right) }\left(t\right)$ and $h^{\left(j\right) }\left(t\right)$.

(9) To compute the CRIs of $\eta$, order $\left\{ \eta^{M_{0}+1}, \eta^{M_{0}+2}, \ldots, \eta^{M}\right\}$ as $\left\{ \eta^{\left[1\right] }, \eta^{\left[2\right] }, \ldots, \eta^{\left[M-M_{0}\right] }\right\}$. Then, the $\left(1-\gamma\right) 100\%$ CRIs of $\eta$ can be given by

4.   Simulation study and comparisons

In this section, through a simulation study using $R$ language, we compare the performances of the MLEs and Bayes estimates of the unknown parameters as well as reliability and hazard rate functions of the DD under SE and LINEX loss functions proposed in the previous sections. Here, we consider the simulation results in the case of $\left(\lambda, \theta, \beta\right) = \left(3.2, 5.9, 1.75\right)$. Then, the true values of $r\left(t\right)$ and $h\left(t\right)$ at time $t = 1.3$ are computed to be $0, 963895$ and $0, 234171$, respectively. The performance of estimators is evaluated in terms of mean squared error (MSE) which computed as

where $\varpi_{1} = \lambda,$ $\varpi_{2} = \theta,$ $\varpi_{3} = \beta,$ $\varpi_{4} = r(t)$ and $\varpi_{5} = h(t)$ for the point estimates (MLE, Lindley, MCMC), also average lengths (ALs) and coverage probability (CPs), which computed as the number of CIs that covered the true values divided by $1000$, for interval estimates (ACIs and CRIs). Bayes estimates and the highest posterior density CRIs are computed based on $M = 12000$ MCMC samples and discard the first values $M_{0} = 2000$ as burn-in. In addition, we assume the informative gamma priors for $\lambda,$ $\theta$ and $\beta$ that is, when the hyper-parameters are $a_{i} = 1.3$ and $b_{i} = 2.2,$ $i = 1, 2, 3$. Moreover, $95\%$ CRIs were computed for each simulated sample. In our study, we consider different sample sizes $n = 5,$ $7,$ $9,$ $12,$ $15,$ $18,$ $20,$ $23$ and $25$. The gamma distribution was chosen as prior distributions of the parameters because it is the most appropriate one that matches the maximum likelihood function. Moreover, they are from the same family. The evidence for this is that two of the full conditional posterior distributions of the parameters resulted in a gamma distribution, which proves the validity of the chose. In addition, choosing another prior distribution or dependent prior will increase the complexity and difficulty of mathematical equations. Gamma distribution is one of the rich distributions, as when changing its parameters (hyper-parameters), we get new data with new information, so it is the focus of attention of most statisticians. Gamma distribution is flexible model with any values of hyper-parameters. A special case: When all hyper-parameters of gamma distribution are zero, we obtain Jaffrey prior. The results of means, MSEs, ALs and CPs of estimates are displayed in Tables 1–5.

From the results in Tables 1–5, we observe the following:

(1) It is clear that from all Tables, as sample size $n$ increases, the MSEs and average interval lengths decrease, also Bayes estimates are performed better than the MLEs of $\lambda$, $\theta$, $\beta$, $r\left(t\right)$ and $h\left(t\right)$ in terms of MSEs and average interval lengths, as expected.

(2) Bayes estimates under MCMC technique perform better than Bayes estimates under Lindley approximation in the sense of having smaller MSEs.

(3) Bayes estimates under Lindley approximation perform better than MLEs in the sense of having smaller MSEs.

(4) Bayes estimates under LINEX loss function perform better than Bayes estimates under SE loss function in the sense of having smaller MSEs.

(5) Bayes estimate under LINEX with $c = 2$ provides better estimates for $\lambda$, $\theta$, $\beta$, $r\left(t\right)$ and $h\left(t\right)$ because of having the smallest MSEs.

(6) Bayes estimates under LINEX for the choice $c = 2$ are performed better than their estimates for the choice $c = -2$ in the sense of having smaller MSEs.

(7) Both MLE and Bayesian methods have very close estimates and their ACIs have quite high CPs (around $0.95$). Also, the Bayesian CRIs have the highest CPs.

(8) In general, we can conclude that the best estimation method is the Bayes method under MCMC technique with LINEX loss function, especially if the prior information about the problem under study is available.

5.   Applications

To illustrate the inferential procedures developed in the preceding sections. we use the data set from Almongy ^[25], see Table 6. A complete sample from the data represents a COVID-19 data belonging to the Netherlands of $30$ days, which recorded from $31$ March to $30$ April $2020$. This data formed of rough mortality rate. The data are as follows

For the goodness of fit test, we compute the Kolmogorov-Smirnov (K-S) distance. It was found that the K-S is $0.072956$ and the associated p-value is $0.9936$. Therefore, according to the p-value we can say that the DD fits quite well to the above data. Empirical, $Q-Q$ and $P-P$ plots are shown in Figure 1, which clear that the DD fits the data very well. Table 7 lists the lower record values from the observed data.

Based on these $12$ lower record values, the MLEs and ACIs for $\lambda$, $\theta$, $\beta$, $r\left(10\right)$ and $h\left(10\right)$ are computed to be as in Tables $6$ and $7$. Moreover, to compute the Bayesian estimates under Lindley's approximation and MCMC technique, the prior distributions of the parameters $\lambda$, $\theta$ and $\beta$ are needed to specify. Since we do not have any prior information available, we assume that the non-informative gamma priors for $\lambda,$ $\theta$ and $\beta$ that is, when the hyperparameters are $a_{i} = 0.0001$ and $b_{i} = 0.0001,$ $i = 1, 2, 3$. Under MCMC technique, the posterior analysis was done across combining M-H algorithm within Gibbs sampler. To conduct the MCMC algorithm, which was described in Subsection $3.2$, the initial values for the parameters $\lambda,$ $\theta$ and $\beta$ were taken to be their MLEs. In addition, we run the chain for $12000$ times and discard the first $2000$ values as 'burn-in' to avoid the effect of the initial values (starting point). The Bayesian estimates as well as $95\%$ CRIs for $\lambda$, $\theta$, $\beta$, $r\left(10\right)$ and $h\left(10\right)$ are presented in Tables 8 and 9.

6.   Conclusions

In this paper, we have estimated the reliability and hazard rate functions of the Dagum distribution based record statistics. To get the best estimators for this purpose, the maximum likelihood and Bayes estimations methods have been utilized. The MLEs and the asymptotic confidence intervals based on the observed Fisher information matrix have been discussed. For Bayesian estimation approach, we assume the gamma priors for the unknown parameters and provide the Bayes estimators under the assumptions of SE and LINEX loss functions. Since the Bayes estimates cannot be obtained in explicit form, we have applied both Lindley's approximation and MCMC technique with Metropolis–Hasting algorithm within Gibbs sampler to compute the approximate Bayes estimates and constructed the credible intervals. The performance of the proposed techniques has been compared via a Monte Carlo simulation study with different sample sizes $n$. It was found that if we have no prior information on the unknown parameters, then it is always better to use the MLEs rather than the Bayes estimators, because the Bayes estimators are computationally more expensive. Finally, for illustrative purpose, the details have been explained using a real-life data set of COVID-19 mortality rates.

Conflict of interest

The authors declare no conflict of interest.