An Overview of Generalized Gamma Mittag–Leffler Model and Its Applications

In model building situations in physical, biological, social and engineering sciences, the usual procedure is to select a probability model from a parametric family of distributions. In many practical problems, it is often found that the selected model is not a good fit for the experimental data, because it requires a model with a thicker or thinner tail than the ones available from the parametric family of distributions. In order to make the tail thicker or thinner, a technique is introduced here by augmenting a series to the original density. Our first step is to construct the thicker or thinner tailed distribution associated with the Mittag–Leffler function, because this function is connected to fractional calculus, the Mittag–Leffler stochastic process, non-Gaussian time series, L$\acute{e}$vi flights and in a limiting process to the topics of Tsallis statistics, superstatistics, as well as to statistical distribution theory.

In reaction rate theory, input-output type situations and reaction-diffusion problems in physics and chemistry, when the integer derivatives are replaced by fractional derivatives, the solutions automatically go in terms of Mittag–Leffler functions and their generalizations; see Haubold and Mathai (2000) [1]. The ordinary and generalized Mittag–Leffler functions interpolate between a purely exponential law and power-law-like behavior of phenomena governed by ordinary kinetic equations and their fractional counterparts; see Kilbas et al. (2004) [2], Kiryakova (2000) [3], Mathai (2010) [4] and Mathai et al. (2010) [5]. This paper examines a new family of statistical distributions associated with Mittag–Leffler functions, which gives an extension to the gamma family, which will then connect to fractional calculus and statistical distribution theory through the theory of special functions. The model investigated in this paper is useful in the study of life testing problems, reliability analysis, in physical situations to describe stable solutions, as well as unstable and chaotic neighborhoods, etc. We will start with the definition of the Mittag–Leffler function.

A two-parameter Mittag–Leffler function is defined as follows: where $\Re (·)$ denotes the real part of $(·)$. Observe that the Mittag–Leffler function is an extension of the exponential function. When $\mathsf{\alpha }=1,\mathsf{\beta }=1$, Equation (1) reduces to ${\mathrm{e}}^x$. In statistical model building, usually, the parameters are real, but since the results to be discussed hold for complex parameters, as well, we will state the corresponding relevant conditions. Various properties, generalizations and applications of the Mittag–Leffler function can be seen from Kilbas et al. (2004) [2].

Consider a probability density of the form: where C is the normalizing constant and $E_{\mathsf{\alpha },\mathsf{\beta }}(-\mathsf{\delta }{x}^{\mathsf{\alpha }})$ is the Mittag–Leffler function. Some interesting special cases of Equation (2) are the following: The density in Equation (2) includes two-parameter gamma, exponential, chi square, noncentral chi square and the like. When $\mathsf{\delta }=0$, Equation (2) reduces to the two-parameter gamma density:

For $\mathsf{\alpha }=\mathsf{\beta }=1$ in Equation (2), we have the exponential density. For $\mathsf{\beta }=\frac{n}{2}$ and $a=\frac{1}{2}$ in $f_1(x)$, we have the chi square density with n degrees of freedom. For $\mathsf{\beta }=p,p=2,3,\cdots$ in Equation (3), we have the Erlang density. For $p=1$ in Erlang density, we have the exponential density. For fixed values of $a,\mathsf{\beta }$ and for various values of δ, we can look at the graphs that give a suitable interpretation to the model in Equation (2).

The above figures show a comparison between gamma density and gamma Mittag–Leffler density for different values of δ. Observe that $\mathsf{\delta }=0$ corresponds to the gamma density. In Figure 1, as the value of δ decreases, the right tail of the new density becomes thicker and thicker compared to that of a gamma density. Similarly the peakedness of the curve slowly decreases. In Figure 2, also, $\mathsf{\delta }=0$ corresponds to the gamma density. When the values of δ increases from $\mathsf{\delta }=0$, the right tail of the new density becomes thinner and thinner compared to that of a gamma density. Similarly, the peakedness of the curve slowly increases. Hence, when we look for a model with a thicker or thinner tail while a gamma density is found to be more or less a proper fit, then a member from the new family of densities introduced here will become quite useful and handy. Observe that the new density is mathematically and computationally easily tractable, just like a gamma density.

The moment-generating function of Equation (2) is given by:

On simplification, we obtain mgfas:

The characteristic function can be obtained if we replace t by $it,i=\sqrt{-1}$. If we put $-t$ instead of t, then we will obtain the Laplace transform of the density. Using the Laplace transform or moment-generating function, we can easily obtain the integer moments by using the following formula:

Thus, the mean value will be obtained as: and:

Variance:

Arbitrary moments can be obtained in terms of the generalized Wright hypergeometric function. That is, which is nothing but the Mellin transform of the function f with $\gamma =s-1$. where ${}_p{\Psi }_q(z)$ is the generalized Wright’s hypergeometric function defined for $z\in C$, complex $a_i,b_j\in C$ and ${\mathsf{\alpha }}_i,{\mathsf{\beta }}_j\in {\Re }_{+}=(0,\infty ),a_i,b_j\ne 0;i=1,2,\dots ,p;j=1,2,\dots ,q$ by the series:

The function in Equation (8) was introduced by Wright and is called the generalized Wright’s hypergeometric function. For convergence conditions, the existence of various contours and other properties, see Wright (1940) [6] or the theory of H-function to be discussed later. If we take γ as integers in Equation (7), then we will obtain integer moments. It may be observed that the distribution function is available in terms of a series of incomplete gamma functions.

In this section, we have given explicit forms of the estimators of the parameters using the method of moments. The motivation for the method of moments comes from the fact that the sample moments are consistent estimators for the corresponding population moments. To start with, let us consider the case for $a=1,\mathsf{\beta }=1$ in Equation (2), then the model will become the Mittag–Leffler extension of the standard exponential distribution and has the density of the form:

The moments can be obtained from Equations (5) and (6) with the parameter value $a=1,\mathsf{\beta }=1$.

The moment estimators of δ and α are given by: and:

Now, consider the Mittag–Leffler extension of standard gamma density. For that, we take $a=1$ in Equation (2); thus, the model has the following form:

Using the same procedure as above, one can obtain the estimators as: and:

If we retain all parameters $\mathsf{\alpha },\mathsf{\beta },\mathsf{\delta }$ and a, then the analytical solution is quite difficult, but the numerical solution can be obtained by using software like MATLAB and MAPLE.

In this section, a generalized density of Equation (2) is considered. Mathai (2005) [7] introduced the pathway model, where the scalar version of the pathway density is given as follows: $\text{where}{c}_1=\frac{\rho {[a(q-1)]}^{\frac{\mathsf{\beta }}{\rho }}\mathrm{\Gamma }(\frac{1}{q-1})}{2\mathrm{\Gamma }(\frac{\beta }{\rho })\mathrm{\Gamma }(\frac{1}{q-1}-\frac{\mathsf{\beta }}{\rho })},\frac{1}{q-1}>\frac{\mathsf{\beta }}{\rho }.$

For $q<1$, writing $q-1=-(1-q)$, we have the form: $\text{where}{c}_2=\frac{\rho {[a(1-q)]}^{\frac{\mathsf{\beta }}{\rho }}\mathrm{\Gamma }(\frac{1}{1-q}+\frac{\mathsf{\beta }}{\rho }+1)}{2\mathrm{\Gamma }(\frac{\mathsf{\beta }}{\rho })\mathrm{\Gamma }(\frac{1}{1-q}+1)}.$

As $q\to 1$, $f_1(x)$ and $f_2(x)$ tend to $f_3(x)$, which we refer to as the extended symmetric generalized gamma distribution, where $f_3(x)$ is given by:

Now, generalize Equation (2) by using Equation (11). Since the normalizing constants obtained in this section are in terms of H-functions, a definition of the H-function is given here. where: where ${\mathsf{\alpha }}_j,j=1,2,...,p,{\mathsf{\beta }}_j,j=1,2,...,q$ are real positive numbers, $a_j,j=1,2,...,p$ and $b_j,j=1,2,...,q$ are complex numbers and L is a contour separating the poles of $\mathrm{\Gamma }(b_j+{\mathsf{\beta }}_{j}s),j=1,2,...,m$ from those of $\mathrm{\Gamma }(1-a_j-{\mathsf{\alpha }}_{j}s),j=1,2,...,n$. The details of convergence conditions, various properties and applications of H-function are available in Mathai et al. (2010) [5].

Consider: where:

The general density in Equation (15), for $q>1$, includes the Mittag–Leffler extension of Type 2 beta density and F density. If $x>0$ is replaced by $|x|,-\infty <x<\infty$, with the appropriate change in the normalizing constant, then we have the Student t density and Cauchy density as special cases, and many more. As $q\to 1$, $f_4(x)$ reduces to $f_5(x)$, where:

This includes the Mittag–Leffler extension of generalized gamma, Weibull, chi-square, Maxwell–Boltzmann and related densities. If $\mathsf{\delta }=0$, Equation (17) reduces to the generalized gamma density. When $\rho =1$, Equation (17) reduces to Equation (2). In particular, if $\rho =\mathsf{\beta }$, the Mittag–Leffler extension of Weibull density is obtained,

For $x>0$ and $q<1$, $f_4(x)$ takes the following form: where:

This general density includes the Mittag–Leffler extension of Type 1 beta density, triangular density, uniform density and many more. As $q\to 1$, Equation (15) reduces to Equation (17). For $\rho =\mathsf{\beta }$, then Equations (15) and (19) will give the two forms of extended q-Weibull density as given below. where the normalizing constants $C_8$ and $C_9$ are obtained from Equations (16) and (20), respectively. In particular, if $\mathsf{\delta }=0$, then $f_8(x)$ and $f_9(x)$ will coincide with q-Weibull density. In particular, for $\mathsf{\beta }=1$ and $\rho =1$, the integral of Equation (19) can be treated as the pathway fractional integral transform of the Mittag–Leffler function with appropriate transformation of the variable. More details about pathway fractional integral transform can be found in Seema Nair ((2009) [8], (2011) [9]).

The model in Equation (15) can be connected to the superstatistics of Beck and Cohen (2003) [10] and Beck (2006) [11] for $q=2,\mathsf{\delta }=0$. There is a vast literature on Beck–Cohen superstatistics in statistical mechanics with many physical interpretations; see, for example, Beck ((2010) [12], (2009) [13]). From a statistical point of view, it can be obtained as an unconditional density from Bayes’ setup, by considering generalized gamma as a conditional density with the assumption that the scale parameter a has a prior gamma density; see Nair and Kattuveettil (2010) [14]. For $\mathsf{\beta }=1,\rho =1,a=1$ and $\mathsf{\delta }=1$, we have Tsallis statistics for $q>1$ and $q<1$, respectively, from the models in Equations (15) and (19). Thus, we can make a connection to nonextensive statistical mechanics; see, for example, Tsallis ((1988) [15], (2004) [16], (2009) [17]) and Mathai et al. (2010) [5]. (It is stated that over 3000 papers have been published on this topic during the past 20 years, and over 5000 people from around the globe are working in this area).

Especially when $a=0$, $\mathsf{\beta }=\mathsf{\alpha }$, such that $0<\mathsf{\alpha }\le 1$ in Equation (2), we have the Mittag–Leffler density, which has the probability density function: and has Laplace transform $\mathsf{\delta }{(s^{\mathsf{\alpha }}+\mathsf{\delta })}^{-1}$, a special case of the general class of Laplace transforms discussed by Mathai et al. (2006) [18]. Similarly, for $a=0$ in (2) and replacing β by $\mathsf{\alpha }\gamma$ and δ by $\frac{1}{\rho }$, one can arrive at the generalized Mittag–Leffler density, which has the probability density function: where $E_{\mathsf{\alpha },\mathsf{\beta }}^{\mathsf{\delta }}(·)$ is the generalized Mittag–Leffler function defined as: where ${(\mathsf{\delta })}_k=\mathsf{\delta }(\mathsf{\delta }+1)(\mathsf{\delta }+2)\cdots (\mathsf{\delta }+k-1)$ and Equation (24) has the Laplace transform of the form:

The Laplace transform in Equation (25) is associated with infinitely-divisible and geometrically infinitely-divisible distributions, α-Laplace and Linnik distributions. The class of Laplace transforms relevant in geometrically infinitely-divisible and α-Laplace distributions is of the form $L(p)=\frac{1}{1+\eta (p)}$, where $\eta (p)$ satisfies the condition $\eta (bp)=b^{\mathsf{\alpha }}\eta (p)$ and $\eta (p)$ is a periodic function for fixed α. In Mathai (2010) [4], it is noted that there is a structural representation in terms of the positive L$\acute{e}$vi random variable by using Equation (25).

A positive L$\acute{e}$vi random variable $u>0$, with parameter α is such that the Laplace transform of the density of $u>0$ is given by ${\mathrm{e}}^{-t^{\mathsf{\alpha }}}$. That is:

Theorem 4.1. Let $y>0$ be a L$\acute{e}$vi random variable with Laplace transform as in (26) and $x\ge 0$ an exponential random variable with parameter δ, and let x and y be independently distributed. Then, $u=y{x}^{\frac{1}{\alpha }}$ is distributed as generalized Mittag–Leffler random variable with Laplace transform (see [4]):

Proof: For proving this result, we will make use of the following lemma of the conditional argument.

Lemma 1.1. For two random variables u and v having a joint distribution, whenever all of the expected values exist, where the inside expectation is taken in the conditional space of u given v and the outside expectation is taken in the marginal space of v.

Now, by applying Equation (27), we have the following: Let the density of u be denoted by $g(u)$. Then, the conditional Laplace transform of g, given x, is given by:

However, the right side of Equation (28) is in the form of a Laplace transform of the density of x with parameter $t^{\mathsf{\alpha }}$. Hence, the expected value of the right side, with respect to x, is available by replacing the Laplace parameter t in ${[1+\rho t]}^{-\gamma }$ (the Laplace transform of a gamma density with the scale parameter ρ and shape parameter γ) by $t^{\mathsf{\alpha }}$, which establishes the result. From Equation (27), one can observe that if we consider an arbitrary random variable y with the Laplace transform of the form: whenever the expected value exists, where $\varphi (t)$ is such that: then from Equation (28), we have:

Now, let x be an arbitrary positive random variable having Laplace transform, denoted by $L_x(t)$ where $L_x(t)=\psi (t)$. Then, from Equations (27) and (30), we have:

In particular, when $\mathsf{\delta }=1$ in Equation (23), we have the standard Mittag–Leffler density introduced by Pillai (1990) [19]. For the properties of Mittag–Leffler distributions connected to the autoregressive process, see the papers of Pillai and Jayakumar (1995) [20] and Jayakumar (2003) [21].

We can look at the model Equation (2) in another way, as well; consider the total integral as: which can be treated as the Laplace transform of the function $x^{\mathsf{\beta }-1}{E}_{\mathsf{\alpha },\mathsf{\beta }}(-\mathsf{\delta }{x}^{\mathsf{\alpha }})$, and hence, $C^{-1}=\frac{a^{\mathsf{\alpha }-\mathsf{\beta }}}{a^{\mathsf{\alpha }}+\mathsf{\delta }}$, where C is the normalizing constant of Equation (2), is nothing but the Laplace transform of the given function. It is shown to be very relevant in fractional reaction-diffusion problems in physics, since it naturally occurs in the derivation of the inverse Laplace transform of the functions of the type $a^{\mathsf{\alpha }}(d+b{a}^{\mathsf{\beta }})$, where a is the Laplace transform parameter and $d\mathrm{and}\mathrm{b}$ are constants.

For more details, see the papers of [8,18,23]. Reaction-diffusion equations are modeling tools for the dynamics presented by a competition between two or more species, activators and inhibitors or production and destruction, that diffuse in a physical medium. As an example, [1] considered the evolution of a star like the Sun, which is governed by a second order system of differential equations, the kinetic equations, describing the rate of change of chemical composition of the star for each species in terms of the reaction rates for destruction and production of that species [24,25,26]. Methods for modeling processes of destruction and production have been developed for bio-chemical reactions and their unstable equilibrium states [27] and for chemical reaction networks with unstable states, oscillations and hysteresis [28].

Consider an arbitrary reaction characterized by a time-dependent quantity $N=N(t)$. It is possible to equate the rate of change $\frac{\mathrm{d}N}{\mathrm{d}t}$ to a balance between the destruction rate d and the production rate p of N, that is $\frac{\mathrm{d}N}{\mathrm{d}t}=-d+p$. In general, through feedback or other interaction mechanisms, destruction and production depend on the quantity N itself: $d=\mathrm{d}(N)$ or $p=p(N)$. This dependence is complicated, since the destruction or production at time t depends not only on $N(t)$, but also on the past history $N(\tau ),\tau <t$, of the variable N. This may be formally represented by: where $N_i$ denotes the function defined by $N_i(t^{*})=N(t-t^{*}),t^{*}>0$. The production and destruction of species is described by kinetic equations governing the change of the number density $N_i$ of species i over time, that is, where $<\sigma \nu {>}_{mn}$ denotes the reaction probability for an interaction involving species mand n, and the summation is taken over all reactions, which either produce or destroy the species i [29]. The first sum in Equation (39) can also be written as: where $a_i$ is the statistically expected number of reactions per unit volume per unit time destroying the species i. It is also a measure of the speed at which the reaction proceeds. In the following, we are assuming that there are $N_j(j=1,\dots i,\dots )$ species j per unit volume and that for a fixed $N_i$, the number of other reacting species that interact with the i-th species is constant in a unit volume. Following the same argument, for the second sum in Equation (39) accordingly, where $N_i{b}_i$ is the statistically expected number of the i-th species produced per unit volume per unit time for a fixed $N_i$. Note that the number density of species i, $N_i=N_i(t)$, is a function of time, while the $<\sigma \nu {>}_{mn}$, containing the thermonuclear functions (see Haubold and Kumar (2008) [30]), are assumed to depend only on the temperature of the gas, but not on the time t and number densities $N_i$. Then, Equation (38) implies that:

For Equation (42), we have three distinct cases, $c_i=a_i-b_i>0$, $c_i<0$ and $c_i=0$, of which the last case says that $N_i$ does not vary over time, which means that the forward and reverse reactions involving species i are in equilibrium; such a value for $N_i$ is called a fixed point and corresponds to a steady-state behavior. The first two cases exhibit that either the destruction $(c_i>0)$ of species i or production $(c_i<0)$ of species i dominates.

For the case $c_i>0$, we have: with the initial condition that $N_i(t=0)=N_0$ is the number density of species i at time $t=0$, and it follows that:

The exponential function in Equation (44) represents the solution of the linear one-dimensional differential Equation (43) in which the rate of destruction of the variable is proportional to the value of the variable. Equation (43) does not exhibit instabilities, oscillations or chaotic dynamics, in striking contrast to its cousin, the logistic finite-difference equation [26,29]. A thorough discussion of Equation (43) and its standard solution in Equation (44) is given in [25].

Let us consider the standard fractional kinetic equation: which is derived by [1]. ${}_0{D}_t^{-\nu }$ is the Riemann–Liouville fractional integral of order ν. It is very easy to find the solution of this equation. On applying Laplace transforms on both sides, with parameter p, one has:

Now, applying inverse Laplace transform,

This solution can be written in terms of the Mittag–Leffler function as:

Theorem 5.1. For $\nu >0$, let the coefficient of $N_0$ of the fractional integral equation be g(t) an arbitrary function, so that the fractional integral becomes: Then, the Laplace transforms of the fractional integral equation will be:

Now, one can take various forms for $G(p)$; all of the special cases are considered in [23,31]. For example:

Let us consider $G(p)=p^{-\mathsf{\beta }}{}_1{F}_1({\gamma }_1,{\mathsf{\beta }}_1,p^{-\mathsf{\alpha }})$ for $\mathsf{\alpha }>0,\mathsf{\beta }>0,\nu >0,c>0$, then Equation (47) becomes:

On applying inverse Laplace transform, we get:

Haubold and Mathai (2000) [1] generalized the standard kinetic Equation (43) to a standard fractional kinetic Equation (45), derived solutions of a fractional kinetic equation that contains the particle reaction rate (or thermonuclear function) as a time constant and provided the analytic technique to further investigate possible modifications of the reaction rate through a kinetic equation. The Riemann–Liouville operator in the fractional kinetic equation introduces a convolution integral with a slowly-decaying power-law kernel, which is typical for memory effects referred to in [32]. This technique may open an avenue to accommodate changes in standard solar model core physics as proposed by [33]. In the solution of the standard fractional kinetic Equation (45), given in Equation (46), the standard exponential decay is recovered for $\nu =1$. However, the Equation (46), for $0<\nu <1$, shows a power-law behavior for $t\to \infty$ and is constant (initial value $N_0$) for $t\to 0$.

In this section, we present an application of generalized gamma Mittag–Leffler density in modeling currency exchange rates. Here, we considered 1152 observations starting from 2004 to 2008, of U.S. dollar-Indian rupee foreign exchange rates, which is available at www.rbi.org. The log returns of the exchange rates are considered. The following Figure 3 is the graph of the log-transformed data embedded with the generalized gamma Mittag–Leffler density.

From the data, the summary statistics for transformed currency exchange rates obtained are as follows:

For convenience, assume that $\alpha =1$; by using the method of moments, estimates of parameters are obtained and are given as:

From the figures, we can observe that the generalized gamma Mittag–Leffler density fits well to the data with the above estimated values of the parameters. Here, we used the chi square statistic to measure the goodness of fit. The calculated chi square value is 2.44, and the corresponding tabled value is 15.507. Hence, we conclude that the model in Equation (2) is a good fit to the dataset considered.

In this section, a multivariate analogue of the density in Equation (2) is considered. The multivariate Mittag–Leffler function will be defined as:

A function close to $E_p(·)$ of Equation (51) is the multi-index Mittag–Leffler function of Kiryakova; see, for example, Kiryakova (2000) [3]. Now, the general density is defined as:

The normalizing constant is available, by proceeding as before, as:

Various properties of one form of a multivariate gamma density can be seen from Griffiths (1984) [34]. The following sections reveal the joint Laplace transform and product moments of the density function.

The Joint Laplace transform is available from Equation (52) by replacing $a_j$ by $a_j+t_j,j=1,\dots ,p$ and then dividing by the normalizing constant. That is:

The arbitrary product moments are obtained in terms of generalized Wright function and are given by:

Conditional density of $x_1$ given $x_2,\dots ,x_p$ is given by: