臨床統計もおもしろいですよ、その3

197 ：卵の名無しさん：2020/05/27(水) 06:40:49.96 ID:uX2xSEBS.net

>>195
モデルは簡単なのでStanやJagsを使わずに完成

WhoInfectedFirst <- function( # 発症間隔から感染源の確率を計算
t=c(1,2), # 発症間隔
k=1e5 # 乱数発生数
){
ln_par1 = 1.434065 # 潜伏期間対数正規分布パラメータ
ln_par2 = 0.6612
n=length(t)+1 # 発症人数

# 潜伏期間
x_incub = matrix(rep(NA,n*k),ncol=n)
for(i in 1:n) x_incub[,i] = rlnorm(k,ln_par1,ln_par2)

# 感染日（一人目の発症日：０）
x_infected = matrix(rep(NA,n*k),ncol=n)
x_infected[,1]= -x_incub[,1]
for(i in 2:n) x_infected[,i] = sum(t[1:(i-1)]) - x_incub[,i]

# i番目の発症者が感染源の確率
fi <- function(i){
mean(apply(x_infected,1,function(x) which.min(x)==i))
}
data.frame(p=round(sapply(1:n,fi),2))
}
WhoInfectedFirst(c(1,2))
WhoInfectedFirst(c(1,0,1,0,0)) # 翌日2人発症 翌々日3人発症