SIR_prob: Transition probabilities of an SIR process

In msuchard/MultiBD: Multivariate Birth-Death Processes

Transition probabilities of an SIR process

Description

Computes the transition pobabilities of an SIR process using the bivariate birth process representation

Usage

SIR_prob(t, alpha, beta, S0, I0, nSI, nIR, direction = c("Forward",
  "Backward"), power = NULL, nblocks = 20, tol = 1e-12, computeMode = 0,
  nThreads = 4)

Arguments

t	time
alpha	removal rate
beta	infection rate
S0	initial susceptible population
I0	initial infectious population
nSI	number of infection events
nIR	number of removal events
direction	direction of the transition probabilities (either Forward or Backward)
power	the power of the general SIR model (see Note for more details)
nblocks	number of blocks
tol	tolerance
computeMode	computation mode
nThreads	number of threads

Value

a matrix of the transition probabilities

Note

The infection rate and the removal rate of a general SIR model are \beta*S^{powS}*I^{powI_inf} and \alpha*I^{powI_rem} respectively. The parameter power is a list of powS, powI_inf, powI_rem. Their default values are powS = powI_inf = pwoI_rem = 1, which correspond to the classic SIR model.

Examples

data(Eyam)

loglik_sir <- function(param, data) {
  alpha <- exp(param[1]) # Rates must be non-negative
  beta  <- exp(param[2])

  if(length(unique(rowSums(data[, c("S", "I", "R")]))) > 1) {
    stop ("Please make sure the data conform with a closed population")
  }

  sum(sapply(1:(nrow(data) - 1), # Sum across all time steps k
             function(k) {
               log(
                 SIR_prob(  # Compute the forward transition probability matrix
                   t  = data$time[k + 1] - data$time[k], # Time increment
                   alpha = alpha, beta = beta,
                   S0 = data$S[k], I0 = data$I[k],       # From: R(t_k), I(t_k)
                   nSI = data$S[k] - data$S[k + 1], nIR = data$R[k + 1] - data$R[k],
                   computeMode = 4, nblocks = 80         # Compute using 4 threads
                 )[data$S[k] - data$S[k + 1] + 1,
                   data$R[k + 1] - data$R[k] + 1]        # To: R(t_(k+1)), I(t_(k+1))
               )
             }))
}

loglik_sir(log(c(3.204, 0.019)), Eyam) # Evaluate at mode

msuchard/MultiBD documentation built on May 19, 2024, 9:30 p.m.