>>750
バグに気づいたので撤回。
infectee==10ではなくてinfector+infectee==10だな。



デバッグ版

R=0.58 # mean of reproductive number
k=0.45 # dispersion parameter
(prob = k/(k+R)) # its probability
Rt=rnbinom(1e5,k,mu=R) # random numbers of negative binomial distribution
hist(Rt,breaks = 'scott',freq=F,ann=F) # show its histgram
sim <- function(n=10){ # simulation
infected=0 # initial value
while(infected!=n){ # while infected is unequal to n
infector=sample(n,1) # prior discrete uniform distribution of infector number
infectee=sum(sample(Rt,infector)) # number of infectee
infected=infectee+infector # number of infected
}
return(infector) # when n infected, 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
summary(spreader)
sum(spreader==1)/length(spreader) # the probability of single super-spreader

まあ、パーティー中に感染させる人数に再生産数を流用していいかは疑問ではあるが、

実効結果は

https://i.imgur.com/ApDSW1J.png

95% CI
> HDInterval::hdi(spreader)[1:2] # 95% credibility interval
lower upper
3 9

中央値、平均値、四分位値
> summary(spreader)
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.000 6.000 7.000 6.844 8.000 10.000

1人のスーパースプレッダーの確率
> sum(spreader==1)/length(spreader) # the probability of single super-spreader
[1] 0.00104