MASARYK UNIVERSITY
FACULTY OF INFORMATICS
}w¡¢£¤¥¦§¨!"#$%&123456789@ACDEFGHIPQRS`ye|
Phase-Type Approximation
Techniques
BACHELOR THESIS
Zuzana Kom´arkov´a
Brno, Spring 2012
Declaration
Hereby I declare, that this thesis is my original authorial work, which I
have worked out by my own. All sources, references and literature used or
excerpted during elaboration of this work are properly cited and listed in
complete reference to the due source.
In Brno, May 18, 2012
Zuzana Kom´arkov´a
Advisor: RNDr. Vojtˇech ˇReh´ak, Ph.D.
iii
Acknowledgement
I would like to thank my advisor, RNDr. Vojtˇech ˇReh´ak, Ph.D., for his help,
support and guidance throughout the work on this thesis. I would also like
to express my gratitude to Luboˇs Korenˇciak for his cooperation on the topic,
his valuable consultations and his never ending patience with me. I would
like to thank Sam Pr´ıvara and Honza Kˇret´ınsk´y for their comments on the
work. Furthermore, I want to thank all the people I had the pleasure to share
the room with for the friendly atmosphere and their help whenever I needed.
Lastly, I would like to thank my family and my friends for their support
during my studies.
v
Abstract
We consider the phase-type approximation problem, i.e. approximating
general distributions by phase-type distributions. We describe the most
common methods to deal with the problem. Further, for each method we
survey the most relevant techniques and the corresponding tool support.
Moreover, we compare the methods and discuss how to choose the most
appropriate method under given circumstances.
vii
Keywords
phase-type distributions, continuous-time stochastic systems, continuoustime
Markov chains, discretization
ix
Contents
1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2 Foundations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.2 Phase-Type Distributions . . . . . . . . . . . . . . . . . . . . . 6
2.3 Demonstration Examples . . . . . . . . . . . . . . . . . . . . 10
3 Moment Matching Method . . . . . . . . . . . . . . . . . . . . . . 13
3.1 Matching Less Than Three Moments . . . . . . . . . . . . . . 13
3.2 Matching Three Moments . . . . . . . . . . . . . . . . . . . . 14
3.3 Matching More Than Three Moments . . . . . . . . . . . . . 16
4 Minimization of a Difference . . . . . . . . . . . . . . . . . . . . . 17
4.1 Kullback-Leibler Divergence . . . . . . . . . . . . . . . . . . . 17
4.2 Other Distance Measures . . . . . . . . . . . . . . . . . . . . . 20
5 Discretization and Discrete Phase-Type Approximation . . . . . 23
5.1 Discrete Phase-Type Distributions . . . . . . . . . . . . . . . 23
5.2 Discretization of Continuous Distributions . . . . . . . . . . 25
5.3 Discrete Phase-Type Approximation . . . . . . . . . . . . . . 27
6 Hybrid Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
6.1 Heavy Tailed Distributions . . . . . . . . . . . . . . . . . . . 30
7 Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
7.1 Accuracy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
7.2 Time Efficiency of the Approximating Algorithm . . . . . . . 33
7.3 Number of Phases . . . . . . . . . . . . . . . . . . . . . . . . . 34
7.4 Generality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
7.5 Demonstration Examples . . . . . . . . . . . . . . . . . . . . 34
7.6 Guideline for the Choice of Suitable Method . . . . . . . . . 35
8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
1
Chapter 1
Introduction
The effort to better understand natural as well as human-made systems
around us has brought significant attention to the modelling and analysis of
complex systems. Deterministic models can appropriately describe systems
based mainly on physical laws, which we can precisely predict. In practice,
totally deterministic systems are unlikely due to influence of unpredictable
factors. Therefore, stochastic (probabilistic) models took their place. To govern
the system dynamics, stochastic models enable one or more variables to
be determined by a random variable.
A lot of common formalisms for stochastic systems are based on
continuous-time Markov chains (CTMC) because of their analytical tractability
[40]. A key assumption in Markov chain analysis is that all the waiting
times have to be distributed according to an exponential distribution. Such
a distribution can well describe many situations, e.g. distances between
mutations on a DNA or service times of agents in queuing systems. Unfortunately,
in many real-life applications, assuming exponential distributions
leads to significant imprecisions, e.g. electronic components failure times
often follow Weibull distribution [1].
A simple extension of CTMC is to allow waiting times to be general,
which leads to the semi-Markov process (SMP), a sequential formalism,
where the target state is first chosen and then the waiting with corresponding
time distribution starts and further, the generalized semi-Markov process
(GSMP), a parallel formalism, where all available waiting times start at once
and after the first one finishes we move to its corresponding state [10]. The
SMP and the GSMP provide better accuracy. On the other hand we have to
deal with much more complicated analysis. A possible approach, how to deal
with such an analysis, is a phase-type approximation (PH-approximation).
The basic idea behind the PH-approximation is to replace general distributions
in a model by phase-type distributions (PH-distributions), which
are defined as the time spent in transient states of a CTMC with one absorbing
state until absorption. This special structure of PH-distributions allows
us to use analytic methods for CTMC in the model. As shown in [32], the
3
1. INTRODUCTION
class of PH-distributions can approximate every non-negative distribution
arbitrarily close.
PH-distributions were first introduced by Neuts in 1975 [33], since then
they have been used in a wide range of stochastic modelling applications
such as queuing theory, finance, biostatistics, survival analysis, and telecommunications.
In 1917, Erlang was the first person who considered a generalization
of exponential distributions by his method of stages (phases)
[17]. He defined a non-negative random variable as the time taken to move
through a fixed number of phases, where the time spending in each phase
is exponentially distributed with rate λ. Such a distribution is said to be an
Erlang distribution. Further, it can be generalized to a Coxian distribution
[11] and both of them are special cases of the family of PH-distributions.
There is a large number of techniques for the PH-approximation. However,
there is no work summarizing them. In this thesis we present and
analyze the main approaches for the PH-approximation and introduce particular
techniques based on them. We focus on their practicability in different
situations. In Chapter 2 we introduce some necessary notions as well as properties
and the proper definition of PH-distributions. In the following four
chapters we show main approaches for the PH-approximation. They are the
following: 1) moment matching method, 2) minimization of a difference,
3) discretization and 4) hybrid methods. We consider conditions on their
applicability, the most recent representatives of such approaches and helpful
tools.
In Chapter 7 we provide a comparison of the main approaches and
characterize situations, for which each of them is suitable. The comparison
is based on the following four criteria: 1) accuracy of the approximation,
2) time efficiency of the approximation algorithm, 3) number of phases in the
resulting PH-distribution and 4) generality of the method, since naturally, it
should work for as broad class of distributions as possible. The third criterion
is important due to the subsequent Markovian analysis of the model.
Here we need to keep the state-space of the underlying Markov chain small
and thus have PH-distributions with very few phases.
4
Chapter 2
Foundations
In this chapter, we introduce the class of continuous PH-distributions and
examine its properties. We also consider some important subclasses of the
class of PH-distributions. At the end of this chapter we introduce distributions,
which we use as examples for demonstration of particular methods in
the following chapters.
2.1 Preliminaries
We need to introduce exponential distributions and continuous-time Markov
chains for the formal definition of PH-distributions.
Definition 1. A continuous random variable X with the probability density function
(pdf)
fX(x) = λe−λx, where x ≥ 0, λ > 0
is said to have an exponential distribution with rate λ.
Rate λ represents the mean number of events per time unit.
The exponential distribution enjoys the status of mainstay in stochastic
modelling for several reasons. One of its most important features is the
memoryless property
Pr(X > x + t|X > x) = Pr(X > t) for t ∈ R ≥ 0,
where Pr(A) stands for the probability of event A. Exponential distributions
are the only continuous distribution that exhibits this property.
Now we can define continuous-time Markov chains.
Definition 2. A Continuous-Time Markov Chain (CTMC) is a tuple (S, R, v0),
where
• S is a finite set of states,
5
2. FOUNDATIONS
• R : S × S → R is a rate matrix, such that R(i, j) ≥ 0 for i = j and
R(i, i) = −
k=i
R(i, k), and
• v0 is an initial probability vector.
We can describe the semantics of CTMC intuitively as follows: the process
starts randomly in a state of S according to the initial probability vector v0
and moves successively from one state to another. Whenever the process
is in a state si waiting times tj are chosen according to the exponential
distributions with rates R(i, j). After the shortest time tj, the process moves
to sj.
Definition 3. A state is called absorbing if once it has been reached it is impossible
to leave it.
Definition 4. A CTMC is called absorbing if it has at least one absorbing state
and from every state it is possible to reach some absorbing state.
Definition 5. In an absorbing CTMC, a state which is not absorbing is called
transient.
2.2 Phase-Type Distributions
The name of phase-type distributions comes from the key idea behind it,
which is to model waiting times as being made up of the time taken to
move through a number exponentially distributed phases. These phases are
equivalent to states of the underlying Markov chain and we will use term
phases when writing about PH-distributions, and states when writing about
CTMC. In this thesis, row vectors are denoted by (x1 x2 ... xn) or x, column
vectors are denoted by (x1 x2 ... xn)T or xT . Denotations 0 and 1T stand
for the row vector of 0’s and the column vector of 1’s. At this moment we
are ready for a formal definition of PH-distributions.
Definition 6 ([33]). Consider an absorbing CTMC with n+1 states, where n ≥ 1,
such that states 1, ..., n are transient and state 0 is an absorbing state. Further, let
the chain have an initial probability of starting in any of the n + 1 phases given by
the probability vector (p0 p). The (continuous) n-phase phase-type distribution
(PH-distribution) is the distribution of time from the above process’s starting until
absorption in the absorbing state.
6
2. FOUNDATIONS
This distribution is described by the initial probability vector (p0 p) and
the transition (n + 1) × (n + 1) rate matrix R of the form
R =
0 0
tT T
,
where T is an n × n matrix. As R is the rate matrix of an absorbing CTMC, it
holds
T1T = −tT and p1T = 1 − p0.
Therefore, we need n2 + n parameters T and p to fully specify an n-phase
PH-distribution. Now, we provide a few basic properties of PH-distributions.
Proposition 1 ( [28]). Assume that X is an n-phase PH-distribution with parameters
T and p. Its cumulative distribution function (cdf) is given by
FX(x) = 1 − p · exp(Tx)1T for x ≥ 0
and its pdf is given by
fX(x) = p · exp(Tx)tT for x > 0,
where the matrix exponential is defined by
exp(A) =
∞
n=0
1
n!
An
.
Proposition 2 ([32]). The set of PH-distributions is dense in the set of non-negative
distributions.
This means that arbitrary non-negative distribution can be arbitrarily
close approximated by a PH-distribution and thus models with PHdistributions
can describe any stochastic process.
One term needed before introducing some examples of PH-distributions
is a parameter called coefficient of variation. This parameter is together
with the second moment, one of the most common measures of dispersion
of a probability distribution. It is often used in relation with the PHapproximation
due to its easier computations than with the second moment.
Definition 7. The coefficient of variation of a random variable X is defined as
cX =
√
E[(X−m1)2]
m1
, i.e. the ratio of the standard deviation and the mean m1.
Here we list some special types of PH-distributions:
7
2. FOUNDATIONS
1 2 n
1
λ λ λ λ
0
Figure 2.1: An n-phase Erlang distribution
Exponential distribution
An exponential distribution is represented as a PH-distribution with
p = (1) and T = (−λ).
Erlang distribution
An Erlang distribution was developed by A. K. Erlang to examine the
number of telephone calls, which might be made at the same time
to the operators of the switching stations. Nowadays, it is used in
the fields of stochastic processes and biomathematics. One important
property, which makes the class of Erlang distributions interesting
in our context, is that an n-phase Erlang distribution has the lowest
squared coefficient of variation among all n-phase continuous PHdistributions,
specifically c2
X = 1
n [14].
The Erlang distribution with n phases and parameter λ is a sequence
of n exponential distributions with the same rate λ, see Figure 2.1. The
Erlang distribution with n phases and parameter λ can be described
as a PH-distribution with
p = (1 0 ... 0) and T =




−λ λ 0 ... 0
0 −λ λ ... 0
: : : : :
0 0 ... 0 −λ



 .
Coxian distribution
Class of Coxian distributions is one of the most studied subclasses
of PH-distributions. Cox in 1955 [11] showed that any distribution
having a rational Laplace-Stieltjes transform (LST) can be represented
by a Coxian distribution. Further, Newton and Reddy [34] proved that
LTS of any distribution function can be approximated arbitrarily close
by a rational function.
A Coxian distribution with n-phases is a sequence of n exponential
distributions with rates λ1, λ2, ..., λn and probabilities a1, a2, ..., an−1
of bypassing the remaining phases, see Figure 2.2. Therefore, to
specify the Coxian distribution we need to get 2n parameters. Its
8
2. FOUNDATIONS
1 2 n
1
(1 − a1)λ1 (1 − a2)λ2 λn
a1 a2
0
Figure 2.2: An n-phase Coxian distribution
PH-distribution representation is
p = (1 0 ... 0) and T =




−λ1 (1 − a1)λ1 0 ... 0
0 −λ2 (1 − a2)λ2 ... 0
: : : : :
0 0 ... 0 −λn



 .
Hyper-Erlang distribution
A hyper-Erlang distribution with m branches is defined as a mixture of
m mutually independent Erlang distributions with initial probabilities
p1, ..., pm. The number of phases of the i-th Erlang distribution is
denoted with ni and its rate with λi, see Figure 2.3. This distribution
can be described by 3m parameters.
The PH-distribution representation of a hyper-Erlang distribution is
given by p = (p1 0 ... p2 0 ... pm 0 ... 0) and
T =













−λ1 λ1 0 ... 0 0 0 ... 0 0
0 −λ1 λ1 ... 0 0 0 ... 0 0
: : : : : : : : : :
0 0 0 ... −λ1 0 0 ... 0 0
0 0 0 ... 0 −λ2 λ2 ... 0 0
: : : : : : : : : :
0 0 0 ... 0 0 0 ... −λm λm
0 0 0 ... 0 0 0 ... 0 −λm













.
The PH-approximation of a distribution Q is an estimation of parameters
p and T of a PH-distribution, so as this PH-distribution is “close” to the
given distribution Q. The measurement of closeness is possible to choose
in many ways and for some of these possibilities see Chapter 5. For easier
computations it would be useful to reduce the number of free parameters.
However, the existence of a minimal canonical representation of n-phase PHdistributions
is an open research problem. Therefore, an intensive research
9
2. FOUNDATIONS
11 21 n1
λ1 λ1 λ1
λ1
0
1m 2m nm
λm λm λm
λm
p1
pm
Figure 2.3: A hyper-Erlang distribution
1 2 n
p1
λ1 λ2 λn−1 λn
p2 pn
0 < λi ≤ λi+1, 1 ≤ i ≤ n − 1
Σpj = 1, 0 ≤ pj ≤ 1, 0 < j ≤ n
0
Figure 2.4: Canonical representation of an n-phase APH-distributions
is targeted at the subclass called acyclic PH-distributions, which provides
this property.
Definition 8. An acyclic PH-distribution (APH-distribution) is a PHdistribution
with an upper triangular rate matrix (under some permutation of
its components). Therefore, R(i, j)=0 holds for all i > j.
There is a minimal canonical representation for any APH-distribution
[12], see Figure 2.4. This canonical representation reduces the number of
parameters to 2n for n-phase APH-distribution compared to n2 + n as we
had for general ones. The restriction to APH-distributions has also other
computational advantage over cyclic ones, because of the upper triangular
shape of the rate matrix of the underlying CTMC.
2.3 Demonstration Examples
Now we introduce a few examples of distributions, which will be used at the
end of the thesis for illustration of fittings by particular PH-approximation
methods.
10
2. FOUNDATIONS
Weibull distribution
Weibull distribution is one of the most widely used distributions as
a time-to-failure of microelectronic devices. It is characterized by two
parameters: scale parameter λ and shape parameter k, such that λ,
k > 0. The failure rate is proportional to a power of time. The parameter
k expresses the power plus one, and so for k < 1, the failure rate
decreases over time. For k = 0 it is constant and for k > 0 it increases
over time.
Weibull distribution is defined as a random variable X with pdf
fX(x) =
kxk−1e−(x
λ )
k
λk
for x > 0.
We will use Weibull distribution as a representant of “well-shaped
functions”, which means functions for which it is easy to find the
satisfactory PH-approximation. We will use the Weibull distribution
with parameters λ = 1 and k = 5, i.e. Weibull(1.0, 5.0).
Uniform distribution
The (continuous) uniform distribution is the distribution where all
subintervals of the same length on the distribution’s support interval
are equally probable. It is characterized by two parameters a, b, which
specify minimum and maximum values of the support interval. Its
pdf is given by
fX(x) =
1
a − b
for a ≤ x ≤ b.
We will use the uniform distribution with parameters a = 0.5 and b =
1.5, i.e. Uniform(0.5, 1.5), as a representant of “bad-shaped functions”.
Data sample of round-trip times
The round-trip time is the length of time it takes for a packet to be sent
plus the length of time it takes for an acknowledgment of that packet
to be received. We will use the data sample consisting of 1158 samples
of round-trip times, where the round-trip time was measured from
muni.cz to seznam.cz. It is the representant of situations, where we
have only sample data, not the formula for probability distribution.
The distribution given by this data sample occurs often in practice and
it is modeled mainly by a shifted gamma distribution.
11
Chapter 3
Moment Matching Method
As the first approach for the PH-approximation we consider the moment
matching method. This method is a member of fitting techniques that utilize
incomplete information of the original distribution. Techniques in this class
use very few parameters to perform approximation by a PH-distribution,
however, all the parameters are matched exactly. The most typically used
parameters are the first k moments of the original distribution. In this chapter,
we focus on the k-moment approximation techniques only and do not
discuss subsidiary techniques approximating e.g. values of complementary
cumulative distribution function at a given set of points [20].
3.1 Matching Less Than Three Moments
Matching only the first moment m1 is simply achieved by the exponential
distribution with rate λ = 1
m1
.
In techniques matching two moments, the coefficient of variation cX is
usually used instead of the second moment, because of consequently easier
representation of results. A technique for distributions with 0 < cX < 1
is shown in [47]. Those distributions can be approximate by a mix of two
Erlang distributions with k and k − 1 phases and same rate λ, so that 1
k ≤
cX ≤ 1
k−1, k ∈ N, k ≤ 3. More specifically, the resulting distribution is with
probability p (resp. 1 − p) the sum of k − 1 (resp. k) independent exponential
distributions with parameters
p =
1
1 + cX
2
kcX
2
− k(1 + cX
2
) − k2
cX
2
1
2
and λ =
k − p
m1
.
Distributions with cX
2 ≥ 0.5 can be approximated by a 2-phase Coxian
distribution with parameters λ1 = 2m1, p1 = 1
2cX
2 and λ2 = λ1 ∗ p1 [30].
Similarly to matching the first moment, matching the first two moments
is not very common for its oversimplification.
13
3. MOMENT MATCHING METHOD
3.2 Matching Three Moments
Designing techniques for matching three moments (similarly to the previous
techniques for one and two moments) is carried out by choosing a convenient
subset of the class of PH-distributions with a small number of parameters
and then fitting those parameters to get a suitable PH-distribution. This
narrows the state-space and thus improves the efficiency of the algorithm.
Matching of the first three moments is the most often used moment
technique. Thus, there have been introduced plenty of algorithms for it. We
would like to choose one that performs “best”. First of all we need to decide
what “best” means in this case. As mentioned in Chapter 1, for comparison
of the PH-approximation methods we use the following criteria: 1) accuracy,
2) time efficiency of the approximating algorithm, 3) minimality of the number
of phases and 4) generality.
We do not have to consider first two criteria in this case. Over all threemoment
techniques accuracy is the same, it matches exactly the first three
moments of the original distribution (for techniques matching the first three
moments and somehow providing better accuracy see Chapter 6 - Hybrid
methods). Also almost all three-moment algorithms have closed-form solution
for the parameters.
Osogami in [37] proposed a technique, which performs very well in the
both remaining criteria. His algorithm works for a whole class of distributions,
which can be approximated by an arbitrary three-moment matching
technique. It is also nearly minimal in the number of phases required. More
precisely, the number of phases required by his technique is at most n + 1,
where n is the minimal number of phases required in an arbitrary threemoment
matching technique using APH-distributions. Instead of moments
he works with normalized moments.
Definition 9 ([37]). The normalized second and third moment of a distribution
X are defined as
mN
2 =
E[X2]
m1
2
and mN
3 =
E[X3]
m1 · E[X2]
, respectively.
Osogami defined subset of PH-distributions called Erlang-Coxian distributions
(EC-distributions), see Figure 3.1. This class of distributions has
only 6 free parameters and it is broad enough to meet properties mentioned
above.
Definition 10. A PH-distribution is said to be an n-phase EC-distribution if
the value is zero with probability 1 − p and an (n − 2)-phase Erlang distribution
followed by a 2-phase Coxian distribution with probability p, for n ∈ N ≥ 2.
14
3. MOMENT MATCHING METHOD
E1 C2
p
λE λE λE (1 − a1)λ1 λ2
(1 − p)
E2 En−2 C1
a1
0
Figure 3.1: An n-phase EC-distribution
Three-moment approximation by a 2-phase Coxian distribution is possible
for any distribution that satisfies mN
2 > 2 and mN
3 > 3
2mN
2 [38]. However
a Coxian distribution requires many more phases for approximating distributions
with low second and third moments. In contrast, the class of Erlang
distributions is known to have the least normalized second moment among
all PH-distributions with a fixed number of phases [14]. By combining them
together Osogami gets distributions, which can fit three moments of distributions
with all ranges of variability and stay small in the number of phases
required [37].
As a reaction to Osogami’s technique Bobbio et al. in [6] proposed a threemoment
matching technique, which improves Osogami’s result by reducing
the number of phases required to exactly the minimal n needed in arbitrary
three-moment matching technique using APH-distributions. Instead
of EC-distributions they use Erlang-Exponential or Exponential-Erlang distributions,
which are defined as (n − 1)-phase Erlang distribution followed
(or preceded) by a single exponential one. Nevertheless, they perform very
similar.
One of the disadvantages of three-moment matching techniques is that
they are not applicable in general. There are distributions which are impossible
to approximate by a PH-distribution so that the first three moments
match exactly [25]. These distributions are characterized by mN
3 = mN
2 = 1.
Nevertheless, the class of distributions, that can be approximated by threemoment
matching techniques is very broad.
Another disadvantage is the number of phases required. Osogami in [39]
defined the set of distributions Sn, which can be approximate by an n-phase
APH-distribution, such that the first three moments match exactly. For characterization
of Sn he used another set of distributions Tn. He showed that
for n ∈ N, n ≥ 2 it holds Sn ⊂ Tn ⊂ Sn+1, where Tn is defined as follows.
Definition 11. For integers n ≥ 2, Tn denotes the set of distributions that satisfy
exactly one of the following two conditions:
15
3. MOMENT MATCHING METHOD
a) mN
2 > n+1
n and mN
3 ≥ n+2
(n+1)mN
2
b) mN
2 = n+1
n and mN
3 = n+2
n .
This shows that the number of phases required can be arbitrary high.
To perform an analysis of the resulting model, we need to keep the size of
the state-space of the corresponding CTMC small. So we may have some
upper bound for the number of phases that we can handle. Therefore, if the
resulting PH-distribution has more phases than this upper bound, there is
no way how to use three-moment matching technique in such situation.
It is well known that two distributions coincide if and only if all of their
moments coincide [13]. Therefore, in the PH-approximation it is desirable to
match as many moments as possible. However, equating large number of
moments makes the analysis enormously difficult. Therefore, matching three
moments is the most common trade-off between accuracy and efficiency of
corresponding algorithms.
3.3 Matching More Than Three Moments
Matching more than three moments has been also considered. Horvath
and Telek in [23] presented an iterative approach to match an arbitrary
number of moments with APH-distributions. An APH-distribution with n
phases can be fully characterized by its first 2n − 1 moments [48], therefore,
Horvath and Telek decided to match first 2n − 1 moments of the original
distribution with n-phase APH-distribution. This is done by one phase
reduction technique. The computational complexity increases exponentially
with n and therefore it can be calculated only up to 8 phases (15 moments)
with a personal computer with 1.5GHz processor and 524MB RAM. The
biggest disadvantage of this technique is that it provides solution only in
a situation, when the given 2n − 1 moments are moments of some n-phase
APH-distribution, which is often not the case.
16
Chapter 4
Minimization of a Difference
As the title indicates, this method is aimed at minimizing a difference between
an original and an approximating distribution. This is done according
to a given difference measure. In contrast to the moment method, in this
approach we first choose the number of phases. Then a technique based
on minimization of a difference finds the best approximation with less or
equal number of phases. There are two important steps in designing such
techniques. The first step comprises the choice of a distance measure of the
two distributions, so as it well captures the behavior important for us. The
second step involves the choice of an appropriate algorithm, which provides
a minimization of such distance.
4.1 Kullback-Leibler Divergence
The Kullback-Leibler divergence, also known as the relative entropy [27], is
the most common considered measure of a difference between two distribution
functions in PH-approximation techniques.
Definition 12. For continuous distributions P, Q with pdf p, q, respectively, the
Kullback-Leibler divergence (KL-divergence) is defined as the integral:
DKL(P||Q) =
∞
−∞
p(x) ln
p(x)
q(x)
.
The relative entropy was introduced by Solomon Kullback and Richard
Leibler in 1951. It is a non-negative measure of a difference between two
probability distributions. It is equal to zero if and only if the distributions
match exactly. This function is not symmetric, therefore it is not a distance
in the strict mathematical sense. The intuition behind the relative entropy
is that it reflects the expected number of additional bits required when encoding
samples from P using a code based on Q. This is also called the code
penalty.
The choice of the relative entropy as a difference measure is common
17
4. MINIMIZATION OF A DIFFERENCE
due to its close relationship with the likelihood theory. If the KL-divergence
is applied to the PH-approximation with empirical samples, the result corresponds
to the maximum likelihood estimation [42].
Definition 13. Suppose we have a sample of independent and identically distributed
observations x1, x2, ..., xn and distributions fθ dependent on parameters θ.
The maximum likelihood estimation (MLE) is a method of estimating parameters
θ, so as the log-likelihood function is maximized, where the log-likelihood function
is defined as:
L(θ|x1, x2, ..., xn) =
n
i=1
ln fθ(xi).
The MLE was vastly popularized by R. A. Fisher in 1920’s. Nowadays, it
is used for a wide range of statistical models, such as linear models, discrete
choice models or hypothesis testing. Due to its common use in statistics,
several techniques for the MLE have been developed.
In 1992, Bobbio and Cumani [3] used the MLE in the PH-approximation.
The MLE is done via non-linear programming in this technique. Due to its
computational complexity, only APH-distributions are considered as a result.
Asmussen et al. in [2] showed a technique for the PH-approximation based
on the MLE, which exploits an expectation-maximization algorithm (EM).
This algorithm can be applied to compute the MLE of any distribution from
the PH-class. Since then, the EM algorithm is the most common choice for
the PH-approximation.
The EM algorithm ([31, 15]) is a general iterative technique for finding
the MLE, where the original distribution can be viewed as incomplete data.
That means we have some observation X, some hidden or missing data Y
of a larger experiment and parameters θ, which are to be estimated. In the
PH-approximation X is the observed set of times to absorption. These observations
provide only incomplete information about the underlying CTMC,
because we do not know anything about the passage through the CTMC.
Therefore the missing data Y is a set of tuples (i1, ..., im−1, s1, ..., sm−1)l for
every xl ∈ X, where m is the number of jumps before absorption, ik are
states of the Markov chain and sk are the sojourn times, ∀k ∈ {1, 2, ..., m−1}.
As the missing data corresponds to the incomplete one, the sojourn times
must satisfy xl = (s1 + ... + sm−1)l. Parameters θ are parameters of PHdistributions,
where the number of phases is chosen by user.
Each iteration of the EM algorithm consists of an expectation step (E-step)
followed by a maximization step (M-step). In the n-th E-step, the missing
data Yn are estimated given the observed data X and current estimate of
18
4. MINIMIZATION OF A DIFFERENCE
parameters θn−1. This is achieved using the conditional expectation. Then
in the n-th M-step, θn are estimated, such that the likelihood function is
maximized under the assumption that the missing data Yn are known from
the n-th E-step. For more information see [2].
Convergence is assured since the algorithm is guaranteed to increase
the likelihood at each iteration [8]. However, the ending stationary point
need not be a global maximum as was shown in [49]. Therefore, running the
procedure from multiple initial parameters is often helpful.
One of the advantages of this algorithm is that it preserves zeros in T and
p. That means, once an element has been estimated to be zero, it remains
zero during the whole computation. Hence if one wants the result in some
subclass of PH-distributions having some elements fixed to zero, one only
needs to start with such zeros in T and p. Therefore, we can limit the algorithm
only on subclasses such as APH-distributions, Erlang distribution, or
hyper-Erlang distribution. The main advantage of using the special structures
is that it reduces the computational complexity and therefore, with the
same amount of time, we can work with a larger number of phases. Another
useful property of Asmussen’s algorithm is that it preserves exactly the
first moment of the original distribution. The main disadvantage is its high
computational complexity in comparison with other techniques, even if we
are restricted to some subclass of PH-distributions. Moreover, it is possible
to end in a saddle point or a local maximum [49].
EMpht1 is a program for the PH-approximation using the Asmussen’s
algorithm [36]. It can be used to fit a PH-distribution to a sample as was
described above. This program can also make the PH-approximation of
continuous distributions. This is done by discretization into a weighted
sample, where the discretization interval is specified by user. The time
of computations, when the number of phases increases above 6, is
several hours on a computer with Intel Core i7-920 2,67GHz and 6GB
RAM.
Asmussen’s algorithm turns out to be extremely costly. Therefore, there
are several results improving it. One of them is algorithm introduced by
Thummler, Buchholz, and Telek [46]. The way to optimize the time complexity
of the EM algorithm for PH-distributions is a restriction to hyper-Erlang
distribution. This subclass is theoretically as powerful as acyclic or general
PH-distributions. However, it still allows the realization of a very efficient
fitting algorithm. In fact its time complexity is in O(mk), where m is the
1. http://home.imf.au.dk/asmus/pspapers.html
19
4. MINIMIZATION OF A DIFFERENCE
number of branches of the hyper-Erlang distribution and k is the number of
samples. Thummler et al. also suggest several optimizations, which can be
used in special cases.
GFIT2 is a program for the PH-approximation based on Thummler, Buchholz,
and Telek’s work [46]. The quality of the fit is in most cases
almost the same as in the general case in EMpht, but it is much faster.
To overcome the problem with saddle points and local optima, there
is a progressive pre-selection in GFIT, which provides best aspirants
for initial value. However, this can make an exponential increase of
computation time as the number of phases increases.
Yet another improvement of the Asmussen’s EM algorithm was proposed
by Okamura et al. [35], where the lower bound on computation cost of the
algorithm is achieved by applying uniformization. It results in general or
acyclic PH-distributions. This approach is implemented in a tool named
PHPACK.
PHPACK3 is another program for the PH-approximation based on the
EM algorithm and it is integrated in mapfit package4. So as the PHapproximation,
this program also provides computation of pdf, cdf and
moments of PH-distributions, which is based on simple uniformization.
This program allows us to use both data samples and probability
distributions.
The last approach based on the EM algorithm we would like to mention is
the one proposed in [26]. This technique uses a generalization of the class of
PH-distributions called Matrix-Exponential (ME) distributions. These distributions
have one important advantage over the PH-distributions, namely,
they have a minimal representation. On the other hand, there is a need for
dealing with results which are not in the class of PH-distributions.
4.2 Other Distance Measures
Naturally, it is also possible to use other difference measures than the KLdivergence.
Some of them are listed bellow, for continuous distributions
with cdf P, Q and with pdf p, q, respectively:
2. http://ls4-www.cs.tu-dortmund.de/home/thummler/
3. http://www.rel.hiroshima-u.ac.jp/okamu/PHPACK/intro.html
4. http://www.rel.hiroshima-u.ac.jp/˜okamu/index.php/software
20
4. MINIMIZATION OF A DIFFERENCE
Pdf area difference is defined as
∞
0
(p(x) − q(x))2
dx.
This measure is used in mathematical procedure named least squares
(LS) for fitting a curve to the observed sample. Some authors use
non-negative square root definition
∞
0
|p(x) − q(x)|dx.
Cdf area difference is defined in similar way as the pdf area difference:
∞
0
(P(x) − Q(x))2
dx.
Weighted area difference is used, when we want to give more stress to
some part of the function, e.g. to capture well the behavior of the tail
of a distribution. Therefore, for the function w(x), which associates
a weight with each x, the weighted area difference is defined as
∞
0
w(x)(p(x) − q(x))2
dx.
Maximum absolute difference between two distributions P, Q, also
known as the total variation distance, is defined as
max
x∈R+
0
|p(x) − q(x)|.
Faddy in [18] and [19] used LS to fit Coxian distributions to real sample
data. By this technique he estimates parameters for models used in drug
kinetics. Also Bux and Herzog in [9] considered minimization of a difference
resulting in Coxian distributions, their chosen distance measure was the
maximum absolute difference. Malhotra and Reibman in [29] suggested LS
technique for the PH-approximation using already developed minimization
algorithms. The most recent algorithm, moreover, with best performance,
is technique introduced by Horvath and Telek [21]. This technique is also
implemented in a tool PhFit [22].
PhFit is a PH-fitting tool designed by Horvath and Telek. Its resultant
distributions are restricted to the class of APH-distributions to narrow
21
4. MINIMIZATION OF A DIFFERENCE
the search space. PhFit allows us to choose from three predefined
distance measures: the relative entropy, the cdf area difference, and the
squared root version of pdf area difference. The minimization is done
via non-linear optimization and simplex method. Initial parameters
are selected from 1000 randomly generated samples.
Use of other difference measures than the KL-divergence is often in
techniques combining the moment method and the distance minimization
method, e.g. Meda [41]. For these techniques see Chapter 6 - Hybrid
methods.
The main advantage of the PH-approximation method based on
minimization of a difference is that, unlike in case of the moment method,
we can choose the number of phases. Therefore, if we have some limitations
on the number of phases in our model, this method allows us to handle such
restriction and finds the solution satisfying it. Other major advantage in
comparison with the moment method is that it can capture the tail behavior
of a function in a better way [29].
The main disadvantage is that this method is computationally more
expensive than the moment method. This is problematic in practice mainly
if we do not want to restrict the result to some subset of PH-distributions.
Further none of these techniques guarantee ending in global optimum.
However, the chance of ending in local optimum and saddle points is low.
22
Chapter 5
Discretization and Discrete Phase-Type Approximation
The term discretization means transformation of a continuous distribution
into a discrete one. A discrete PH-approximation captures a discrete probability
distribution by a discrete PH-distribution. Therefore, to use the discrete
PH-approximation and exploit the properties of discrete PH-distributions in
a model with continuous probability distributions, we first have to discretize
such continuous distributions. After that the PH-approximation by discrete
PH-distributions can be done.
5.1 Discrete Phase-Type Distributions
To define discrete PH-distributions we need to introduce discrete-time
Markov chains.
Definition 14. A discrete-time Markov chain (DTMC) is a tuple (S, P, v0),
where
• S is a finite set of states,
• P : S × S → [0, 1] is a matrix of transition probabilities, such that ∀i ∈
{0, ..., |S| − 1} :
|S|−1
j=0
pij = 1 and
• v0 is an initial probability vector.
The intuitive semantics is as follows. The process starts in a state of S
with probability given by the initial vector v0. Then it moves successively
from one state to another. If the chain is currently in si, then it moves to
sj with probability P(i, j) from the transition matrix P. Similarly to the
continuous case we can define absorbing DTMC.
Definition 15. A state is called absorbing if it is impossible to leave it. A DTMC
is called absorbing if it has at least one absorbing state and it is possible to reach
23
5. DISCRETIZATION AND DISCRETE PHASE-TYPE APPROXIMATION
some absorbing state from every state. In an absorbing DTMC, a state which is not
absorbing is called transient.
Now we have all notions needed for the definition of discrete PH-
distributions.
Definition 16. A discrete n-phase PH-distribution (DPH-distribution) is
defined as the time until absorption in a DTMC with n + 1 states, where n ≥ 1,
such that states 1, ..., n are transient and state 0 is an absorbing state and the initial
vector is given by (p0 p).
Similarly to the continuous case, this distribution is described by the
initial probability vector (p0 p) and the probability transition (n+1)×(n+1)
matrix
P =
1 0
tT T
,
where T is an n × n matrix. As P is the matrix of transition probabilities of
an absorbing DTMC, it holds that
T1T = 1T − tT and p1T = 1 − p0.
Therefore the DPH-distribution is fully specified by T and p. For its probability
mass function (pmf) and cdf, the following proposition holds.
Proposition 3 ([4]). Assume that X is an n-phase DPH-distribution with parameters
T and p. Its cdf is given by
F(k) =
k−1
i=0
pTi
tT
= 1 − pTk
1T
for k ∈ N0
and its pmf is given by
f(k) = pTk−1tT for k ∈ N.
DPH-distributions as well as continuous PH-distributions were introduced
and formalized in [33]. However, they have received little attention
in applied stochastic modelling since then, because the main research activity
was oriented towards continuous models and therefore towards the
continuous PH-approximation. In recent years, the discrete models have
drawn considerable attention. Since DPH-distributions have some useful
properties in contrast to continuous PH-distributions, there have been introduced
techniques, that profit from DPH-distributions in continuous cases. In
[4], authors provide detailed study of the class of DPH-distributions and its
24
5. DISCRETIZATION AND DISCRETE PHASE-TYPE APPROXIMATION
properties. DPH-distributions inherit many properties from the continuous
PH-distributions. Properties, in which DPH-distributions are different, are
the most important for us.
One important well-known result for the continuous PH-distributions is
that the squared coefficient of variation c2
X cannot be less than 1
n , where n
is the number of phases [14]. The equation is reached by an n-phase Erlang
distribution and this limiting case is independent to the mean. It is easy to
show that this relation does not hold for the class of DPH-distributions, since
there is always at least one n-phase DPH-distribution with cX = 0 (this is
achieved by deterministic timeouts). Unlike the continuous PH-distributions,
the minimal coefficient of variation for DPH-distributions depends on the
mean.
Proposition 4 ([44]). The minimal squared coefficient of variation c2
X of a DPHdistribution
X with n phases and mean m1 is:
c2
X =
R(m1)(1−R(m1))
m2
1
if m1 ≤ n
1
n − 1
m1
if m1 > n,
where R(x) = x − x and x stands for the integer part of x.
Another useful property of DPH-distributions is that they can exactly
represent deterministic values and accurately approximate finite-support distribution
functions such as uniform distributions, which CPH-distributions
cannot.
Similarly to the continuous case, there is a subclass called acyclic DPH-
distributions.
Definition 17. A DPH-distribution is called acyclic (ADPH-distribution) if
its states can be ordered in such a way that matrix T is an upper triangular matrix.
Moreover, as shown in [4], this subclass has also a unique minimal
canonical representation, see Figure 5.1.
5.2 Discretization of Continuous Distributions
Naturally, we can use the discrete PH-approximation in models where we
already have a discrete sample of data points. We can also profit from different
properties of DPH-distributions in models with continuous probability
distributions. In such cases, first of all, we need to discretize the continuous
distribution using a given discretization step and then generate discrete
25
5. DISCRETIZATION AND DISCRETE PHASE-TYPE APPROXIMATION
1 2 n
p1
q1 q2 qn−1 qn
p2 pn
1 − q1 1 − q2 1 − qn
Σpj = 1, 0 ≤ pj ≤ 1, 0 < j ≤ n
0 < qi ≤ qi+1 ≤ 1, 1 ≤ i ≤ n − 1
0
Figure 5.1: Canonical representation of an n-phase DPH-distribution
samples with the associated mass probabilities. After that we are ready for
the discrete PH-approximation.
Discretization introduces errors mostly related to the size of chosen
discretization step. Therefore, it is desirable to pay some attention to its
choice. In general, the smaller discretization interval is the better accuracy is
achieved. On the other hand, the discretization interval changes the scale
of the representation. Therefore, to provide a good fit, it should hold that
for as many i as possible, E(Xi) is close to δiE(Xi
δ) for the discretization
interval δ, discretized distribution Xδ, which is expressed in δ units, and
the original distribution X [4]. So the discretization changes the mean, but
almost preserves the coefficient of variation. As mentioned in the previous
section, the minimal coefficient of variation of an n-phase DPH-distribution
is a function of the mean and therefore the choice of discretization step has
important effect on the variability of possible DPH-distributions. Due to the
choice of a smaller δ, more phases are needed to achieve the same spread of
a distribution.
Authors in [4] provided lower and upper bounds on a good discretization
interval for the PH-approximation resulting in n-phase ADPH-distributions.
These bounds mainly come from the dependence of the minimal coefficient
of variation cX and mean m1. The upper bound is expressed by m1
n−1. For the
lower bound it holds that
• if cX > 1
n, then any small value of δ is suitable,
• for cX ≤ 1
n , in order to allow the ADPH-distribution reach such small
cX, δ should be greater than ( 1
n − c2
X)m1.
After the discretization step was chosen, a discretization can be performed.
The most natural and common one is when the elements are integer multiples
of the discretization step δ. One possible rule for discretization of
26
5. DISCRETIZATION AND DISCRETE PHASE-TYPE APPROXIMATION
a continuous function X with cdf FX over the discrete set S = {x1, x2, ...} is
then the following:
pi = FX
xi + xi+1
2
− FX
xi−1 + xi
2
for i > 1 and p1 = FX
x1 + x2
2
where pi is the probability associated with each xi.
5.3 Discrete Phase-Type Approximation
Techniques for the discrete PH-approximation correspond to the ones for the
continuous PH-approximation. Three-moment matching method restricted
to 2-phase ADPH-distributions was introduced in [45]. Further, the bounds
for the moments of a distribution have been given, such that it is possible
to approximate the distribution by a 2-phase ADPH-distribution. In [4] the
EM-algorithm for the discrete PH-approximation was introduced. In order
to reduce the complexity, the technique is restricted to the class of ADPH-
distributions.
In [24] general technique for the discrete PH-approximation based on
minimization of a difference was presented. The minimization is viewed as
optimization problem and authors consider several algorithms for it.
PhFit [22], as mentioned in the previous chapter, is a tool for the PHapproximation
based on the minimization of a difference. This tool
also provides the discrete PH-approximation in the same way as the
continuous PH-approximation. The predefined difference measures
for DPH-distributions are the cdf area difference, a discrete version of
the pdf (pmf) area difference, and the relative entropy.
Comparative results on the use of DPH-distributions and continuous PHdistributions
in models are provided in [5]. A more suitable class of distributions
is chosen based on the optimization of a parameter called scale factor.
This parameter describes the time associated with each step in a DTMC.
Therefore, a DPH-distribution converges to a continuous PH-distribution
with the same number of phases as the scale factor tends to zero.
The advantages of the PH-approximation method based on discretization
are mainly the following:
• better fitting of distributions with small coefficient of variation,
• easier combination with constant durations, which are often needed in
models, and
27
5. DISCRETIZATION AND DISCRETE PHASE-TYPE APPROXIMATION
• easier computation, since this method works with sums instead of
integrals.
We also list the main disadvantages, which were considered in [4]. They are
• bounds on moments for n-phase ADPH-distributions,
• limited number of waves of the pmf, more specifically not greater than
n
2 , where n is the number of phases and
• no possibility to capture sharp changes of the pmf at high ( n) time
instances, since the sharpest possible change of the pmf at time k
(k > n) is achieved by a discrete Erlang distribution.
Using DPH-distributions in a continuous model has another very important
disadvantage, that is the need of discretization of the entire model. This
can lead to a problem with coincident events. Furthermore, it is necessary
to unify discretization intervals of different distributions in the model or
to transform the underlying DTMC of resulting DPH-distributions in such
a way, that they have the same scale. The unification of discretization intervals
can lead to significant imprecision, since the role of the discretization
interval on the discrete PH-approximation is crucial. On the other hand,
transformation of resultant DPH-distributions can lead to a high increase of
the state-space.
28
Chapter 6
Hybrid Methods
Hybrid methods combine previous techniques for the PH-approximation
to achieve better accuracy or to overcome limitations of particular methods.
Smickler in 1992 [41] was one of the first, who considered a hybrid method.
His program Meda (Mixed Erlang Distributions for Approximation) for
a given set of samples provides a PH-distribution with the first three moments
matched exactly. Moreover, Meda minimizes the cdf area difference
via a non-linear programming and the result of this method is restricted
to the class of hyper-Erlang distributions. This method achieves little bit
better accuracy, on the other hand it still have the main disadvantages of
both methods. For the moment matching method it is the fixed number of
phases and bounds on moments. For the method based on the minimization
of the cdf area difference it is a higher time complexity.
A similar method was introduced by Reibman and Malhotra in [29].
It combines the moment matching method and the LS. It is implemented
for lognormal distribution in extension to the modelling toolkit SHARPE,
called GSHARPE. Results are restricted to a hyper-Erlang distribution with
2 branches. This method has the same disadvantages as Meda.
A different combination of the moment matching and minimization of
a difference is considered in [46]. The authors introduced moment matching
technique for hyper-Erlang distributions with m branches, so as matching
first three moments requires adjusting only two of the m Erlang branches.
Therefore, there are m(m−1)
2 possibilities for the choice of the two branches.
Therefore, to choose the best one, one can use some specified difference
measure. In this case a likelihood was used. The authors also provided some
other strategies for combining the moment matching and the MLE using
their tool GFIT.
Yet a another approach in hybrid techniques is introduced in [16]. The
authors suggested a method for the PH-approximation based on matching
of moments. This method overcomes bounds on moments in the moment
matching method by permitting deviations from the specified k moments.
29
6. HYBRID METHODS
In other words, it matches the first k moments as close as possible. This is
done via mathematical programming and it can be viewed also as a method
minimizing a difference, for the difference measure chosen as the sum of
weighted differences in the first k moments. This method allows to match
moments even if the moment matching method is unable to provide a result.
6.1 Heavy Tailed Distributions
The attention to the tail behavior of a distribution is devoted due to a common
occurrence of heavy tailed distributions and their importance in recent
telecommunication systems.
Definition 18. A random variable X with cdf FX is said to have a heavy tail if
for all λ > 0 :
lim
x→∞
eλx
(1 − FX(x)) = ∞
Intuitively, a distribution is heavy tailed if its tail decays slower than exponentially.
However, common problem with PH-approximation techniques
is that they either provide a good fit of the body of the distribution and do
not pay enough attention to the tail (e.g. EM methods), or they prioritize the
tail over the body (e.g. minimization of weighted difference measures). To
deal with such limitations, it is natural to consider a hybrid method, which
splits the distribution into the body and the tail and then uses different
techniques for each of them.
Feldmann and Whitt in [20] proposed very effective method for
monotony decreasing functions, which is well applicable for fitting heavy
tails. A heavy tailed distribution is approximated by a finite mixture of
exponentials using a recursive algorithm. The tail of all PH-distribution
decays exponentially [43]. Therefore, of course the tail is exponential also in
the result of this method, but it allows to match longer part of the original
decreasing tail.
This result was exploited by Horvath and Telek in the tool PhFit [22]. In
[21], they introduced hybrid method which minimizes an arbitrary distance
measure for the body approximation and Feldmann and Whitt’s method for
the approximation of the tail. The point of split is chosen by user, but it has
to satisfy the condition on a monotony decreasing tail.
Hybrid techniques provide interesting ideas how to improve previous
methods and overcome their limitations. However, in practice the only really
interesting hybrid method is the last mentioned method implemented in
PhFit, which provides better fit of the tail behaviour of a distribution.
30
Chapter 7
Comparison
It is not possible to choose the best approach for the PH-approximation in
general. For the comparison of PH-approximation methods we use four
criteria:
1. Accuracy
This is probably the most natural and also most often considered
criterion. The goal of the PH-approximation is to find a PH-distribution
that can replace an original general distribution in a model. In order to
preserve parameters of the model as much as possible, it is desirable
to get the PH-distribution, which is as close as possible to the original
distribution.
2. Time efficiency of the approximating algorithm
To apply the PH-approximation in a model it is preferable that the
running time of the approximating algorithm is as short as possible.
Furthermore, the computational time of some techniques is affected
by the number of phases required. As a result, some techniques with
higher computational time complexity cannot provide result, whenever
the number of phases required is high, see Section 7.2.
3. Number of phases
The PH-approximation is used to allow Markovian analysis of the
resulting model. The complexity of such analysis depends on the size
of the state-space of the underlying Markov chain, therefore we need to
keep the number of states low. Therefore, the influence of the number
of phases used in the PH-distribution on the size of the state-space is
crucial. On the one hand, in sequential formalism such as SMP, the
state is first chosen and then the corresponding waiting time starts
running. Therefore the PH-approximation increases the SMP statespace
linearly. On the other hand in parallel formalism such as GSMP,
all the waiting times corresponding to the available states from the
current state start and we move to the state with the first finished
31
7. COMPARISON
waiting time. Therefore, in cases containing two or more outgoing
non-exponential distributions from one state in a GSMP generating
only PH-distributions is not sufficient. To capture the behavior of the
original GSMP, it is needed to generate cross-products of matching
PH-distributions and this can lead to an exponential increase of the
state-space, so-called state-space explosion. Thus, mainly in such cases,
it is desirable to generate matching PH-distribution with very few
phases.
4. Generality
Ideally the algorithm should work for as many models as possible.
7.1 Accuracy
Comparison of PH-approximation approaches by accuracy of their fittings
is not easy to do. To say which PH-distribution is closer to the original
distribution, we have to define some measure of proximity. We can use all
measures mentioned in Chapter 4, but we can define it also as the sum of differences
in moments or e.g. some combination of all of them. To standardize
the question of accuracy in the PH-approximation, there was a benchmark
developed at the workshop on “Fitting phase type distribution”, organized
by Asmussen in Aalborg, 1991. A reformulated version of the benchmark
was introduced in [7]. It provides representative distributions for fitting and
relevant measures, which are
• relative errors in the first three moments,
• pdf area difference, and
• minus the cross entropy, which is defined as
∞
0
log p(x)dQ(x), where
p(x) is pdf of the approximating PH-distribution and Q(x) is cdf of the
original distribution.
This benchmark was adopted by many authors in the field of PHapproximations
to compare their algorithm with others. The results are
not surprising. The moment matching method performs best in relative
errors in moments, on the other hand the techniques based on the MLE
perform best in the cross entropy measure, since it is equal up to the constant
to the likelihood. Techniques based on the minimization of cdf or pdf area
difference often perform best in the pdf area difference measure. Therefore, it
is not possible to say in general, which one perform the best. That is heavily
32
7. COMPARISON
dependent on features we want to compute from the model afterwards.
Another way to indicate quality of the PH-approximation is by density
plots. To provide a good fit, an algorithm should perform very well in most
of the measures from the benchmark and this can be roughly seen from
the plot. Results are often very similar for “well-shaped functions”, where
all of the PH-approximation methods perform very well and require small
number of phases. Distributions with small coefficient of variation and sharp
changes need high number of phases to be satisfactory approximated. The
only method, which is very different in visual inspection, is the method
based on the discrete PH-approximation, see Section 7.5.
However, there are specific situations when it is important to choose a
particular method to reach satisfactory accuracy such as models with heavy
tailed distributions. The hybrid technique for such distributions, which is
also implemented in PhFit, has already been mentioned. This technique
allows to better approximate the tail. But even if we do not desire good
approximation of the tail, we need to be careful with the choice of the right
technique. In particular, the moment method is heavily dependent on the
tail of a distribution and for heavy tailed distributions the results can be
unsatisfactory in the body part as well as in the tail part. To prevent this, one
can first preprocess the input data by cutting of the heavy tail.
Another situation, where it is very helpful to choose a particular method
is in models, where we can exploit different properties of DPH-distributions.
That is in models with distributions with small coefficient of variation and
sharp changes, finite-support distributions and deterministic values.
7.2 Time Efficiency of the Approximating Algorithm
With respect to this criterion, the moment matching techniques perform best,
since they mostly have closed-form solutions. This is also the main disadvantage
of the method based on minimization of a difference. Empht tool for the
PH-approximation based on the minimization of a difference on personal
computer (Intel Core i7-920 2,67GHz and 6GB RAM) could not provide
results for more than 10 phases within a day. More recent tool GFIT could
not provide results for more than 200 phases. The 200-phases continuous
PH-approximation and 400-phases discrete PH-approximation took several
hours in PhFit. Therefore, it is not recommended to use them in situations,
when good approximations need high number of phases (and we can handle
so many phases in the model). In such situations it is recommended to use
moment matching techniques and this is also the only situation, when we
recommend them.
33
7. COMPARISON
7.3 Number of Phases
Unlike the algorithm efficiency, number of phases is the main advantage
of other methods over the moment matching approach. In the moment
matching method, the result is guaranteed to be minimal in the number
of phases, such that the first k moments are equal. However, this minimal
number can be arbitrarily high. In most cases we need to keep this number
small to get a satisfactory size of the state-space of the underlying Markov
chain. In such cases it is often impossible to find a solution in the moment
matching method. In other approaches, we can first choose the number of
phases and then search for the most accurate result satisfying the choice.
However, for a small number of phases, the resulting approximation can
still be unsatisfactory.
7.4 Generality
The only restriction on the approximated distribution appears in the moment
matching method. To provide a solution, the original distribution has
to satisfy bounds for its moments. For the three-moment technique it is
mN
3 > mN
2 > 1, where mN
3 and mN
2 stand for the third and the second normalized
moments, respectively. Other approaches work for arbitrary general
distribution.
Yet another restriction is on the form of such distribution, since the
method based on the MLE requires data sample as input. Therefore, there is
a need to create the data sample from the original continuous distribution
and bring thus some other errors to the approximation.
The last limitation is for the discrete PH-approximation in a model, since
there is a need of possibly unwanted discretization of the entire model.
7.5 Demonstration Examples
Figures 7.1, 7.2, and 7.3 show density plots of our demonstration examples
introduced in Chapter 2 and the resultant distributions of several PHapproximation
methods. They are the following:
1) Osogami’s technique, which represents the three-moment matching
method,
2) results provided by tool GFIT based on the EM algorithm,
34
7. COMPARISON
3) results provided by tool PhFit for APH-distributions, which implements
the method based on minimization of an arbitrary distance measure, and
4) PhFit results using ADPH-distributions.
As a distance measure in PhFit, the relative entropy is chosen, since it provides
the best results. The number of phases used in (a) is 6 and these plots
show, how particular methods work with very few phases. There is no possibility
to reduce the number of phases required in Osogami’s technique and
therefore, in our case there are no results of Osogami’s technique using 6
phases. In (b) the number of phases is determined by Osogami’s technique
and so we can compare the moment matching method with the others.
Figure 7.1 compares PH-approximations of the Weibull(1.0, 5.0) distribution.
Osogami’s technique results in PH-distribution with 22 phases and
discretization intervals for the discrete PH-approximation are δ = 0.1 in (a)
and δ = 0.07 in (b). Weibull(1.0, 5.0) distribution is “well PH-approximable”
function and, as can be seen from the plots, for more phases all of the mentioned
methods perform very well.
Figure 7.2 shows PH-approximations on the Uniform(0.5, 1.5) distribution.
Osogami’s technique results in 16 phases and for the discrete PHapproximation
the discretization interval is in both cases δ = 0.1. Except for
the discrete PH-approximation, none of methods is able to capture sharp
changes of the distribution.
Figure 7.3 shows PH-approximations on the data sample of round trip
times with 1158 samples. Osogami’s technique results in 33 phases and discretization
intervals for the discrete PH-approximation are δ = 1.0 in (a) and
δ = 0.2 in (b). Similarly to the approximation of the uniform distribution,
the discrete PH-approximation performs the best and all the other methods
perform almost similar in the visual inspection.
7.6 Guideline for the Choice of Suitable Method
In Figure 7.4 we provide a decision diagram, which should work as a guideline
for the choice of a suitable PH-approximation method. To keep the
diagram small this guideline is given only for the most usual situations.
There are situations, when we would like to get as good approximation as
possible, but we do not have strictly specified requirements on the proximity
measure. In other words, it does not cover specific situations, when we want
to, e.g. match exactly first three moments, or provide good approximation of
a heavy tail. These specific situations are discussed in detail in the previous
chapters.
35
7. COMPARISON
(a) 6 phases (b) 22 phases
Figure 7.1: Weibull(1.0, 5.0)
(a) 6 phases (b) 16 phases
Figure 7.2: Uniform(0.5, 1.5)
(a) 6 phases (b) 33 phases
Figure 7.3: Data sample of round-trip times
36
7. COMPARISON
Can we
discretize the
entire model?
Do we need a
lot of phases
and can we
handle them?
Are there only
“well-shaped”
functions in
the model?
Do we need a
lot of phases
and can we
handle them?
Minimization
of a difference
Moment
method
Discrete PH-
approximation
Minimization
of a difference
Moment
method
no yes
yesno no yes
yesno
Figure 7.4: Guideline for the choice of suitable method
The guideline is based on observations we made, while trying different
PH-approximation techniques. First of all we need to decide if we can handle
discretization of the entire model. If the answer is no, we do not have to
consider the discrete methods and therefore, the only remaining decision is
whether to use the moment method or the method based on minimization
of a difference.
From our experience, Osogami’s technique—and therefore also all well
working moment techniques—provides results with plots very similar to
the plots of techniques based on minimization of a difference with the same
number of phases. For smaller number of phases the moment method cannot
provide a result and it is possible to use the method minimizing a difference.
For higher number of phases, the method minimizing a difference provides
37
7. COMPARISON
better result than the moment method and therefore, it is again preferable
to use it. Because of that, the use of the moment method is recommended
only in the case, when the moment method needs “a lot of phases” and
we can handle so many phases. “A lot of phases” means such a number of
phases that the method based on minimization of a difference cannot provide
a result for it and therefore the only available solution for this number of
phases is the one of the moment method. The number of phases, say n,
needed in Osogami’s technique is given by
n =



pmN
2
pmN
2 −1
if pmN
3 = 2pmN
2 − 1, and pmN
2 ≤ 2,
pmN
2
pmN
2 −1
+ 1 otherwise,
where p is given by
p =



(mN
2 )2+2mN
2 −1
2(mN
2 )2 l if mN
3 > 2mN
2 − 1, and 1
mN
2 −1
is an integer,
1
2mN
2 −mN
3
if mN
3 < 2mN
2 − 1,
1 otherwise.
Naturally, we also have to check the bounds on moments in the moment
matching method. If it cannot provide the result, then we should use the
method based on minimization of a difference for the number of phases as
high as possible.
If the first answer is yes and we can thus use the discrete method, then
we have to differentiate the moment method, the method based on minimization
of a difference and the discrete PH-approximation. Therefore, there
is one more decision in such a situation. If the distribution is “well-shaped”
then there is no reason to use other methods than the method based on
minimization of a difference. This is due to the fact, that the moment method
in such situations need a small number of phases, therefore it is better to use
the method minimizing a difference and to have the possibility to regulate
the number of phases. Further, there is no need to use the discrete method
and thus to complicate the situation with discretization of the entire model,
since the method minimizing a difference provides very similar and satisfactory
results.
Otherwise, when the distribution has small coefficient of variation, sharp
changes, finite support or if we have deterministic timeouts in the model,
then we can exploit different properties of the discrete method. Again, we
have to consider the suitability of the moment method. This is done in a way
similar to the to decision between the moment method and the method minimizing
a difference, however “a lot of phases” is in this case such a number,
38
7. COMPARISON
that the discrete method cannot provide result for it. And for smaller number
of phases the discrete method cannot provide satisfactory result.
39
Chapter 8
Conclusion
The aim of the thesis was to provide a summary of the PH-approximation
techniques and to describe their suitability for particular models. We provided
a detailed description of the main approaches. In particular, we discussed
the moment method, the method based on minimization of a difference,
the method based on discretization and hybrid methods. We gave
an overview of most relevant techniques and tools, which can be used for
such methods and illustrated results of approximations by these methods
on three examples. This thesis also includes a guideline for the choice of
suitable method, which is based on our experience and observations, while
trying different PH-approximation techniques.
For the comparison we considered four criteria, namely accuracy, time
efficiency of the approximating algorithm, number of phases, and generality.
We realized that we are not able to decide in general which one is the
best, since the results are dependent on the particular situation and each
of the methods suits best for some of them. However, the method based
on discretization performed very well with respect to all the criteria and
provided the best results considering the accuracy of the approximation of
a probability distribution. Its main disadvantage is the need of discretization
of the entire model. More detailed description of problems related with the
consequential discretization of the model and characterization of models,
where the discretization is appropriate, is left as the future work.
41
Bibliography
[1] R.B. Abernethy, J.E. Breneman, C.H. Medlin, and G.L. Reinman. Weibull
Analysis Handbook. Aero Propulsion Laboratory, Wright-Patterson
AFB Ohio, 1983.
[2] S. Asmussen, O. Nerman, and M. Olsson. Fitting phase-type distributions
via the EM algorithm. Scandinavian Journal of Statistics, pages
419–441, 1996.
[3] A. Bobbio and A. Cumani. ML estimation of the parameters of a
PH distribution in triangular canonical form. Computer Performance
Evaluation, 22:33–46, 1992.
[4] A. Bobbio, A. Horv´ath, M. Scarpa, and M. Telek. Acyclic discrete phase
type distributions: properties and a parameter estimation algorithm.
Performance Evaluation, 54(1):1–32, 2003.
[5] A. Bobbio, A. Horv´ath, and M. Telek. The scale factor: a new degree
of freedom in phase-type approximation. Performance Evaluation,
56(1):121–144, 2004.
[6] A. Bobbio, A. Horv´ath, and M. Telek. Matching three moments with
minimal acyclic phase type distributions. Stochastic models, 21(2-3):303–
326, 2005.
[7] A. Bobbio and M. Telek. A benchmark for PH estimation algorithms:
results for Acyclic-PH. Stochastic models, 10(3):661–677, 1994.
[8] S. Borman. The Expectation Maximization Algorithm-A short tutorial.
2004.
[9] W. Bux and U. Herzog. The phase concept: Approximation of measured
data and performance analysis. Computer Performance, pages 23–38,
1977.
[10] C.G. Cassandras and S. Lafortune. Introduction to discrete event systems,
volume 11. Kluwer academic publishers, 1999.
43
8. CONCLUSION
[11] D.R. Cox. A use of complex probabilities in the theory of stochastic
processes. In Mathematical Proceedings of the Cambridge Philosophical
Society, volume 51, pages 313–319. Cambridge Univ Press, 1955.
[12] A. Cumani. On the canonical representation of homogeneous Markov
processes modelling failure-time distributions. Microelectronics and
reliability, 22(3):583–602, 1982.
[13] J.H. Curtiss. A note on the theory of moment generating functions. The
Annals of Mathematical Statistics, 13(4):430–433, 1942.
[14] A. David and S. Larry. The least variable phase type distribution is
Erlang. Stochastic Models, 3(3):467–473, 1987.
[15] A.P. Dempster, N.M. Laird, and D.B. Rubin. Maximum likelihood from
incomplete data via the EM algorithm. Journal of the Royal Statistical
Society. Series B (Methodological), pages 1–38, 1977.
[16] S.E. Elmaghraby, R. Benmansour, A. Artiba, and H. Allaoui. On The
Approximation of Arbitrary Distributions by Phase-Type Distributions.
[17] A.K. Erlang. Solution of some problems in the theory of probabilities
of significance in automatic telephone exchanges. The Post Office
Electrical Engineers Journal, 10:189–197, 1917.
[18] M.J. Faddy. Compartmental models with phase-type residence-time
distributions. Applied stochastic models and data analysis, 6(2):121–
127, 1990.
[19] M.J. Faddy. A structured compartmental model for drug kinetics. Biometrics,
pages 243–248, 1993.
[20] A. Feldmann and W. Whitt. Fitting mixtures of exponentials to long-tail
distributions to analyze network performance models. Performance
evaluation, 31(3-4):245–279, 1998.
[21] A. Horvath and M. Telek. Approximating heavy tailed behaviour with
phase type distributions. 2000.
[22] A. Horv´ath and M. Telek. Phfit: A general phase-type fitting tool.
Computer Performance Evaluation: Modelling Techniques and Tools,
pages 1–14, 2002.
[23] A. Horv´ath and M. Telek. Matching more than three moments with
acyclic phase type distributions. Stochastic models, 23(2):167–194, 2007.
44
8. CONCLUSION
[24] C. Isensee and G. Horton. Approximation of discrete phase-type distributions.
In Proceedings of the 38th annual Symposium on Simulation,
pages 99–106. IEEE Computer Society, 2005.
[25] M.A. Johnson and M.R. Taaffe. Matching moments to phase distributions:
Mixtures of Erlang distributions of common order. Stochastic
Models, 5(4):711–743, 1989.
[26] K. Kim and N. Thomas. A fitting method with generalized Erlang
distributions. Simulation Modelling Practice and Theory, 2011.
[27] S. Kullback. Information theory and statistics. Dover Publications,
1997.
[28] G. Latouche and V. Ramaswami. Introduction to matrix analytic methods
in stochastic modeling. ASA-SIAM, Philadelphia, 1999.
[29] M. Malhotra and A. Reibman. Selecting and implementing phase
approximations for semi-Markov models. Stochastic Models, 9(4):473–
506, 1993.
[30] R. Marie. Calculating equilibrium probabilities for λ(n)/Ck/1/N
queues. In ACM Sigmetrics Performance Evaluation Review, volume 9,
pages 117–125. ACM, 1980.
[31] G.J. McLachlan and T. Krishnan. The EM algorithm and extensions,
volume 274. Wiley New York, 1997.
[32] R. Nelson. Probability, stochastic processes, and queueing theory: the
mathematics of computer performance modelling. Springer, 1995.
[33] M.F. Neuts. Probability distributions of phase type. Liber Amicorum
Prof. Emeritus H. Florin, 173:206, 1975.
[34] D.J. Newman and A.R. Reddy. Rational approximation to ex and to
related functions. Journal of Approximation Theory, 25(1):21–30, 1979.
[35] H. Okamura, T. Dohi, and K.S. Trivedi. A refined EM algorithm for PH
distributions. Performance Evaluation, 68(10):938, 2011.
[36] M. Olsson. The empht-programme. Relat´orio t´ecnico, Department
of Mathematics, Chalmers University of Technology, and G¨oteborg
University, Sweden, 1998.
45
8. CONCLUSION
[37] T. Osogami and M. Harchol-Balter. A closed-form solution for mapping
general distributions to minimal PH distributions. Computer Performance
Evaluation. Modelling Techniques and Tools, pages 200–217,
2003.
[38] T. Osogami and M. Harchol-Balter. Necessary and sufficient conditions
for representing general distributions by Coxians. Computer Performance
Evaluation. Modelling Techniques and Tools, pages 182–199,
2003.
[39] T. Osogami and M. Harchol-Balter. Closed form solutions for mapping
general distributions to quasi-minimal PH distributions. Performance
Evaluation, 63(6):524–552, 2006.
[40] A. Reibman, R. Smith, and K. Trivedi. Markov and Markov reward
model transient analysis: An overview of numerical approaches. European
Journal of Operational Research, 40(2):257–267, 1989.
[41] L. Schmickler. Meda: mixed Erlang distributions as phase-type representations
of empirical distribution functions. Stochastic Models,
8(1):131–156, 1992.
[42] J. Shore. On a relation between maximum likelihood classification
and minimum relative-entropy classification (Corresp.). Information
Theory, IEEE Transactions on, 30(6):851–854, 1984.
[43] Y. Takahashi. Asymptotic exponentiality of the tail of the waiting-time
distribution in a PH/PH/c queue. Advances in Applied Probability,
pages 619–630, 1981.
[44] M. Telek. The minimal coefficient of variation of discrete phase type
distributions. 2000.
[45] M. Telek and A. Heindl. Matching moments for acyclic discrete and
continuous phase-type distributions of second order. 2002.
[46] A. Thummler, P. Buchholz, and M. Telek. A novel approach for phasetype
fitting with the EM algorithm. Dependable and Secure Computing,
IEEE Transactions on, 3(3):245–258, 2006.
[47] H.C. Tijms. Stochastic models: an algorithmic approach. John Wiley &
Sons New York, 1994.
46
8. CONCLUSION
[48] A. Van de Liefvoort. The moment problem for continuous distributions.
Unpublished technical report, University of Missouri, WP-CM-1990-02,
Kansas City, 1990.
[49] CF Wu. On the convergence properties of the EM algorithm. The
Annals of Statistics, 11(1):95–103, 1983.
47