- 1 名前：１３２人目の素数さん [2020/02/29(土) 02:18:41 ID:twdO677Q.net]
- 東大数学科卒の元官僚はこう分析してるが、お前らはどうなると思う？
www.zakzak.co.jp/soc/news/200220/dom2002200003-n2.html
中国国外感染者の中国国内との比率をみると、
１月２０日の数字公表以降は、０・８～２・６％で比較的安定している。
これは、新型肺炎の感染者のほとんどは中国国内、それも湖北省に集中しているからだ。
ちなみに中国国外での感染者数は、中国国内の１・１％だ（２月１６日現在）。
本コラムで紹介したが、現時点では、最終的な中国国内の感染者数は２０万人超と筆者は推計している。
となると、中国国外の感染者は数千人程度になるだろう。
中国国外のうち日本の比率は１割弱なので、日本の感染者数は数百人程度であろう。
その場合、死者も数人から１０人程度になるだろう。
こうした推計をすると、今の感染者は氷山の一角だと思われるが、今後の増加ペースはどうなるだろうか。
新型コロナウイルスの検査は簡単に行えるので、今後、日本での感染者数は増えていくだろう。
ある時点ではそれがネズミ算的に増えるかのように思える局面もあるだろうが、
筆者の推計が正しければ、現時点ではせいぜい数百人が一つのメドだ。
- 119 名前：１３２人目の素数さん mailto:sage [2020/03/24(火) 11:57:10 ID:TnHQvRcs.net]
- 上記の準備をして以下で実行
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％で有病率を推定
- 120 名前：１３２人目の素数さん mailto:sage [2020/03/24(火) 13:15:20 ID:/QqkwKRd.net]
- >>99
期待値というのは、無次元量ではない。観測値とか物理量と同じように単位をつけて議論できる量。
従って「期待値が10%増える」等という言葉があれば、期待値が1.1倍になるのだろうと感じるのが普通。
そのような性質を持つ期待値に対し、「10%増える」と表現し、
「期待値の値そのものが、0.1増えることを意味している」
と説明しなければならないならば、やはり誤解を招きやすい表現だと思う。
今回の期待値は比率であり、無次元量であったから、「10%」と言うのが、
どちらの意味としても、通用したため発生したとは言えるが、読み手の立場に立った表現を望む。
似た議論に、選挙時の投票率がある。前回の投票率が40%。今回の投票率が50%だとする。
「前回に比べ、今回は10%増えました」
「前回に比べ、今回は25%増えました」
どちらも、言い得る表現。聞き手の混乱を避けるため、前者の意味で使う場合、
「10%ポイント増えました」とコメントするのを最近聞くようになった。
私にはよい傾向と感じるが、中には、違いは何かとか、混乱の源の存在さえ意識していない人もいるようだ。
「３割増も４割強増も大した差ではない」には、「式が違っても結果が誤差範囲なら問題ない」
という考えが背景に見える。そのような方が、混乱を引き起こしかねない表現を用いた。
だから、補足した。果たして本当に杞憂だったのだろうか？
- 121 名前：１３２人目の素数さん mailto:sage [2020/03/24(火) 14:42:23 ID:TnHQvRcs.net]
- 富山では62人PCR検査して陽性0人（3月22日までの集計)有病率を推定とその信頼区間を推定したい。
www.pref.toyama.jp/cms_pfile/00021629/01366377.pdf
PCR検査の感度は最頻値0.6標準偏差0.１、特異度は最頻値0.9標準偏差0.05のベータ分布（正規分布は負になったり1を超えるので不適）、
有病率は一様分布として、推定される有病率の期待値と95％を計算せよ。
図示するとこんな感じ。
https://i.imgur.com/Ip6gSCa.png
stanのモデルのスクリプトはこれ
sn,spはβ分布のパラメータ、その計算法は既述
data{
int n; // sample size
int x; // positive test result
real<lower=0,upper=1> ul; // uniform(0,ul)
real<lower=0> sn[2]; // sen ~ beta(sn[1],sn[2])
real<lower=0> sp[2]; // spc ~ beta(sp[1],sp[2])
}
parameters{
real<lower=0,upper=1> prev; // prevalence
real<lower=0,upper=1> sen; // sensitivity
real<lower=0,upper=1> spc; // specificity
}
transformed parameters{
real<lower=0,upper=1> p;
p = prev*sen + (1-prev)*(1-spc) ; // probability of positive test result
}
model{
x ~ binomial(n,p);
prev ~ uniform(0,ul);
sen ~ beta(sn[1],sn[2]);
spc ~ beta(sp[1],sp[2]);
}
- 122 名前：１３２人目の素数さん mailto:sage [2020/03/24(火) 14:54:23 ID:TnHQvRcs.net]
- >>116
ここで問題
感度特異度の分布はそのベータ分布として
何人陰性が続けば95%信頼区間の上限が0.05を下回るか？
- 123 名前：１３２人目の素数さん mailto:sage [2020/03/24(火) 15:14:26 ID:mBslr8ul.net]
- >>117
感度、特異度の分布？？？
- 124 名前：１３２人目の素数さん mailto:sage [2020/03/24(火) 15:56:48 ID:TnHQvRcs.net]
- >>118
何でも確率変数にするのがベイズ推計。
ｐ値の分布すら考えるぞ。
- 125 名前：１３２人目の素数さん mailto:sage [2020/03/24(火) 16:10:56 ID:TnHQvRcs.net]
- PCR検査の感度は最頻値0.6標準偏差0.１、特異度は最頻値0.9標準偏差0.05のベータ分布を事前分布にしたけど、
事後分布はstanによるMCMCで
感度は期待値0.57 95％信頼区間は[0.37,0.77]
特異度は期待値0.96 95%信頼区間は[0.91,0.99]
とコンピュータが計算してくれる。
- 126 名前：１３２人目の素数さん mailto:sage [2020/03/24(火) 17:38:05 ID:TnHQvRcs.net]
- >116のように弱情報事前分布を設定することで事後分布は次のように描ける。
https://i.imgur.com/J1Xqdfj.png
- 127 名前：１３２人目の素数さん mailto:sage [2020/03/24(火) 17:43:56 ID:TnHQvRcs.net]
- >>54
いや、特異度の事前分布を設定することで事後分布をMCMCで求めることができる。
>116の設定での結果が>121
- 128 名前：１３２人目の素数さん [2020/03/24(火) 18:10:37.68 ID:EUfp1x4d.net]
- >>122
結局事前分布の設定次第ってことはないの？
- 129 名前：１３２人目の素数さん mailto:sage [2020/03/24(火) 18:45:58.22 ID:TnHQvRcs.net]
- >>123
日本人の平均身長を推測するのにその値は１～２ｍの間であるという弱情報事前分布は合理的。
感度特異度の分布に正規分布を使うのはアホ。
負になったり、１を超えたりするから。
- 130 名前：１３２人目の素数さん mailto:sage [2020/03/24(火) 19:08:27.01 ID:TnHQvRcs.net]
- >>123
感度を0.4-0.8の一様分布、特異度を0.8-1.0の一様分布にしても有病率の推定値は
> round(re$mci,5)
mean lower upper
0.02827 0.00000 0.08592
であまり変わらないね。
- 131 名前：１３２人目の素数さん mailto:sage [2020/03/24(火) 19:20:40 ID:TnHQvRcs.net]
- sensitivity ~ N(m=0.6,sd=0.1) specificity ~ N(m=0.9, sd=0.05)
にしても推測有病率は平均3％弱で 95%CIは0-8％とあまり分布の形にはよらないね。
mean lower upper
0.026841384 0.000000153 0.081071379
確率だと定義域が0-1で計算しやすいのでβ分布を使うことが多い。
- 132 名前：１３２人目の素数さん mailto:sage [2020/03/24(火) 21:43:16 ID:mBslr8ul.net]
- >>119
そう？
統計の推定の理論で推計する母数は確率変数ではないと習ったけど？
- 133 名前：１３２人目の素数さん mailto:sage [2020/03/25(水) 05:45:40 ID:jmNOx22O.net]
- >>127
時代は頻度主義統計からベイズ統計だよ。
- 134 名前：１３２人目の素数さん mailto:sage [2020/03/25(水) 06:04:40 ID:jmNOx22O.net]
- 頻度主義統計でも最尤推定では
データを固定してパラメータを動かすだろ。
- 135 名前：１３２人目の素数さん mailto:sage [2020/03/25(水) 06:30:36 ID:jmNOx22O.net]
- >>127
階層ベイズモデルを扱ったことないの？
>112は簡単な実例。
- 136 名前：１３２人目の素数さん mailto:sage [2020/03/25(水) 07:18:32 ID:yWXBkNWD.net]
- >>128 >>130
何を持ってベイズ統計っていってんのか知らん。
pcr検査の感度とは被験者が感染者である場合の検査結果が陽性となる条件付き確率でしょ？
条件付き確率の分布ってどういうことよ？
確率がまた確率になるってなんの話してんの？
変数Xの平均とか分散とかは統計学においては推定すべき定数であって確定値。
それの分散なんて数学的に意味不明。
一体どこの統計学の教科書にそんなデタラメ書いてあんの？
- 137 名前：１３２人目の素数さん mailto:sage [2020/03/25(水) 07:23:50 ID:r1V62jxn.net]
- まちがえた。
確率の平均がまた確率変数になるってどういうことよ、ね。
式でかけば確率変数Xの平均E(X)の分散ってなんの話ってことになる。
確率変数Xはある標本空間上の関数だけどE(X)は実数だよ？
- 138 名前：１３２人目の素数さん mailto:sage [2020/03/25(水) 09:57:10 ID:2o2
]
- [ここ壊れてます]
- 139 名前：7M3ww.net mailto: >>131
ベイズ階層モデルも組めない奴とは議論にならんね。
分散の事前分布に逆ガンマ分布でなく半コーシーを使う方がいいとかいう議論も理解できんだろ。 [] - [ここ壊れてます]
- 140 名前：１３２人目の素数さん mailto:sage [2020/03/25(水) 10:06:13 ID:2o27M3ww.net]
- >>131
ベータ分布は定義域が[0,1]で二項分布の確率の確率密度関数としてベイズ階層モデルでは頻用されるよ。
ベイズ階層モデルを使わずにこの計算できるならやってみてくれ。
020/3/24 11:00時点で検査人数での陽性率は171/2013であるという。
新型コロナ肺炎のPCR検査の感度は５～７割、特異度は９割前後らしい。幅をもたせた値を使って検査をうけたグループの有病率を計算せよ。
- 141 名前：１３２人目の素数さん mailto:sage [2020/03/25(水) 11:27:50 ID:82yASlvk.net]
- >>133
まぁ言わんとする事はもちろんわかるし伝わるけど、疫学だから数学やってる人間がなんとなく伝わるではダメだろ？
数学だけの話ではなく、疫学は実社会とキチンと繋がってるんだから？
統計学ではあくまで検定する母数は定数。
それは確率モデルでは実数値であり、定数。
そして統計データを確率変数に割り当てる。
当然それらの確率変数は一つの測度空間の一つしかない確率変数であり、平均も分散もひとつしかない定数値。
それらをいっぱい考えてどうこう言ってるんだろうとは思うけどそんなの統計学や疫学の一般的な考えにはない。
何故なら現実世界はひとつしかなく、確率変数に対応している統計量も一個しかない。
もちろん母数がめちゃめちゃ大きい統計量で例えば10000個のデータを100こずつ切って100個の統計量を100の世界からとってきたなんて考えが無理クリできなくはないが、そんな考え方は普通しない。
それはあくまで100個ずつに区切られた10000個の一つの世界の確率変数としか扱わない。
そういうオリジナルな考えで捉えたいならそれは勝手だけど、それならそれで話の中で明示しないとダメ。
数学の世界なら言わずもがなの話は言わなくてもエスパーしてもらえても、疫学、統計学の世界では実社会とつながる話だからダメ。
- 142 名前：１３２人目の素数さん mailto:sage [2020/03/25(水) 12:30:50.29 ID:jmNOx22O.net]
- >>135
能書きいいから、
ベイズ階層モデルを使わずにこの計算できるならやってみてくれ。
020/3/24 11:00時点で検査人数での陽性率は171/2013であるという。
新型コロナ肺炎のPCR検査の感度は５～７割、特異度は９割前後らしい。幅をもたせた値を使って検査をうけたグループの有病率を計算せよ。
- 143 名前：１３２人目の素数さん mailto:sage [2020/03/25(水) 12:39:07.64 ID:jmNOx22O.net]
- >>136
こういう判断が現実には必要。
検査特性を無視して単純な割り算だと検査を受けた人の有病率は8.5%弱になるけどこれは過大評価か過小評価か？
- 144 名前：１３２人目の素数さん mailto:sage [2020/03/25(水) 14:59:13.04 ID:jmNOx22O.net]
- 検査感度が5-7割、特異度が9割前後なら
検査陽性率=有病率とすると常に過大評価かどうか気になったので陽性数を変化させて計算してみた。
検査感度はmode=0.6,sd=0.1 特異度はmode=0.9,sd=0.05のベータ分布に設定してJAGSでベイズ階層モデルをたてて計算。
https://i.imgur.com/zTdxRrb.png
陽性率が20%未満のときは過大評価、それ以上のときは過小評価である、という結論になった。
ベイズ統計を理解できている人の検証希望。
- 145 名前：１３２人目の素数さん mailto:sage [2020/03/25(水) 17:30:41 ID:jmNOx22O.net]
- >>138
プログラムの練習がてらに、
MCMCのアルゴリズムの異なるstanでベイズ階層モデルを組んで検証。
当然ながら、同様の結果。 検査陽性率が２０％を境に過大評価と過小評価が入れ替わる。
https://i.imgur.com/ItSNWdD.png
- 146 名前：１３２人目の素数さん mailto:sage [2020/03/25(水) 21:21:08.00 ID:jmNOx22O.net]
- >>136(自己レス)
今日の都の発表で(171+41)/(2013+89) に検査陽性率が増えたので再計算。
https://i.imgur.com/THdYDqT.png
- 147 名前：１３２人目の素数さん mailto:sage [2020/03/25(水) 21:33:22.10 ID:jmNOx22O.net]
- ");
//]]>-->
- 148 名前：8" rel="noopener noreferrer" target="_blank" class="reply_link">>>138
サンプリング回数を増やしてグラフを完成。
https://i.imgur.com/kLjCD2y.png []- [ここ壊れてます]
- 149 名前：１３２人目の素数さん mailto:sage [2020/03/26(木) 16:25:58 ID:+rQz06p5.net]
- >>140
89は検査数で検査人数は74という。
計算し直すと
> PCRj2(N,r,SEN=0.6,SD1=0.1,SPC=0.9,SD2=0.05,N.ITER=5e5)
|**************************************************| 100%
mean lower upper
0.05720165 0.00000015 0.1332385
- 150 名前：１３２人目の素数さん mailto:sage [2020/03/26(木) 16:34:28 ID:+rQz06p5.net]
- 41/74の推測有病率は
mean lower upper
0.8121975 0.5957315 0.9999992
- 151 名前：１３２人目の素数さん mailto:sage [2020/03/27(金) 11:07:27 ID:sdGiAEI7.net]
- オリンピック延期発表後の検査陽性率は88/169で52%だが、
PCR検査の感度と特異度がはっきりしないので、検査陽性率をこの集団の有病率とするのは正しくない。
88/169のときの感度・特異度と推定有病率の関係をグラフにしてみた。
https://i.imgur.com/iQC88tZ.png
感度0.6、特異度0.9のときの推定有病率は85%で陽性率からの憶測は過小評価といえる。
- 152 名前：１３２人目の素数さん [2020/03/27(金) 18:36:04.97 ID:8rq7DP6B.net]
- 検査陽性率が小さいときには、実際の有病率より過大評価してるし、
検査陽性率が高いときは、過小評価してるだろうってことでしょ。
そのくらいは定性的に理解できる。
- 153 名前：１３２人目の素数さん mailto:sage [2020/03/27(金) 20:59:08.67 ID:sdGiAEI7.net]
- >>145
どこが境目かは直感じゃわからんね。
- 154 名前：１３２人目の素数さん [2020/03/27(金) 22:15:19.84 ID:8rq7DP6B.net]
- そりゃ感度や特異度次第だからな。
まあ、数％と数十％では違うんだということがわかればいいんじゃね？
境目なんかどうでもいいでしょ。
- 155 名前：１３２人目の素数さん mailto:sage [2020/03/27(金) 22:22:39.69 ID:sdGiAEI7.net]
- 陽性率が15%でこれを有病率の推測値に使うのは過大評価なのか過小評価がわからんのはまずいね。
- 156 名前：１３２人目の素数さん mailto:sage [2020/03/27(金) 23:24:27 ID:sdGiAEI7.net]
- オリンピック延期決定以後の検査数と陽性数
subjects=c(74,95,87)
positives=c(41,47,40)
PCRs3(subjects,positives,iter=10000,warmup=1000)
として、
感度・特異度を考慮した推定有病率は
mean lower upper
0.77417 0.56756 0.99944
>
日々の陽性数が二項分布に従うとして計算。
- 157 名前：１３２人目の素数さん [2020/03/28(土) 03:24:30.89 ID:NK6wIjWT.net]
- 志村けんみたいな有名人がコロナに感染してることから日本全体のコロナ感染者数を推定してみる。
まず、日本の有名人が1000人いるとしよう。
つぎに、日本でコロナに感染していない確率をxとしよう。
すると、有名人1000人が一人も感染していない確率は、xの1000乗となる。これをyとおこう。すると、有名人が一人でも感染している確率は(1-y)となる、これをzとおこう。
まとめると以下の関係がなりたつ。
・コロナに感染しない確率:x
・有名人が一人もコロナに感染しない確率:y=x^1000
・有名人が一人でもコロナに感染している確率:z = 1-y
- 158 名前：１３２人目の素数さん [2020/03/28(土) 03:25:09.78 ID:NK6wIjWT.net]
- 志村は感染したわけなので、以下、2つのケースにわける
ケース1: zが10%のとき
z=0.1, 故にy = 1-0.1=0.9
故にx = y^0.001よりx=0.9^0.001=0.999894
これがコロナに感染していない確率なので、
コロナ感染確率は、1-0.999894=0.000106
よって日本のコロナ感染者数は推定
120,000,000*0.000106=12,720人
ケース2:zが50%のとき
ケース1と同様の計算で、
日本のコロナ感染者数は推定
120,000,000*0.000693=83,160人
- 159 名前：１３２人目の素数さん [2020/03/28(土) 09:37:19.12 ID:uwBdnirU.net]
- 検査が少ないから感染者増が緩やか？数学的に検証してみた
agora-web.jp/archives/2045047.html
主な関係国について、新型コロナ感染者数の片対数グラフがある。
agora-web.jp/cms/wp-content/uploads/2020/03/WS000876.jpg
FT.comより
感染者数の伸びが日本は緩やかと解釈するのが普通だが、検査が少ないからとする解釈もある。本当はどうなのか計算してみる。
結論を先に書くと、検査が多いか少ないかは関係ない。
- 160 名前：１３２人目の素数さん mailto:sage [2020/03/28(土) 09:40:23.67 ID:QZo3p56d.net]
- 対数をとると係数（感染者の発見率）は定数項になり、今回の片対数グラフ�
- 161 名前：ﾌ整理法の前提としてキャンセルされる
日本が展開しているのは患者認定の精度上昇であり、医療リソースの効果を最大化して死者数を低く抑えている要因の一つといえる [] - [ここ壊れてます]
- 162 名前：１３２人目の素数さん mailto:sage [2020/03/28(土) 10:39:35.99 ID:BJlezchp.net]
- キャバクラ客100人から無作為に5人から検体を採取してこの検体を混合攪拌してコロナ検査したところ陽性であった。
（１）100人のキャバクラ客の陽性数の期待値と95％信頼区間を求めよ。
（２）PCR検査の感度0.6、特異度0.9として100人のキャバクラ客の感染数の期待値と95%信頼区間を求めよ。
- 163 名前：１３２人目の素数さん mailto:sage [2020/03/28(土) 11:58:45.95 ID:BJlezchp.net]
- >>151
> m=1000 # 有名人の人数
> n=1.268e5 # 日本の人口
> x=0:n # 感染者数:x, 非感染数:n-x
> pmf=1- chooseZ(n-x,m)/chooseZ(n,m) # 1 - (ｍ人全員非感染の確率)
> pdf=pmf/sum(pmf) # 確率密度関数化して
> (E=sum(x*pdf)) # 期待値を計算
Big Rational ('bigq') :
[1] 63590201/1002
> as.numeric(E)
[1] 63463.27
6万3000人と計算された。
- 164 名前：１３２人目の素数さん [2020/03/28(土) 12:46:07.99 ID:NK6wIjWT.net]
- >>155
良く分からんが、ありがとう。
こちとら高校レベルの確率の知識しかないもんで。
- 165 名前：１３２人目の素数さん mailto:sage [2020/03/28(土) 15:42:06.57 ID:BJlezchp.net]
- >>156
n(=10)人の中にi人の感染者がいるとき無作為にm(=2)人を選ぶ。
選ばれた2人の中に少なくとも一人の感染者がいる確率をP[x]として、
n個からr個選ぶ組み合わせの数をChoose(n,r)で表すと
P[xi]=1- choose(10-x,2)/choose(10,2)
xを0から10まで変化させて、
Σx*P[x]/(ΣP[x])で
期待値が求まる。
- 166 名前：１３２人目の素数さん mailto:sage [2020/03/28(土) 15:42:43.27 ID:BJlezchp.net]
- タイプミス修正
P[x]=1- choose(10-x,2)/choose(10,2)
- 167 名前：１３２人目の素数さん mailto:sage [2020/03/28(土) 16:07:51.35 ID:qsSYTF8t.net]
- 何このアホスレ？
- 168 名前：１３２人目の素数さん mailto:sage [2020/03/28(土) 16:53:08.32 ID:BJlezchp.net]
- 有名人の数を増やしてみても同様の結果になった。
> # 有名人が感染
> library(gmp)
> m=18200 # 有名人の数(桜を見る会参加人数)
> n=1.268e5 # 日本の人口
> x=0:n # 感染者数:x, 非感染数:n-x
> pmf=1- chooseZ(n-x,m)/chooseZ(n,m) # 1 - (ｍ人全員非感染の確率)
> pdf=pmf/sum(pmf) # 確率密度関数化して
> (E=sum(x*pdf)) # 期待値を計算
Big Rational ('bigq') :
[1] 1154070201/18202
> as.numeric(E) # E=63463.27 (m=1000) , E=1154070201/18202=63403.48(m=1.268e5)
[1] 63403.48
- 169 名前：１３２人目の素数さん [2020/03/28(土) 19:42:56.69 ID:NK6wIjWT.net]
- >>160
なんだってー。直感に反するな
- 170 名前：１３２人目の素数さん mailto:sage [2020/03/29(日) 09:23:06.20 ID:WogCQeQk.net]
- >>161
総人口100人として有名人の数を1～100人まで変化させて、有名人に感染者がいたときの100人中の感染者の数をグラフにすると
https://i.imgur.com/SMFnNNl.png
有名人の数を変化さえても期待値にさほどの変化はない。
- 171 名前：１３２人目の素数さん [2020/03/29(日) 10:18:20.74 ID:2PsxdXJm.net]
- >>162
感染者が1名以上という条件だと、
有名人の割合が一定以上になると飽和するんだな。
- 172 名前：１３２人目の素数さん mailto:sage [2020/03/29(日) 10:39:47.55 ID:WogCQeQk.net]
- Ax: x人の感染者がいる(x=0～n)という事象
B:最低一人の感染陽性判定という事象
Pr[Ax|B]=Pr[B|Ax]Pr[Ax]/Pr[B]
Pr[Ax]:事前確率
Pr[B|Ax]:尤度
Pr[B]:周辺尤度(規格化定数)
求めたい期待値Eは
Σ(x*Pr[Ax|B])/ΣPr[Ax|B] = Σ(x*Pr[B|Ax]Pr[Ax])/Σ(Pr[B|Ax]Pr[Ax])
Pr[Ax]がxにかかわらず定数であれば
E=Σ(x*Pr[B|Ax])/Σ(Pr[B|Ax])
事前確率分布を一様分布と仮定しての計算
つまり、感染者が1人の確率も50人の確率も100人の確率,....も一定という前提での計算。
- 173 名前：１３２人目の素数さん mailto:sage [2020/03/29(日) 10:47:28.57 ID:WogCQeQk.net]
- >>163
そうみたいですね。
> data.frame(有名人=1:10,期待値=sapply(1:10,function(x) fn(100,x)$mean))
有名人 期待値
1 1 67.00000
2 2 62.75000
3 3 60.20000
4 4 58.50000
5 5 57.28571
6 6 56.37500
7 7 55.66667
8 8 55.10000
9 9 54.63636
10 10 54.25000
> data.frame(有名人=1:10*10,期待値=sapply(1:10*10,function(x) fn(100,x)$mean))
有名人 期待値
1 10 54.25000
2 20 52.31818
3 30 51.59375
4 40 51.21429
5 50 50.98077
6 60 50.82258
7 70 50.70833
8 80 50.62195
9 90 50.55435
10 100 50.50000
- 174 名前：１３２人目の素数さん [2020/03/29(日) 10:58:30.70 ID:1Oo79tY3.net]
- 「有名人」を「wikに載ってる人」と定義し
その数を１００００人としてそのうち４人（志村、藤浪、長坂、伊藤隼人）
感染したとしても結果は変わらない
- 175 名前：１３２人目の素数さん mailto:sage [2020/03/29(日) 10:58:36.48 ID:WogCQeQk.net]
- 昨日の東京のコロナ陽性者は87人検査して63人陽性であったという。
検査の感度0.6 特異度0.9と仮定して、87人中に感染者は何人と推定されるか？
真陽性率＝感度=0.6
偽陽性率＝１－特異度=0.1
87人中の感染者数をxとすると
陽性者数= 感染者数*真陽性率 + 非感染者数*偽陽性率
63=x*0.6+(87-x)*0.1
これを解くとあり得ない答になる。
- 176 名前：１３２人目の素数さん mailto:sage [2020/03/29(日) 11:48:31.42 ID:WogCQeQk.net]
- >>166
総人口ｎ人、有名人ｍ人、そのうち感染者ｋ人とすると
ｎ人中の感染者の期待値は
x = 0 ～ nとして 、xCkはｘ人からｋ人選ぶ組み合わせの数を表す
Σ(x*(xCk/nCm))/Σ(xCk/nCm) = =Σ(x*(xCk))/Σ(xCk)
となるのでｍの値には依存しない。
ｎ
- 177 名前：１３２人目の素数さん [2020/03/29(日) 14:27:34.63 ID:2PsxdXJm.net]
- >>168
するとこの計算で出てくる推定感染者数6万人って値は意味ない感じですか？
- 178 名前：１３２人目の素数さん mailto:sage [2020/03/29(日) 14:33:09.96 ID:WogCQeQk.net]
- >>167
陽性者数が87人中63人になるような感度と特異度を最小二乗法で求めると。
> (opt=optim(c(0.6,0.9,63),nazo,method='CG'))
$par
[1] 0.916014625 0.779617519 63.002729987
- 179 名前：１３２人目の素数さん [2020/03/29(日) 14:48:03.55 ID:0jXKnAa1.net]
- 学術の巨大掲示板群 - アルファ・ラボ
ttp://x0000.net
数学 物理学 化学 生物学 天文学 地理地学
IT 電子 工学 言語学 国語 方言 など
- 180 名前：１３２人目の素数さん mailto:sage [2020/03/29(日) 15:31:10 ID:WogCQeQk.net]
- >>170
初期値に依存するから意味のないスクリプトであると判明したので撤回します。
- 181 名前：１３２人目の素数さん mailto:sage [2020/03/29(日) 15:31:33 ID:WogCQeQk.net]
- >>169
単なる数字の遊びだろうね。
- 182 名前：１３２人目の素数さん mailto:sage [2020/03/29(日) 15:37:58 ID:WogCQeQk.net]
- >>169
前提となっているのが、
日本人1億2680万人いるとして
日本人の感染者数が１人である確率も1億人である確率も同じと、一様分布を仮定しているのが現実離れしている。
よって現実的には意味がない。
- 183 名前：１３２人目の素数さん mailto:sage [2020/03/31(火) 03:21:38.60 ID:5/cy/U/F.net]
- https://youtu.be/WUMN_71p3Js?t=56
専門家会議がモデルを出したから議論してくれ
- 184 名前：１３２人目の素数さん mailto:sage [2020/03/31(火) 06:08:43.61 ID:2llZ2I8j.net]
- >>175
Reed Frost モデルかな？
何を使ったかには言及がなかった。
- 185 名前：１３２人目の素数さん mailto:sage [2020/03/31(火) 06:12:02.74 ID:2llZ2I8j.net]
- Reed -Frostはパラメータが１個ですむから推定しやすいんだろう。
- 186 名前：１３２人目の素数さん mailto:sage [2020/03/31(火) 08:54:47.69 ID:2llZ2I8j.net]
- >>76
54119人という値になった。
計算プログラムは以下の通り。
# width of 99% confidence interval when 1000 subjects are examined
p2w <- function(
prevalence,
subjects=1000,
sensitivity=0.6,
specificity=0.9,
conf.level=0.99){
# prevalence -> width of 99% confidence interval
n=subjects
p=prevalence*sensitivity+(1-prevalence)*(1-specificity) # positive rate=prev*TP+(1-prev)*FP
q=1-p
2*qnorm(1-(1-conf.level))*sqrt(p*q/n) # width of 99%CI
}
p2w=Vectorize(p2w)
prevalence=seq
- 187 名前：(0,1,by=0.01)
plot(prevalence,p2w(prevalence),bty='l',type='l',lwd=2,ylab='99%CI width',
main='subjects:1000\nsensitivity:0.6\nspecificity:0.9')
optimize(p2w,c(0,1),maximum=TRUE)
#
sj2w <- function(subjects){ # subjects -> maximum 99%CI width & its prevalence
optimize(function(prev) p2w(prev,subjects),c(0,1),maximum = TRUE)
}
# at how many subjects 99%ci width equals 0.01
uniroot(function(x,u0=0.01) sj2w(x)$objective-u0,c(1000,100000)) [] - [ここ壊れてます]
- 188 名前：１３２人目の素数さん [2020/03/31(火) 09:55:37.96 ID:cpD4Fk2x.net]
- 上って、灘校東大理IIIの超秀才のはずなのに、なんで
あんなに頭の悪い発言ばかりしてんの？
変な宗教にでも取り憑かれて理性が狂わされてるのかな？
- 189 名前：１３２人目の素数さん mailto:sage [2020/03/31(火) 10:07:35.24 ID:2llZ2I8j.net]
- 日本人1億2680万人からX人を無作為に抽出してPCR検査して、感染者数(≠検査陽性者数)を信頼区間99%誤差±1%で検定したい。
PCR検査は感度0.6，特異度0.9とする。
何人を抽出すれば十分といえるか？
54000人程度になったけど、あってる？
- 190 名前：１３２人目の素数さん mailto:sage [2020/03/31(火) 14:43:06 ID:2llZ2I8j.net]
- >>179
超秀才は理Iに行くんじゃないの？
- 191 名前：１３２人目の素数さん mailto:sage [2020/03/31(火) 14:50:29 ID:ncBHjUEo.net]
- >>180
感染率の程度、感度・特異度の値の精度の言及無しに出された結論に、ほとんど説得力は無い。
- 192 名前：１３２人目の素数さん mailto:sage [2020/03/31(火) 15:19:09 ID:2llZ2I8j.net]
- >>182
感度 beta(13.6991,9.4661)でmode 0.6 sd=0.1
特異 beta(36.172,4.908) でmode 0.9 sd=0.05
でベイズの階層モデルを組んでみるかな。
- 193 名前：１３２人目の素数さん mailto:sage [2020/03/31(火) 15:45:31.45 ID:2llZ2I8j.net]
- >>183
そのβ分布を弱情報事前分布に設定して、乱数発生させて計算すると
54000人で99%信頼区間の幅の分布は
> summary(s2w(54000))
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.008144 0.009912 0.009981 0.009927 0.010005 0.010011
となるから、まあ、概ねあっていると思うな。
- 194 名前：１３２人目の素数さん mailto:sage [2020/03/31(火) 17:50:11.35 ID:ncBHjUEo.net]
- 最も重要なファクターは事前感染率。
事前感染率はいくらに設定したの？
偽陽性が調査対象の10パーセント程含まれる。
医者が怪しいと判断した場合にのみ検査をする場合は、真陽性が調査対象の数十パーセントが期待できる。
このような場合は、真陽性は偽陽性より多数であることが期待でき、検査対象の正確な感染率は把握できるが、
「日本人1億2680万人からX人を無作為に抽出」のような方法だと、感染率0.01%(←現在確認できている感染者の
7倍程度が実際の感染者数に相当)辺りが妥当だと思われるが、この場合、五万人調査して、真陽性5人、偽陽性5000人
のような数字が出てくる。感染率0.02%だったとすると、真陽性10人、偽陽性5000人だ。
中央値のみで判断すると、例えば、5005人の陽性が出ると、0.01%で、5010人の陽性者が出ると0.02% のような
データが出てくる。誤差との見極めは困難。
このような数字から、信頼できる感染率が出せるのか？
- 195 名前：１３２人目の素数さん mailto:sage [2020/04/01(水) 07:44:43.76 ID:xwYPMdxl.net]
- >>185
一様分布
- 196 名前：１３２人目の素数さん mailto:sage [2020/04/01(水) 07:48:29.51 ID:xwYPMdxl.net]
- 確率の分布を考えずにスポットで考える思考のやつとは議論にならんな。
ベイズ階層モデルやったことないの？
- 197 名前：１３２人目の素数さん mailto:sage [2020/04/01(水) 09:12:32 ID:bZbNlxPT.net]
- 0%～100% までの一様分布のようだな。
つまり、事前確率全く不明だから、1/2教の経典に従い、0.5=50%でやったということ。
医者が検査を行った方がよいと判断した集団でも、なかなか有病率50%はいかない。
そのような結果は、無作為抽出で必要なの調査人数はどれくらいか等という議論では使えない。
全住民を対象にした無作為抽出なら、十
- 198 名前：万人に一人 以上いる（いた）のは確実だった一方、
百人に一人 という程たくさんはいないだろう と見積もれる。0.001%～1% 辺りで行うべき。
ちょっと考えれば判ることを指摘しているに過ぎない。
調査対象の有病率0.01以下の集団に対し、特異度90%の性能の機器で調査しても、ほとんどがエラー。
せめて 有病率 は、 1-特異度 と同じオーダーか、1-特異度 より大きくないと、扱えない。
特異度99.99%の機器を用意するか、でなければ、有病率を１０パーセント程度以上に煮詰めてからやれというお話。 [] - [ここ壊れてます]
- 199 名前：１３２人目の素数さん mailto:sage [2020/04/01(水) 09:19:12 ID:deMoC1lt.net]
- >>188
東京都の行政検査では陽性率が50%を越える日があるぞ。
- 200 名前：１３２人目の素数さん mailto:sage [2020/04/01(水) 09:26:31 ID:deMoC1lt.net]
- 有病率の事前分布を一様分布として
日々の陽性数は二項分布に従うとして
オリンピック延期決定後の検査を受けた集団での有病率をMCMC出だすと
（感度特異度は既述のβ分布を仮定）
> subjects=c(74,95,87,143,244,330)
> positives=c(17,41,47,40,63,68)
> PCRs3(subjects,positives,iter=10000,warmup=1000)
mean lower upper
0.37288732 0.09822213 0.63719043
- 201 名前：１３２人目の素数さん mailto:sage [2020/04/01(水) 09:31:13 ID:deMoC1lt.net]
- >>188
別に有病率を(0,0.1)の一様分布にしても計算できるけど
都の行政検査も陽性率が50%を越える日もあったから一様分布でいいと思うね。行政検査に回った集団の話だけど。
感度・特異度も弱情報事前分布が設定できる。
- 202 名前：１３２人目の素数さん mailto:sage [2020/04/01(水) 09:33:18 ID:deMoC1lt.net]
- 一変数のポイント確率しか計算できない奴との議論は不毛だね。
- 203 名前：１３２人目の素数さん mailto:sage [2020/04/01(水) 09:43:53 ID:HHJL1yTu.net]
- 結局なんの疫学データにも基づかない、疫学データで追試することもできない、なんの理論的根拠もない統計仮説下のお話なんて統計学、疫学できないな意味なんかないんだよな。
計算機で遊んでる以上の意味なんかない。
- 204 名前：１３２人目の素数さん mailto:sage [2020/04/01(水) 09:46:36 ID:bZbNlxPT.net]
- >>189
だからきちんと「なかなかいかない」と書きました。
>>191
目的が「日本人1億2680万人からX人を無作為に抽出してPCR検査して、感染者数(≠検査陽性者数)
を信頼区間99%誤差±1%で検定したい。 」なのだから、あなたの主張は前提を無視ししている。
- 205 名前：１３２人目の素数さん mailto:sage [2020/04/01(水) 12:55:21.15 ID:xwYPMdxl.net]
- >>193
計算機で遊ぶこともできずに電卓で計算して必死で書いていて虚しくない？
CTの診断能を検討した論文。
https://doi.org/10.1148/radiol.2020200823
誰でも鑑別できるのか疑問に思った
このペーパのTable 3に3人の読影医の結果が載っている。
TP FP TN FN sen spc PPV NPV accuracy
1 158 13 192 61 0.72 0.94 0.92 0.76 0.83
2 157 24 181 62 0.72 0.88 0.87 0.74 0.80
3 206 156 49 13 0.94 0.24 0.57 0.79 0.60
陽性尤度比、陰性尤度比、Diagnostic Odd Ratio(陽性尤度比/陰性尤度比)を計算して加えると
TP FP TN FN sen spc PPV NPV acc PLR NLR DOR
1 158 13 192 61 0.72 0.94 0.92 0.76 0.83 11.4 0.30 38
2 157 24 181 62 0.72 0.88 0.87 0.74 0.80 6.1 0.32 19
3 206 156 49 13 0.94 0.24 0.57 0.79 0.60 1.2 0.25 5
PPV,accuracy,DORから読影医3が劣っているようにみえる。
PPVで三者を検定してみる。多重比較になるので一番厳しいBonferri法で補正
Pairwise comparisons using Pairwise comparison of proportions
data: TP out of TP + FP
1 2
2 0.4 -
3 1e-15 2e-11
明らかに３が劣っている。
- 206 名前：１３２人目の素数さん mailto:sage [2020/04/01(水) 12:55:59.40 ID:xwYPMdxl.net]
- 読影医1,2を加算して計算すると
感度72% [67-76]
特異度91% [88-94]
という結果が得られた。
しかし、現実には何でもコロナと診断する傾向のある読影医3も紛れこむからこういう読影医も加算して計算しないと現実的でないね。
問題
3人を統合したときの感度・特異度とその95%信頼区間を述べよ。
- 207 名前：１３２人目の素数さん mailto:sage [2020/04/01(水) 12:59:59.63 ID:YULTPcko.net]
- 昔パソコンは習うより慣れろ、理屈なんかわからなくても使ってたらわかるってのがあったけど、まさに正反対の方向にダメダメだな。
学問に対するなんの畏敬の念もない。
- 208 名前：１３２人目の素数さん mailto:sage [2020/04/01(水) 13:08:31.65 ID:xwYPMdxl.net]
- >>197
＞学問に対するなんの畏敬の念
ひょっとしてアホなの？
- 209 名前：１３２人目の素数さん mailto:sage [2020/04/01(水) 13:14:52.37 ID:xwYPMdxl.net]
- Housefield数の計算原理がわからなくても
この画像が新型コロナ肺炎かどうか、診断できる方が有用なんだよな。
https://pubs.rsna.org/na101/home/literatum/publisher/rsna/journals/content/radiology/0/radiology.ahead-of-print/radiol.2020200823/20200309/images/medium/radiol.2020200823.fig3d.gif
中心極限定理の証明できなくても、学問への畏敬とかなくても、二項分布を正規分布で近似して計算できる。
- 210 名前：１３２人目の素数さん [2020/04/01(水) 18:34:14.58 ID:zMY/D89k.net]
- >>168 他皆様
有名人の感染者が増えてきましたが
市中感染率に影響はないという県警でよろしいのでしょうか？
- 211 名前：１３２人目の素数さん mailto:sage [2020/04/02(木) 06:13:09 ID:+vJJzaTC.net]
- >>200
サンプルサイズは期待値の信頼区間幅に影響するけど期待値そのものに影響しないってことでは？
- 212 名前：１３２人目の素数さん mailto:sage [2020/04/02(木) 09:26:20.10 ID:mzm7EAoV.net]
- 市中感染率が増加の時はもちろんそうだが、一定、あるいは、減少傾向であっても、
経過日数が多くなれば、感染者数は多くなる。
例えば、十日に一人有名人の感染が報告されるというのが継続されていたなら、感染率は一定と
考えられるが、それが、一週間に一人 → 五日に一人 → 三日に一人 → ほぼ毎日 →．．．
のように、報告されるペースに変化があると、感染率も変化していると考えられる。
- 213 名前：１３２人目の素数さん mailto:sage [2020/04/02(木) 09:30:08.05 ID:mzm7EAoV.net]
- 補足だが、あまりにも、有名人感染の報告頻度が多くなると、ニュースとしての価値が低くなり、
以前であったら報告されていたであろうケースが報告されなくなるということもあるので、
その辺も考慮して考える必要はある。
- 214 名前：１３２人目の素数さん mailto:sage [2020/04/03(金) 11:52:46 ID:cch/ocoF.net]
- 横浜市立大学データサイエンス学部佐藤彰洋教授のCOVID-19（新型肺炎）の感染拡大抑止に関する研究・検討資料内容を共有するページ
https://www.fttsus.jp/covinfo/#Tokyo
矢原 徹一：九州大学理学研究院教授の試算
https://jbpress.ismedia.jp/articles/-/59942
- 215 名前：１３２人目の素数さん mailto:sage [2020/04/04(土) 11:37:34 ID:ZFu90Xbq.net]
- SEIR MODEL
dS(t)/dt = mu*(N-S) - b*S(t)*I(t)/N - nu*S(t)
dE(t)/dt = b*S(t)I(t)/N - (mu+sig)*E(t)
dI(t)/dt = sig*E(t) - (mu+g)*I(t)
dR(t)/dt = g*I(t) - mu*R + nu*S(t)
mu:自然死亡率 b:感染率(S->I)
nu:ワクチン有効率(S->R) sig:発症率(E->I),g:回復率(I->R)
の微分方程式の数値解を使ってシミュレーション
対策しない（外出を控えず、マスクもしない）方が患者や死者は増えるけど早く収束するな。
contacｔ_rate と trannsmission_probabilityを変化させてグラフにしてみた。
https://i.imgur.com/6OgJkDb.png
- 216 名前：１３２人目の素数さん mailto:sage [2020/04/04(土) 15:28:24 ID:zerwqPau.net]
- 一次�
- 217 名前：Y業ごと消滅していいならそうかもな []
- [ここ壊れてます]
- 218 名前：１３２人目の素数さん mailto:sage [2020/04/05(日) 09:54:53.42 ID:fV/kgtmE.net]
- オリンピック延期決定以後の東京都の行政PCR検査での陽性率をグラフにすると
https://i.imgur.com/XMKErLV.png
（陽性数より検査件数の公表は２～３日遅れる）
PCR検査は感度60％、特異度90％くらいなので検査を受けた集団の有病率はもっと多いはず。
感度(最頻値0.6 標準偏差0.1)、特異度（最頻値0.9 標準偏差0.05)のベータ分布に設定、有病率は(0,1)の一様分布でMCMCしたみた。
https://i.imgur.com/1ZkOk33.png
有病率40％くらいありそうだな。
- 219 名前：１３２人目の素数さん mailto:sage [2020/04/05(日) 23:57:53 ID:fV/kgtmE.net]
- 新型コロナ肺炎に再感染があるとして流行具合をシミュレーションしてみた。
赤が感染者
上：再感染率０％
中：再感染率１％
下：再感染１％に治癒確率を5倍にする治療薬がある場合
https://i.imgur.com/nhCe0aZ.png
|
 |
