# R - Fitting a Coxian phase-type distribution to data

Using R, I would like to fit a Coxian phase-type distribution to a vector of waiting times ($$t_i$$) where the $$\mu$$'s and $$\lambda$$'s are unknown i.e. I want to find the most optimal values of the $$\mu$$'s and $$\lambda$$'s. For reference, I want to minimize:

$$$$L=\sum_{i=1}^n \log(\textbf{p}\,\exp\{\textbf{Q}t_i\}\,\textbf{q}),$$$$

where: $$$$\textbf{p}=(1 0 0 \dots 0 0),$$$$

$$$$\textbf{q}=(\mu_1 \mu_2 \dots \mu_n)^T,$$$$

and:

$$$$\mathbf{Q}= \begin{bmatrix} -(\lambda_1 + \mu_1) & \lambda_1 & 0 & \dots & 0 & 0 \\ 0 & -(\lambda_2 + \mu_2) & \lambda_2 & \dots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \dots & -(\lambda_{n-1} + \mu_{n-1}) & \lambda_{n-1} \\ 0 & 0 & 0 & \dots & 0 & -\mu_n \\ \end{bmatrix}.$$$$

I've attempted using several optimization functions - optim() seems to be the most appropriate - but I keep running into errors. So I tried breaking it down into smaller problems. For instance, when I know the exact values of $$\mu_1$$, $$\mu_2$$ and $$\lambda_1$$ of a two-phase Coxian distribution, I can easily obtain the value of $$L$$. I used the following code to do this:

library(neldermead)
library(expm)

t <- ed$Total.Time x1 <- 0.1 # $$x_1$$ is $$\mu_1$$ x2 <- 0.1 # $$x_2$$ is $$\mu_2$$ x3 <- 0.1 # $$x_3$$ is $$\lambda_1$$ p <- matrix(c(1,0), nrow=1, ncol=2, byrow=TRUE) q <- matrix(c(x1,x2), nrow=2, ncol=1, byrow=TRUE) likely <- vector() for (i in 1:nrow(ed)) { likely[i] <- p%*%expm(matrix(c(-(x3+x1)*t[i],x3*t[i],0,-x2*t[i]), nrow=2, ncol=2, byrow=TRUE))%*%q } log.likely <- vector() for (i in 1:nrow(ed)) { log.likely[i] <- log(likely[i]) } L = -(sum(log.likely))  So that was fine. However, the errors occur when I try to run this code with unknown values of $$x$$, as shown below: p2 <- matrix(c(1,0), nrow=1, ncol=2, byrow=TRUE) q2 <- function(x) {matrix(c(x[1],x[2]), nrow=2, ncol=1, byrow=TRUE)} likely2 <- vector() for (i in 1:nrow(ed)) { likely2[i] <- function(x) {p2%*%expm(matrix(c(-(x[3]+x[1])*t[i],x[3]*t[i],0,-x[2]*t[i]), nrow=2, ncol=2, byrow=TRUE))%*%q2} }  The error message I get is: "Error in likely2[i] <- function(x) { : incompatible types (from closure to logical) in subassignment type fix"  Had I not received an error message, I would have continued with the following code: log.likely2 <- vector() for (i in 1:nrow(ed)) { log.likely2[i] <- log(likely2[i]) } x0 <- c(0.1,0.1,0.1) L <- function(x){-(sum(log.likely2))} optim(x0, L, method="Nelder-Mead")  If anyone has any suggestions, that would be greatly appreciated! Also, I don't think the Nelder-Mead method within optim() can take lower and upper bounds of the parameters. This is problematic because obviously the $$\mu$$'s and $$\lambda$$'s have to be between $$0$$ and $$1$$. • What is ed$Total.time and more generally what is ed? – JimB Feb 22 at 3:42
• @JimB ed is a data frame containing information on all patients who have been through a hospital's Emergency Department over the last few years. ed\$Total.time is the total amount of time a patient spent in ED from arrival to departure. – Jordan Feb 22 at 10:31
• I was too subtle. If you gave an explicit data frame of ed (even fake values in ed), there would be a minimal working example and we wouldn't need to make up some input data and you'd have more folks willing and able to help. But upon further inspection it appears that you are trying to assign a function definition to an element in a vector (which you can't do in R): likely2[i] <- function(x) {p2%*%expm(matrix(c(-(x[3]+x[1])*t[i],x[3]*t[i],0,-x[2]*t[i]), nrow=2, ncol=2, byrow=TRUE))%*%q2}. This is a programming issue rather than a statistical issue and is best asked elsewhere. – JimB Feb 22 at 16:33