Mathematical Problems in Engineering

Volume 2014 (2014), Article ID 610907, 13 pages

http://dx.doi.org/10.1155/2014/610907

## Constructing Matrix Exponential Distributions by Moments and Behavior around Zero

^{1}Dipartimento di Informatica, Università di Torino, Corso Svizzera 185, 10149 Torino, Italy^{2}Faculty of Computing and Information Technology in Rabigh, King Abdulaziz University, Rabigh, P.O. Box 344, Rabigh 21911, Saudi Arabia

Received 1 September 2014; Accepted 5 December 2014; Published 24 December 2014

Academic Editor: Bo Shen

Copyright © 2014 Alessio Angius et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

This paper deals with moment matching of matrix exponential (ME) distributions used to approximate general probability density functions (pdf). A simple and elegant approach to this problem is applying Padé approximation to the moment generating function of the ME distribution. This approach may, however, fail if the resulting ME function is not a proper probability density function; that is, it assumes negative values. As there is no known, numerically stable method to check the nonnegativity of general ME functions, the applicability of Padé approximation is limited to low-order ME distributions or special cases. In this paper, we show that the Padé approximation can be extended to capture the behavior of the original pdf around zero and this can help to avoid representations with negative values and to have a better approximation of the shape of the original pdf. We show that there exist cases when this extension leads to ME function whose nonnegativity can be verified, while the classical approach results in improper pdf. We apply the ME distributions resulting from the proposed approach in stochastic models and show that they can yield more accurate results.

#### 1. Introduction

Probability distributions that can be expressed as the composition of exponential stages have gained widespread acceptance in recent years, specifically as ways of modeling nonexponential durations maintaining the Markov property of the underlying stochastic process. Distributions with these characteristics are represented by the family of phase type (PH) distributions which are given by the distribution of time to absorption in Markov chains (see, e.g., Chapter 2 in [1]). In the literature, mainly continuous-time PH distributions have been studied, but more recently some works that deal with the discrete-time version have also been proposed [2]. Phase type distributions are practically interesting since they allow studying with Markov models problems that are characterized by distributions with coefficient of variation either smaller or larger than one.

Matrix exponential (ME) distributions [3] are an extension of the PH class. They can provide more compact approximation than the PH class [4] and can capture cases in which the density functions that we want to represent may be zero on positive real numbers (i.e., they may exhibit multimodal shapes). ME distributions have the same algebraic form as PH distributions, but they do not enjoy a simple stochastic interpretation (a stochastic interpretation of ME distributions is given in [5], but this is much more complicated and less practical than that possessed by PH distributions).

The main reason to apply PH and ME distributions in stochastic modeling is that these distributions can be easily used as building blocks of more complex models. This fact is well-known for PH distributions because it is straightforward that if we are given a system in which all sojourn times are according to PH distributions and the next state distribution is Markov then the overall system behavior can be described by a Markov chain. Starting from [6], the construction of the Markov chain, often referred to as the expanded Markov chain, was proposed in the literature for several modeling formalisms, such as stochastic Petri nets (SPN), stochastic process algebras, or stochastic automata networks [7–9]. The downside of dealing with the expanded chain is that if the model has many states and/or it describes many activities performed in parallel, then the chain can have a huge number of states. To alleviate this problem several authors proposed techniques for the compact representation of the expanded chain. Such techniques are based on either Kronecker algebra (among the first works see, e.g., [10]) or, more recently, decision diagrams techniques (e.g., [11]).

The necessary theoretical background to use ME distributions as building blocks was developed instead only recently. In this case the overall model is not a Markov chain, but the transient system behavior can still be described by a set of ordinary differential equations. In the context of SPN, this was shown for a subclass of ME distributions in [12], while in [4] the result was extended to the whole family. In the context of quasi-birth-and-death processes the possibility of using ME distributions was investigated in [13].

In order to apply PH or ME distributions in stochastic modeling one needs methods to create distributions that capture real measured durations. Most of the existing techniques fall into two categories: those that are based on the maximum likelihood (ML) principle and those that aim to match moments.

One of the first works on ML estimation considered acyclic PH distributions [14] (i.e., PH distributions whose underlying Markov chain is acyclic), while an approach for the whole family, based on the expectation-maximization method, is proposed in [15]. Since these early papers, many new methods and improvements have been suggested for the whole PH family and for its subclasses (see, e.g., [16, 17]). Much less research tackled ML based fitting of ME distributions because of the lack of a practical stochastic interpretation. One such method, based on semi-infinite programming, is described however in Chapter 9 of [18] where the computational complexity of the problem is discussed and an algorithm is devised.

For what concerns moment matching methods the following results are available. For low-order (≤3) PH and ME distributions (the order of PH and ME distributions is the size of their generator matrices) either moment bounds and moment matching formulas are known in an explicit manner or there exist iterative numerical methods to check if given moments are possible to capture [19–21]. For higher orders there exist matching algorithms, but these often result in improper density functions (by proper pdf we mean a pdf which is nonnegative and normalized) and the validity check is a nontrivial problem [18, 22]. In [23] a simple method is provided that constructs a minimal order acyclic PH distribution given three moments. Characterization of moments of PH and ME distributions is discussed in [24]. Moreover, tool support is available for the construction of PH and ME distributions. Specifically, ML based fitting is implemented in PhFit [25] and a set of moment matching functions is provided in BuTools [26].

In this paper we consider moment matching of ME distributions by Padé approximation [18]. This method provides an order- ME distribution given moments. The problem is that the probability density function (pdf) associated with the distribution can be improper (i.e., it may assume negative values). We show that Padé approximation can be extended in such a way that it considers not only moments, but also the behavior of the original pdf around zero (called also zero-behavior in the sequel). We present examples where this extension leads to ME distributions with proper pdf for cases where the original approach results in improper pdf. Moreover, it can give a better approximation of the shape of the original pdf.

The paper is organized as follows. In Section 2, we provide the necessary theoretical background. The extended Padé approximation is introduced in Section 3. The relation between the zero-behavior of ME distribution and its moments is discussed in Section 4. Numerical illustration of the proposed approach is provided in Section 5. Application of ME distributions resulting from the extended Padé approximation in stochastic models is discussed in Section 6. The paper is concluded then in Section 7.

#### 2. Background

The cumulative distribution function of a matrix exponential random variable, , is of the form
where is a row vector, referred to as the* initial vector*, is a square matrix, referred to as the* generator*, and is a column vector of ones. When the cardinality of the vector is , then the distribution is called* order-** matrix exponential distribution* (ME()). The vector-matrix pair is called the representation of the distribution. Throughout the paper we assume that the distribution is without mass at time zero, that is, , but the extension to the case with mass at zero is straightforward. Note that the entries of may be negative; that is, it is not necessarily a probability vector. If is a probability vector and is the infinitesimal generator of a transient Markov chain, then the distribution belongs to the class of PH distributions.

The pdf of ME distribution can be computed as
and the moment generator function can be expressed as
where is the identity matrix. The last expression results in a rational function of the form
where if there is no probability mass at zero (i.e., if as we have assumed before). If the ME distribution is* nonredundant*, then the number of poles of (4) is equal to the order of the matrix representation of the ME distribution [27]. In this paper, we consider only nonredundant distributions.

The moments can be obtained from the moment generating function and the th moment is It can be seen from (4) that ME() distribution without probability mass at zero is determined by parameters and, consequently, moments can be matched by ME() if the moments are inside the moment bounds of the ME() family (i.e., if the moments are such that they can be realized by ME() distribution).

There exist methods to construct ME distributions given moments. The applicability of these methods is limited by the facts that(i)explicit moment bounds of the ME family are known only for low-order () ME distributions,(ii)the validity check of the pdf for higher order () ME distribution is an open research problem. In the following we report those results that are relevant to our paper.

For order-2 distributions the PH and the ME families are equivalent and the moment bounds are provided in [19] in terms of the squared coefficient of variation defined as . The bounds are meaning that if three moments satisfy the above conditions, then there exists a proper ME distribution that realizes the moments.

For order-3 ME distributions the moment bounds are not known explicitly. In [21] an algorithm is proposed which, given ME density, checks if it is proper. In some cases the algorithm requires the numerical solution of a transcendent equation. This means that, given 5 moments, one can create ME distribution and then check its validity by the method provided in [21]. This check is implemented in the BuTools package [26].

#### 3. Padé Approximation

Padé approximation for ME distributions exploits the fact that for any distribution the following relation holds between the moment generating function and the moments of the distribution: where we assume having which implies that the pdf is normalized. Accordingly, one may consider approximating the moment generating function by a rational function as where the parameters of the moment generating function of the ME distribution, that is, , and , , can be determined based on the first moments of the distribution whose moments we aim to match. More precisely, we look for such parameters with which the first coefficients of the polynomial are equal to the first coefficients of the polynomial The above corresponds to solving a linear system of equations given as with if , , , and if .

The above-described procedure results in a moment generating function that exactly matches the first moments of the original distribution. It may, however, fail as there is no guarantee for the nonnegativity of the corresponding pdf.

The procedure can be extended to capture the behavior of the original pdf around zero. This is done by adding new equations to those provided by (11) and removing those referring to the highest order moments. For example, the pdf of ME() at zero assumes simply the value while the derivative at zero can be obtained as Matching means that the number of moments that can be captured is decreased to , while matching both and means that we can capture moments.

In general, the value of the th derivative of the pdf at zero can be obtained by applying two properties of the moment generating functions. First, the moment generating function of the th derivative of is given by Second, the value at zero can be obtained as Equations (14) and (15) provide a recursive calculation procedure to calculate as function of the parameters of the generating function. For example, in case of , we have

It can be seen from (12) that matching directly determines . Having determined , (13) provides a linear equation that must be satisfied if one aims to capture . In general, capturing for adds simple linear equations to those that are used to match the moments.

#### 4. Relation of Derivatives at Zero and Moments by Hankel Matrices

In this section, we discuss the relation between the zero-behavior of ME distribution and its moments. Let us consider an order ME distribution with representation and introduce the Hankel matrix where . The entries of have the following interpretation: (i)for we have that is related to the th moment of the distribution through ,(ii),(iii)for we have that is related to the derivative at zero of the cdf of the ME distribution because

The rank of the equals [27] which implies that is nonzero for and zero for . As a consequence, given the quantities , one can determine based on the equation For example, given 3 moments of ME distribution, the fourth moment can be calculated by considering or, given and 2 moments, the third moment can be determined based on

The above relation of the derivatives and the moments can be exploited to transform moment matching procedures into zero-behavior matching procedures or into procedures matching part of the zero-behavior and some moments. Moreover, existing moment bounds can be transformed into bounds regarding zero-behavior.

For example, the set of inequalities provided in (6) can be transformed into inequalities regarding , and by solving (21) for and using the result in (6). This leads to meaning that if , and satisfy the above conditions, then there exists a proper ME distribution that realizes them.

#### 5. Numerical Examples

In this section we provide numerical examples of considering the zero-behavior in constructing ME distributions. In Section 5.1, we show that there are cases in which considering the moments of the distribution leads to ME distribution with improper pdf, while considering also the behavior around zero results in a proper distribution. In Section 5.2, we provide an example for which moment matching gives a proper pdf, but the shape of the pdf can be improved by considering the value of the original pdf at zero instead of considering the highest order moment. In Section 5.3, we show that even if the original characteristics around zero are not possible to capture, it is possible to set the characteristics of the ME distribution around zero in such a way that the fitting of the shape of the pdf improves considerably.

##### 5.1. Proper Distribution by Behavior around Zero

Consider the lognormal distribution whose pdf is given by with parameters

Using our proposed approach, we derived two order-3 ME distributions to match the characteristics of the above pdf. The first one matches 5 moments and has an improper pdf as it assumes negative values around 0. Its representation is given by the vector-matrix pair The second one matches three moments and also and and results in a proper distribution. Its representation is The validity of the second ME distribution was checked by using BuTools [26]. In Figure 1, we depicted the original pdf and the two matching ME distributions. In Table 1, we provide the first five moments of the original and the second matching distribution (the moments of the first one are identical to those of the original distribution). For the second ME distribution, even if the 4th and 5th moments are not matched, they are close to those of the lognormal distribution and also the shape of the pdf is similar to that of the lognormal distribution.