数学 統計に詳しい人が語るコロナウイルス

1 名前：１３２人目の素数さん [2020/02/29(土) 02:18:41 ID:twdO677Q.net] 東大数学科卒の元官僚はこう分析してるが、お前らはどうなると思う？ www.zakzak.co.jp/soc/news/200220/dom2002200003-n2.html 中国国外感染者の中国国内との比率をみると、 １月２０日の数字公表以降は、０・８〜２・６％で比較的安定している。 これは、新型肺炎の感染者のほとんどは中国国内、それも湖北省に集中しているからだ。 ちなみに中国国外での感染者数は、中国国内の１・１％だ（２月１６日現在）。 本コラムで紹介したが、現時点では、最終的な中国国内の感染者数は２０万人超と筆者は推計している。 となると、中国国外の感染者は数千人程度になるだろう。 中国国外のうち日本の比率は１割弱なので、日本の感染者数は数百人程度であろう。 その場合、死者も数人から１０人程度になるだろう。 こうした推計をすると、今の感染者は氷山の一角だと思われるが、今後の増加ペースはどうなるだろうか。 新型コロナウイルスの検査は簡単に行えるので、今後、日本での感染者数は増えていくだろう。 ある時点ではそれがネズミ算的に増えるかのように思える局面もあるだろうが、 筆者の推計が正しければ、現時点ではせいぜい数百人が一つのメドだ。 552 名前：１３２人目の素数さん mailto:sage [2020/05/16(土) 08:36:22 ID:/d9XIHEO.net] >>520 アヒル本の著者だね。俺も読んだ。 553 名前：１３２人目の素数さん [2020/05/16(土) 12:36:46 ID:cGFgpNwX.net] 4/2の時点で感染者数を６０００くらいと見積もってるね。＞アヒル本の人 実数の三倍程度。いい線かもしれない。 554 名前：１３２人目の素数さん [2020/05/16(土) 13:39:32 ID:XSllT5Kn.net] https://github.com/contactmodel/COVID19-Japan-Reff/blob/master/nishiura_Rt%E4%BC%9A%E8%AD%B0_12May2020.pdf 555 名前：１３２人目の素数さん mailto:sage [2020/05/16(土) 14:29:29 ID:BxGcLzV+.net] int<lower = K> upper_bound; ↓ int　upper_bound; にしてもエラーがでる。 stan_dataで upper_bound = 147 にすると動くけど、何をやってんのか自分でもよくわからん。 556 名前：１３２人目の素数さん mailto:sage [2020/05/16(土) 16:47:14 ID:UBjBTXVe.net] アヒル本ってなんですか？ 557 名前：１３２人目の素数さん mailto:sage [2020/05/16(土) 17:55:00 ID:BxGcLzV+.net] >>525 名著として名高いStanの入門書 StanとRでベイズ統計モデリング (Wonderful R) (日本語) 単行本 ? 2016/10/25 表紙の色がアヒルのくちばしの色に似ているかららしい。 558 名前：１３２人目の素数さん mailto:sage [2020/05/16(土) 17:56:04 ID:BxGcLzV+.net] >>524 自己解決 select ↓ dplyr::select でちゃんと動作した 559 名前：１３２人目の素数さん mailto:sage [2020/05/16(土) 19:16:41 ID:5/7mBSAW.net] >>526 thx 560 名前：１３２人目の素数さん mailto:sage [2020/05/16(土) 19:28:52 ID:BxGcLzV+.net] R に　EpiEstimというパッケージがあって、再生産数を算出する関数が搭載されている。 結局、infecterとinfecteeが発症するまでの期間serial intervalの分布をどう設定するかで結果が変わるみたいだなぁ。 Rのヘルプファイルを解読中。 Rのヘルプファイルは不親切設計で有名（理解できている人の備忘録みたいな性格だから）。 561 名前：１３２人目の素数さん mailto:sage [2020/05/16(土) 20:49:50 ID:5T7nYzW3.net] ここに居る人達って何者なの 学部の知識超えてるよね 統計でご飯食べてる人たち？ 562 名前：１３２人目の素数さん [2020/05/16(土) 21:03:16 ID:YPW1s8+S.net] しかし、なんでも揃ってるなRって 563 名前：１３２人目の素数さん mailto:sage [2020/05/16(土) 22:57:43 ID:BxGcLzV+.net] >>529 Stanでの西浦モデルではinfecterとinfecteeが発症するまでの期間serial intervalの分布に ## Serial interval [Nishiura et al 2020 - only certain cases] param1_SI = 2.305, param2_SI = 5.452, // serial interval vector[K] gt = pweibull(param1_SI, param2_SI, K); として使われているので、平均値などを出してみた。 乱数発生と理論値 > x=rweibull(1e5,param1_SI,param2_SI) > mean(x) ; param2_SI*gamma(1+1/param1_SI) [1] 4.829273 [1] 4.830129 > var(x) ; param2_SI^2*(gamma(1+2/param1_SI)-(gamma(1+1/param1_SI))^2) [1] 4.907755 [1] 4.940682 > median(x) ; param2_SI*(log(2)^(1/param1_SI)) [1] 4.655777 [1] 4.6505 > density(x)$x[which. 564 名前：max(density(x)$y)] ; param2_SI*(1-1/param1_SI)^(1/param1_SI) [1] 4.116837 [1] 4.259624 > optimise(function(x) dweibull(x,param1_SI,param2_SI),c(0,10),maximum = T)$max [1] 4.259623 グラフにしてみた。 https://i.imgur.com/9vvCJuZ.png 正規分布で近似してもよさそうな感じだな。 [] [ここ壊れてます] 565 名前：１３２人目の素数さん mailto:sage [2020/05/17(日) 00:04:44 ID:u5AEq3c8.net] >>532 ワイブル分布の平均 4.830129と標準偏差2.222765をそのまま正規分布のパラメータに使って、グラフを重ねてみる。 https://i.imgur.com/TnzGwWx.png ワイブル分布で発生させた乱数をワイブルでフィットさせてAICを出してみた Goodness-of-fit criteria 1-mle-weibull Akaike's Information Criterion 438377.2 Bayesian Information Criterion 438396.2 ワイブル分布で発生させた乱数を正規分布でフィットさせてAICを出してみた。 Goodness-of-fit criteria 1-mle-norm Akaike's Information Criterion 444280.9 Bayesian Information Criterion 444299.9 まぁ、許容範囲。 これで、 library(EpiEstim)の例にある、　mean_si std_siが求まった ## Estimate R with assumptions on serial interval res <- estimate_R(incid, method = "parametric_si", config = make_config(list( mean_si = 4.83, std_si = 2.22))) domestic ,　imported, unobserved の分類がよくわからんが、全部足してグラフを描いてみた。 https://i.imgur.com/rKBeWgq.png 566 名前：１３２人目の素数さん mailto:sage [2020/05/17(日) 00:18:59 ID:u5AEq3c8.net] 別の論文だと対数正規分布がフィットすると西浦氏は記載している。 serila interval : infector と infecteeの発症間隔 https://www.ijidonline.com/article/S1201-9712(20)30119-3/pdf その分布が平均4.7　標準偏差2.9の対数正規分布が最もフィットするのはいいんだが、 その分布を与えるパラメータの記述がほしい。 最小二乗法で求めてみた。　 $par [1] 1.3862713 0.5679836 ワイブル分布にも似るとか書いてあるがパラメータ記載なし この対数正規分布をワイブル分布で近似してみる。 Fitting of the distribution ' weibull ' by maximum likelihood Parameters: estimate Std. Error shape 1.757488 0.00392072 scale 5.316986 0.01014664 https://i.imgur.com/Uzg6u84.png 点線が2項分布で実線がワイブル分布 567 名前：１３２人目の素数さん [2020/05/17(日) 07:09:05 ID:/ApxcPC4.net] “頭脳王”東大生・河野玄斗 基本的な数学でコロナウイルス検査を全員にしても意味がないことを証明してみた https://www.youtube.com/watch?v=jMIScCb04qs 568 名前：１３２人目の素数さん mailto:sage [2020/05/17(日) 07:28:25 ID:mpH4uuPx.net] 実効再生産数を感染者数の推移から推定する数理的原理をキチンと解説した本、または論文誰か知りませんか？ 勉強してみたい。 569 名前：１３２人目の素数さん mailto:sage [2020/05/17(日) 07:33:35 ID:u5AEq3c8.net] >>536 RのEpiEstimの著書の論文なんかどうでしょう？ A New Framework and Software to Estimate Time-Varying Reproduction Numbers During Epidemics https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3816335/ 570 名前：１３２人目の素数さん mailto:sage [2020/05/17(日) 07:55:41 ID:u5AEq3c8.net] Rで引数なしで関数を実行させようとすると、ソースが表示される。 その関数が呼び出した関数のソースを次々に辿っていけば内部計算がわかるので、そこから原理がつかめるかも（俺には無理なのでしないけど） 内部関数のときは:が3つとか、パッケージ名:::関数、パッケージ名:::関数.default（例、ｔ検定のソース表示は　stats:::t.test.default　）で表示される。 EpiEstim::estimate_R EpiEstim:::process_si_data EpiEstim:::process_config_si_from_data coarseDataTools::dic.fit.mcmc 571 名前：１３２人目の素数さん mailto:sage [2020/05/17(日) 08:10:41 ID:u5AEq3c8.net] >>529 感染者に０が続くと再生産数の信頼区間幅がどんどん広くなってくる。 まあ、疫病用のソフトウェアと理解しておこう。 https://i.imgur.com/QbwNydN.png infected=c(0,1,1,1,0,0,0,2,0,3,2,3,1,1,1,1,3,0,1,0,0,0,0,0,1,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0) ## Estimate R with assumptions on serial interval res <- estimate_R(infected, method = "parametric_si", config = make_config(list( mean_si = 4.83, std_si = 2.22))) 572 名前：１３２人目の素数さん mailto:sage [2020/05/17(日) 09:00:46 ID:u5AEq3c8.net] >>539 ちゃんと、記載されていた　：P the precision of these estimates depends directly on the number of incident cases in the time window [t ? τ + 1; t]. This allows us to control the precision by adjusting the window size. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3816335/pdf/kwt133.pdf 573 名前：１３２人目の素数さん mailto:sage [2020/05/17(日) 09:34:29 ID:u5AEq3c8.net] 基本コンセプトはこれだろうな。 Therefore, in practice, we apply our method to data consisting of daily counts of onset of symptoms where the infectivity profile ws is approximated by the distribution of the serial interval. 公表された感染者数で計算するために発症から診断/公表までの時間も考慮した点が西浦モデルの優れた点ではないかと感じている。 574 名前：１３２人目の素数さん mailto:sage [2020/05/17(日) 12:07:27 ID:S4iLXC97.net] >>537 ありがとう。 読んでみます。 でも"a new frame work"っていうならもっと大本のこの道の研究者なら「Rtの推定はコレ」みたいな標準的なものがあるんですかね？ どなたかご存知ないですか？ 575 名前：１３２人目の素数さん mailto:sage [2020/05/17(日) 12:11:01 ID:S4iLXC97.net] >>540 おぉ、こっちがRで採用されてる推定法なんですね。 読んでみます。 576 名前：１３２人目の素数さん mailto:sage [2020/05/17(日) 12:16:27 ID:S4iLXC97.net] と思ったら>>537と>>540は同じかorz 577 名前：１３２人目の素数さん mailto:sage [2020/05/17(日) 15:31:13 ID:u5AEq3c8.net] 結局、これが核心部分 Rt =I t (number of new infections generated at time t) / Σ[s=1,t] I t-s * Ws ( = total infectiousness of infected individuals at time t) Ws : an infectivity profile given by a probability distribution ws, dependent on time since infection of the case, s, but independent of calendar time, t. E[I t] = Rt Σ[s=1,t] I t-s * Ws Σ[s=1,t] I t-s * Wsの部分が畳み込み積分で Ws ∝　serial interval　　ガンマ分布で近似するのが定石らしい。 578 名前：１３２人目の素数さん mailto:sage [2020/05/17(日) 20:41:31 ID:/XIaXEJI.net] >>545 それはそのモデルでのRtの定義ですね、 ではなく例えば4/1〜4/30までのデータx1,x2,‥,x30までが与えられた時、各日付のCIをどう定義してるのかがわからないんです。 統計量Xが一つあってその観測値xが一つある時そのレベルλのレベルλのCIがIであるとは P(X<x|θ)>(1-λ)/2 & P(X>x|θ)>(1-λ)/2 ⇔ θ∈I (θがIに入らないときはxが小さすぎてそんな観測値が得られる確率が(1-λ)/2以下になるか、xが大きすぎてそんな事が起こる確率が(1-λ)/2以下になってほとんど起こり得ない) となりますが、統計量が複数になるとこの大きすぎて、小さすぎてと二つのハズレ領域を考えるだけでは済まなくなります。 ハズレ領域の設定の仕方は色々考えられるけど統計の問題なので自分で俺様流の領域を設定するわけにもいかないのでpublicにはどう処理してるのだろうと。 しかし疫学の教科書わざわざ買うほどには興味もないしw 知っ 579 名前：てる人いたら教えてもらおうと。 まだ論文読んでないので書いてあるかもですけど。 [] [ここ壊れてます] 580 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 05:56:32 ID:hVD/4gbI.net] >>546 >523氏が挙げてくれた 実効再生産数とその周辺　に　記述があったが、publicと呼べるのかどうかは門外漢なので知らない。 https://github.com/contactmodel/COVID19-Japan-Reff/blob/master/nishiura_Rt%E4%BC%9A%E8%AD%B0_12May2020.pdf こんな記述があった >> Several non mathematical definitions. A: あるカレンダー時刻 t で起こっている 2 次感染数の 1 人あたり平均値 B: あるカレンダー時刻 t で感染した 1 人がその後の経過で生み出す 2 次感染者数の平均値 C: 罹患率 有病割合比などから推定される 1 人あたりの 2 次感染者数（ actual reproduction number とか） D: 予防接種など流行対策下での 1 人当たりが生み出す 2次感染者数 << と定義は一義的ではないみたいだが、 西浦モデルでのR1(t) R2(t)は　>540の論文では　Rt　Rct　（ｃはcohortの頭文字)として言及されているから、まあ、理論疫学者の間ではcommonな定義なんだろうと思う。 581 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 06:16:59 ID:hVD/4gbI.net] >>546 >4/1〜4/30までのデータx1,x2,‥,x30 x1,...,x30が感染者数なら非負整数で与えられるから、CIを考える必要があるかな？ 集計ミスで19人が117人とかあったらしいから、信頼区間を考える必要があるのかもしれないとは思うけど。 確かに、発症日に関してはファジーな部分があるとは思う。 いつから熱がありましたか？味がわからなくなりましたか？と問われても１〜３日位の幅はでるだろう。 RのEpiEstim::estimate_Rをヘルプファイルの実例を使って走らせてみた。 serial intervalの分布をデータから推定させるのに 　infecterの発症日の下限・上限 　infecteeの発症日の下限・上限 を設定する項目が(EL,ER,SL,SR)があった。 このデータに合致する分布関数(ガンマ、ワイブル、対数正規分布が指定可能）のパラメータを算出して計算させているみたい。 西浦モデルではワイブル分布で固定。 潜伏期間にも変動があるから、誰がinfecterで誰がinfecteeかを決定するのも難しいだろうとは思う。 後から発症した人間が感染源というのもありうるし。 582 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 11:15:46 ID:iMUOaQs/.net] >>548 統計データから真の罹患日を推定する方法もあるようですが そこではないんです。 しりたいのはCIのハズレ領域の設定です。 1変数の場合、母数θに対して分布Fθが定まっている場合、レベルλに対して[0,1]の部分集合J(λ)が決まって、観測値xに対する信頼区間I(λ,x)は θ∈I(λ,x)⇔Fθ(x)∈J(λ) を満たす区間として定まります。 上の方でI(λ)が上下対称に取らないのはなぜという話題がありましたが、コレがその理由です。 J(λ)の方は上下の(1-λ)/2を削って((1-λ)/2,1-(1-λ)/2)をとり、上下対称に“ハズレ領域”をとりますが、それをもとに計算されるI(λ,x)は対称とはならないからです。 問題は観測値が2変数以上ある場合“ハズレ領域”をどう設定するものかわからないのです。 私が大学で勉強した時はそこまでやらなかったので。 普通に考えればI^3の中の立方体の体積がλになるように真ん中にとるんだろうなぁと思うんですけど。 583 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 13:21:08 ID:1P7V5xJn.net] 職場でも最初に発症した人が感染源のように扱われるけど 潜伏期間の分布を考えたら断定はできない。 COVID19の潜伏期間の論文 https://www.nejm.org/doi/full/10.1056/NEJMoa2001316 から、 潜伏期間は対数正規分布で近似できてそのパラメータは #--- incubation period --- # from Li et al NEJM 2020 # log 584 名前：normal mean = 5.2 ln.par1 = 1.434065 ln.par2 = 0.6612 という。 ある人物Xが新型コロナ肺炎に罹患したとする。 行動調査によって発症前にキャバクラに行っており接客したキャバ嬢がX発症の２日後に発症していたことがわかった。 Xがキャバ嬢から移された確率を求めよ。 [] [ここ壊れてます] 585 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 13:25:20 ID:1P7V5xJn.net] >>549 Highest Density Probability Intervalを求めればいいんじゃないの？ 586 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 14:04:29 ID:jxOnYm2h.net] >>551 何ですかそれは？ 587 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 14:08:09 ID:1P7V5xJn.net] >>550 正確にはキャバ嬢がXより先に感染していた確率だな。 588 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 14:08:59 ID:1P7V5xJn.net] >>552 確率分布が対称じゃないときの信頼区間 589 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 14:13:34 ID:1P7V5xJn.net] >>554 こんな感じ。https://i.imgur.com/C8jOPlx.jpg 590 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 14:26:48 ID:jxOnYm2h.net] >>554 ただ、それだとそもそもモードに近いとこをやってます。 信頼区間は密度関数を横に切るのではなく両裾を縦に切ってハズレ部分が1-λになるようにするので少しイメージが違うしきがします。 モードなのかメジアンなのかの違いです。 いずれにせよ、こうやればいいという拡張のための俺様ルールを設定するのはいくらでもできますが、統計の話なのでそんな俺様ルールについて語っても意味ありません。 591 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 15:34:14 ID:1P7V5xJn.net] >>556 単峰性の場合、信頼区間幅が最小になるのがHighest Density Interval >550なら HDI > HDInterval::hdi(x)[1:2] lower upper 0.5822687 12.5635525 分位数 > quantile(x,c(0.025,0.975)) 2.5% 97.5% 1.148711 15.334698 HDIの方が幅が小さい。 592 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 15:36:39 ID:jxOnYm2h.net] >>557 ？？？ 593 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 15:44:16.03 ID:jxOnYm2h.net] ああ、わかった。HDIやCIの意味を誤解してませんか？ HDIでググって調べたらコレ↓ですよ。 https://rindalog.blogspot.com/2015/10/hdi-highest-density-interval.html?m=1 594 名前：１３２人目の素数さん [2020/05/18(月) 19:32:28.41 ID:4BbElJo8.net] >>179 単にお勉強ができただけ。 頭が良くないのさ。 自分の頭で物事を考えるってことができない。 595 名前：１３２人目の素数さん [2020/05/18(月) 19:34:11.31 ID:4BbElJo8.net] >>181 その通り。 具体的には理学部の数学科と物理学科。 工学部にも時々もの凄いのがいる。 596 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 19:49:54 ID:1P7V5xJn.net] >>559 いや、Rのパッケージ　HDIntervalで計算しているから誤解していないと思う。 内部の処理コードをみると信頼区間幅が最小になるのを最尤法で出しているね。 597 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 19:51:25 ID:1P7V5xJn.net] >>561 ほんとうに頭のいい人は医学部でなく理学部か工学部に行く。 ほんとうに頭の悪い奴は底辺シリツ医大に行く。 シリツ医大には手先の器用な人もいるが、頭が器用な奴をみたことがない。 598 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 20:05:42 ID:1P7V5xJn.net] >>550(自答) # # 人物ｄが発症してdelay日後に濃厚接触したキャバ嬢ｃが発症 # ｃの感染がｄより先行していた確率は？ rm(list=ls()) stancode= " data{ real onset_delay; real ln_par1; real ln_par2; } parameters{ real <lower=0> d_incubation; real <lower=0> c_incubation; } transformed parameters{ real infection_delay = onset_delay + d_incubation - c_incubation; } model{ d_incubation ~ lognormal(ln_par1,ln_par2); c_incubation ~ lognormal(ln_par1,ln_par2); } " model=stan_model(model_code = stancode) #--- incubation period --- # from Li et al NEJM 2020 # lognormal mean = 5.2 ln_par1 = 1.434065 ln_par2 = 0.6612 fn.stan <- function(delay, print=FALSE, ...){ dataList=list(onset_delay=delay,ln_par1=ln_par1,ln_par2=ln_par2) fit=sampling(model,data=dataList, ...) ms=rstan::extract(fit) if(print) BEST::plotPost(ms$infection_delay,compVal=0,xlab='infection delay') mean(ms$infection_delay < 0) } fn.stan(2,print=T,iter=5000,warmup=1000) onset_delays=0:20 y=sapply(onset_delays,fn.stan) plot(onset_delays,y, ylab='Pr[ Infected Later ])') ２日後の発症だと > fn.stan(2,print=T,iter=5000,warmup=1000) [1] 0.2945 ３割くらいはあとから発症した方が先に感染していた可能性がある。 599 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 20:48:18.33 ID:1P7V5xJn.net] Temporal dynamics in viral shedding and transmissibility of COVID-19 https://www.nature.com/articles/s41591-020-0869-5 のRのコード https://github.com/ehylau/COVID-19/blob/master/Fig1c_Rscript.R と 西浦モデルのコード https://nbviewer.jupyter.org/github/contactmodel/COVID19-Japan-Reff/blob/master/scripts/C.%20Calculating%20the%20Rt%20in%20Stan.ipynb から発症間時間(serial interval)の分布を重ねてみた。 https://i.imgur.com/vrnra5F.png 600 名前：１３２人目の素数さん [2020/05/18(月) 21:01:41.77 ID:4BbElJo8.net] >>563 一方的で申し訳ないが私立大医学部は金持ちのバカ息子が行くイメージ。 601 名前：１３２人目の素数さん [2020/05/18(月) 21:02:26.04 ID:4BbElJo8.net] 西浦さんはさんざん適当なことを言って世論を煽ってどう責任を取るのかな？ 602 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 21:18:52.42 ID:AArRB0Ix.net] >>567 少なくとも西浦氏は算出コードを公開しているだけでも好感が持てる。 603 名前：１３２人目の素数さん mailto:sage [2020/05/18(月) 21:20:19.41 ID:AArRB0Ix.net] >>566 それは正しい認識。 凄いのはド底辺シリツ医の馬鹿さ加減だよ。 裏口バカと呼ばれるのがよくわかる。 imagizer.imageshack.com/img923/2715/RosCsf.jpg i.imgur.com/XBFnEcU.jpg 馬鹿だという自覚がないので救いようがない。 ICU Bookの最終章の冒頭で著者がこう書いている。 In clinical matters, ignorance can be dangerous, but ignorance of ignorance can be fatal. 「叱られないと勉強しない」の対偶を「勉強すると叱られる」 と答えるのはignorance can be dangerousの範疇だが、 ドヤ顔で >対偶をとれば意味が逆になる例文。 というのは、まさに ignorance of ignorance can be fatal. 604 名前：１３２人目の素数さん mailto:sage [2020/05/21(木) 11:39:47 ID:hAZkHjNF.net] 西浦モデルの再生算数の事前分布を変化させてグラフにしてみた。 西浦モデルでのデフォルト https://i.imgur.com/G1wVYgI.png 事前分布を一様分布(0,10)近似の正規分布で近似させた場合 https://i.imgur.com/doS5LEu.png 再生算数の平均０、標準偏差１の場合 https://i.imgur.com/doS5LEu.png 印象としては、西浦モデルは頑強性robustnessのあるモデルとは言えない。 事前分布に大きく影響されるモデルだと思う。 605 名前：１３２人目の素数さん mailto:sage [2020/05/21(木) 11:42:59 ID:hAZkHjNF.net] （url訂正) 西浦モデルの再生算数の事前分布を変化させてグラフにしてみた。 西浦モデルでのデフォルト https://i.imgur.com/G1wVYgI.png 事前分布を一様分布(0,10)近似の正規分布で近似させた場合 https://i.imgur.com/doS5LEu.png 再生算数の平均０、標準偏差１の場合 https://i.imgur.com/0J1RpDa.png 印象としては、西浦モデルは頑強性robustnessのあるモデルとは言えない。 事前分布に大きく影響されるモデルだと思う。 606 名前：１３２人目の素数さん mailto:sage [2020/05/21(木) 11:59:00 ID:hAZkHjNF.net] >>571 誤解されるのは不本意なの追記するけど、ソースを公開する西浦先生の姿勢は高く評価している。 隠蔽・改竄・破棄する安倍とは大違い。 607 名前：１３２人目の素数さん mailto:sage [2020/05/21(木) 12:04:41 ID:RwIzsMl5.net] だからベイズ推定が統計学の世界でメジャーにならんのだろ？ 論理的根拠のない“事前分布”なるもので答えがひょいひょい変わるのでは社会的な影響が大きい防疫政策の決定には使えない。 普通の統計学の検定なら理論的に根拠のある数字、推論しか使わない。 計算は大変だけど。 608 名前：１３２人目の素数さん mailto:sage [2020/05/21(木) 12:36:51 ID:hAZkHjNF.net] >>573 成人の平均身長を１〜２ｍに事前分布にするのは納得できるし、 生まれる子供が男子である 609 名前：確率は０．４〜０．６というのも俺は納得できる。 PCRの感度が３０〜７０％として計算するのも納得できるからその設定で階層ベイズモデルを組むことには異論はないな。 [] [ここ壊れてます] 610 名前：１３２人目の素数さん mailto:sage [2020/05/21(木) 12:43:16 ID:hAZkHjNF.net] こういう問題 あるタクシー会社のタクシーには１から通し番号がふられている。 タクシー会社の規模から保有タクシー台数は１００台以下とわかっている（弱情報事前分布）。 この会社のタクシーを５台みかけた。最大の番号が６０であった。 この会社の保有するタクシー台数の期待値と９５％信用区間(信頼区間)を求めよ。 をベイズで解くときは、 ６０台〜１００台である確率を一様分布として処理しているから これに異論があるのは理解できるけど 日本人成人の平均身長を推定に１〜２ｍを事前分布に想定するのには俺は異論はないね。 一様分布ではなくてガンマ分布にすべきだというのは議論になるとは思うけど。 611 名前：１３２人目の素数さん mailto:sage [2020/05/21(木) 13:00:03 ID:hAZkHjNF.net] >>573 ベイズ信奉者から、ベイズ論者を採用したGAFA(Google等）が成長したと教わった。 迷惑メールのフィルタリングとか雑音の除去とか日常的に役立っているというぞ。 612 名前：１３２人目の素数さん mailto:sage [2020/05/21(木) 13:04:49 ID:hAZkHjNF.net] >普通の統計学の検定なら理論的に根拠のある数字、推論しか使わない。 >計算は大変だけど。 普通の統計学こそ、無理やり既知の分布に当てはめようとするんだよね。 MCMCの方の方が応用が広い。 613 名前：１３２人目の素数さん mailto:sage [2020/05/21(木) 13:12:59.01 ID:hAZkHjNF.net] プロ野球選手の打率は？と問われたら選手次第で異なる、と誰でもわかるのに 確定不能の平均値が存在していると妄想して計算を始めるのが古典主義統計。 つまり、値は存在するけど確定できないという信仰の世界。 昨今の新コロナでいえば、PCRの検査の感度・特異度が一定と考えるのが古典主義。 プロ野球選手の打率と同じでそんなのは場面場面で変化するよ、と考えて計算するのがベイズ。 614 名前：１３２人目の素数さん mailto:sage [2020/05/21(木) 13:16:04.83 ID:EdpRqUGW.net] >>576 もちろん推定の有力な方法であるにせよ、元の仮定に何の根拠もないわけだからそれから得られる結論には論理的根拠はない、ないが、数学的に伝統的な手法で与えるた結論と大差ない事がなんらかの保証があるなら、有力になる。 それが論理的に“大差ない結論が得られる”事が示されてるなら単なる計算手法に過ぎないし、示されていなくても経験的に“よい結果ぎ得られる事が多い”ならそのジャンルではそこそこ信頼するに値するんだろう。 しかしなんらかのモデルでは答えが一意に定まらず、事前分布の選び方により大きく答えが違ってしまう場合があっても不思議はないし、そのような場合ではやはり、“ではどう計算するのが正しいかのか”の論証を待たなければ信頼するのは危険になる。 615 名前：１３２人目の素数さん mailto:sage [2020/05/21(木) 13:19:12.39 ID:hAZkHjNF.net] 古典的（頻度主義）統計信者って、この計算はどうやるんだ？ 俺は乱数発生させて計算できるけど、 そればベイズなのかどうかは知らんが、条件付き確率なのでベイズなんだろうな。 （開業医スレに投稿したけど、回答できるやつは０） COVID19の潜伏期間の論文 https://www.nejm.org/doi/full/10.1056/NEJMoa2001316 結論は #--- incubation period --- # from Li et al NEJM 2020 # lognormal mean = 5.2 ln.par1 = 1.434065 ln.par2 = 0.6612 ある開業医が新型コロナ肺炎に罹患したとする。 行動調査によって発症前にキャバクラに行っており接客したキャバ嬢が開業医発症の２日後に発症していたことがわかった。 キャバ嬢は開業医から移されたと主張して１億円の賠償を求めている。 潜伏期間には幅がありキャバ嬢から移された可能性もあると主張してその確率を計算して賠償金を値切りたい。 いくら値切れるか計算せよ。 616 名前：１３２人目の素数さん mailto:sage [2020/05/21(木) 13:31:57 ID:3Y6HuNza.net] >>580 潜伏期間なんて数学的に決定できるハズないやん？ 数学がやるのは例えば感染日と実際 617 名前：発症した日が確定できるような症例がある程度以上あって、それが従うと病理学的に信頼できる分布の族があって、その中からデータに最も“沿う”分布を選び出すことしかできない。 それは何件も実際の患者のウイルス量の統計をとったり、ウイルスが体内でどのように増えていくのかの病理学的研究データがあって初めて可能になる。 [] [ここ壊れてます] 618 名前：１３２人目の素数さん mailto:sage [2020/05/21(木) 17:20:41 ID:Du+rgnHD.net] >>571 よく理解できていないので質問ですけど 事前分布とは具体的に何の分布ですか？ 基本再生算数の推定値の分布？ 実行再生算数の推定値の分布？ 実行再生算数の事前分布は基本再生算数の分布としたらいいのかなと思いますけど 619 名前：１３２人目の素数さん mailto:sage [2020/05/21(木) 17:26:10 ID:BjgAHZWs.net] まぁこのスレは用語がめちゃくちゃだからなぁ。 620 名前：１３２人目の素数さん mailto:sage [2020/05/21(木) 17:58:56 ID:hAZkHjNF.net] モデル前提での計算できないアホ発見！ 621 名前：１３２人目の素数さん mailto:sage [2020/05/21(木) 18:22:49 ID:hAZkHjNF.net] 潜伏期は病原体と罹患者のパワーバランスで決まるだろうから 定数でなくてばらつきはあると思うね。 622 名前：１３２人目の素数さん [2020/05/21(木) 18:25:39 ID:8+K/mYXQ.net] >>572 >隠蔽・改竄・破棄する安倍とは大違い。 安倍ちゃんがやってるわけじゃないだろ。 指示もされてねーのに、官僚がやってんだよ。 そんなこともわからんようでは、あんたの解析もあてにならんな。 623 名前：１３２人目の素数さん mailto:sage [2020/05/21(木) 18:26:26 ID:hAZkHjNF.net] >>582 西浦のソースだと　Rt　~ normal(2.4 ,2) 624 名前：１３２人目の素数さん mailto:sage [2020/05/21(木) 18:28:21.45 ID:hAZkHjNF.net] >>586 安倍じゃなければ官僚もまともだったんじゃねえの？ 安倍らの集団に組み込まれると東大卒も馬鹿になるようだぞ。 625 名前：１３２人目の素数さん [2020/05/21(木) 18:31:24.22 ID:8+K/mYXQ.net] >>581 その通り。 結局、そういう本質的なデータや理論抜きでは、ベイズ推定やったって 限界があるし、結果の説得力もない。 まあ、適用限界ってものがあるのは何やっても同じだけど。 626 名前：１３２人目の素数さん [2020/05/21(木) 18:32:55.80 ID:8+K/mYXQ.net] >>588 >安倍らの集団に組み込まれると東大卒も馬鹿になるようだぞ。 統計的に確認してみれば？ベイズ推定でもっともらしい結果が出るんじゃね？ｗ 627 名前：１３２人目の素数さん mailto:sage [2020/05/21(木) 18:33:59.08 ID:hAZkHjNF.net] >>590 今は2/2かなｗ 加藤と西村 628 名前：１３２人目の素数さん mailto:sage [2020/05/21(木) 18:40:02 ID:hAZkHjNF.net] >>590 馬鹿になる事後確率分布はこんな結果になりました。 https://i.imgur.com/ZhpWJdo.png 629 名前：１３２人目の素数さん mailto:sage [2020/05/21(木) 18:41:26 ID:hAZkHjNF.net] 95%　CI は　0.5005265 1 630 名前：１３２人目の素数さん mailto:sage [2020/05/21(木) 18:46:18 ID:hAZkHjNF.net] >>592 ちなみに事前分布は一様分布に設定。 631 名前：１３２人目の素数さん mailto:sage [2020/05/21(木) 18:53:16.25 ID:hAZkHjNF.net] 結局、統計ってこういう道楽なんだよなぁ。 632 名前：１３２人目の素数さん mailto:sage [2020/05/21(木) 19:01:50.90 ID:hAZkHjNF.net] >>581 いやウイルスの振る舞いも素粒子の振る舞いもばらつきがあるんじゃないの？ 存在も位置も確率的にしかわからない、という議論になると思う。 633 名前：１３２人目の素数さん mailto:sage [2020/05/21(木) 21:14:17 ID:hAZkHjNF.net] >>571 最初の２つでは頑強と言えそう。 平均身長のたとえだと 事前分布に1〜２ｍを選んでも１〜１０ｍを選んでも同様の結果だが、０〜１ｍを選ぶと現実離れした結果が返ってくるという感じだな。 634 名前：１３２人目の素数さん mailto:sage [2020/05/21(木) 21:27:14 ID:hAZkHjNF.net] 次のおもちゃ しばらく、これで遊べそう。 臨床所見からロジスティック回帰でCOVID19の確率を出すペーパーがでるだろうなと思っていた。 Real-time tracking of self-reported symptoms to predict potential COVID-19 https://www.nature.com/articles/s41591-020-0916-2 635 名前：１３２人目の素数さん [2020/05/22(金) 01:26:24 ID:MSJJjK3u.net] >>591 馬鹿の定義は？単なるお前の主観じゃねーかｗ ベイズ統計にのめり込むと馬鹿になる事後確率分布でも求めてろ。 636 名前：１３２人目の素数さん mailto:sage [2020/05/22(金) 04:53:19 ID:C0tPqEF8.net] こういうのも興味ある人多い？感染してからの日数とPCR陰性になる確率の関係。 https://twitter.com/AdamJKucharski/status/1260839061318705152 (deleted an unsolicited ad) 637 名前：１３２人目の素数さん mailto:sage [2020/05/22(金) 05:39:46 ID:s9txG4+B.net] 各国のロックダウンの度合いを数値化してるところ。色々と分析に使えるかも。 https://ourworldindata.org/grapher/covid-stringency-index?tab=chart&year=2020-05-07&country=JPN+NOR+SWE+USA 638 名前：１３２人目の素数さん mailto:sage [2020/05/22(金) 08:35:38 ID:DQesskhT.net] >>599 あれ3/3になったのか。 それだと、95％CrIは0.5559329 1 事前分布にはJeffereysで計算 639 名前：１３２人目の素数さん mailto:sage [2020/05/22(金) 09:06:15 ID:DQesskhT.net] >>600 ありがとうございます。 おねーちゃんと濃厚接触したあとは、何日目に検査したらいいかの参考になりました。 640 名前：１３２人目の素数さん mailto:sage [2020/05/22(金) 09:49:33 ID:DQesskhT.net] >>600 論文を検索してみた、この論文からのグラフのよう。 Variation in False-Negative Rate of Reverse Transcriptase PolymeraseChain Reaction?Based SARS-CoV-2 Tests by Time Since Exposure https://www.acpjournals.org/doi/pdf/10.7326/M20-1495 641 名前：１３２人目の素数さん mailto:sage [2020/05/22(金) 10:18:18 ID:DQesskhT.net] Measurements: A Bayesian hierarchical model was fitted to estimate the false-negative rate by day since exposure and symptom onset. とのことで、 こんな、ありふれたベイズ階層モデル x[j,t] ~ Binomial(n[j,t],p[j,t]) logit(p[j,t]) = βj + β*1log(t) + β2*log(t)^2 + β3*(t)^3 βj ~ Normal(β0,σ2) 642 名前：１３２人目の素数さん mailto:sage [2020/05/22(金) 10:36:46 ID:DQesskhT.net] Rとstanのコードが公表されているので、遊べそう。 https://github.com/HopkinsIDD/covidRTPCR 643 名前：１３２人目の素数さん [2020/05/22(金) 17:38:12 ID:MSJJjK3u.net] >>602 1/1だよ。あんただけｗ 644 名前：１３２人目の素数さん mailto:sage [2020/05/22(金) 19:01:13 ID:DQesskhT.net] >>607 信頼区間出せない馬鹿発見！ 645 名前：１３２人目の素数さん mailto:sage [2020/05/22(金) 19:09:59 ID:xdbfxKAO.net] このスレででてる信頼区間なんかほとんど統計の検定試験では出てこない俺様信頼区間だけどな。 信頼区間の定義ちゃんと言える奴の方が少ないだろ。 646 名前：１３２人目の素数さん mailto:sage [2020/05/22(金) 19:11:40 ID:DQesskhT.net] 1/1でも信頼区間を出せちゃうのがベイズ 647 名前：１３２人目の素数さん mailto:sage [2020/05/22(金) 19:13:01 ID:DQesskhT.net] >>609 ベイズでは信用区間と呼んで信頼区間と区別する人もいるね。 CIでなくてCrIという記載も目にする。 648 名前：１３２人目の素数さん mailto:sage [2020/05/22(金) 19:20:12 ID:DQesskhT.net] ベイズではconfidence interval でなくてcredibility interval 649 名前：１３２人目の素数さん [2020/05/22(金) 23:34:07 ID:MSJJjK3u.net] >>612 あんたの独占スレになったな。1/1達成、おめでとうｗ 650 名前：１３２人目の素数さん mailto:sage [2020/05/23(土) 06:44:49 ID:COZ69QMb.net] bootstrap使うと大抵の数値の信頼区間は出せるけど1/1には無理だな。 651 名前：１３２人目の素数さん mailto:sage [2020/05/23(土) 15:08:38 ID:COZ69QMb.net] JAGSを使って1/1のときの95％ Credibility Intervalを求めるスクリプト 事前分布はJeffereys library(rjags) data=list(x=1,n=1) cat(' model{ x ~ dbin(p,n) # binomial distribution p ~ dbeta(0.5,0.5)　# Jeffereys prior }',file='tmp.txt') jagsModel=jags.model('tmp.txt',data=data,n.chains=4) update(jagsModel) codaSamples=coda.samples(jagsModel,n.iter=1e5,var='p') gelman.plot(codaSamples) # 収束を確認 js=as.data.frame(as.matrix(codaSamples)) HDInterval::hdi(js$p) # 95% Highest Densty Interval # 既存のパッケージで確認。 HDInterval::hdi(qbeta,shape1=0.5+1,shape2=0.5) binom::binom.bayes(1,1) 652 名前：１３２人目の素数さん mailto:sage [2020/05/23(土) 17:34:11 ID:COZ69QMb.net] エクセルのツールをCDCが公開しているね

---

---

[ 続きを読む ] / [ 携帯版 ]

read.cgi ver5.27 [feat.BBS2 +1.6] / e.0.2 (02/09/03) / eucaly.net products.
担当:undef