[過去ログ]
数学 統計に詳しい人が語るコロナウイルス (1002レス)
数学 統計に詳しい人が語るコロナウイルス http://rio2016.5ch.net/test/read.cgi/math/1582910321/
上
下
前
次
1-
新
通常表示
512バイト分割
レス栞
このスレッドは過去ログ倉庫に格納されています。
次スレ検索
歴削→次スレ
栞削→次スレ
過去ログメニュー
114: 132人目の素数さん [sage] 2020/03/24(火) 11:57:10.20 ID:TnHQvRcs 上記の準備をして以下で実行 PCRj2 <- function( N,X, UL=1, SEN=0.7, SPC=0.9, SD=0.05, print=TRUE){ # UL:upper limit of dunif(0,UL) library(rjags) library(BEST) sn=Mv2ab(SEN,SD^2) sp=Mv2ab(SPC,SD^2) modelstring=paste0(' model { x ~ dbin(p,n) p <- prev*sen + (1-prev)*(1-spc) sen ~ dbeta(sn[1],sn[2]) spc ~ dbeta(sp[1],sp[2]) prev ~ dunif(0,ul) } ') writeLines(modelstring,'TEMPmodelj.txt') dataList=list(n=N,x=X,ul=UL,sen=SEN,spc=SPC,sn=sn,sp=sp) jagsModel = jags.model( file="TEMPmodelj.txt" ,data=dataList, quiet=TRUE) update(jagsModel) codaSamples = coda.samples( jagsModel , variable=c("prev","p","sen","spc"), n.iter=1e5, thin=5) js=as.matrix(codaSamples) if(print){ BEST::plotPost(js[,'prev'],xlab='prevalence',showMode = TRUE) lines(density(js[,'prev']),col='skyblue')} re=c(mean=mean(js[,'prev']),HDInterval::hdi(js[,'prev'])[1:2]) return(re) } options(digits = 5) options(scipen = 5) PCRj2(1000,10) # 陽性率1%で有病率を推定 PCRj2(1000,300) # 陽性率30%で有病率を推定 PCRj2(1000,600) # 陽性率60%で有病率を推定 http://rio2016.5ch.net/test/read.cgi/math/1582910321/114
メモ帳
(0/65535文字)
上
下
前
次
1-
新
書
関
写
板
覧
索
設
栞
歴
あと 888 レスあります
スレ情報
赤レス抽出
画像レス抽出
歴の未読スレ
AAサムネイル
Google検索
Wikipedia
ぬこの手
ぬこTOP
0.080s*