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:

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

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

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

and:

$$\begin{equation} \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}. \end{equation}$$

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:

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,x), nrow=2, ncol=1, byrow=TRUE)} likely2 <- vector() for (i in 1:nrow(ed)) { likely2[i] <- function(x) {p2%*%expm(matrix(c(-(x+x)*t[i],x*t[i],0,-x*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