Skip to content

Commit afad90e

Browse files
authored
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.
1 parent 025b19e commit afad90e

File tree

1 file changed

+9
-3
lines changed

1 file changed

+9
-3
lines changed

R/ctmcProbabilistic.R

Lines changed: 9 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -61,12 +61,17 @@ rctmc <- function(n, ctmc, initDist = numeric(), T = 0, include.T0 = TRUE,
6161
i <- 1
6262
while (i <= n){
6363
idx <- which(ctmc@states == state)
64-
t <- t + rexp(1, -ctmc@generator[idx, idx])
65-
state <- ctmc@states[sample(1:dim(ctmc), 1, prob = trans[idx, ])]
64+
if (ctmc@generator[idx, idx] == 0) {
65+
# absorbing state; stay here forever
66+
t <- Inf
67+
} else {
68+
t <- t + rexp(1, -ctmc@generator[idx, idx])
69+
}
6670

67-
if(T > 0 & t > T)
71+
if((T > 0 & t > T) | (is.infinite(t)))
6872
break
6973

74+
state <- ctmc@states[sample(1:dim(ctmc), 1, prob = trans[idx, ])]
7075
states <- c(states, state)
7176
time <- c(time, t)
7277
i <- i + 1
@@ -84,6 +89,7 @@ rctmc <- function(n, ctmc, initDist = numeric(), T = 0, include.T0 = TRUE,
8489
stop("Not a valid output type")
8590
}
8691

92+
8793
#' @title Return the generator matrix for a corresponding transition matrix
8894
#'
8995
#' @description Calculate the generator matrix for a

0 commit comments

Comments
 (0)