Statistical moments of phase-type distributions

Computing (factorial) moments of a given order of phase-type distributed random variables.

moments(object, i, all = FALSE)

Arguments

object	either a continuous phase-type distributed object of class contphasetype or a discrete phase-type distributed object of class discphasetype.
i	a positive number stating the order of the desired moment
all	a logical value indicating whether the function should compute all moments up to the given order. The default is equal to FALSE.

Source

Mogens Bladt and Bo Friis Nielsen (2017): Matrix-Exponential Distributions in Applied Probability. Probability Theory and Stochastic Modelling (Springer), Volume 81.

Value

For all = FALSE, the function either returns the \(i\)'th-order moment (if the object is continuous phase-type distributed) or the \(i\)'th factorial moment (if the object is discrete phase-type distributed). In both cases, the length of the output is one. For all = TRUE, the function computes all (factorial) moments up to the given order, hence the output is a vector of length \(i\).

Details

In the discrete case \(( \tau ~ DPH(initDist,P) )\), the factorial moments are given by $$E[\tau(\tau-1) ··· (\tau-i+1)] = i! initDist P^(i-1) (I-P)^(-i) e,$$ where initDist is the initial distribution and \(P\) is the sub-transition probability matrix. For \(\tau ~ PH(initDist, T)\), the \(i\)'th-order moment is defined as $$E[\tau^i] = i! initDist (-T)^(-i) e,$$ where initDist is again the initial distribution and \(T\) is the sub-intensity rate matrix. In both cases, \(e\) is a vector with one in each entry.

Examples

## Using the function moments() to compute the mean
## and variance of a phase-type distribution

## For n=4, the time to the most recent common ancestor is
## phase-type distributed with initial distribution
initDist <- c(1,0,0)
## and sub-intensity rate matrix
T_Mat <- matrix(c(-6,6,0,
                   0,-3,3,
                   0,0,-1), nrow = 3, ncol = 3, byrow = TRUE)
## Defining an object of type "contphasetype"
TMRCA <- contphasetype(initDist, T_Mat)
## Computing all moments up to order 2
m <- moments(TMRCA, i=2, all = TRUE)
## We get the desired numbers
m[1]
#>   1 
#> 1.5 
phmean(TMRCA)
#> [1] 1.5

m[2] - m[1]^2
#>        2 
#> 1.138889 
phvar(TMRCA)
#> [1] 1.138889

## For theta=2, the number of segregating sites plus one is
## discrete phase-type distributed with initial distribution
initDist <- c(1,0,0,0)
## and sub-transition probability matrix
P_Mat <- matrix(c(0.4, 0.3, 4/30, 2/30,
                   0, 0.5, 2/9, 1/9,
                   0, 0, 2/3, 0,
                   0, 0, 0, 2/3), nrow = 4, ncol = 4, byrow = TRUE)
## Defining an object of type "discphasetype"
S_Total <- discphasetype(initDist, P_Mat)
## Computing all moments up to order 2
m <- moments(S_Total, i=2, all = TRUE)
## We get the desired numbers
m[1]
#>        1 
#> 4.666667 
phmean(S_Total)
#> [1] 4.666667

m[2] + m[1] - m[1]^2
#>        2 
#> 9.111111 
phvar(S_Total)
#> [1] 9.111111