moments: Statistical moments of phase-type distributions

In aumath-advancedr2019/PhaseTypeGenetics: The Phase-Type Distribution with Applications in Population Genetics

Description

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

Usage

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.

Details

In the discrete case ( τ ~ DPH(initDist,P) ), the factorial moments are given by

E[τ(τ-1) ··· (τ-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 τ ~ PH(initDist, T), the i'th-order moment is defined as

E[τ^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.

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.

Source

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

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]
phmean(TMRCA)

m[2] - m[1]^2
phvar(TMRCA)

## 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]
phmean(S_Total)

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

aumath-advancedr2019/PhaseTypeGenetics documentation built on Dec. 3, 2019, 7:16 a.m.