3. Univariate phase-type distributions
3.1. The continuous univariate phase-type
Let \(\{X_t\}_{t\ge0}\) be a Markov jump process with \(p\) transient states and one absorbing state. If \(\tau\) is defined as the time when \(\{X_t\}_{t\ge0}\) reaches the absorbing state, then \(\tau\) follows a continuous phase-type distribution \(\tau \sim PH(, \boldsymbol{\pi},\boldsymbol{T})\), where \(\boldsymbol{T}\) is the sub-intensity matrix of the corresponding Markov jump process, and \(\boldsymbol{\pi}\) is the vector of initial probabilities.
As specified before, a continuous phase-type distribution object can
be created in PhaseTypeR using the PH()
function:
subintensity_matrix <- matrix(c(-1.5, 0, 0,
1.5, -1, 0,
0, 1, -0.5), ncol = 3)
initial_probabilities <- c(0.9, 0.1, 0)
ph <- PH(subintensity_matrix, initial_probabilities)
ph
#> $subint_mat
#> [,1] [,2] [,3]
#> [1,] -1.5 1.5 0.0
#> [2,] 0.0 -1.0 1.0
#> [3,] 0.0 0.0 -0.5
#>
#> $init_probs
#> [,1] [,2] [,3]
#> [1,] 0.9 0.1 0
#>
#> $defect
#> [1] 0
#>
#> attr(,"class")
#> [1] "cont_phase_type"We can easily plot the phase-type distribution as a network by first converting the phase-type object into a graph, and then plotting it:
library(igraph)
#>
#> Attaching package: 'igraph'
#> The following objects are masked from 'package:stats':
#>
#> decompose, spectrum
#> The following object is masked from 'package:base':
#>
#> union
net <- phase_type_to_network(ph)
set.seed(8)
plot(net, edge.curved=.1, edge.label = E(net)$weight,
edge.color = ifelse(as_data_frame(net, what="edges")$from == 'V0',
'purple', 'grey'),
layout = layout_with_fr(net, weights = rep(1, length(E(net)))))
In the above graph the initial probabilities are represented as purple arrows, while rates are represented by gray arrows.
In the continuous case, we can also plot the transition probability matrix at time t as a graph:
oldpar <- par()[c('mfrow', 'mar')]par(mfrow=c(2,3), mar=c(2,2,2,2))
set.seed(6)
for (i in c(0, 1, 3, 5, 10, 20)) {
net <- phase_type_to_network(ph, i)
plot(net, edge.arrow.size=0.5, edge.curved=.1, edge.label = E(net)$weight, main=paste0('time = ', i),
edge.color = ifelse(as_data_frame(net, what="edges")$from == 'V0',
'purple', 'grey'),
layout = layout_with_fr(net, weights = rep(1, length(E(net)))))
}
par(oldpar)The mean of a continuous phase-type distribution is provided in Theorem 3.1.16 and Corollary 3.1.18 in (Bladt and Nielsen 2017, 81:p.135)
\[ \mathbb{E} (\tau) = \boldsymbol{\pi} (-\boldsymbol{T})^{-1} \boldsymbol{e}, \] and from Corollary 3.1.18 (Bladt and Nielsen 2017, 81:p.135), the higher order moments are given by
\[ \mathbb{E} (\tau^n) = n!\boldsymbol{\pi} (-\boldsymbol{T})^{-n} \boldsymbol{e}, \quad n=1,2,\ldots. \]
From this formula we find the variance
\[ \mathbb{V} (\tau) = \mathbb{E} (\tau^2)-\mathbb{E} (\tau)^2= (2\boldsymbol{\pi} (-\boldsymbol{T})^{-2} \boldsymbol{e}) -(\boldsymbol{\pi} (-\boldsymbol{T})^{-1} \boldsymbol{e})^2. \]
Calculating the mean and variance of a phase-type distributed random
variable is straightforward in PhaseTypeR. We can do so by
using the generic functions mean() and var()
on a cont_phase_type object. For example:
cat('\n', 'Mean: ',mean(ph),'\n')
#>
#> Mean: 3.6
cat(' Variance: ',var(ph),'\n \n')
#> Variance: 5.44
#> These two quantities are also printed when using the
summary() function:
summary(ph)
#>
#> Subintensity matrix:
#> [,1] [,2] [,3]
#> [1,] -1.5 1.5 0.0
#> [2,] 0.0 -1.0 1.0
#> [3,] 0.0 0.0 -0.5
#>
#> Initial probabilities:
#> [,1] [,2] [,3]
#> [1,] 0.9 0.1 0
#>
#> Defect:
#> [1] 0
#>
#> Mean: 3.6
#>
#> Variance: 5.44The probability density function (PDF) of a continuous phase-type distribution is provided in Theorem 3.1.7 (Bladt and Nielsen 2017, 81:p.~131)
\[ f(u) = \boldsymbol{\pi}e^{\boldsymbol{T}u}\boldsymbol{t}, \quad u\ge0, \] where \(\boldsymbol{t} = -\boldsymbol{T}\boldsymbol{e}\) is the exit rate vector of the sub-intensity matrix (the vector of rates from state \(i\) to the absorbing state). The cumulative distribution function (CDF) is
\[ F(u) = 1-\boldsymbol{\pi}e^{\boldsymbol{T}u}\boldsymbol{e}, \quad u\ge0. \]
We calculate the PDF and the CDF with dPH() and
pPH(), respectively,
x <- seq(0,8,length.out = 100)
pdf <- dPH(x, ph)
{plot(x, pdf, xlab = "Time to absorption", ylab = "PDF", col = "blue",
type = 'l', lwd=2)
title('Probability density function (PDF)')}
x <- seq(0,8,length.out = 100)
cdf <- pPH(x, ph)
{plot(x, cdf, xlab = "Time to absorption", ylab = "CDF", col = "orange",
type = 'l', lwd=2)
title('Cumulative density function (CDF)')
}
Additionally, we can compute the quantile function using
qPH():
x <- seq(0.05, 0.95, 0.01)
{plot(x, qPH(x, ph), col = "green", type="l", lwd=2, xlim=c(0,1),
ylab = "Time to absorbtion", xlab = "Quantile")
title('Quantile function')}
We can also make random draws from the phase-type distribution using
rPH():
cat('10 random samples: \n', rPH(5, ph), '\n', rPH(5, ph))
#> 10 random samples:
#> 3.162411 1.325381 2.67847 0.8816998 4.542533
#> 1.044944 9.829242 5.278639 1.063361 1.689269Here is a histogram of 200 samples from the phase-type distribution:
hist(rPH(500, ph), breaks = 10)
set.seed(6)
tab_FullPH <- rFullPH(ph)
x <- c(tab_FullPH$time, sum(tab_FullPH$time)/length(tab_FullPH$time))
x <- cumsum(x)
y <- c(tab_FullPH$state, 4)
plot(x,y,type="n", xlim=c(0,x[length(x)]))
segments(c(0,x[-length(x)]),y,x,y)
points(c(0,x[-length(x)]),y,pch=16)
points(x[-length(x)],y[-length(x)],pch=1)
points(x[length(x)],y[length(x)],pch='>')
3.2. The discrete univariate phase-type
Instead of using a Markov Jump Process (also sometimes called a continuous time Markov chain), it is also possible to consider a discrete Markov chain (MC). Let’s consider \(\tau \sim DPH(\boldsymbol{\pi},\boldsymbol{T})\). Here \(\tau\) is the number of state transitions before entering the absorbing state. The sub-intensity matrix \(\boldsymbol{T}\) contains the probabilities of jumping from one state to another.
Discrete phase-type distributions can be created in
PhaseTypeR using the DPH() function by
specifying a sub-intensity matrix that meets the requirements for a
discrete phase-type distribution:
subintensity_matrix <- matrix(c(0, 0.2, 0.8,
0.5, 0.5, 0,
0, 0, 0.4), ncol = 3, byrow = T)
initial_probabilities <- c(0.7, 0.3, 0)
dph <- DPH(subintensity_matrix, initial_probabilities)
dph
#> $subint_mat
#> [,1] [,2] [,3]
#> [1,] 0.0 0.2 0.8
#> [2,] 0.5 0.5 0.0
#> [3,] 0.0 0.0 0.4
#>
#> $init_probs
#> [,1] [,2] [,3]
#> [1,] 0.7 0.3 0
#>
#> $defect
#> [1] 0
#>
#> attr(,"class")
#> [1] "disc_phase_type"We can also plot discrete phase-type functions the same way as continuous ones:
net <- phase_type_to_network(dph)
set.seed(8)
plot(net, edge.curved=.1, edge.label = E(net)$weight,
layout = layout_with_fr(net, weights = rep(1, length(E(net)))))
The moments are given by Proposition 2.7.iii (Campillo Navarro 2019, p.12)
\[ \mathbb{E} (\tau(\tau - 1) \cdots (\tau-(n-1))) = n!\boldsymbol{\pi T}^{n-1} (\boldsymbol{I}-\boldsymbol{T})^{-n} \boldsymbol{e}, \quad n=1,2,\ldots. \]
This leads to a mean and variance given by
\[ \mathbb{E} (\tau) = \boldsymbol{\pi}(\boldsymbol{I}-\boldsymbol{T})^{-1}\boldsymbol{e}, \quad \] \[ \mathbb{E} (\tau(\tau-1))= 2\boldsymbol{\pi}\boldsymbol{T}(\boldsymbol{I}-\boldsymbol{T})^{-2}\boldsymbol{e}, \quad \] \[ \mathbb{V} (\tau) = \mathbb{E} (\tau(\tau-1)) + \mathbb{E} (\tau) - \mathbb{E} (\tau)^2. \]
In PhaseTypeR the mean and variance of a discrete
phase-type distribution can be obtained using the generic functions
mean() and var():
cat('\n', 'Mean: ',mean(dph),'\n')
#>
#> Mean: 4.016667
cat(' Variance: ',var(dph),'\n \n')
#> Variance: 5.863611
#> Theorem 1.2.58 and 1.2.59 in (Bladt and Nielsen 2017, 81:p.30) provides the probability distribution of a discrete phase-type distribution
\[ \mathbb{P}(\tau=n) = \boldsymbol{\pi T}^{n-1}\boldsymbol{t}, \quad n=1,2,\ldots, \]
and the cumulative distribution function is
\[ F_\tau(n) = 1-\boldsymbol{\pi T}^{n}\boldsymbol{e}, \quad n=1,2,\ldots. \]
As for the continuous case, it is possible to investigate the
probability distribution and cumulative distribution by using a similar
set of functions, namely dDPH and pDPH:
x <- seq(0,10,by=1)
pdf <- dDPH(x, dph)
{plot(x, pdf, xlab = "Time to absorption", ylab = "PDF and CDF", col = "blue",
type = 'l', lwd=2)
title('Probability function (PDF)')}
x <- seq(0,10,by=1)
cdf <- pDPH(x, dph)
{plot(x, cdf, xlab = "Time to absorption", ylab = "CDF", col = "orange",
type = 'l', lwd=2)
title('Probability function (CDF)')}
Similarly, we have a function for the quantile
(qDPH)
x <- seq(0.05, 0.95, 0.01)
{plot(x, qDPH(x, dph), col = "green",
ylab = "Time to absorption",
xlab = "Quantile",xlim=c(0,1),ylim=c(0,10),type="l")
title('Quantile function')}
Finally, here is a histogram of 200 samples from the phase-type
distribution drawn from `rDPH``:
hist(rDPH(500, dph), breaks = 10)
