>>752
レスありがとうございます。

パラメータの分布を勝手に設定して階層モデルもどきでシミレーションしてみた。数字のお遊びでしょうけど。

From the empirical offspring distribution and
fitted negative binomial distribution shown in Figure 2B,
we estimated an observed reproductive number (R) of 0.58 (95% CI: 0.45 - 0.71)
and dispersion parameter (k) of 0.45 (95% CI: 0.31 - 0.76).

Rは正規分布、kはガンマ分布として、期待値と95%CIが上記に適合するようなパラメータを
ニュートン・ラフソン法で算出すると
正規分布の標準偏差sd＝0.0663277494
> pnorm(0.71,0.58,0.0663277494)-pnorm(0.45,0.58,0.0663277494)
[1] 0.95
ガンマ分布の形状パラメータsh＝24.38466,尺度パラメータ=0.01845
> pgamma(0.76,sh,scale=sc)-pgamma(0.31,sh,scale=sc)
[1] 0.95
が得られた。

各パラメータが上記の分布に従うとして
R ~ normal(mean=0.58,sd=sd)
k ~ gamma(shape=sh,scale=sc)
Rt ~ nbinom(dispersion=k,mean=R)
のモデルでシミレーション

既述のコードに
RR=rnorm(1e5,R,sd) # R ~ dnorm(R,sd)
kk=rgamma(1e5,sh,scale=sc) # k ~ dgamma(shape=sh,scale=sc)
RRt <- function(x){ # Rt ~ dnbinom(k,mu=R)
k=sample(kk,1)
R=sample(RR,1)
rnbinom(x,k,mu=R) # pick x random Rt numbers
}
を追加して
while loop内のinfecteeの数を
infectee=sum(RRt(infector)) # number of infectee
に変更。

最初の10人に何人が感染していたかの95％推定
> HDInterval::hdi(spreader)[1:2] # 95% credibility interval
lower upper
3 9
> BEST::plotPost(spreader,xlim=c(1,10)) # graph with 95%CI & mean
https://i.imgur.com/LKu4hKD.png
> summary(spreader)
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.000000 6.000000 7.000000 6.894006 8.000000 10.000000
> sum(spreader==1)/length(spreader) # the probability of single super-spreader
[1] 0.001023

予想通り、ほとんど、結果は変わらず。乖離した方がおかしい。