negBinExtinction.R

invisible('
Extinction probability as a function of 
- Reproduction number R
- Dispersion k
- Number imported infections

Assumes negative binomial offspring distribution
') 

# negative binomial probability generating function 
#
# @parms R mean (reproduction number)
# @param k overdispersion
g <- function(x, R, k ){
	sigma2 <- R * (1 + R / k )
	mu <- R
	p <- (sigma2 - mu) / sigma2
	r <- mu^2 / ( sigma2 - mu)
	(  (1-p) / ( 1 - p*x ) )^r
}

# extinction probability as function of R , k 
ext1 <- function( R, k, u = 1e-9, iter = 1e3 )
{
	for ( i in 1:iter) 
		u <- g( u, R, k )
	u
}

# extinction probability as function of number of imports, R, k,  ... 
ext <- function( nimports , R, k , ... )
{
	e1 <- ext1( R, k, ... )
	e1^nimports
}

# probability of epidemic = 1 - extinction probability 
pepidemic <- function( nimports, R, k, ... )
{
	1 - ext( nimports, R, k, ... )
}


n <- 1:10
p <- pepidemic( n, 2.2, 0.16 )

png('ext0.png', width = 480, height = 300)
plot( p , xlab = 'Number of undetected imported cases' , ylab = 'Probability of epidemic', type = 'b')
dev.off()