数学統計に詳しい人が語るコロナウイルス

[過去ﾛｸﾞ] 数学統計に詳しい人が語るコロナウイルス (1002ﾚｽ)
上下前次1-新
通常表示 512ﾊﾞｲﾄ分割ﾚｽ栞
抽出解除ﾚｽ栞

このｽﾚｯﾄﾞは過去ﾛｸﾞ倉庫に格納されています｡
次ｽﾚ検索歴削→次ｽﾚ栞削→次ｽﾚ過去ﾛｸﾞﾒﾆｭｰ

ﾘﾛｰﾄﾞ規制です｡10分ほどで解除するので､他のﾌﾞﾗｳｻﾞへ避難してください｡

750: １３２人目の素数さん [sage] 2020/07/05(日) 13:11:55 ID:rpUuAKzr

>672のデータ
https://www.researchsquare.com/article/rs-29548/v1

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).

を使って乱数発生させてシミュレーションしてみた。

mu=0.58
size=0.45
(prob = size/(size+mu))
Rt=rnbinom(1e5,size=size,mu=mu) # reproductive number
hist(Rt,breaks = 'scott',freq=F,ann=F)
sim <- function(n=10){
infector=sample(1:n,1) #  how many infectors?
infectee=sum(sample(Rt,infector)) # how many infectees?
if(infectee==n) return(infector)  # is the number of infectees n?
else return(NA)
}
s=replicate(1e6,sim())
spreader=s[!is.na(s)]
hist(spreader,freq=F,ylab='',axes=F,breaks='scott') ; axis(1)
BEST::plotPost(spreader,xlim=c(0,10))
mean(spreader)
HDInterval::hdi(spreader)[1:2]
sum(spreader==1)/length(spreader)

期待値は
> mean(spreader)
[1] 8.032594
95％CIは
> HDInterval::hdi(spreader)[1:2]
lower upper
4    10
1人のスーパースプレッダーから広がった確率は
> sum(spreader==1)/length(spreader)
[1] 0.001580299

http://rio2016.5ch.net/test/read.cgi/math/1582910321/750

751: １３２人目の素数さん [sage] 2020/07/05(日) 20:02:50 ID:rpUuAKzr

>>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

http://rio2016.5ch.net/test/read.cgi/math/1582910321/751

上下前次1-新書関写板覧索設栞歴

ｽﾚ情報赤ﾚｽ抽出画像ﾚｽ抽出歴の未読ｽﾚ AAｻﾑﾈｲﾙ

ぬこの手ぬこTOP 0.034s