Multivariate Discrete Phase-Type Distributions with Their

Applications ^{∗}

Matt Goff Department of Mathematics Washington State University Pullman, WA 99164

March 15, 2001

Abstract

A dynamic random environment usually induces statistical dependence among components of a system operating in such an environment thereby making explicit computation of the system performance measures very difficult and possibly even intractable. The multivariate phase-type distribution is a powerful tool for unifying a variety of stochastic models and for constructing new models that yield to algorithmic analysis in such applications. In this project, we introduce the multivariate discrete phasetype distribution and use it to study the dependence structures of multi-component systems operating in dynamic random environments. The advantages of multivariate discrete phase-type distributions stem from the facts that such distributions can be easily computed and they have close ties with multivariate continuous phase-type distributions.

In this project, detailed properties of multivariate discrete phase-type distributions will be established, and computational algorithms, both matrix algorithms and simulation procedures, will be developed. Using stochastic comparison methods and computational experimentation, the dependence structures of multivariate discrete phase-type distributions will be throughly examined.

The methods and results to be obtained in this project will also be applied to analyze various operating systems, including assemble-to-order systems, synchronized communication networks, and reliability systems with simultaneous subsystem failures. The objectives are to develop good computable bounds, efficient heuristics, and important characterizations of system performance measures whose explicit computations are either intractable or computationally intensive. The insights gained from the proposed research will be summarized in a set of design and control principles for the efficient and effective management of synchronized operating systems.

∗Supported in part by the NSF grant DMI 9812994

1 Introduction

Many stochastic models that have matrix-geometric solutions involve in one way or another (univariate or multivariate) probability distributions of phase-type. However, the closed form solutions for such stochastic models with multivariate phase-type distributions are in general computationally intensive. In this project, we study multivariate phase-type distributions and related stochastic models using stochastic comparison methods and computational experimentations in order to better understand their dependence nature. The objectives are to develop good computable bounds, efficient heuristics, and important characterizations of the multivariate phase-type distributions and system performance measures of related stochastic models.

A nonnegative random variable T (or its distribution function) is said to be of phase-type (PH) if T is the time until absorption in a finite-state continuous-time Markov chain (Neuts 1981). The theory of univariate PH distributions has been developed by several researchers and has been found to have various applications in queueing theory and reliability theory (see Neuts 1981, Assaf, Langberg, Savits and Shaked 1984 and references therein). Univariate PH distributions that they (and their densities, Laplace transforms and all moments) can be written in closed form in terms of the parameters of the underlying Markov chains. Furthermore the set of univariate PH distributions is dense in the set of all distributions with support in [0, ∞).

Several classes of multivariate phase-type (PH) distributions have also been introduced and studied in the literature (Assaf, Langberg, Savits and Shaked 1984, and Kulkarni 1989). To explain the most interesting class of multivariate PH distributions, consider a continuous-time and right-continuous Markov chain X = {X(t), t ≥ 0} on a finite state space E with generator Q. Let E_{i}, i = 1, …, m, be m nonempty stochastically closed subsets of E such that ∩_{i=1}^{m}E_{i} is a proper subset of E (A subset of the state space is said to be stochastically closed if once the process X enters it, X never leaves). We assume that absorption into ∩_{i=1}^{m}E_{i} is certain. Since we are interested in the process only until it is absorbed into ∩_{i=1}^{m}E_{i}, we may assume, without loss of generality, that ∩_{i=1}^{m}E_{i} consists of one state, which we shall denote by Δ. Thus, without loss of generality, we may write E = {1, 2, …, s} ∪ (∪_{i=1}^{m}E_{i}) for some integer s ≥ 1. Let β be an initial probability vector on E such that β(Δ) = 0. Then we can write β = (α, 0). Define

T_{i} = inf{t ≥ 0 : X(t) ∈ E_{i}}, i = 1, …, m.  (1.1)

As in Assaf et al. (1984), for simplicity, we shall assume that P (T_{1} > 0, T_{2} > 0, …, T_{m} > 0) = 1. The joint distribution of (T_{1}, T_{2}, …, T_{m}) is called a multivariate phase-type distribution with representation (α, Q, E_{i}, i = 1, …, m), and (T_{1}, T_{2}, …, T_{m}) is called a (multivariate) phase-type random vector. Let MP H denote the class of all the multivariate phase-type distributions. As in the univariate case, those multivariate PH distributions (and their densities, Laplace transforms and moments) can be written in a closed form. The set of m-dimensional PH distributions is dense in the set of all distributions on [0, ∞)^{m}.

An example of multivariate PH distributions is the Marshall-Olkin distribution (Marshall and Olkin 1967); one of the most widely discussed multivariate life distributions in reliability theory (see, for example, Barlow and Proschan 1981). The multivariate Marshall-Olkin distribution has also been used recently to model certain correlated queueing systems, such as assemble-to-order systems (see Li and Xu (1999) and the references therein). The multivariate PH distributions defined by (1.1), their properties, and some related applications in reliability theory were discussed in detail in Assaf et al. (1984), and Li and Xu (2000). Another class of multivariate PH distributions, based on the total accumulated reward until absorption in a continuous-time Markov chain, was introduced in Kulkarni (1989).

The approach used in this project is different from that in the literature (see, for example, Li (1999)). We first introduce multivariate discrete phase-type distributions and decompose the multivariate continuous phase-type distribution into the multivariate discrete phase-type distribution and Erlang distributions. The multivariate phase-type distributions and related stochastic models can be studied by using this relation and the fact that the discrete phase-type distribution is much easier to analyze.

2 Multivariate Discrete Phase-Type Distributions

Let {X_{n}, n ≥ 0} be a discrete-time Markov chain with a finite state space E. Let Γ_{1}, …, Γ_{m} be non-empty stochastically closed subsets such that ∩_{k=1}^{m}Γ_{k} ≠ ∅. Without loss of generality, we assume that ∩_{k=1}^{m}Γ_{k} consists of a single state Δ. We also assume that absorption in ∩_{k=1}^{m}Γ_{k} is certain. Let P = {p_{ij}} be the transition matrix of the process. Hence, p_{ij} = 0 if i ∈ Γ_{k} and j ∈ Γ_{k}^{c} for some k.

Let α be an initial probability vector on E with α(Δ) = 0. Define

T_{j} = inf{n ≥ 0 : X_{n} ∈ Γ_{j}}, j = 1, …, m.  (2.1)

Definition 2.1 The joint distribution of (T_{1}, …, T_{m}) defined by (2.1) is called a multivariate discrete phasetype distribution.

Example 2.2 Consider an m-component system that is subject to multiple types of shocks. Assume that a Bernoulli process N(t) with parameter p^{K}, K ⊆ E, governs the occurrence of shocks fatal to the system components, where 0 ≤ p^{K} ≤ 1 and ∑_{K⊆E} p^{K} = 1. At each arrival epoch, with probability p^{K}, the shock destroys all the components in subset K ⊆ E simultaneously. With probability p^{∅}, all the components survive at each arrival epoch. Let T_{j}, j = 1, …, m, be the lifetime of component j. The joint distribution of (T_{1}, …, T_{m}) is called the multivariate discrete Marshall-Olkin distribution.

To see that this distribution is a special case of multivariate discrete phase-type distributions, we let {X_{n}, n ≥ 0} be a Markov chain with state space {K : K ⊆ E}, and starting at X_{0} = ∅. State K ⊆ E means that all the components in K have failed and the others are still operational. The transition probabilities of the chain is given by

p_{KL} = ∑_{K^{c}∩L⊆J⊆L} p^{J}, K ⊆ L.

Let Γ_{j} = {K : j ∈ K ⊆ E}, j = 1, …, m be the set of failure states of component j. Clearly,

T_{j} = inf{n : X_{n} ∈ Γ_{j}}, j = 1, …, m.

Thus, the joint distribution of (T_{1}, …, T_{m}) is a multivariate discrete phase-type distribution.

Example 2.3 Consider a system of m components operating in a random environment. Let T_{i} denote the lifetime of component i, i = 1, …, m. A component fails when it receives a fatal shock from the random environment. As long as all the components are functioning, the shock arrival process to the system is a Bernoulli processes with parameters α_{j}^{∅}, j ∈ E, where 0 ≤ α_{j}^{∅} ≤ 1 and ∑_{j} α_{j}^{∅} = 1. With probability α_{j}^{∅}, an arriving shock destroys component j, j = 1, …, m, but no simultaneous failure occurs at any arrival epoch. Upon the first component failure, say component j_{1}, the shock arrival process to the functioning components j, j ∈ E, j ≠ j_{1}, change to a Bernoulli process with parameters α_{j}^{j_{1}}. Suppose component index j_{i} corresponds to the ith component failure, i ∈ E. Upon the failure of the lth component, with the order of failures given by j_{1} … j_{l}, the shock arrival process to the remaining functioning components change to a Bernoulli process with parameters α_{j}^{j_{1}...j_{l}}, j ∈ E, j ≠ j_{i}, i = 1, …, l. With probability α_{j}^{j_{1}...j_{l}}, an arriving shock destroys component j, j ∈ E, j ≠ j_{i}, i = 1, …, l, and no simultaneous failure occurs at any arrival epoch. Clearly, the lifetime vector (T_{1}, …, T_{m}) has dependent components. When m = 2, the joint distribution of the lifetime vector is a discrete version of the bivariate extension of the exponential distribution introduced by Freund (1961).

To see that this model is a special case of the multivariate discrete phase-type distributions we introduced above, let the state space E contain all the permutations of each subset of E. The initial state of the Markov chain {X_{n}, n ≥ 0} is ∅. Its transition probability matrix is given as follows, for any state K = j_{1} ⋯ j_{l}, l ∈ E,

p_{KL} = α_{j}^{K}, where L = j_{1} … j_{l}j,

and zero otherwise. Let the set of failure states of component j be {K : j ∈ K} where j ∈ K means that j belongs to the set of indices of permutation K. Clearly

T_{j} = inf{n : X_{n} ∈ Γ_{j}}, j = 1, …, m.

Thus, the joint distribution of (T_{1}, …, T_{m}) is a multivariate discrete phase-type distribution.

Theorem 2.4 Let (T_{1}, …, T_{m}) be a multivariate discrete phase-type distribution with the underlying Markov chain X = {X_{n}, n ≥ 0}. Assume that X has the initial probability vector α = (β, 0) and the following transition matrix,

P = [ \begin{matrix}A & B \\ 0 & 1 \end{matrix} ] ,

where A is a (|E| − 1) × (|E| − 1) matrix, and B is a (|E| − 1) × 1 matrix. Then for any n_{1} ≥ n_{2} ≥ … ≥ n_{m}, we have

P (T_{1} > n_{1}, …, T_{m} > n_{m}) = βA^{n_{m}}g_{m}A^{n_{m−1}−n_{m}}g_{m−1} … A^{n_{1}−n_{2}}g_{1}e,

where g_{k} is a diagonal (|E| − 1) × (|E| − 1) matrix whose the ith diagonal element is 1 if i ∈ Γ_{k}^{c} and zero otherwise, and e denotes the vector of 1′s.

Several properties of multivariate discrete phase-type distributions are in order.

Theorem 2.5 Let (T_{1}, …, T_{m}) and (S_{1}, …, S_{n}) be independent multivariate discrete phase-type vectors. Then the conjunction (T_{1}, …, T_{m}, S_{1}, …, S_{n}) is a multivariate discrete phase-type vector.

Theorem 2.6 The class of multivariate discrete phase-type distributions is closed under finite mixture.

Theorem 2.7 All marginals of a multivariate discrete phase-type distribution are multivariate discrete phase-type distributions.

Theorem 2.8 Let (S_{1}, …, S_{m}) be a random vector with the multivariate phase-type distribution as defined in Assaf et al. (1984) (see (1.1)). Then we have

S_{j} = ∑_{i=1}^{T_{j}} E_{i}, j = 1, …, m,  (2.2)

where (T_{1}, …, T_{m}) is defined as in Definition 2.1 and E_{i}, i = 1, 2, … are independent and identically distributed with exponential distributions.

The dependence structures (dependence properties and bounds) of (S_{1}, …, S_{m}) and related models can be analyzed by using this relation and the structural properties of (T_{1}, …, T_{m}).

We believe this approach will be effective because the multivariate discrete phase-type distributions can easily be analyzed both theoretically (using the stochastic comparison methods (see, for example Joe (1997)) and numerically (using matrix methods and Markov chain simulations). The knowledge obtained from this research will help fill important voids in understanding the dependence nature of the multivariate phase-type distributions, and can be also applied to analyze various practical operating systems, including assemble-to-order systems, synchronized communication networks, and reliability systems with simultaneous subsystem failures.

3 Computational Issues

Although the multivariate phase-type distributions can be expressed in the closed form in terms of the initial probability distribution and transition kernel of the underlying Markov chain, their joint distributions and higher moments are in general computationally intensive. Kulkarni (1989) developed the methods using partial differential equations and Laplace Stieltjes transforms to compute the joint distributions of a general class of multivariate phase-type distributions. However, as mentioned in Kulkarni (1989), these two methods both face numerical difficulties. In this research, we propose a simulation method for multivariate discrete phase-type distributions, and then use the relation (2.2) to simulate the joint distribution of multivariate phase-type distributions.

Let (T_{1}, …, T_{m}) have a multivariate discrete phase-type distribution with the underlying Markov chain X = {X_{n}, n ≥ 0}. Let α be the initial distribution of X, and P = {p_{ij}} be the transition matrix of X. Discrete Algorithm:

1. Initially, sample an i_{0} from distribution α, and set k = 0.

2. Let X_{k} = i_{k}. Sample an i_{k+1} from the conditional distribution (p_{i_{k}1}, p_{i_{k}2}, …). If i_{k+1} ∈ Γ_{j} for the first time, then T_{j} = k + 1.

3. If i_{k+1} ∈ ∩_{j=1}^{m}Γ_{j} then stop. Otherwise, set k = k + 1 and go to (2).

Use the discrete algorithm to generate n independent copies (T_{1}^{(j)}, …, T_{m}^{(j)}), j = 1, …, n, of (T_{1}, …, T_{m}). Then for any function f : R^{m} → R, we have

\frac{∑_{j=1}^{n} f(T_{1}^{(j)}, …, T_{m}^{(j)})}{n} → Ef(T_{1}, …, T_{m}),

as n → ∞. The joint distribution and higher moments of a multivariate phase-type distribution can be computed by using this discrete algorithm and the relation (2.2). An error estimate for this method and comparison of this method with the matrix method are also proposed for study.

Since higher joint moments of multivariate phase-type distributions are of particular interest, we also propose to find an algorithm for computing these higher joint moments of discrete phase-type distributions recursively, using the idea from Laplace Stieltjes transform method discussed in Kulkarni (1989).

4 Dependence Structures of Multivariate Phase-Type Distribu-

tions

Dependence structures of multivariate phase-type distributions can be examined using stochastic comparison methods and the numerical methods we have proposed. Numerical experimentations on multivariate phasetype distributions will shed new light in analyzing their dependence structure in addition to the properties which can be established through various comparison techniques.

Many different notions of dependence have been introduced and studied extensively in the literature (see, for example, Tong 1980, and Joe 1998). Here we only discuss certain notions of dependence orders between two random vectors and dependence among the components of a random vector that are most relevant to this research.

To compare the dependence structures between two random vectors X = (X_{1}, …, X_{m}) and Y = (Y_{1}, …, Y_{m}) in the weakest sense, one may compare their covariance matrices element-wise. However, some stronger notions of dependence orders have been developed in the literature. These are summarized in the following definition.

Definition 4.1 Let X = (X_{1}, …, X_{m}) and Y = (Y_{1}, …, Y_{m}) be two R^{m}-valued random vectors.

1. X is said to be larger (smaller) than Y in the upper (lower) orthant order, denoted by X ≥_{uo} (≤_{lo}) Y, if P (X > (≤) x) ≥ P (Y > (≤) x), for all x ∈ R^{m}.

2. X is said to be more upper (lower) orthant dependent than Y if X ≥_{uo} (≤_{lo}) Y, and X_{i} =_{st} Y_{i} for all i. (‘=_{st}’ denotes the equality in distribution).

3. X is said to be larger than Y in the usual stochastic order, denoted by X ≥_{st} Y, if P (X ∈ U) ≥ P (Y ∈ U) for all upper sets U ⊆ R^{m} (U is said to be upper if x ∈ U and x ≤ y imply that y ∈ U).

The following facts are easy to verify (see, for example, Tong 1980):

X ≥_{st} Y ⟹ X ≥_{uo} Y and X ≥_{lo} Y

Note that the orthant orders, coupled with identical marginals, emphasize the comparisons of dependence strengths of the two vectors by separating the marginals from consideration.

To express the nature of dependence among the components of a random vector, one can of course use its covariance matrix. X is said to be positively (negatively) pairwise-correlated if the covariance matrix Cov(X_{i}, X_{j}) ≥ (≤) 0, for all i ≠ j. However, the following notions of dependence are stronger and frequently used in the literature. Definition 4.2 Let X = (X_{1}, …, X_{m}) be an R^{m}-valued random vector. Let X^{I} = (X_{1}^{I}, …, X_{m}^{I}) denote a vector of real random variables such that X_{j} =_{st} X_{j}^{I} for each j and X_{1}^{I}, …, X_{m}^{I} are independent.

1. X is said to be positively upper (lower) orthant dependent (PUOD, PLOD) if X ≥_{uo} (≤_{lo}) X^{I}. X is said to be negatively upper (lower) orthant dependent (NUOD, NLOD) if X ≤_{uo} (≥_{lo}) X^{I}.

2. X is said to be associated if Cov(f(X), g(X)) ≥ 0 whenever f and g are nondecreasing. X is said to be negatively associated if for every subset K ⊆ {1, …, m}, Cov(f(X_{i}, i ∈ K), g(X_{j}, j ∈ K^{c})) ≤ 0 whenever f and g are nondecreasing.

It is known that (Esary, Proschan and Walkup 1967, Szekli 1995)

X is (positively) associated ⟹ X is PUOD and PLOD,  (4.1)
X is negatively associated ⟹ X is NUOD and NLOD.  (4.2)

The PLOD (PUOD, NLOD, NUOD)-property of a random vector means that its joint distribution or survival function can be bounded below or above by the products of its marginal distributions or survival functions. Various properties of these classifications have been discussed in Tong 1980, and Szekli 1995.

The following problems are proposed for study.

1. Relation (2.2) naturally leads to the question of how many dependence properties multivariate phasetype distributions inherit from the underlying multivariate discrete phase-type distributions? We will study the relationship of dependence properties of multivariate continuous and discrete phase-type distributions.

2. It is known from Li (1999) that the multivariate phase-type distribution is associated if the underlying Markov chain is stochastically monotone. But this result depends on a result from Harris (1977), which does not hold for the discrete-time Markov chain. We propose to extend this result to the discrete case by using the theory of association of Markov chains from Lindqvist (1987).

3. It is known that multivariate Freund distributions with certain structures are negatively correlated. We will apply the discrete algorithm to multivariate discrete Freund distributions with the goal of identifying the relationship between parameters which leads to negative dependence.

4. Study the comparison results on dependence structures of multivariate discrete Marshall-Olkin distributions using orthant comparison techniques.

5 Applications

Since any nonnegative random vector can be approximated by a sequence of random vectors with multivariate phase-type distributions, the multivariate phase-type distribution should be a powerful tool for stochastic modeling of various applications. In fact, some simple multivariate phase-type distributions, such as Marshall-Olkin distribution and Freund distribution, have found many applications in reliability and queueing modeling (see, for example, Li (1999) and Li and Xu (2000) and the references therein). However, the numerical difficulties and relatively unknown dependence structure limit the widespread use of multivariate phase-type distributions in applications. The multivariate discrete phase-type distribution shares a similar conceptual framework as the continuous version, but it is relatively straightforward to compute numerically. It should therefore be an effective tool in applications that involve discrete-time stochastic models. We propose the following three applications for further study.

1. Cumulative Shock Models with Simultaneous Damages. Li and Xu (2001) studied the dependence structure of a general class of continuous-time multi-component shock models with simultaneous failures. Here we propose to study a multi-component cumulative shock model where the shocks arrive according to a discrete-time process and any arrival shock inflicts random damage on certain subset of components simultaneously. The dependence properties of the system components are of interest as well as the explicit expression of higher joint moments of component lifetimes. These can be obtained by using the method we proposed in Section 3.

2. Shock Models with Monotone Markovian Arrivals. Consider a shock model in which the shocks arrive according to a discrete-time process. Any arrival shock, with certain probability, destroys some subset of components simultaneously. This failure pattern could be changed at any shock arrival epoch. Such a shock model combines the features from both Marshall-Olkin distribution and Freund distribution. Our goal is to obtain some computable bounds for dependent component lifelengths using the idea of association.

3. Discrete-Time Queueing Systems with Correlated Arrivals. A discrete-time m parallel GI/M/1 queueing system with correlated arrivals has many applications in assemble-to-order systems and synchronized communication networks. Our goal is to calculate joint system performance measures such as the covariance matrix of the queue length vector.

References

[1] Assaf, D., Langberg, N., Savits, T., and Shaked, M. (1984), Multivariate phase-type distributions, Operations Research 32, 688–701.

[2] Barlow, R. E. and Proschan, F. (1981), Statistical Theory of Reliability and Life Testing, To Begin With, Silver Spring, MD.

[3] Esary, J. D. and F. Proschan (1970), A reliability bound for systems of maintained, independent components, Journal of American Statistical Association, 65, 329-338.

[4] Freund, J. E. (1961), A bivariate extension of the exponential distribution, J. Am. Statist. Asso. 56, 971–977.

[5] Harris, T. E. (1977). A correlation inequality for Markov processes in partially ordered state spaces. Ann. Prob. 5, 451-454.

[6] Joe, H. (1998), Multivariate Models and Dependence Concepts,

[7] Kulkarni, V. G. (1989), A new class of multivariate phase type distributions, Operations Research 37, 151–158.

[8] Li, H. (1999), Association of multivariate phase-type distributions with applications to shock models. Preprint.

[9] Li, H. and Xu, S. H. (2000) Directionally convex comparison of correlated first passage times. Preprint.

[10] Li, H. and Xu, S. H. (1999), On the dependence structure and bounds of correlated parallel queues and their applications to synchronized stochastic system. Journal of Applied Probability, to appear.

[11] Lindqvist, H. (1988), Monotone and Associated Markov Chains, with Applications to Reliability Theory, Journal Appl. Prob., 24, 679-695.

[12] Neuts, M. F. (1981), Matrix-Geometric Solutions in Stochastic Models: An Algorithmic Approach, The Johns Hopkins University Press, Baltimore.

[13] Marshall, A. W. and Olkin, I. (1967). A multivariate exponential distribution. J. Amer. Statist. Assoc. 62, 30-44.

[14] Szekli, R. (1995), Stochastic Ordering and Dependence in Applied Probability, Springer-Verlag.

[15] Tong, Y.L. (1980), Probability Inequalities in Multivariate Distributions, Academic Press, New York.