SIR_prob: Transition probabilities of an SIR process

View source: R/SIR_prob.R

SIR_probR Documentation

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.