FITTING PHASE–TYPE SCALE MIXTURES TO HEAVY–TAILED DATA AND
DISTRIBUTIONS

arXiv:1705.04357v1 [math.ST] 11 May 2017

MOGENS BLADT AND LEONARDO ROJAS-NANDAYAPA
A BSTRACT. We consider the fitting of heavy tailed data and distribution with a special attention to
distributions with a non–standard shape in the “body” of the distribution. To this end we consider a
dense class of heavy tailed distributions introduced in [3], employing an EM algorithm for the the
maximum likelihood estimates of its parameters. We present methods for fitting to observed data,
histograms, censored data, as well as to theoretical distributions. Numerical examples are provided
with simulated data and a benchmark reinsurance dataset. We empirically demonstrate that our model
can provide excellent fits to heavy–tailed data/distributions with minimal assumptions.
Keywords. Statistical inference, heavy–tailed, phase–type, scale mixtures, approximating distributions, EM algorithm.

1. I NTRODUCTION
In this paper we consider the maximum likelihood estimation for a dense class of nonnegative
heavy–tailed tailed distributions, referred to as scale mixtures of phase–type distribuitions (NPH),
which was defined in [3]. Distributions in the NPH class allow for the simultaneous modelling of
the “body” and the “tail” of general distributions which are assumed to be absolutely continuous
and nonnegative while their “tails” are assumed to belong to some general class of heavy tailed
distributions, like for instance Regularly Varying or Weibullian. These very general assumptions
allows us to fit heavy–tailed distributions which may look distinctively different from distributions
usually found in catalogues.
Apart from providing an adequate description of data, distributions from NPH can also be seen
as infinite–dimensional phase–type distributions, which to all intents and purposes are as tractable
as their finite–dimensional counterparts. Much of the machinery available for finite–dimensional
phase–type distributions is also applicable to the extended class. Algorithms are for example available for the exact calculations of properties related to renewal theory, random walks (ladder processes) and ruin probabilities (see [3] for details).
The maximum likelihood estimation will be carried out employing an EM algorithm similarly as
for finite–dimensional phase–type distributions [2]. The main challenge we face is the algorithmic
implementation resulting from the extension to infinite dimensions since we cannot make a pre–
fixed cut–off in the number of dimensions as this would be equivalent to a light–tailed estimation.
We shall see, however, that it is possible to reduce the formulas to simple (one–dimensional) infinite
series which can be numerically evaluated to a specified degree of precision. We are also able to
simplify certain expressions of the original paper ([2]), obtaining explicit expressions which allow
for an improved numerical performance.
The rest of the paper is organized as follows. In Section 3 we develop the main algorithm for
estimating independent and identically distributed data sampled from an NPH distribution. We
show that the algorithm essentially uses the empirical cumulative distribution function, and that
a significant increase in speed may be obtained by grouping the data into intervals on parts of the
support and considering the corresponding histograms. In Section 4 we adjust the algorithm to cope
with the presence of left-, right- or interval censored data. Section 5 provides an EM–algorithm for
adjusting an NPH distribution to a theoretical distribution function F, which in turn is equivalent to
finding the distribution in the NPH class with a given order which minimizes the Kullback–Leibler
1

2

M. BLADT AND L. ROJAS

distance to F. In section 6 we provide some numerical examples from both simulated data, real
data and a fits to theoretical distribution. In there, we highlight some minimal assumptions made
for efficiently fitting general Regularly Varying and Weibullian distributions.
2. P HASE – TYPE DISTRIBUTIONS AND THE NPH CLASS
Before proceeding with a more detailed account, we first provide some background on phase–
type distributions and the extended class NPH of scale mixtures of phase–type distributions.
The class PH of phase–type distributions consists of distributions of (random) times until a finite
state Markov jump process exits a set of transient states. This can be made precise by letting
E = {1, 2, ..., p, p + 1} denote the state–space of a Markov jump process {Xt }t≥0 , where states
1, 2, ..., p are transient and p + 1 is absorbing. The intensity matrix for {Xt }t≥0 can then be written
on the form


T t
,
Λ=
0 0
where T is a p× p sub–intensity matrix, t is a p–dimensional column vector and 0 the p–dimensional
row vector of zeroes. We follow the convention that matrices are written in capital bold and their
elements with the corresponding minuscule (like A = {ai j }), bold minuscule Greek letters (like
α ) are row vectors while bold minuscule Roman letters (like t ) are column vectors. Their dimensions are usually clear from the context and left unspecified unless necessary. Let αi = P(X0 = i),
p
∑i=1 αi = 1 and α = (α1 , ..., α p ). Then we say that
τ = inf{s > 0 : Xs = p + 1}
α , T ). Since t = −T
T e , where e denotes the
has a phase–type distribution with representation PH p (α
column vector of ones, the distribution of τ is fully specified in terms of α and T . For further
background on phase–type distributions we refer e.g. to [?], [4], [7] or [8].
The class of phase–type distributions is widely used in the area of Applied Probability, where
they may often provide exact (or even explicit) solutions in complex stochastic models. This is for
example the case for ruin probabilities in risk theory or waiting time distributions for queues. Any
distribution with support on the positive real numbers may be approximated arbitrarily close by a
phase–type distribution. In spite of this denseness property, phase–type distributions are all light
tailed and consequently any approximating (finite–dimensional) phase–type distribution will not be
able to capture a possible heavy tailed behaviour.
In [3] the dense class, NPH, of genuinely heavy tailed distributions is proposed in terms of
infinite–dimensional phase–type distributions with a finite number of parameters. The idea is very
simple and goes a follows. Let N be a discrete random variable with support on {si : i ∈ N, si > 0}
and distribution π = {πi }i≥1 where πi = P(N = si ). If N, for example, is a discretization of a
continuous distribution G at step length ∆, then si = i∆ and πi = G(i∆) − G((i − 1)∆), i = 1, 2, . . . .
α , T ) be independent of N. When sampling from Y , we first draw an index i (referred
Let τ ∼ PH p (α
α , T /si ). Hence, the
to as the level) from π and then a phase–type random variable from PH p (α
random variable Y = Nτ may be seen as a scale mixture of phase–type distributions, so we often
refer to N as the scaling random variable and its distribution π as the scaling distribution of the
NPH. The density of Y can written as
∞

fY (y) = ∑ πi α eT y/si t /si ,

y > 0.

i=1

The distribution of Y can also be seen as an infinite dimensional phase–type distribution since we
can write
π ⊗ α ) eΓ y γ ,
fY (y) = (π

FITTING PHASE–TYPE SCALE MIXTURES TO HEAVY–TAILED DATA AND DISTRIBUTIONS

3

where



Γ=



T1 0 0 0
0 T2 0 0
0 0 T3 0
0 0 0 T4
..
..
..
..
.
.
.
.

···
···
···
···
..
.





 
 
=
 
 

T /s1
0
0
0
0
T /s2
0
0
0
0
T /s3
0
0
0
0
T /s4
..
..
..
..
.
.
.
.

···
···
···
···
..
.








Γe . Here e is now the infinite–dimensional column vector of ones and we notice that the
and γ = −Γ
exponential of the infinite dimensional matrix Γ is well defined since it is a bounded operator (as
T i e . We shall write
the sequence {si : si > 0} is bounded away from zero). We also let t i = −T
π , α , T ).
Y ∼ NPH p (π
While the support of π is important, we will not denote it explicitly in the parametrisation of Y .
A key feature of the NPH class is that it contains a rich variety of heavy–tailed distributions. For
instance, if π has unbounded support, then Y has a heavy–tailed distribution [?, cf.]]RojasXie. In
the regularly varying case, the tail of Y greatly follows that of N. Breiman’s lemma implies that if
N has a regularly varying with tail index −α < 0, then the tail of Y will also be regularly varying
with the same index. More precisely,
P(Y > t) = E[τ α ]P(N > t)(1 + o(1)),

t → ∞.

A similar feature occurs for the the class of Weibullian distributions [1], which is defined as the
collection of nonnegative distributions having survival function
(1)

p

F(x) = xδ e−(λ x) (C + o(1)),

x > 0, λ , p > 0, δ ∈ R.

If the parameter p ∈ (0, 1), then the distribution is heavy–tailed (and light–tailed otherwise). Since
a PH distribution is Weibullian with parameter p = 1, then Lemma 2.1 of [1] implies that if the
π , α , T ) is Weibullian with
scaling distribution π is Weibullian with parameter p1 then Y ∼ NPH(π
parameter p1 /(1 + p1 ).
3. E STIMATION
In this section we address the problem of fitting an NPH distribution to a data set. We assume
that y1 , y2 , ..., yM forms an i.i.d. data set sampled from NPH p (π(θθ ), α , T ), where π (θθ ) is some
parametric distribution describing the distribution of N. We assume that the support for π (θθ ) does
not depend on θ . We shall estimate the parameters θ , α and T .
Hence, with probability π (θθ )i , yn is the realization of the i–th level phase–type distribution
α , T /si ), but both the level i and the actual realization of the underlying Markov jump process
PH p (α
are not observable. Thus, we resort to the EM algorithm for approximating the maximum likelihood
estimators. To this end we first attend the corresponding estimation problem for complete data.
Assume that apart from y1 , ..., yM we have also observed all levels and all underlying Markov
jump process. Let In denote the level of the phase–type distribution of the n’th path, and let
M

Li =

∑ 1{In = i}

n=1

denote the number of i–level processes in the data. Next consider the Markov jump process
(n)
{Ju }u≥0 underlying the n’th phase–type distribution (which generates the data yn ) and let
M

Bik

=

(n)

∑ 1{J0

n=1

= k, In = i}

4

M. BLADT AND L. ROJAS

be the number of i–level processes that are initiated in state k. Define
M Z yn

Zki

=

∑

n=1 0

(n)

1{Ju = k, In = i}du,

i
which is the total time all underlying i–level Markov jump processes spend in state k and let Nk`
denote the total number of jumps from state k to ` within all i–level Markov jump processes. Finally,
let Nki be the number of i–level processes that exit to the absorbing state from state k.
Then the complete data likelihood is easily seen to be (see e.g. [2] or [?] for further comments
on this)


 i

  i

p
p 
∞
tk` Nk`
tk` i p tk Nk
tk i
Bik
Li
Lc (θθ , α , T ) = ∏ πi (θθ ) ∏ αk ∏
exp − Zk ∏
exp − Zk
si
si
si
k,`=1
i=1
k=1
k=1 si
`6=k

with corresponding log–likelihood
∞

`c (θθ , α , T ) =

∞

i

∑ L log πi(θθ ) + ∑ ∑

i=1

Bik log αk +

i=1 k=1


tk` i ∞ p i
tk
− ∑ ∑ Zk + ∑ ∑ Nk log
si
i=1 k=1
i=1 k,`=1 si
∞

(2)

p

p

∞

p



∑∑

i
Nk`
log

−∑

∑ ski Zki .

i=1 k,`=1
`6=k
 ∞ p

tk`
si



t

i=1 k=1

`6=k

α , T̂
T ), which is the point that maximizes
In order to calculate the maximum likelihood estimator (θ̂θ , α̂
the likelihood (or log–likelihood) function, we calculate first order partial derivatives of tk` , tk , while
we use the Lagrange multiplier methods for αi and πi (θθ ) due to the constraints on them summing
up to one.
Consider the parameter tk` , k 6= `. Then
∞
∞ Zi
∂ `c
i 1
= ∑ Nk`
−∑ k =0
∂tk` i=1
tk` i=1 si

implying
tˆk` =

i
∑∞
i=1 Nk`
.
i
∑∞
i=1 Zk /si

Similarly,
i
∑∞
i=1 Nk
.
i
∑∞
i=1 Zk /si
The diagonal terms tˆkk = − ∑`6=k tˆk` − tˆk . Regarding α , consider the Lagrange function

tˆk =

∞

α ) = ∑ Bik log αk + µ(1 − ∑ αk ),
M(α
i=1

k

where µ is a Lagrange multiplier. Then
∞ Bi
∂M
= ∑ k − µ = 0,
∂ αk i=1 αk

which result in
∞

αk µ = ∑ Bik .
i=1

FITTING PHASE–TYPE SCALE MIXTURES TO HEAVY–TAILED DATA AND DISTRIBUTIONS

5

Summing over k yields
∞

µ=∑

p

∑ Bik = M

i=1 k=1

so
α̂k =

1 ∞ i
∑ Bk .
M i=1

Concerning π (θθ ), it depends on the particular form of the discrete distribution whether it can be
estimated explicitly or numerically. We shall consider some particular examples.
We now consider the case of incomplete data. We shall employ the EM–algorithm, which optimizes the incomplete data likelihood L,
M

L(θθ , α , T ; y ) = ∏ fY (yk ; θ , α , T ),
k=1

using the complete likelihood Lc (or `c ) in the following way. Here fY (y; θ , α , T ) denotes the
density function of Y ∼ NPH p (π(θθ ), α , T ).
0: Initialize with some “arbitrary” (θθ 0 , α 0 , T 0 ) and let n = 0.
1: (E–step) Calculate the function
Y = y) .
h : (θθ , α , T ) → E(θθ n ,αα n ,TT n ) (`c (θθ , α , T )|Y
2: (M–step) Let (θθ n+1 , α n+1 , T n+1 ) := argmax(θθ ,αα ,TT ) h(θθ , α , T ).
3: n=n+1; GOTO 1.
In each iteration, the incomplete likelihood is increased (i.e. L(θθ n , α n , T n ) ≤ L(θθ n+1 , α n+1 , T n+1 ))
and the procedure hence converges (possibly to a local maximum or saddlepoint though).
We notice that the actual calculations in the EM–algorithm only involve the complete data likelihood. From the actual form of the log–likelihood (2) we see that it is a linear function of the
sufficient statistics, and the conditional expected value of the log–likelihood given the data will
then be the log–likelihood function with the sufficient statistics replaced by their conditional expectations given the data. Hence the M–step is trivial since we only have to plug in the conditional
expectations given the data instead of the sufficient statistics which are not available
Concerning the E–step we proceed as follows. First we consider one (generic) data point (M = 1)
and let y = y1 . We need to calculate the conditional expected values of the sufficient statistics given
Y = y. All distributions and expectations are under P = PΨ , Ψ = (θθ , α , T ), but we will omit the
index in order to ease the exposition.
Concerning Li , it equals one if i is the chosen level and zero otherwise. Let I denote the random
variable indicating the chosen level. Then
E(Li |Y = y) = P(I = i|Y = y)
P(Y ∈ dy|I = i)P(I = i)
=
P(Y ∈ dy)
α exp(T
T i y)tt i
πi (θθ )α
=
.
fY (y; θ , α , T )

(3)
Similarly, for Bik we have

E(Bik |Y = y) = E(1{J0 = k, I = i}|Y = y)
P(Y ∈ dy|J0 = k, I = i)P(I = i)P(J0 = k)
=
P(Y ∈ dy)
0
T i y)tt i
πi (θθ )αk e k exp(T
=
fY (y; θ , α , T )

6

M. BLADT AND L. ROJAS

so

!

∞

∑ Bik Y = y

E

=

i=1

Regarding

θ )αk e 0k exp(T
T i y)tt i
∑∞
i=1 πi (θ
.
fY (y; θ , α , T )

Zki ,
E

Zki |Y

Z y


=y = E
1{Ju = k, I = i}du Y = y
0


Z y
= E 1{I = i}
1{Ju = k}du Y = y
0
Z y

= E
1{Ju = k}du Y = y, I = i P(I = i|Y = y)
0

Z y

=
=
=
=
=

P(Y ∈ dy|I = i)P(I = i)
P(Y ∈ dy)
0
Z y
P(Y ∈ dy, Ju = k|I = i) P(Y ∈ dy|I = i)P(I = i)
du
P(Y ∈ dy|I = i)
P(Y ∈ dy)
0
Z y
P(I = i)
P(Y ∈ dy, Ju = k|I = i)du
P(Y ∈ dy) 0
Z y
P(I = i)
P(Y ∈ dy|Ju = k, I = i)P(Ju = k|I = i)du
P(Y ∈ dy) 0
Z y
πi (θθ )
e 0 eT i (y−u)t i α eT i u e k du.
fY (y; θ , α , T ) 0 k
P(Ju = k|Y = y, I = i)du

i |Y = y) is derived in a similar way, resulting in
The formula for E(Nk`
i
E(Nk`
|Y

πi (θθ )
tk`
= y) =
si fY (y; θ , α , T )

Z y
0

e 0` eT i (y−u)t i α eT i u e k du.

Finally,
tk
πi (θθ )
α eT i y e k .
si fY (y; θ , α , T )
The θ needs to be treated separately on a case–by–case basis and it generally involves a numerical
procedure for finding the next iteration θ (n+1) . We need to maximize
E(Nki |Y = y) =

∞

θ → ∑ E(Li |Y = y) log πi (θθ )
i=1

subject to

θ ) = 1,
∑∞
i=1 πi (θ

where (see (3))

πi (θθ )
T i y)tt i .
α exp (T
fY (y; θ , α , T )
In general, for M > 1 datapoints we simply sum the previous formulas with arguments y j , j =
1, ..., M.
E(Li |Y = y) =

Remark 3.1. We see that the formulas in the E–step involve both matrix–exponentials and integrals
thereof, and by defining the generic integral
Z y

J (y) := J (y; α , T ) =

eT (y−u)t α eT u du,

0

we have that (see [12])

(4)

exp

T tα
0 T

   Ty

e
J (y)
y =
.
0 eT y

FITTING PHASE–TYPE SCALE MIXTURES TO HEAVY–TAILED DATA AND DISTRIBUTIONS

7

T y) and J (y) is by calculating
Thus a simple (and numerically efficient) way of obtaining both exp(T
the matrix exponential on the left hand side.
The EM–algorithm can now be stated as follows.
Theorem 3.2 (EM–algorithm).
0: Initialize with some “arbitrary” (θθ , α , T ).
1: (E–step) Compute
M

Y = y ) = πi (θθ ) ∑
E(Li |Y

j=1

T i y j )tt i
α exp(T
fY (y j ; θ , α , T )

T i y j )tt i
αk e 0k exp(T
j=1 fY (y j ; θ , α , T )
M

Y = y ) = πi (θθ ) ∑
E(Bik |Y

M J (y /s ; α , T )

j i
kk
Y = y = πi (θθ ) ∑
E Zki |Y
f
(y
;
θ
,
α
,
T
)
j=1 Y j
M

i
Y = y) =
E(Nk`
|Y

πi (θθ )
si

J (y j /si ; α , T )`ktk`
j=1 fY (y j ; θ , α , T )

πi (θθ )
si

M

Y = y) =
E(Nki |Y

∑
∑

j=1

α eT i y j e ktk
fY (y j ; θ , α , T )

2: (M–step) Maximize
∞

Y = y ) log πi (θθ )
θ → ∑ E(Li |Y
i=1

subject to ∑i πi (θθ ) = 1 and let θ̂θ be the argument at which the maximum is attained. Let
α̂k =

1 ∞
∑ E(Bik |YY = y )
M i=1

tˆk` =

i Y = y)
∑∞
i=1 E(Nk` |Y
1
i Y = y)
∑∞
i=1 si E(Zk |Y

tˆk =

i Y = y)
∑∞
i=1 E(Nk |Y
1
i Y = y)
∑∞
i=1 si E(Zk |Y

T = {tˆk` }k,`=1,...,p .
and assign diagonal values tˆkk = − ∑l6=k tˆk` − tˆk so T̂
3: Reassign values to initial parameters
θ := θ̂θ
α
α := α̂
T
T := T̂
4: GOTO 1.
Remark 3.3. Consider the E–step of Theorem 3.2. Suppose that there are repeated values among
the data points such that among the data y1 , ..., yM there are ỹ1 , ..., ỹD different values and that ỹi
appears ki ≥ 1 times in the original data. Then k1 + · · · kD = M and the sums over j = 1, ..., M in
expected values can then be reduced to weighted sums of fewer terms instead. For example,
D

J (ỹ j /si ; α , T )kk
Y = y = πi (θθ ) ∑ k j
E Zki |Y
.
fY (ỹ j ; θ , α , T )
j=1

8

M. BLADT AND L. ROJAS

This shows that the EM–algorithm of Theorem 3.2 can also be used to estimate an NPH distribution
when data are represented by a histogram rather than the raw data. This is particularly suitable when
the amount of data is large since the histogram may here adequately represent the target distribution.
Remark 3.4. In order to reduce the computational burden, we may replace the raw data by a
histogram in certain regions of the support where data gather rather densely. Heavy–tailed data
with support on [0, ∞) will typically concentrate in an interval [0, T ) (the ”body” of the distribution)
while it will be more scarce in [T, ∞) (the ”tail”).
Now assume that M is large and that the support for the data may be split into [0, T ) and [T, ∞)
for some T > 0 such that the concentration of data points in [0, T ] is large. Then divide [0, T ] into K
disjoint sub–intervals [Ti , Ti+1 ), i = 0, 1, . . . , K − 1 where 0 = T0 < T1 < · · · < TK = T and ki denote
the number of data points falling into [Ti , Ti+1 ). Then we may use the repeated data reduction of
i+1
Remark 3.3 for the interval [0, T ) by treating Ti +T
, i = 0, ..., K − 1 as data points with counts ki .
2
In fact we might choose any point in [Ti , Ti+1 ) as a representative, including any data point falling
in this interval, however, if we choose the left points of the intervals, Ti , as representatives, we have
to make special arrangements regarding the first interval in order not to provoke an atom at zero.
Concerning the tail, [T, ∞), we continue to use the raw data since the scattering in the tail is
typically too diffuse in order to be represented by a histogram without too much loss of information.
Example 3.5. Here we present an important special case for the choice of the scaling distribution {πi (θθ )} of N for the Regularly Varying case. It features an explicit solution to the M–step
maximization problem of the function
∞

θ → ∑ log πi (θ )wi ,
i=1

Y = y ). Define
where wi = E(Li |Y
πi (θ ) = P(N = si ) = F(si+1 ) − F(si ) =

1
1
− θ ,
θ
si
si+1

si = e(i−1)c .

= (e−θ c )(i−1) (1 − e−θ c ),

i = 1, 2, . . .

This distribution corresponds to a the discretization of a Pareto distribution F and the resulting discrete distribution is supported by {eic : i =∈ N}. For c > 0, N with parameter e−θ c . The advantage
of representing the scaling Pareto distribution on this form is that its argument which maximizes
∞

θ → ∑ log πi (θ )wi ,
i=1

Y = y ), and it is not difficult to see it has an explicit solution given by
where wi = E(Li |Y


∞
∑ wi 

1


i=1
θb = − log 1 − ∞
.
c


∑ i wi
i=1

The selected support of N is also convenient for implementation purposes. In practice, the infinite
series defining the estimators can only be computed up to a finite number of terms. The sequence
π(θ ) converges faster to 0 so a significant reduction of terms to be computed is required in order
to attain a certain specified precision.
Since the distance between consecutive points in the support of N is unbounded, then the distribution of N is not regularly varying and Breiman’s lemma no longer applies. Nevertheless, the tail
probability of Y oscillates between two regularly varying functions, so this model can provide an
accurate approximation to any regularly varying distribution [11].

FITTING PHASE–TYPE SCALE MIXTURES TO HEAVY–TAILED DATA AND DISTRIBUTIONS

9

4. C ENSORING
In certain situations some data may be censored. Again we consider first the case of a single data
point y taken from a realization of Y ∼ NPH p (π(θθ ), α , T ). The data point is right censored at t if
the only knowledge about Y is that Y > t, left censored at t if Y ≤ t and interval censored at (s,t]
if Y ∈ (s,t]. Left censoring is a special case of interval censoring with s = 0 while right censoring
can be obtained by fixing s and letting t → ∞. Formulas for right censoring will, however, appear
as a part of the derivation of interval censoring.
The EM algorithm works entirely the same way as for uncensored data with the only difference
that we are no longer observing a data point Y = y but Y ∈ (s,t]. This will only change the E–step
where we now have to calculate the following conditional expectations:
i
E(Li |Y ∈ (s,t]), E(Bik |Y ∈ (s,t]), E(Zki |Y ∈ (s,t]), E(Nk`
|Y ∈ (s,t]) and E(Nki |Y ∈ (s,t]).

Concerning Li , notice that
E(Li 1{Y ≤ t}) = P(Y ≤ t|I = i)P(I = i)
= πi (θθ )(1 − α eT it e )
and
α eT it e .
E(Li 1{Y > t}) = πi (θθ )α
Thus
E(Li 1{Y ∈ (s,t]})
P(Y ∈ (s,t])
E(Li ) − E(Li 1{Y ≤ s}) − E(Li 1{Y > t})
=
P(Y ∈ (s,t])
α eT it e
πi (θθ ) − πi (θθ )(1 − α eT i s e ) − πi (θθ )α
=
P(Y ∈ (s,t])
T
s
α e i e − α eT it e )
πi (θθ )(α
=
.
α )(α
α eT i s e − α eT i t e )
∑∞
i=1 πi (α

E(Li |Y ∈ (s,t]) =

For right censored data we get
α eT it e
α eT it e
πi (θθ )α
πi (θθ )α
= ∞
.
α eT it e
P(Y > t)
∑i=1 πi (θθ )α
The rest of the formulas are derived as for censored phase–type distributions (see [9]) conditionally on the level Li , with parameters α , T i , which happens with probability πi (θθ ). Thus we get
that

πi (θθ ) αk e 0k eT i s e − αk e 0k eT i e e
i
E(Bk |Y ∈ (s,t]) =
P(Y ∈ (s,t])
Z t

Z t
T iu
0 T i (t−u)
T iu
πi (θθ )tk` /si
α e e k du − e ` e
e α e e k du
s
s
i
E(Nk` |Y ∈ (s,t]) =
P(Y ∈ (s,t])
Z t

Z t
T iu
T iu
0 T i (t−u)
πi (θθ )
α e e k du − e k e
e α e e k du
s
s
i
E(Zk |Y ∈ (s,t]) =
P(Y ∈ (s,t])
E(Li |Y > t) =

Z t

E(Nki |Y ∈ (s,t]) =

πi (θθ )tk /si

s

α eT i u e k du

P(Y ∈ (s,t])

10

M. BLADT AND L. ROJAS

with similar and obvious formulas for the right censored case. The integrals are calculated similarly
as in the uncensored case, namely
Z t

T is
α eT i u e k du = α T −1
− eT it )eek
i (e

s

= si α T −1 (eT i s − eT it )eek
Z t
s

e 0` eT i (t−u) e α eT i u e k du = J i (t)`k − J i (s)`k ,

where


 
T i eα
J (y) = exp
y .
0 Ti
For more than one data point, the data are split into a group of uncensored data and into other
groups of different types of censored data. The conditional expectations are then calculated for all
data points subject to their group classification, and all conditional expectations of the same kind
(jumps, occupation times etc.) are then summed over all data. This amounts to the E–step in an
EM algorithm, the rest of which is identical to Theorem 3.2.
i

Remark 4.1. In Remark 3.4 we suggested the use of a histogram for representing the “body” of
the distribution if the amount of data is large. Interval censoring provides a feasible alternative in
this same direction.
5. F ITTING TO A KNOWN DISTRIBUTION
It may be of interest to approximate a given (heavy tailed) distribution H by a distribution
G ∈ NPH if for example methods for calculation the ruin probability based on the claim size distribution H is not known. In [2] it was shown how the EM algorithm can be modified in order
to approximate phase–type distributions to a given distribution F with non–negative support. The
EM algorithm then converges to a limit which minimizes the Kullback–Leibler distance between
phase–type distributions of a specified order and the distribution F.
The idea is to let the number of data points M → ∞ such that the distribution H is seen as an
empirical distribution function for a data set with M = +∞. If we let M → ∞ then (see Theorem
3.2)
Z ∞
∞
∞
T i y j )tt i
T i y)tt i
αk e0k exp(T
1 M αk e0k exp(T
α̂k = ∑ πi (θθ ) ∑
→ ∑ πi (θθ )
dH(y).
M j=1 fY (y j ; θ , α , T )
fY (y; θ , α , T )
0
i=1
i=1
Similarly it is not difficult to see that
tˆk` =

1
i Y = y)
i Y = y)
∑∞
∑∞
i=1 M E(Nk` |Y
i=1 E(Nk` |Y
=
1
1 1
i Y = y)
i Y = y)
∑∞
∑∞
i=1 si E(Zk |Y
i=1 si M E(Zk |Y
∞

→

πi (θθ )
∑ si
i=1
∞

Z ∞
J (y/si ; α , T )`ktk`

fY (y; θ , α , T )

0

πi (θθ )
∑ si
i=1

Z ∞
J (y/si ; α , T )kk

fY (y; θ , α , T )

0

dH(y)

dH(y)

and
∞

tˆk →

πi (θθ )
∑ si
i=1
∞

∑

i=1

Z ∞
0

α eT i y e ktk
dH(y)
fY (y; θ , α , T )

Z
πi (θθ ) ∞ J (y/si ; α , T )kk

si

0

fY (y; θ , α , T )

dH(y)

.

FITTING PHASE–TYPE SCALE MIXTURES TO HEAVY–TAILED DATA AND DISTRIBUTIONS

11

Concerning θ , maximizing
∞

Y = y ) log πi (θθ )
θ → ∑ E(Li |Y
i=1

is equivalent to maximizing
∞

1
Y = y) log πi (θθ ).
E(Li |Y
M
i=1

θ→∑
As M → ∞ the latter converges to
∞ Z
θ →∑
i=1

0

∞


α exp(T
T i y)tt i
πi (θθ )α
dH(y) log πi (θθ ),
fY (y; θ , α , T )

which is then the function to be maximized. Let θ̂θ denote the argument which maximizes this
function.
In general none of the integrals will have explicit solutions, and approximations (e.g. numerical
integration, Quasi Monte Carlo methods or more sophisticated variants) will have to be employed.
The EM algorithm now works as follows:
0: initiate with some parameters (θθ , α , T ).
α , T̂
T and θ̂θ .
1: calculate α̂
α , T̂
T ).
2: assign (θθ , α , T ) = (θ̂θ , α̂
3: GOTO 1.
6. E XAMPLES
In this section we provide five examples. In the first example we consider the simplest case of
simulated data from a scale mixture of exponential distributions, where the scaling distribution has
a Pareto type of tail. In the second example we fit NPH distributions to a real data set (Danish
reinsurance data of fire claims), while in the last three examples we consider the fitting of NPH
distributions to the theoretical distributions log–gamma, Weibull and log–normal respectively.
Example 6.1 (Erlang distributions). A q dimensional Erlang distribution, ERq (λ ), is a phase–type
distribution with canonical representation




−λ λ
0 ··· 0
0
 0 −λ λ · · · 0 
 0 




0
0 −λ · · · 0  and t =  0  .
α = (1, 0, 0, ..., 0), T = 
 .
 . 
..
..
......
.. 
 .. 
 ..
.
.
...
. 
0
0
0 · · · −λ
λ
Hence, the correspoding NPH distribution has density
 q q−1
λ
y
e−yλ /si ,
fY (y; θ , λ ) = ∑ πi (θθ )h(y/si ; q, λ ) = ∑ πi (θθ )
s
(q
−
1)!
i
i=1
i=1
∞

(5)

∞

and the maximum likelihood estimator for λ is easily calculated to be
1 ∞ Zi
b
λ −1 =
∑ si ,
qM i=1

N

Zi =

∑ y j · 1(I j = i),

j=1

while the maximum likelihood estimate for θ has the general form
∞

θb = arg max ∑ Li log πi (θ ),
θ

i=1

N

Li =

∑ 1(I j = i).

j=1

12

M. BLADT AND L. ROJAS

The E–step is similar as for the general model. In the case of the sufficient statistic Z i we have
Y = y) =
E(Z i |Y

M

M

∑ y j P(I j = i|Y = y j ) =

∑ yj

j=1

j=1
M

=

∑ yj

j=1

P(I j = i)P(Y ∈ dy j |I = i)
P(Y ∈ dy j )
πi (θθ )h(y j ; q, λ /si )
.
f (y j ; θ , λ )

As for the sufficient statistic Li , its expected value is analogue for the general model in the
previous section and found to be equal to
M


Y =y =
E Li |Y

πi (θθ )h(y j /si ; q, λ )
.
f (y j ; θ , λ )
j=1

∑

We present a small simulation study for the case of q = 1 (corresponding to infinite dimensional
hyperexponential distribution), si = i and scaling distribution
πi (θ ) =

i−θ
, i = 1, 2, ...,
ζ (θ )

where

∞

ζ (θ ) = ∑ i−θ
i=1

is the Riemann Zeta function with parameter θ . This distribution is also known as the Riemann
Zeta or discrete Pareto distribution since its tail resemble that of a Pareto distribution.
In order to find the EM–estimate of θ (n+1) , we have to maximize the function
∞

i

∞

Y = y)θ log(i) − ∑ E(Li |Y
Y = y) log(ζ (θ )).
h(θ ) = − ∑ E(L |Y
i=1

i=1

Differentiating with respect to θ and rearranging implies that we have to solve the following equation w.r.t. θ ,

1 ∞
ζ 0 (θ )
Y = y log(i),
= − ∑ E Li |Y
(6)
ζ (θ )
M i=1
which is done numerically using a simple Newton–Raphson procedure.
Results of the simulation study, which is shown in Table 1, reveals that the EM algorithm is able
to recover the underlying structure of the data, and that the estimation, as expected, improves with
larger sample sizes.
parameters
λ
θ
100
500
1.0 2.0 1.04 1.88 0.95 2.09
1.0 2.5 0.85 2.67 0.90 2.52
1.0 5.0 0.92 7.1 1.09 7.26

sample size
1000
5000
10000
1.14 1.93 0.99 2.00 1.01 2.01
1.01 2.63 0.94 2.64 1.00 2.52
1.05 4.65 1.00 4.85 1.01 5.03

Table 1: EM–estimates of (b
λ , θb) infinite dimensional hyper–exponential distributions with varying
parameters and sample size
Example 6.2. We consider 2167 reinsurance data for Danish fire insurance claims above 1 million
DKR for the period 1980–1993. These data have been widely studied in Extreme Value Theory
[6, 10]. The data corresponds to claims in millions of Danish Kroner for the period 1980–1993, the
amounts being adjusted for inflation to prices of 1985. We subtracted 1 (million) from all data in

FITTING PHASE–TYPE SCALE MIXTURES TO HEAVY–TAILED DATA AND DISTRIBUTIONS

13

order to shift the support to [0, ∞) which is the natural support for a phase–type distributions. Of
the 2167 data, 519 are repeated values so only 1648 have different values.
We propose an NPH model for the data by mixing a five dimensional phase–type distribution
with an N that is assumed to follow the discretized Pareto distribution of Example 3.5.
Just above 90% of the data is below 5, while less than 10% falls between 5 and the maximum
observation of 262.2504. Thus we select the “body part” of the distribution as the densely populated
region [0, 5] while the much more diffuse region [5, ∞) is then considered the “tail”.
First we test the difference in performance between the EM–algorithm using raw data and one
using the hybrid method proposed in Remark 3.4.
We divide [0, 5] into 100 subintervals of the same size, thus treating losses within 50,000 DKR
as all being the same. We choose the centre of intervals as data points. In the tail we have 186 data
points of which 8 are repeated values so a total of 178 different values. This amounts to a total of
278 distinct data points. The result of the experiment is following. While the actual discretization
of [0, 5] does not appear to have any effect at all on the estimation, as compared with the EM–
algorithm of the raw data, the speed was increased by a factor 5.75 which is more or less consistent
with the execution times depending linearly on the number of data points, in which case we should
have expected an increase in speed by a factor 1648/278 = 5.92.
In all the EM algorithms which follows we have used the hybrid method of Remark 3.4 with
[0, 5] divided into 100 equally sized subintervals.
We now investigate the effect the choice of the parameter c (see Example 3.5) will have on the
estimation. Intuitively, the smaller we choose the discretization steps of the geometric progression,
the better the estimation. Fixing the dimension of the phase–type distribution to five, we select
three different values for c, 1, 1/2 and 1/4. To obtain the EM estimates we randomly selected
various sets of initial and iterated the EM until the relative change in the loglikelihood was smaller
than 10−8 . We kept the model with the largest likelihood. The estimates are as follows
c = 1:
θ̂ = 1.2743
α = (0.6415, 0.0099, 0.0055, 0.2115, 0.1316)
α̂


−2.7430 1.3565
0
0
0

 0.0003 −3.0398
0
0
0



T = 
0
0
−2.6313 1.1226
0.2167 
T̂


0
0
0.7223 −1.3953 0.4089 
0
0
1.5129
0.8388 −2.4779
c = 1/2:
θ̂ = 1.3136
α = (0.7692, 0.0149, 0.1154, 0.0148, 0.0857)
α̂


−3.0713 0.1234
0.2319
0.0655
1.6457
 0.2239 −1.9576 0.3787
0.9651
0.3899 



T =  0.1824
0.7410 −3.0705 0.7978
0.4776 
T̂

 0.3099
0.6150
0.4751 −1.5588 0.1588 
0.0034
0.0446
0.2263
0.0446 −3.4154
c = 1/4:
θ̂ = 1.3230
α = (0.0267, 0.4563, 0.0010, 0.2603, 0.2556)
α̂

14

M. BLADT AND L. ROJAS



−0.7896 0.2531
0.0385
0.2652
0.0265
 0.0810 −3.7022 1.7356
1.1626
0.3934 



T = 
0.5838
0.0013
−3.7890
0.0576
0.0302
T̂


 0.3386
0.0606
1.0320 −3.7735 0.1957 
0.6694
0.0292
0.2376
0.3920 −3.4129
α , T ) look distinct is because phase–type representations are by no means
That the estimates for (α
unique, however, the estimated densities which are plotted in Figure 1 against the histogram of the
data are almost identical (actually indistinguishable in the plotted range of [0, 6]). It is clear that the
different phases do not have a physical interpretation but are merely dummy (or black box) states
used for the only purpose of obtaining a proper adjustment of the distribution.
In order to inspect the tail behaviour we plotted the log–survival function in Figure 2. We have
provided two plots of the survival function in different ranges in order to be able to inspect the
tail behaviour well beyond the largest data value (which is 262.2504). The different values for θ
(obtained as a consequence of the different choices of c) does only result in a slightly different tail
behaviour for very large claim sizes. We conclude that the estimation does not depend significantly
on the choice of c within the possible choices of c = 1, 1/2 or 1/4, and can probably be extrapolated
to concluding robustness against choices in c within any “reasonable” range for c.

0.8

Data
c=1
c=1/2
c=1/4

0.6
0.4
0.2
0
0

1

2

3

4

5

6

F IGURE 1. Estimated densities against histogram of the NPH models for Danish
reinsurance data.
Finally we compare the fit for c = 1/4 to one obtained by [6] and [10]. The first reference fits
a Generalized Pareto Distribution via maximum likelihood estimation while the second employs
the Hill estimator. The implementation of the methods described above involve the selection of an
appropriate threshold which should be chosen by the modeler. According to [10] the values in the
interval [1.40, 1.46] produce excellent fits, the value recommend by the same author being 1.45.
We finally compare our estimated model with c = 1/4 to a NPH model where we fix θ = 1.45
as recommended by [10]. This can be done by running the EM algorithm in the usual way but
avoiding any adjustments in θ . The results of the estimation are given next
θ = 1.45: Fixed tail.
α = (0.1114, 0.0080, 0.3366, 0.3600, 0.1840)
α̂


−3.0292 0.8810
0.0655
0.1600
0.1538
 0.1204 −0.6765 0.3295
0.0829
0.0827 



T =  1.4308
0.2854 −4.2611 1.1133
1.4316 
T̂

 0.4463
0.2717
0.3082 −3.6472 1.3508 
0.0958
0.0739
0.0230
0.0265 −3.2740

FITTING PHASE–TYPE SCALE MIXTURES TO HEAVY–TAILED DATA AND DISTRIBUTIONS

0

15

c=1
c=1/2
c=1/4

-1

-2
0

5

10

15

20

0

c=1
c=1/2
c=1/4

-2

-4
0

100

200

300

400

F IGURE 2. Estimated survival functions of the NPH models for Danish reinsurance
data corresponding to θ equal to 1.2748, 1.3115, 1.3387.

Again this representation does not resemble any of the previous estimations (due to non–uniqueness
of representations) but from figures 3 and 4 it is clear that there is no distinguishable difference
between densities for the fixed and EM adjusted tails in the range [0, 6] where the main body of the
distribution is situated. The tail behaviour also looks quite similar though the EM fitted tail seems
to be slight heavier than the fixed tail.

0.8

Data
EM Adjusted Tail
Fixed Tail

0.6
0.4
0.2
0
0

2

4

6

F IGURE 3. Estimated densities of the NPH models with EM adjusted (θ̂ = 1.3387)
and fixed (θ̂ = 1.45) tails against histogram of the Danish reinsurance data.
Example 6.3 (Fitting to a theoretical distribution: Loggamma). Next we consider the problem of
approximating a theoretical distribution via an NPH model. In our first example we take as target a
log–gamma distribution with shape parameter α, scale parameter β and shifted one unit to the left

16

M. BLADT AND L. ROJAS

0

EM Adjusted Tail
Fixed Tail

-1

-2
0

5

10

15

20

0

EM Adjusted Tail
Fixed Tail

-2

-4
0

100

200

300

400

F IGURE 4. Estimated survival functions of the NPH models with EM adjusted (θ̂ =
1.3387) and fixed (θ̂ = 1.45) tails.

so its support is [0, ∞), which is the natural support for the class NPH. Its density function is then
given by
β α logα−1 (x + 1)
f (x) =
,
x ≥ 0.
Γ(α)(x + 1)β +1
This distribution is regularly varying with parameter β . For this example we choose the parameters
α = β = 2, for the target distribution with the purposes of analysing a distribution with a mode away
from 0 and a moderately heavy–tail. We consider an NPH model where the underlying phase–type
distribution has five phases and the scaling distribution is that of Example 3.5 with c = 1. With this
example, we want to test if general Regularly Varying distributions can be correctly fitted with the
general model suggested in Example 3.5.
We employed Quasi Monte Carlo ideas to approximate the integrals in the EM steps and iterated
the algorithm until the relative error was smaller than 10−9 . The results are given below
θ̂ = 1.6031
α = (0.5717, 0.0330, 0.0000, 0.3954, 0.0000)
α̂


−1.9634 0.0609
0.5025
0.1249
1.2751
 0.0616 −0.3372 0.0775
0.0382
0.1428 



T =  0.7529
0.1178 −2.2723 0.4797
0.0068 
T̂

 0.7278
0.3060
1.1458 −4.8966 2.7170 
0.8923
0.0317
0.0482
0.2021 −3.4321
The densities of the log-gamma distribution and its NPH approximation are plotted in Figure 5.
The densities of the log-gamma and its NPH approximation are almost indistinguishable from each
other in the body region. The tail behavior is correctly captured by the NPH model as seen in Figure
6. The shape of the tail of the NPH is very close to that of the Loggamma distribution although the
NPH estimate has a heavier tail (the estimated regularly varying parameter was θ̂ = 1.6031).

FITTING PHASE–TYPE SCALE MIXTURES TO HEAVY–TAILED DATA AND DISTRIBUTIONS

17

EM
LG(2,2)

0.4
0.2
0
0

1

2

3

4

F IGURE 5. Density of the Loggamma distribution with parameters (2, 2) and an
approximation via an NPH model with EM adjustment.

0

EM
LG(2,2)

-1
-2
0

10

20

30

40

50

0

EM
LG(2,2)

-2

-4
0

100

200

300

400

F IGURE 6. Survival functions of the Loggamma distribution LG(2, 2) and an approximation via an NPH model with EM adjustment.
Example 6.4 (Fitting to a theoretical distribution: Weibull). Next we move away from the Regularly Varying case and consider instead the Weibullian case. As a target distribution
we consider a
√
√
−
x
classical two–parameter Weibull with λ = 1 and p = 1/2 so its density is e
/2 x, x ≥ 0. For
adjusting this model we consider an NPH family of distributions where the phase–type part has five
phases and the scaling distribution is supported over a geometric progression ec , e2c , e3c , . . . with
c = 1 and taken as a discretization of a classical two–parameter Weibull distribution with δ = p − 1
More precisely, its density is given by
p

f (x) = pλ p x p−1 e−(λ x) ,

x > 0.

Notice, that the target distribution is in the same two–parameter family of Weibull distributions.
The results are given below.
λ̂ = 0.6181, p̂ = 1.1673
α = (0.2370, 0.3349, 0.0901, 0.3380, 0)
α̂

18

M. BLADT AND L. ROJAS




−3.0858
0.6964
0.2188
0.3355
0.5572
 51.0253 −207.2799 18.3961 18.5851 46.0724 



T = 
0.4369
0.2511
−0.9839
0.0122
0.0487
T̂


 1.1297
0.2448
0.9388 −11.0548 0.7637 
0.7583
1.0994
0.9844
0.4882 −3.3303
The agreement between the Weibull distribution and its NPH approximation is very good in the
body of the distribution. The approximation of the tail is also particularly good for values going
up to 100 which correspond to probabilities of order 10−4 . The NPH adjusted model is Weibullian
with parameter pb(1 + pb)−1 ≈ 1.1673/2.1673 ≈ 0.5386. Recall that the parameter of the target
distribution is 0.5.

15

EM
Weibull(1,0.5)

10
5
0
0

0.2

0.4

0.6

0.8

1

F IGURE 7. Density of the Weibull distribution with parameters (1, 0.5) and an approximation via NPH models with EM adjustment.

0

EM
Weibull(1,0.5)

-1
-2
-3
0

10

20

30

40

50

0

EM
Weibull(1,0.5)

-2
-4
-6
-8
0

100

200

300

400

F IGURE 8. Survival functions of the Weibull(1, 0.5) and an approximation via an
NDPH model with EM adjustment. The EM model is Weibullian with shape parameter p = pb(1 + pb) ≈ 0.5386.

FITTING PHASE–TYPE SCALE MIXTURES TO HEAVY–TAILED DATA AND DISTRIBUTIONS

19

Example 6.5 (Fitting to a theoretical distribution: Lognormal). Finally we consider a Lognormal
distribution with location parameter µ = 0 and dispersion parameter σ 2 = 1. The lognormal distribution is heavy–tailed but its survival function ultimately decays faster than a Regularly Varying
function (but slower than a Weibullian function). We consider the lognormal case more difficult
because a scaling distribution π (θθ ) for which the NPH model has the same tail behavior as the
lognormal distribution is unknown.
Therefore, we consider two alternative NPH models to adjust the data. For both we have selected
a phase–type part has eight phases and the scaling distribution is chosen as a discretization of a
Lognormal distribution LN(µ, σ 2 ) and supported over the set {si = ei : i = 0, 1, . . . }. The difference
between the two models is that for the first one we let the EM algorithm to estimate the values of
the parameters µ and σ , while for the second one we take these values to be fixed and equal to 0
and 1 respectively. The results for the first model are given below
µ̂ = 0.1909, σ̂ = 0.5979
α = (0.0000, 0.0000, 0.0000, 0, 0.0000, 0.0439, 0.1360, 0.8201)
α̂


−7.1741 0.1557
0.3338
6.2239
0.0195
0.1985
0.1061
0.1366
 0.2056 −1.1622 0.0528
0.2434
0.1133
0.0519
0.0328
0.0782 


0.0891 −6.2221 4.7969
0.1873
0.2546
0.1238
0.0176 
 0.7528


0.1133
0.1303 −9.4035 0.1845
0.0654
0.2161
0.1597 
 0.3840
T = 
T̂

0.4274
1.5430
0.9113 −5.2223 0.2926
1.0875
0.0056 
 0.9550
 0.2341
0.0842
2.3467
0.0000
1.3606 −4.6381 0.0621
0.5504 


 1.8359
0.0083
1.0178
0.0000
0.0211
0.2123 −3.3825 0.2870 
4.9757
0.0000
0.9689
0.0000
3.4702
0.2689
0.3086 −9.9923

EM
EM fixed
LN(0,1)

0.6
0.4
0.2
0
0

1

2

3

4

F IGURE 9. Density of the Lognormal distribution with parameters (0, 1) and an
approximation via NPH models with EM adjustment.
From Figure 9 we can observe that the adjustment of the body of the distributions of the first
model is excellent, but the tail probabilities differ. In the lognormal case, the heaviness of the
distribution is mostly determined by the parameter σ . The larger the parameter σ the more heavier
its tail. In our estimations we obtained an estimate σb = 0.5979 which suggest a lighter tail, and this
is confirmed in the last panel of Figure 10.
For our second model we obtained a very poor fit, but this is somewhat expected since the tail
distribution of the NPH model is very different from the tail behavior of its scaling distribution
(as in the case of Regularly Varying or Weibullian cases). In fact, [11] demonstrate that the tail
probability of the NPH model is significantly heavier and different to the scaling distribution. More

20

M. BLADT AND L. ROJAS

0

EM
EM fixed
LN(0,1)

-2
-4
0

10

20

30

40

50

0

EM
EM Fixed
LN(0,1)

-5

-10
0

100

200

300

400

F IGURE 10. Survival functions of the Lognormal distribution LN(0, 1) and an approximation via an NPH model with EM adjustment.
precisely, it is shown that if N ∼ LN(0, 1), then
lim

t→∞

P(τN > t)
= ∞.
P(N > t)

7. C ONCLUSIONS
In this paper we have introduced a new methodology for adjusting phase–type scale mixture
distributions to heavy–tailed data. While the model only requires the specification of the dimension
of the underlying phase–type distribution and a parametric family of discrete scaling distributions,
the suggested algorithm simultaneously fit both the body and the tail of general distributions.
Since the class of NPH distributions are generally genuinely heavy tailed (if N has unbounded
support) and dense in the class of heavy tailed distributions with support on R+ , we may in principle approximate any heavy tailed distribution (data) arbitrarily close by a NPH distribution. In
particular, for the case of regular varying and Weibullian distributions the aforementioned approximation is not only in the limit (denseness) but can be effectively carried out in praxis.
R EFERENCES
[1] Marek Arendarczyk and Krzysztof Dȩbicki. Asymptotics of supremum distribution of a gaussian process over a
weibullian time. Bernoulli, 17(1):194–210, 2011.
[2] S. Asmussen, O. Nerman, and M. Olsson. Fitting Phase-Type Distributions via the EM Algorithm. Scandinavian
Journal of Statistics, 23:419–441, 1996.
[3] M Bladt, B. F. Nielsen, and G. Samorodnitsky. Calculation of ruin probabilities for a dense class of heavy tailed
distributions. Scandinavian Actuarial Journal, pages 573–591, 2015.
[4] Mogens Bladt. A review on phase-type distributions and their use in risk theory. Astin Bulletin, 35(1):145–161,
2005.
[5] Arnoldo Frigessi, Ola Haug, and Håvard Rue. A dynamic mixture model for unsupervised tail estimation without
threshold selection. Extremes, 5(3):219–235, 2002.
[6] McNeil Alexander J. Estimating the Tails of Loss Severity Distributions Using Extreme Value Theory. ASTIN
Bulletin, 27(1):117–137, 1997.

FITTING PHASE–TYPE SCALE MIXTURES TO HEAVY–TAILED DATA AND DISTRIBUTIONS

21

[7] Guy Latouche and Vaidyanathan Ramaswami. Introduction to matrix analytic methods in stochastic modeling,
volume 5. Society for Industrial Mathematics, 1987.
[8] Marcel F Neuts. Matrix-geometric Solutions in Stochastic Models. An Algorithmic Approach. Courier Dover
Publications, 1981.
[9] M Olsson. Estimation of phase-type distributions from censored data. Scandinavian Journal of Statistics, pages
443–460, 1996.
[10] Sidney I. Resnick. Discussion of the Danish Data on Large Fire Insurance Losses. ASTIN Bulletin, 27(1):139–151,
1997.
[11] Leonardo Rojas-Nandayapa and Wangyue Xiei. Asymptotic tail behavior of phase-type scale mixture distributions. Submitted.
[12] Charles Van Loan. Computing integrals involving the matrix exponential. IEEE Transactions on Automatic Control, 23(3):395–404, 1978.
I NSTITUTE FOR A PPLIED M ATHEMATICS AND S YSTEMS , NATIONAL U NIVERSITY OF M EXICO , A.P. 20-726,
01000 M EXICO , D.F., M EXICO
E-mail address: bladt@sigma.iimas.unam.mx
S CHOOL OF MATHEMATICS AND PHYSICS , T HE U NIVERSITY OF Q UEENSLAND , S T. L UCIA 4072, B RISBANE ,
AUSTRALIA
E-mail address: l.rojas@uq.edu.au

