Multivariate birth-death processes

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

1 2 |

`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 |

`nblocks` |
number of blocks |

`tol` |
tolerance |

`computeMode` |
computation mode |

`nThreads` |
number of threads |

a matrix of the transition probabilities

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
``` |

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.

Please suggest features or report bugs with the GitHub issue tracker.

All documentation is copyright its authors; we didn't write any of that.