From afad90e6b88ad542973c3543c24aade872596c5c Mon Sep 17 00:00:00 2001 From: ntguardian <cgmil@msn.com> Date: Wed, 21 Aug 2024 16:18:04 -0400 Subject: [PATCH] Corrected bug for continuous time Markov chain simulation with absorbing states rtcmc originally could not handle continuous time Markov chains with absorbing states because rexp cannot work with rate zero parameters, which would be input for continuous time Markov chains. Added fix that checks for transition rate of 0 and makes transition time Inf if so. Chain will terminate early if it reaches an absorbing state. --- R/ctmcProbabilistic.R | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/R/ctmcProbabilistic.R b/R/ctmcProbabilistic.R index 526aedd..cf8ddfe 100644 --- a/R/ctmcProbabilistic.R +++ b/R/ctmcProbabilistic.R @@ -61,12 +61,17 @@ rctmc <- function(n, ctmc, initDist = numeric(), T = 0, include.T0 = TRUE, i <- 1 while (i <= n){ idx <- which(ctmc@states == state) - t <- t + rexp(1, -ctmc@generator[idx, idx]) - state <- ctmc@states[sample(1:dim(ctmc), 1, prob = trans[idx, ])] + if (ctmc@generator[idx, idx] == 0) { + # absorbing state; stay here forever + t <- Inf + } else { + t <- t + rexp(1, -ctmc@generator[idx, idx]) + } - if(T > 0 & t > T) + if((T > 0 & t > T) | (is.infinite(t))) break + state <- ctmc@states[sample(1:dim(ctmc), 1, prob = trans[idx, ])] states <- c(states, state) time <- c(time, t) i <- i + 1 @@ -84,6 +89,7 @@ rctmc <- function(n, ctmc, initDist = numeric(), T = 0, include.T0 = TRUE, stop("Not a valid output type") } + #' @title Return the generator matrix for a corresponding transition matrix #' #' @description Calculate the generator matrix for a