On the estimation of partially observed continuous-time Markov chains

Abstract

Motivated by the increasing use of discrete-state Markov processes across applied disciplines, a Metropolis–Hastings sampling algorithm is proposed for a partially observed process. Current approaches, both classical and Bayesian, have relied on imputing the missing parts of the process and working with a complete likelihood. However, from a Bayesian perspective, the use of latent variables is not necessary and exploiting the observed likelihood function, combined with a suitable Markov chain Monte Carlo method, results in an accurate and efficient approach. A comprehensive comparison with simulated and real data sets demonstrate our approach when compared with alternatives available in the literature.

Appendices

Appendix A: Metropolis-Hastings algorithm (5) details

The Metropolis–Hastings moves used in our work are the following:

• For $\lambda \in \Lambda = \left\{ \lambda _1,\ldots ,\lambda _m\right\}$ let $g\,:\, {\mathbb {R}}\rightarrow {\mathbb {R}}^+$ be given by $g(x)=e^x$ and propose moves from $\lambda$ to ${\tilde{\lambda }}$ by

  $$\begin{aligned} {\tilde{\lambda }}= g( g^{-1}(\lambda ) +Z ) \end{aligned}$$

  with $Z\sim \text {Norma}(0,1)$ so ${\tilde{\lambda }}$ has a LogNormal distribution. The corresponding transition kernel is given by $q({\tilde{\lambda }}\,|\,\lambda )\propto e^{-0.5(\log (\lambda )-\log ({\tilde{\lambda }}))^2}{\tilde{\lambda }}^{-1}$.

• For transition matrix P with no zeros in off-diagonal entries and

  $\pmb {p}\in \left\{ \left( p_{j,1},\ldots ,p_{j,j-1},p_{j,j+1},\ldots ,p_{j,m} \right) \, : \, 1 \le j \le m \right\}$, let $c>0$ and propose

  $$\begin{aligned} \tilde{\pmb {p}}|\pmb {p}\sim \text {Dirichlet}(\pmb {1}+c\pmb {p}) \end{aligned}$$

  with $\pmb {1}=(1,\ldots ,1)$ so $\tilde{\pmb {p}}$ has mode $\pmb {p}$.

• For transition matrix P with zeros in off-diagonal entries and

  $p\in \left\{ p_{j,k} \, : \, 1 \le j\le m,\,1 \le k\le m,\, j\ne k\right\}$, let $g\,:\, {\mathbb {R}}\rightarrow {\mathbb {R}}^+$ be a standard logistic function $g(x)=e^x/(1+e^x)$ and propose moves from p to ${\tilde{p}}$ by

  $$\begin{aligned} {\tilde{p}}= g( g^{-1}(p) +Z ) \end{aligned}$$

  with $Z\sim \text {Normal}(0,1)$ so ${\tilde{p}}$ has a LogitNormal distribution. Denote $\pmb {w}$ as the vector $\pmb {p}$ with entry value p changed to ${\tilde{p}}$. Finally we do a normalization

  $$\begin{aligned} \tilde{\pmb {p}} = \pmb {w}/\sum _{k=1}^m w. \end{aligned}$$

  The transformation $f(x_1,\ldots ,x_m) = \left( \frac{x_1}{\sum _{i=1}^m x_i} , \ldots ,\frac{x_{m-1}}{\sum _{i=1}^m x_i},\sum _{i=1}^m x_i \right)$ is known to have Jacobian $\left( \sum _{i=1}^m x_i \right) ^{-m+1}$. So the corresponding transition kernel is given by $q(\tilde{\pmb {p}}\,|\,\pmb {p})\propto \exp \left( -0.5(\text {logit}(p_k)-\text {logit}({\tilde{p}}))^2\right) \frac{ \left( {\tilde{p}} + \sum _{i\ne k}p_i \right) ^{m-1}}{{\tilde{p}}(1-{\tilde{p}})}$. Note that we do not marginalize the random variable $S=\sum _{i=1}^m w_i$ so we have to draw simulations of the auxiliary variable ${\tilde{p}}$ to evaluate the above conditional density.

Appendix B: Second simulation study details

For a generator matrix

$$\begin{aligned} G = \begin{pmatrix} -\lambda _1 & \lambda _1 \\ \lambda _2 & - \lambda _2 \end{pmatrix} \end{aligned}$$

with $\lambda _1,\lambda _2>0$, we have explicitly that

$$\begin{aligned} Q(t) = \exp (tG) = \begin{pmatrix} \frac{\lambda _2}{\lambda _1 + \lambda _2} + \frac{\lambda _1}{\lambda _1 + \lambda _2}e^{-(\lambda _1 + \lambda _2) t} & \frac{\lambda _1}{\lambda _1 + \lambda _2} - \frac{\lambda _1}{\lambda _1 + \lambda _2}e^{-(\lambda _1 + \lambda _2) t} \\ & & \\ \frac{\lambda _2}{\lambda _1 + \lambda _2} - \frac{\lambda _2}{\lambda _1 + \lambda _2}e^{-(\lambda _1 + \lambda _2) t} & \frac{\lambda _1}{\lambda _1 + \lambda _2} + \frac{\lambda _2}{\lambda _1 + \lambda _2}e^{-(\lambda _1 + \lambda _2) t} \end{pmatrix}. \end{aligned}$$

For discrete observation over a time mesh $\Delta k$, $\Delta \in {\mathbb {R}}^+$, and $k \in \{1,\ldots , n\}$ for some $n\in {\mathbb {N}}$ of the Markov chain associated G; let $n^{eq}_i$ be the number of times the state remains equal for transitions starting in state i and $n^{ch}_i$ the number of times the state changes for transitions starting in state i along the time mesh, with $i\in \{1,2\}$. The likelihood is readily seen to be

$$\begin{aligned} {\mathcal {L}}(\lambda )&= \left( \frac{\lambda _2}{\lambda _1 + \lambda _2} + \frac{\lambda _1}{\lambda _1 + \lambda _2}e^{-(\lambda _1 + \lambda _2) \Delta } \right) ^{n^{eq}_1} \left( \frac{\lambda _1}{\lambda _1 + \lambda _2} + \frac{\lambda _2}{\lambda _1 + \lambda _2}e^{-(\lambda _1 + \lambda _2) \Delta } \right) ^{n^{eq}_2} \\&\times \left( \frac{\lambda _1}{\lambda _1 + \lambda _2} - \frac{\lambda _1}{\lambda _1 + \lambda _2}e^{-(\lambda _1 + \lambda _2) \Delta } \right) ^{n^{ch}_1} \left( \frac{\lambda _2}{\lambda _1 + \lambda _2} - \frac{\lambda _2}{\lambda _1 + \lambda _2}e^{-(\lambda _1 + \lambda _2) \Delta } \right) ^{n^{ch}_2}. \end{aligned}$$

Simulations for the above CTMC with $(\lambda _1,\lambda _2)=(1,1)$ and $(\lambda _1,\lambda _2=(2,1)$ were drawn. The prior distributions were $\lambda \sim \text {Gamma}(1,1)$ and $\left\{ p_{j,k} \right\} _{j\ne k}\sim \text {Dirichlet}(1)$ for all the experiments. A CTMC was generated until 1000 transitions were obtained and we considered the discrete observations given by times $\Delta k$ with $\Delta =1$ and $1\le k \le \min \left\{ n \, ;\ X(n\Delta )\text { was fully observed}\right\}$.

In Figs. 6 and 7 we compare the true posterior distributions given the simulated data with the Gibbs sampler 1 realizations of $\lambda _1$ and $\lambda _2$; whereas in Figures 8 and 9 we do the same experiment for the Gibbs sampler 2 and in 10, 11 for the Metropolis–Hastings sampler. The choice of priors was taken as indicated above for all the three samplers. The variance in the LogNormal proposals for $\lambda _i$ in the Metropolis–Hastings were set to 0.0025. We ran 50000 iterations of each sampler after having diagnosed stationarity via the potential scale reduction factor of Gelman and Rubin (1992). We highlight that the run time for the Metropolis–Hastings algorithm is significantly faster than the Gibbs samplers as showed in Tables 1 and 2. This is due to the computational cost of simulating the conditioned trajectories between discretely observed times for the CTMC.

Fig. 6

Marginal and bivariate histograms of $\pmb {\lambda }$ for the Gibbs sampler 1 ran for 50000 iterations after diagnosing stationarity compared with the the true posterior distribution. corresponding to $\pmb {\lambda }=(1,1)$

Fig. 7

Marginal and bivariate histograms of $\pmb {\lambda }$ for the Gibbs sampler 1 ran for 50000 iterations after diagnosing stationarity compared with the the true posterior distribution. corresponding to $\pmb {\lambda }=(2,1)$

Fig. 8

Marginal and bivariate histograms of $\pmb {\lambda }$ for the Gibbs sampler 2 ran for 50000 iterations after diagnosing stationarity compared with the the true posterior distribution. corresponding to $\pmb {\lambda }=(1,1)$

Fig. 9

Marginal and bivariate histograms of $\pmb {\lambda }$ for the Gibbs sampler 2 ran for 50000 iterations after diagnosing stationarity compared with the the true posterior distribution. corresponding to $\pmb {\lambda }=(2,1)$

Fig. 10

Marginal and bivariate histograms of $\pmb {\lambda }$ for the Metropolis–Hastings sampler ran for 50000 iterations after diagnosing stationarity compared with the the true posterior distribution. corresponding to $\pmb {\lambda }=(1,1)$

Fig. 11

Marginal and bivariate histograms of $\pmb {\lambda }$ for the Metropolis–Hastings sampler ran for 50,000 iterations after diagnosing stationarity compared with the the true posterior distribution. corresponding to $\pmb {\lambda }=(2,1)$

Appendix C: Further comparisons for the third simulation study

In this appendix we present further mean generator comparisons between the Gibbs and Metropolis–Hastings samplers for discretely observed CTMCs. Also we present comparisons of generators fitted with the EM algorithm and our Metropolis–Hastings approach.

1.1 C.1: Gibbs and Metropolis–Hastings comparisons

In Figs. 13, 15 and we present the mean generator comparisons for the Gibbs and Metropolis–Hastings samplers with simulation studies determined, respectively, by $\Delta = 0.5,2$ (Figs. 12, 13, 14, 15, 16, 17, 18, 19, 20).

As mentioned in the main document the entry-wise distance of the Gibbs mean generator with the true generator tends to be greater than the one with respect to the Metropolis–Hasting mean generator; this is particularly illustrated in the mean generators for $\Delta =0.5$ as shown in Table 1 of the main document in terms of Frobenius distances.

Fig. 12

Comparison of estimated generator matrices for Metropolis–Hastings and Gibbs chains for $\Delta = 4$ with $\lambda = 0.25$ (first row), $\lambda =0.5$ (second row) and $\lambda = 1$ (third row)

Fig. 13

Comparison of estimated generator matrices for Metropolis–Hastings and Gibbs chains for $\Delta = 2$ with $\lambda = 2$ (first row), $\lambda =1$ (second row) and $\lambda = 0.5$ (third row)

Fig. 14

Comparison of estimated generator matrices for Metropolis–Hastings and Gibbs chains for $\Delta = 0.5$ with $\lambda = 8$ (first row), $\lambda =4$ (second row) and $\lambda = 2$ (third row)

Fig. 15

Comparison of estimated generator matrices for Metropolis–Hastings and Gibbs chains for $\Delta = 0.25$ with $\lambda = 4$ (first row), $\lambda =8$ (second row) and $\lambda = 16$ (third row)

1.2 C.2: EM and Metropolis–Hastings comparison

We use the EM algorithm of the ctmcd R package, see Pfeuffer (2017), to compare with the proposed Metropolis–Hastings algorithm in the context of the third simulation study.

Fig. 16

Comparison of estimated generator matrices for Metropolis–Hastings and EM algorithms for $\Delta = 4$ with $\lambda = 1$ (first row), $\lambda =0.5$ (second row) and $\lambda = 0.25$ (third row)

Fig. 17

Comparison of estimated generator matrices for Metropolis–Hastings and EM algorithms for $\Delta = 2$ with $\lambda = 2$ (first row), $\lambda =1$ (second row) and $\lambda = 0.5$ (third row)

Fig. 18

Comparison of estimated generator matrices for Metropolis–Hastings and EM algorithms for $\Delta = 1$ with $\lambda = 4$ (first row), $\lambda =2$ (second row) and $\lambda = 1$ (third row)

Fig. 19

Comparison of estimated generator matrices for Metropolis–Hastings and EM algorithms for $\Delta = 0.5$ with $\lambda = 8$ (first row), $\lambda =4$ (second row) and $\lambda = 2$ (third row)

Fig. 20

Comparison of estimated generator matrices for Metropolis–Hastings and EM algorithms for $\Delta = 0.25$ with $\lambda = 16$ (first row), $\lambda =8$ (second row) and $\lambda = 4$ (third row)