Nothing
## ----setup, include=TRUE, echo=FALSE, cache=FALSE---------------------------------------
library(knitr)
# set global chunk options
opts_chunk$set(fig.path='figure/minimal-', fig.align='center', fig.show='hold', cache=FALSE)
options(formatR.arrow=TRUE,width=90,tidy=FALSE)
## ----ODE, echo=TRUE, warning=FALSE------------------------------------------------------
library(matrixStats)
# solutions to phi1 and phi2, the generating function ODES
# phi1 is generating function starting with 1 suscep
# phi2 is starting with 1 infected
phi1 <- function(t, s1, s2, beta, gam, I0){
return( 1 + exp(-gam*t)*beta*I0*(s2-1)/(beta*I0-gam) +
exp(-beta*I0*t)*(s1 - 1 - beta*I0*(s2-1)/(beta*I0-gam) ) )
}
phi2 <- function(t, s2, gam){
return(1 + (s2-1)*exp(-gam*t))
}
#solve over a grid: s1.seq and s2.seq will be vectors of complex numbers
#along unit circle
phi1.grid <- function(t, s1.seq, s2.seq, beta, gam, I0) {
#this will store the grid
grid <- matrix(nrow = length(s1.seq), ncol = length(s2.seq))
for(i in 1:length(s1.seq)){
for(j in 1:length(s2.seq)){
grid[i,j] <- phi1(t, s1.seq[i], s2.seq[j], beta, gam, I0)
}
}
return(grid)
}
phi2.grid <- function(t, s2.seq, gam){
#this will store the grid
grid <- matrix(nrow = length(s2.seq), ncol = length(s2.seq))
for(j in 1:length(s2.seq)){
grid[,j] <- phi2(t, s2.seq[j], gam) #columns will be constant
}
return(grid)
}
#exponentiate grid solutions appropriately and fast fourier transform
#to get trans probs
getTransProbsODE <- function(t, gridLength, beta, gam, S0, I0){
s1.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
s2.seq <- exp(2*pi*1i*seq(from = 0, to = (gridLength-1))/gridLength)
jointGrid <- phi1.grid(t, s1.seq, s2.seq, beta, gam, I0)^S0 *
phi2.grid(t, s2.seq, gam)^I0
fourierGrid <- fft(jointGrid)/length(jointGrid)
return(Re(fourierGrid))
}
## ----BranchingFunctions-----------------------------------------------------------------
# expression for the phi1 cross-partials
A <- function(t, m, n, k, j, beta, gam){
if(j > m-k){return(0)} else{
return( factorial(m) * exp(-k*beta*n*t) *
( 1 - beta*n*exp(-gam*t)/(beta*n-gam) - exp(-beta*n*t) *
(1 - beta*n/(beta*n-gam) ) )^(m-k-j) *
( beta*n*(exp(-gam*t) - exp(-beta*n*t)) /
(beta*n - gam) )^(j) / factorial(m-k-j) )
}
}
# expression for the phi2 partials
B <- function(t, n, j, gam){
if(j>n){return(0)} else{
return( factorial(n)* (1-exp(-gam*t))^(n-j) *
exp(-gam*j*t)/factorial(n-j) )
}
}
# remove if statements from A(), B(), so that vector arguments work
# input j must be greater than m-k
A.vectorized <- function(t, m, n, k, j, beta, gam){
return( factorial(m) * exp(-k*beta*n*t) *
( 1 - beta*n*exp(-gam*t)/(beta*n-gam) - exp(-beta*n*t) *
(1 - beta*n/(beta*n-gam) ) )^(m-k-j) *
( beta*n*(exp(-gam*t) - exp(-beta*n*t)) /
(beta*n - gam) )^(j) / factorial(m-k-j) )
}
#log version of previous, to handle larger numbers
A.vec.log <- function(t, m, n, k, j, beta, gam){
return( lgamma(m+1) - lgamma(m-k-j+1) - (k*beta*n*t) +
(m-k-j)* log( ( 1 - beta*n*exp(-gam*t)/(beta*n-gam) -
exp(-beta*n*t) *
(1 - beta*n/(beta*n-gam) ) ) )
+ j * log( ( beta*n*(exp(-gam*t) - exp(-beta*n*t)) /
(beta*n - gam) ) )
)
}
B.vectorized <- function(t, n, j, gam){
return( factorial(n)* (1-exp(-gam*t))^(n-j) * exp(-gam*j*t) /
factorial(n-j) )
}
B.vec.log <- function(t, n, j, gam){
return( lgamma(n+1) - lgamma(n-j+1) + (n-j)*log(1 - exp(-gam*t))
- gam*j*t )
}
# loops over appropriate indices in a Leibniz-rule sum
# for partial derivative solution.
# uses expressions A(), B(), and returns P_{(m,n),(k,l)}(t) given
# beta and gamma.
# slow version with for loop: included here since it's more transparent
# to read
TransProb_mnkl_old <- function(t, m, n, k, l, beta, gam){
AA <- BB <- rep(0,l+1)
c <- choose(l,seq(0,l))
for(i in 1:(l+1)){
AA[i] <- A(t, m, n, k, l-i+1, beta, gam)
BB[i] <- B(t, n, i-1, gam)
}
# print(AA)
# print(BB)
return( sum(AA*BB* c) / (factorial(k)*factorial(l) ) )
}
# note: of course k can never be greater than m
# uses log versions to handle large populations
# put k and l first for use with outer()
TransProb_mnkl <- function(k,l, t, m, n, beta, gam){
if(k>m){return(0)}
logAA <- logBB <- rep(0,l+1)
aj <- rev(seq(0,min(m-k,l)))
bj <- seq(0, min(n,l) )
c <- lgamma(l+1) - lgamma(seq(0,l) + 1) - lgamma( l+1 - seq(0,l))
logAA[( l+1-min(m-k,l) ):(l+1)] <- A.vec.log(t, m, n, k, aj, beta, gam)
logBB[ 1: (min(n,l)+1) ] <- B.vec.log(t, n, bj, gam)
# print(AA)
# print(BB)
# print(c)
term <- c + logAA + logBB
return( exp( logSumExp(term) - lgamma(k+1) - lgamma(l+1) ) )
}
#vectorize the previous for use with outer
vectorized <- Vectorize(TransProb_mnkl, vectorize.args = c('k','l'))
# S_0, I_0 are m,n: fixes those and returns a grid of transitions
# over different end states for comparison with getTransProbs()
#vectorized version
getTransProbsClosed <- function(t, gridLength, beta, gam, S0, I0){
return( outer(0:gridLength, 0:gridLength, vectorized,
t = t, beta = beta, gam = gam, m = S0, n = I0))
}
# simple simulation of true SIR model over a time interval
simulateSIR <- function(t.end, S, I, beta, gam, maxEvents = 99999999){
t.cur <- 0
for(i in 1:maxEvents){
if( S<0 || I <0 ){
print("Negative population? Error")
return(-99) #error code
}
if( S == 0 || I ==0){
#print("S or I is zero, end epidemic")
return( c(S,I) ) #end epidemic
}
infectRate <- S*I*beta
recovRate <- I*gam
rates <- c(infectRate, recovRate)
t.next <- rexp(1, sum(rates)) #time until next event
t.cur <- t.cur+t.next
if(t.cur > t.end){ #end of simulation period
return( c(S,I) )
}
#sample the type of next event
decision <- rbinom(1, 1, infectRate/sum(rates))
if(decision == 1) { #infection
S <- S - 1; I <- I + 1
} else { #recovery
I <- I - 1
}
}
return(-99) #error code for testing
}
#run simulation once with error catch
sim.once <- function(t.end, S, I, beta, gam, maxEvents = 99999999){
res = -99 # error catch
while(res[1] == -99){
res <- simulateSIR(t.end, S, I, beta, gam, maxEvents) }
return(res)
}
getTrans.MC <- function(N, t.end, S, I, beta, gam){
result <- replicate(N, sim.once(t.end, S, I, beta, gam))
#make big enough to account for all events: count end states
trans.count <- matrix(0, S+I,S+I)
for(i in 1:N){
id <- result[,i]+1
#indices in the resulting transition count: ie if you end at (1,1),
# you add a count to the (2,2) entry of the count matrix, etc
trans.count[id[1], id[2]] = trans.count[id[1], id[2]] + 1
}
tpm <- trans.count/sum(trans.count)
return(tpm)
}
## ---------------------------------------------------------------------------------------
#Choose initial S and I population here
S <- 140
I <- 10
beta = .5/(S+I)
gamma = .1
N <- 1000 #number of MC realizations: increase N in practice
t.end <- .5 #time interval length
## ----MonteCarlo, eval=FALSE-------------------------------------------------------------
# #monte carlo estimate and standard error
# tpm.MC <- getTrans.MC(N, t.end, S, I, beta, gamma)
# sd.MC <- sqrt( (tpm.MC)*(1 - tpm.MC)/N )
## ----doMonteCarlo, echo=FALSE, eval=TRUE, include=FALSE---------------------------------
# Provide manual caching because knitr's caching
# is not working in my hands
file <- system.file("vignetteCache", "monte_carlo.RData", package="MultiBD")
if (file.exists(file)) {
load(file)
} else {
#monte carlo estimate and standard error
tpm.MC <- getTrans.MC(N, t.end, S, I, beta, gamma)
sd.MC <- sqrt( (tpm.MC)*(1 - tpm.MC)/N )
save(tpm.MC, sd.MC, file = "../inst/vignetteCache/monte_carlo.RData")
}
## ----GenFunc, eval=FALSE----------------------------------------------------------------
# # now, calculate probabilities using generating functions
# gridLength = 256
# t1 <- system.time(
# tpm1 <- getTransProbsODE(t.end, gridLength,
# beta, gamma, S, I)[1:(S+I),1:(S+I)])
# t2 <- system.time(
# tpm2 <- getTransProbsClosed(t.end, gridLength,
# beta, gamma, S, I)[1:(S+I),1:(S+I)])
## ----doGenFunc, echo=FALSE, eval=TRUE, include=FALSE------------------------------------
# Provide manual caching because knitr's caching
# is not working in my hands
file <- system.file("vignetteCache", "gen_func.RData", package="MultiBD")
if (file.exists(file)) {
load(file)
} else {
# now, calculate probabilities using generating functions
gridLength = 256
t1 <- system.time(
tpm1 <- getTransProbsODE(t.end, gridLength,
beta, gamma, S, I)[1:(S+I),1:(S+I)])
t2 <- system.time(
tpm2 <- getTransProbsClosed(t.end, gridLength,
beta, gamma, S, I)[1:(S+I),1:(S+I)])
save(tpm1, tpm2, t1, t2,
file = "../inst/vignetteCache/gen_func.RData",
compression_level = 9)
}
## ---------------------------------------------------------------------------------------
t1
t2
## ---------------------------------------------------------------------------------------
# total errors:
sum(abs(tpm1 - tpm2)) # compare the two methods
sum(abs(tpm1 - tpm.MC)) #compare to true model Monte Carlo probabilities
## ----marginalizedSIRplots---------------------------------------------------------------
#marginalize over susceptibles: c
infectiveProbs.MC <- colSums(tpm.MC)
infectiveSD.MC <- sqrt( (infectiveProbs.MC)*(1 - infectiveProbs.MC)/(N) )
#check that this SD is correct...
lower <- infectiveProbs.MC - infectiveSD.MC*1.96
upper <- infectiveProbs.MC + infectiveSD.MC*1.96
infectiveProbs.FFT <- round(colSums(tpm2),4)
require(plotrix)
plot(seq(S + I), xlim = c(0,50), infectiveProbs.FFT, pch = 3,
col = 'purple', main = "Probabilities of ending with x infectives")
plotCI(seq(S + I), xlim = c(0,50), infectiveProbs.MC, pch = 16,
col = 4, ui = upper, li = lower, add=TRUE)
#marginalize other way:
suscepProbs.MC <- rowSums(tpm.MC)
suscepSD.MC <- sqrt( (suscepProbs.MC)*(1 - suscepProbs.MC)/(N) )
lower <- suscepProbs.MC - suscepSD.MC*1.96
upper <- suscepProbs.MC + suscepSD.MC*1.96
suscepProbs.FFT <- round(rowSums(tpm2),4)
plot(seq(S+I), xlim = c(110,160), suscepProbs.FFT, pch = 3,
col = 'purple', main = "Probabilities of ending with x susceptibles")
plotCI(seq(S+I), xlim = c(110,160), suscepProbs.MC, pch = 16,
col = 4, ui = upper, li = lower, add=TRUE)
## ----multiBDsetup-----------------------------------------------------------------------
library(MultiBD)
tList <- c(.1, .2, .25, .3 ,.35, .4, .5, .6, .7, .8, .9, 1)
gridLength = 128
a0 = 110 # S_0
b0 = 15 # I_0
A = 0
B = gridLength - 1
alpha = 3.2 #3.2 #this is death rate
beta = .025 #.019 #this is transition or infection rates
nSim = 4000 #number of MC simulations
brates1=function(a,b){0}
drates1=function(a,b){0}
brates2=function(a,b){0}
drates2=function(a,b){alpha*b}
trans=function(a,b){beta*a*b}
## ----MakeTPMs, cache=FALSE, eval=FALSE--------------------------------------------------
# #indexed by time, type of computation, and dimensions of the tpm
# tpmArray <- array(NA, dim = c(length(tList), 3, 52, 25 )
# ) #store a subset of transition probabilities
# time1 <- rep(0, length(tList)); time2 <- rep(0, length(tList))
#
# for(i in 1:length(tList)){
# t.end <- tList[i]
# time1[i] <- system.time(
# tpm.Closed <- getTransProbsClosed(t.end, gridLength,
# beta, alpha, a0, b0) )
# tpm1 = tpm.Closed[1:(a0+1),] #using 2-type branching approximation
#
# #using continued fractions via MultiBD
# time2[i] <- system.time(
# tpm2 <- dbd_prob(t.end, a0, b0, drates1, brates2, drates2, trans,
# a=A, B))#, computeMode=2))
# #MC simulation "ground truth"
# tpm.MC <- getTrans.MC(nSim, t.end, a0, b0, beta, alpha)
# tpm3 <- tpm.MC[1:(a0+1), ]
#
# #store subset of matrices containing about 99 percent of the mass:
# tpmArray[i,1,,] <- tpm1[60:(a0+1),1:25]
# tpmArray[i,2,,] <- tpm2[60:(a0+1),1:25]
# tpmArray[i,3,,] <- tpm3[60:(a0+1),1:25]
# }
## ----doMakeTPMs, echo=FALSE, eval=TRUE, include=FALSE-----------------------------------
# Provide manual caching because knitr's caching
# is not working in my hands
file <- system.file("vignetteCache", "tpm_array.RData", package="MultiBD")
if (file.exists(file)) {
load(file)
} else {
#indexed by time, type of computation, and dimensions of the tpm
tpmArray <- array(NA, dim = c(length(tList), 3, 52, 25 )
) #store a subset of transition probabilities
time1 <- rep(0, length(tList)); time2 <- rep(0, length(tList))
for(i in 1:length(tList)){
t.end <- tList[i]
time1[i] <- system.time(
tpm.Closed <- getTransProbsClosed(t.end, gridLength,
beta, alpha, a0, b0) )
tpm1 = tpm.Closed[1:(a0+1),] #using 2-type branching approximation
#using continued fractions via MultiBD
time2[i] <- system.time(
tpm2 <- dbd_prob(t.end, a0, b0, drates1, brates2, drates2, trans,
a=A, B))#, computeMode=2))
#MC simulation "ground truth"
tpm.MC <- getTrans.MC(nSim, t.end, a0, b0, beta, alpha)
tpm3 <- tpm.MC[1:(a0+1), ]
#store subset of matrices containing about 99 percent of the mass:
tpmArray[i,1,,] <- tpm1[60:(a0+1),1:25]
tpmArray[i,2,,] <- tpm2[60:(a0+1),1:25]
tpmArray[i,3,,] <- tpm3[60:(a0+1),1:25]
}
save(tpmArray, time1, time2, file = "../inst/vignetteCache/tpm_array.RData")
}
## ---------------------------------------------------------------------------------------
#for example, look at the ones with t.end = .5
small1 <- tpmArray[5,1,,]
small2 <- tpmArray[5,2,,]
small3 <- tpmArray[5,3,,]
#they comprise most of transition probability mass:
sum(small1); sum(small2); sum(small3)
# mean errors
mean(abs(small1- small3 ) ) #2-type vs MC
mean(abs(small2 - small3) ) #Continued Frac vs MC
# scaled heatmap images to compare tpm visually
par(mfrow=(c(3,1)))
image(small1, main = "Two-type branching approximation")
image(small2, main = "Continued Fraction expansion")
image(small3, main = "Monte Carlo estimates")
## ---------------------------------------------------------------------------------------
par(mfrow=(c(3,1)))
image(tpmArray[12,1,,], main = "Two-type branching approximation")
image(tpmArray[12,2,,], main = "Continued Fraction expansion")
image(tpmArray[12,3,,], main = "Monte Carlo estimates")
## ----transitionCompare------------------------------------------------------------------
library(plotrix)
inds <- t(which(tpmArray[7,2,,] >= sort(tpmArray[7,2,,], decreasing=T)[16],
arr.ind=TRUE))
#ind1 <- sample(52,25, replace=T); ind2 <- sample(25,25,replace=T)
par(mfrow = c(4,4), oma = c(5,4,0,0) + 0.1, mar = c(0,0,1,1) + 0.1)
for(i in 1:16){
plot(tList, tpmArray[,2,inds[1,i], inds[2,i] ], pch = 17, col = 'red',
ylim = c(0,max(tpmArray[,,inds[1,i], inds[2,i]])),
yaxt = 'n', xlab = "dt")
MCp <- tpmArray[,3,inds[1,i], inds[2,i] ] #MC prob
plotCI(tList, MCp, pch = 4, col = 'green',
ui = MCp+1.96*sqrt(MCp*(1-MCp)/nSim),
li <- MCp-1.96*sqrt(MCp*(1-MCp)/nSim), add = TRUE)
points(tList, tpmArray[,1,inds[1,i], inds[2,i] ],
col='purple', pch = 16)
}
## ----LargerPopulation,eval=FALSE--------------------------------------------------------
# tList <- c( .5, 1)
# gridLength = 256
# a0 = 235 # S_0
# b0 = 15 # I_0
# A = 0
# B = gridLength - 1
# alpha = 3.2 #3.2 #this is death rate
# beta = .025 #.019 #this is transition or infection rates
# nSim <- 10000
# #store a subset of transition probabilities
# tpmArray <- array(NA, dim= c(length(tList),2, (a0+1), 240 ))
# time1 <- rep(0, length(tList))
#
# for(i in 1:length(tList)){
# t.end <- tList[i]
# time1 <- system.time(
# tpm2 <- dbd_prob(t.end, a0, b0, drates1, brates2, drates2, trans,
# a=A, B))#, computeMode=2))
# #MC simulation "ground truth"
# tpm.MC <- getTrans.MC(nSim, t.end, a0, b0, beta, alpha)
# tpm3 <- tpm.MC[1:(a0+1), ]
#
# #store subset of matrices containing about 99 percent of the mass:
# tpmArray[i,1,,] <- tpm2[1:(a0+1),1:60]
# tpmArray[i,2,,] <- tpm3[1:(a0+1),1:60]
# }
## ----doLargerPopulation, echo=FALSE, eval=TRUE, include=FALSE---------------------------
# Provide manual caching because knitr's caching
# is not working in my hands
file <- system.file("vignetteCache", "larger_pop.RData", package="MultiBD")
if (file.exists(file)) {
load(file)
} else {
tList <- c( .5, 1)
gridLength = 256
a0 = 235 # S_0
b0 = 15 # I_0
A = 0
B = gridLength - 1
alpha = 3.2 #3.2 #this is death rate
beta = .025 #.019 #this is transition or infection rates
nSim <- 10000
#store a subset of transition probabilities
tpmArray <- array(NA, dim= c(length(tList),2, (a0+1), 240 ))
time1 <- rep(0, length(tList))
for(i in 1:length(tList)){
t.end <- tList[i]
time1 <- system.time(
tpm2 <- dbd_prob(t.end, a0, b0, drates1, brates2, drates2, trans,
a=A, B))#, computeMode=2))
#MC simulation "ground truth"
tpm.MC <- getTrans.MC(nSim, t.end, a0, b0, beta, alpha)
tpm3 <- tpm.MC[1:(a0+1), ]
#store subset of matrices containing about 99 percent of the mass:
tpmArray[i,1,,] <- tpm2[1:(a0+1),1:60]
tpmArray[i,2,,] <- tpm3[1:(a0+1),1:60]
}
save(tpmArray, time1,
file = "../inst/vignetteCache/larger_pop.RData",
compression_level = 9)
}
## ---------------------------------------------------------------------------------------
par(mfrow = c(2,1))
image(tpmArray[1,1,,1:60], main = "Continued Fraction approximation, t=.5")
image(tpmArray[1,2,,1:60], main = "Monte Carlo estimates")
par(mfrow = c(2,1))
image(tpmArray[2,1,,1:60], main = "Continued Fraction approximation, t=1")
image(tpmArray[2,2,,1:60], main = "Monte Carlo estimates")
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.