Transition probabilities of an SIR process

Share:

Description

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

Usage

1
2
SIR_prob(t, alpha, beta, S0, I0, nSI, nIR, direction = c("Forward",
  "Backward"), 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)

nblocks

number of blocks

tol

tolerance

computeMode

computation mode

nThreads

number of threads

Value

a matrix of the transition probabilities

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
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