# EMpht.jl

A Julia port of the EMpht.c program, used for fitting Phase-Type distributions via an EM algorithm.

The original C — which is available on Søren Asmussen's website — is well documented and has a decent performance for phase-type distributions with a small or moderate number of phases. However it is quite slow for when the number of phases is large (>= 20), and the UX is very old-school Unix.

This port is much simpler and faster. See the examples given below and the test cases.

## Examples

To fit a phase-type distribution to some data:

```
using Distributions
using EMpht
data = rand(Exponential(1/10), 1_000) # Generate some data to fit
sample = EMpht.Sample(obs=data) # Create an EMpht Sample object with this data
ph = empht(sample, p=5) # Fit the data using p=5 phases
xGrid = range(0, 8, length=1_00) # Create a grid to evaluate the density function over
fitPDFs = pdf.(ph, xGrid) # The probability density function of the fitted phase-type
```

The default structure of the phase-type is "Coxian" (see below for details).
For large values of `p`

the "CanonicalForm1" is recommended.
To impose no structure on the phase-type, use "General", though the results degrade quickly with `p > 5`

or so.
Another available option is "GeneralisedCoxian".

```
phGen = empht(sample, p=20, ph_structure="General")
phCox = empht(sample, p=20, ph_structure="Coxian")
phCF1 = empht(sample, p=50, ph_structure="CanonicalForm1")
```

If the data is not fully observed, i.e. the data is binned (interval-censored), then the Sample object is updated like so:

```
# The intervals
int = [1.5 2.0; 2.0 2.5; 2.5 3.0; 3.0 3.5; 3.5 4.0; 4.0 4.5; 4.5 5.0; 5.0 5.5;
5.5 6.0; 6.0 6.5; 6.5 7.0; 7.0 7.5]
# The number of observations falling into each interval
intweight = [4.0, 34.0, 107.0, 170.0, 202.0, 222.0, 140.0, 77.0, 24.0, 14.0,
4.0, 2.0]
# Create an EMpht Sample object with this data
sInt = EMpht.Sample(int=int, intweight=intweight)
# Fitting the interval-censored data
phCF1 = empht(sInt, p=100, ph_structure="CanonicalForm1")
xGrid = range(0, 8, length=1_000)
fitPDFs = pdf.(phCF1, xGrid)
```

To choose the algorithm used to fit the data (see papers below for details):

```
phunif = empht(sample, p=5, method=:unif) # Fit using the uniformization technique (default)
phode = empht(sample, p=5, method=:ode) # Fit using the more traditional ODE solving technique
```

EMpht.jl can read all necessary information from a JSON file (the number of phases to fit, the special structure of the phase-type, the sample to fit). For example, if you download the Coxian100.json file inside the test directory, the following will launch a fit based on those parameters:

`ph100 = empht("Coxian100.json")`

## Resources

The relevant papers for the algorithms are:

- S. Asmussen, O. Nerman & M. Olsson,
*Fitting phase-type distribution via the EM algorithm*, Scandinavian Journal of Statistics 23, 419-441 (1996), - M. Olsson,
*Estimation of phase-type distributions from censored data*, Scandinavian Journal of Statistics 23, 443-460 (1996). - H. Okamura, T. Dohi, K.S. Trivedi,
*A refined EM algorithm for PH distributions*, Performance Evaluation 68, 938-954 (2011) - H. Okamura, T. Dohi, K.S. Trivedi,
*Improvement of expectation-maximization algorithm for phase-type distributions with grouped and truncated data*, Appl. Stochastic Models Bus. Ind. 29, 141-156 (2013)

Some case studies using this package are:

- S. Asmussen, P.J. Laub, H. Yang,
*Phase-type models in life insurance: fitting and valuation of equity-linked benefits*, Risks 7(1), 17 pages (2019) - A. Vuorinen,
*The blockchain propagation process: a machine learning and matrix analytic approach*, University of Melbourne Masters Thesis (2019), see website or thesis.