MultiBD (version 1.0.0)

SIR_prob: 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
)

Value

a matrix of the transition probabilities

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

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