>>17
数値を変化させてグラフ化できるように関数化
事前分布は一様分布でなくJeffereysを採用
必要に応じて指定

"
PCR検査の感度を0.7、特異度を0.9とする。
広島県で第一号の感染発見例は
県の検査で１回陰性、病院の検査で２回陽性、症状軽快した現時点で陰性（何回やったか報道がないので１回陰性とする）であるという。

検査前の感染確率の分布が一様分布であると仮定して、
現在患者が感染している確率とその95%CIを計算してみた。
"
PCR2prob <- function(
sn=0.7, # sensitivity
sp=0.9, # specificity
plus=2, # how many positive result?
minus=2, # how many negative result?
n=1e5,
p0=rbeta(n,0.5,0.5), # prior : Jeffreys' distribution
print=TRUE) # how large the simulation
{
oz0=p0/(1-p0) # prob -> odds
pLR=sn/(1-sp) # TP/FP
nLR=(1-sn)/sp # FN/TN
oz1=oz0*pLR^plus*nLR^minus # Bayesian formula
p1=oz1/(1+oz1) # odds -> prob
if(print & length(p0)>1){ # p0 ~ some distribution
BEST::plotPost(p1,showMode =T) # show mode instead of mean
print(HDInterval::hdi(p1)) # Highest Density Interval
print(quantile(p1,c(.025,0.5,.975))) # 95%CI by quantile
print(summary(p1)) # mean, median
dens = density(p1) # print mode
mode_i = which.max(dens$y)
print(c(Mode=dens$x[mode_i]))
}
if(length(p0)==1) print(p1) # when p0 is point-designated
invisible(p1)
}

PCR2prob()
PCR2prob(p0=rbeta(1e5,1,1)) # p0 ~ uniform distiribution
PCR2prob(p0=0.5) # point probability
PCR2prob(minus=3) # one more negative result
PCR2prob(minus=4) # two more negative result