rm(list=ls())

mu=0.58 # mean of reproductive number
size=0.45 # dispersion parameter
(prob = size/(size+mu)) # its probability
Rt=rnbinom(1e5,size=size,mu=mu) # random numbers of negative binomial distribution
hist(Rt,breaks = 'scott',freq=F,ann=F) # show its histgram
sim <- function(n=10){ # simulation
infectee=0 # initial value
while(infectee!=n){ # while infectee is unequal to n
infector=sample(n,1) # prior discrete uniform distribution of infector number
infectee=sum(sample(Rt,infector)) # total number of infectee
}
return(infector) # when n infectee, return infector number
}
spreader=replicate(1e5,sim()) # simulation & calculation
hist(spreader,freq=F,ylab='',axes=F,breaks='scott') ; axis(1)
HDInterval::hdi(spreader)[1:2] # 95% credibility interval
BEST::plotPost(spreader) # graph with 95%CI & mean
sum(spreader==1)/length(spreader) # the probability of single super-spreader