- 98 名前：１３２人目の素数さん mailto:sage [2020/03/23(月) 03:28:01.22 ID:uvHIelYA.net]
- 日本人の平均身長を推測するのにその値は１〜２ｍの間であるという弱情報事前分布は合理的。
現時点での新型コロナの有病率は0.1未満の一様分布という弱情報事前分布として
【富山県最強伝説】新型コロナウイルスPCR検査件数 54人 陽性0人
ある集団から54人を無作為に選んでPCR検査したら陽性0であった。感度0.7 特異度0.9としてこの集団の有病率の期待値と9信頼区間を推測する。
事前分布のパラメータを変えるとstanだとコンパイルが必要になるのでjagsでプログラムを組んでみた。
# 感度SEN, 特異度SPCの検査でN人中X人が陽性であったときの推定有病率prevalence
# 弱情報事前分布はprevalence ~ dunif(0,UL) , UL:上限
library(rjags)
PCRj <- function(N,X,UL=1,SEN=0.7,SPC=0.9,verbose=TRUE){ # UL:upper limit of dunif(0,UL)
modelstring=paste0('
model
{
x ~ dbin(p,n)
p <- prev*sen + (1-prev)*(1-spc)
prev ~ dunif(0,',UL,')
}
')
if(verbose & UL!=1) cat(modelstring)
writeLines(modelstring,'TEMPmodel.txt')
dataList=list(n=N,x=X,sen=SEN,spc=SPC)
jagsModel = jags.model( file="TEMPmodel.txt" ,data=dataList,quiet=TRUE)
update(jagsModel)
codaSamples = coda.samples( jagsModel ,
variable=c("prev","p"), n.iter=1e6, thin=10)
js=as.matrix(codaSamples)
BEST::plotPost(js[,'prev'],xlab='prevalence',showMode = TRUE) ; lines(density(js[,'prev']),col='skyblue')
round(c(mean=mean(js[,'prev']),HDInterval::hdi(js[,'prev'])[1:2]),10)
}
実行結果
> PCRj(54,0,UL=0.1)
model
{
x ~ dbin(p,n)
p <- prev*sen + (1-prev)*(1-spc)
prev ~ dunif(0,0.1)
}
|**************************************************| 100%
mean lower upper
0.0245104429 0.0000003782 0.0703606657
|
 |
