数学 統計に詳しい人が語るコロナウイルス
■ このスレッドは過去ログ倉庫に格納されています
東大数学科卒の元官僚はこう分析してるが、お前らはどうなると思う? http://www.zakzak.co.jp/soc/news/200220/dom2002200003-n2.html 中国国外感染者の中国国内との比率をみると、 1月20日の数字公表以降は、0・8〜2・6%で比較的安定している。 これは、新型肺炎の感染者のほとんどは中国国内、それも湖北省に集中しているからだ。 ちなみに中国国外での感染者数は、中国国内の1・1%だ(2月16日現在)。 本コラムで紹介したが、現時点では、最終的な中国国内の感染者数は20万人超と筆者は推計している。 となると、中国国外の感染者は数千人程度になるだろう。 中国国外のうち日本の比率は1割弱なので、日本の感染者数は数百人程度であろう。 その場合、死者も数人から10人程度になるだろう。 こうした推計をすると、今の感染者は氷山の一角だと思われるが、今後の増加ペースはどうなるだろうか。 新型コロナウイルスの検査は簡単に行えるので、今後、日本での感染者数は増えていくだろう。 ある時点ではそれがネズミ算的に増えるかのように思える局面もあるだろうが、 筆者の推計が正しければ、現時点ではせいぜい数百人が一つのメドだ。 >>179 単にお勉強ができただけ。 頭が良くないのさ。 自分の頭で物事を考えるってことができない。 >>181 その通り。 具体的には理学部の数学科と物理学科。 工学部にも時々もの凄いのがいる。 >>559 いや、Rのパッケージ HDIntervalで計算しているから誤解していないと思う。 内部の処理コードをみると信頼区間幅が最小になるのを最尤法で出しているね。 >>561 ほんとうに頭のいい人は医学部でなく理学部か工学部に行く。 ほんとうに頭の悪い奴は底辺シリツ医大に行く。 シリツ医大には手先の器用な人もいるが、頭が器用な奴をみたことがない。 >>550 (自答) # # 人物dが発症してdelay日後に濃厚接触したキャバ嬢cが発症 # cの感染がdより先行していた確率は? 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 ])') 2日後の発症だと > fn.stan(2,print=T,iter=5000,warmup=1000) [1] 0.2945 3割くらいはあとから発症した方が先に感染していた可能性がある。 >>563 一方的で申し訳ないが私立大医学部は金持ちのバカ息子が行くイメージ。 西浦さんはさんざん適当なことを言って世論を煽ってどう責任を取るのかな? >>567 少なくとも西浦氏は算出コードを公開しているだけでも好感が持てる。 >>566 それは正しい認識。 凄いのはド底辺シリツ医の馬鹿さ加減だよ。 裏口バカと呼ばれるのがよくわかる。 http://imagizer.imageshack.com/img923/2715/RosCsf.jpg http://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. 西浦モデルの再生算数の事前分布を変化させてグラフにしてみた。 西浦モデルでのデフォルト https://i.imgur.com/G1wVYgI.png 事前分布を一様分布(0,10)近似の正規分布で近似させた場合 https://i.imgur.com/doS5LEu.png 再生算数の平均0、標準偏差1の場合 https://i.imgur.com/doS5LEu.png 印象としては、西浦モデルは頑強性robustnessのあるモデルとは言えない。 事前分布に大きく影響されるモデルだと思う。 (url訂正) 西浦モデルの再生算数の事前分布を変化させてグラフにしてみた。 西浦モデルでのデフォルト https://i.imgur.com/G1wVYgI.png 事前分布を一様分布(0,10)近似の正規分布で近似させた場合 https://i.imgur.com/doS5LEu.png 再生算数の平均0、標準偏差1の場合 https://i.imgur.com/0J1RpDa.png 印象としては、西浦モデルは頑強性robustnessのあるモデルとは言えない。 事前分布に大きく影響されるモデルだと思う。 >>571 誤解されるのは不本意なの追記するけど、ソースを公開する西浦先生の姿勢は高く評価している。 隠蔽・改竄・破棄する安倍とは大違い。 だからベイズ推定が統計学の世界でメジャーにならんのだろ? 論理的根拠のない“事前分布”なるもので答えがひょいひょい変わるのでは社会的な影響が大きい防疫政策の決定には使えない。 普通の統計学の検定なら理論的に根拠のある数字、推論しか使わない。 計算は大変だけど。 >>573 成人の平均身長を1〜2mに事前分布にするのは納得できるし、 生まれる子供が男子である確率は0.4〜0.6というのも俺は納得できる。 PCRの感度が30〜70%として計算するのも納得できるからその設定で階層ベイズモデルを組むことには異論はないな。 こういう問題 あるタクシー会社のタクシーには1から通し番号がふられている。 タクシー会社の規模から保有タクシー台数は100台以下とわかっている(弱情報事前分布)。 この会社のタクシーを5台みかけた。最大の番号が60であった。 この会社の保有するタクシー台数の期待値と95%信用区間(信頼区間)を求めよ。 をベイズで解くときは、 60台〜100台である確率を一様分布として処理しているから これに異論があるのは理解できるけど 日本人成人の平均身長を推定に1〜2mを事前分布に想定するのには俺は異論はないね。 一様分布ではなくてガンマ分布にすべきだというのは議論になるとは思うけど。 >>573 ベイズ信奉者から、ベイズ論者を採用したGAFA(Google等)が成長したと教わった。 迷惑メールのフィルタリングとか雑音の除去とか日常的に役立っているというぞ。 >普通の統計学の検定なら理論的に根拠のある数字、推論しか使わない。 >計算は大変だけど。 普通の統計学こそ、無理やり既知の分布に当てはめようとするんだよね。 MCMCの方の方が応用が広い。 プロ野球選手の打率は?と問われたら選手次第で異なる、と誰でもわかるのに 確定不能の平均値が存在していると妄想して計算を始めるのが古典主義統計。 つまり、値は存在するけど確定できないという信仰の世界。 昨今の新コロナでいえば、PCRの検査の感度・特異度が一定と考えるのが古典主義。 プロ野球選手の打率と同じでそんなのは場面場面で変化するよ、と考えて計算するのがベイズ。 >>576 もちろん推定の有力な方法であるにせよ、元の仮定に何の根拠もないわけだからそれから得られる結論には論理的根拠はない、ないが、数学的に伝統的な手法で与えるた結論と大差ない事がなんらかの保証があるなら、有力になる。 それが論理的に“大差ない結論が得られる”事が示されてるなら単なる計算手法に過ぎないし、示されていなくても経験的に“よい結果ぎ得られる事が多い”ならそのジャンルではそこそこ信頼するに値するんだろう。 しかしなんらかのモデルでは答えが一意に定まらず、事前分布の選び方により大きく答えが違ってしまう場合があっても不思議はないし、そのような場合ではやはり、“ではどう計算するのが正しいかのか”の論証を待たなければ信頼するのは危険になる。 古典的(頻度主義)統計信者って、この計算はどうやるんだ? 俺は乱数発生させて計算できるけど、 そればベイズなのかどうかは知らんが、条件付き確率なのでベイズなんだろうな。 (開業医スレに投稿したけど、回答できるやつは0) 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 ある開業医が新型コロナ肺炎に罹患したとする。 行動調査によって発症前にキャバクラに行っており接客したキャバ嬢が開業医発症の2日後に発症していたことがわかった。 キャバ嬢は開業医から移されたと主張して1億円の賠償を求めている。 潜伏期間には幅がありキャバ嬢から移された可能性もあると主張してその確率を計算して賠償金を値切りたい。 いくら値切れるか計算せよ。 >>580 潜伏期間なんて数学的に決定できるハズないやん? 数学がやるのは例えば感染日と実際発症した日が確定できるような症例がある程度以上あって、それが従うと病理学的に信頼できる分布の族があって、その中からデータに最も“沿う”分布を選び出すことしかできない。 それは何件も実際の患者のウイルス量の統計をとったり、ウイルスが体内でどのように増えていくのかの病理学的研究データがあって初めて可能になる。 >>571 よく理解できていないので質問ですけど 事前分布とは具体的に何の分布ですか? 基本再生算数の推定値の分布? 実行再生算数の推定値の分布? 実行再生算数の事前分布は基本再生算数の分布としたらいいのかなと思いますけど 潜伏期は病原体と罹患者のパワーバランスで決まるだろうから 定数でなくてばらつきはあると思うね。 >>572 >隠蔽・改竄・破棄する安倍とは大違い。 安倍ちゃんがやってるわけじゃないだろ。 指示もされてねーのに、官僚がやってんだよ。 そんなこともわからんようでは、あんたの解析もあてにならんな。 >>582 西浦のソースだと Rt ~ normal(2.4 ,2) >>586 安倍じゃなければ官僚もまともだったんじゃねえの? 安倍らの集団に組み込まれると東大卒も馬鹿になるようだぞ。 >>581 その通り。 結局、そういう本質的なデータや理論抜きでは、ベイズ推定やったって 限界があるし、結果の説得力もない。 まあ、適用限界ってものがあるのは何やっても同じだけど。 >>588 >安倍らの集団に組み込まれると東大卒も馬鹿になるようだぞ。 統計的に確認してみれば?ベイズ推定でもっともらしい結果が出るんじゃね?w >>581 いやウイルスの振る舞いも素粒子の振る舞いもばらつきがあるんじゃないの? 存在も位置も確率的にしかわからない、という議論になると思う。 >>571 最初の2つでは頑強と言えそう。 平均身長のたとえだと 事前分布に1〜2mを選んでも1〜10mを選んでも同様の結果だが、0〜1mを選ぶと現実離れした結果が返ってくるという感じだな。 次のおもちゃ しばらく、これで遊べそう。 臨床所見からロジスティック回帰でCOVID19の確率を出すペーパーがでるだろうなと思っていた。 Real-time tracking of self-reported symptoms to predict potential COVID-19 https://www.nature.com/articles/s41591-020-0916-2 >>591 馬鹿の定義は?単なるお前の主観じゃねーかw ベイズ統計にのめり込むと馬鹿になる事後確率分布でも求めてろ。 各国のロックダウンの度合いを数値化してるところ。色々と分析に使えるかも。 https://ourworldindata.org/grapher/covid-stringency-index?tab=chart& ;year=2020-05-07&country=JPN+NOR+SWE+USA >>599 あれ3/3になったのか。 それだと、95%CrIは0.5559329 1 事前分布にはJeffereysで計算 >>600 ありがとうございます。 おねーちゃんと濃厚接触したあとは、何日目に検査したらいいかの参考になりました。 >>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 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) このスレででてる信頼区間なんかほとんど統計の検定試験では出てこない俺様信頼区間だけどな。 信頼区間の定義ちゃんと言える奴の方が少ないだろ。 >>609 ベイズでは信用区間と呼んで信頼区間と区別する人もいるね。 CIでなくてCrIという記載も目にする。 ベイズではconfidence interval でなくてcredibility interval >>612 あんたの独占スレになったな。1/1達成、おめでとうw bootstrap使うと大抵の数値の信頼区間は出せるけど1/1には無理だな。 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) エクセルのツールをCDCが公開しているね。 https://www.cdc.gov/coronavirus/2019-ncov/hcp/COVIDSurge.html COVID-19Surge is a spreadsheet-based tool that hospital administrators and public health officials can use to estimate the surge in demand for hospital-based services during the COVID-19 pandemic. A user of COVID-19Surge can produce estimates of the number of COVID-19 patients that need to be hospitalized, the number requiring ICU care, and the number requiring ventilator support. The user can then compare those estimates with hospital capacity, using either existing capacity or estimates of expanded capacity. With COVID-19Surge, users define the population in the hospital “catchment area” or local jurisdiction. The user also enters the number of cases to date, and the available hospital resources (non-ICU beds, ICU beds, and mechanical ventilators). Users can assess up to three community mitigation strategies simultaneously and compare the impacts on hospital resources versus a “no intervention” scenario. @kazu_fujisawa 西浦先生が一番の功労者だというのは僕も同意見だな。あと、細かいこと言うと、 あの限られたデータしかない中で、8割削減なら1ヶ月で収束、接触6-7割削減なら2ヶ月で収束、という予測精度もお見事としか言いようがない。 わかった後からだと 「簡単じゃん!」とか言えるけど 先んじて発表した西浦先生は すごくまともな研究者だと思います。 中国・武漢全域で行われた新型コロナウイルスのPCR検査で、189人が無症状の感染者と確認されました。 武漢では35日ぶりに新たな感染者が確認されたため、14日から「10日間大戦争」と銘打ち、延べ900万人にPCR検査を行いました。中国メディアによりますと、そのうち657万人の結果が出て、189人が無症状の感染者でした。 10万人のうち2.87人が無症状感染者の計算になるとしています。中国のSNSには、検査人数のあまりの多さに看護師が泣き叫ぶ動画が投稿されていて、 一日の検査能力である最大10万件を大幅に超えた検査の精度を疑問視する声もあります。10日間ですべての検査が終わらなかったので、検査は26日午後も続いています。 [2020/05/26 13:22] https://news.tv-asahi.co.jp/news_international/articles/000184795.html https://news.tv-asahi.co.jp/articles_img/000184795_640.jpg 189/6570000は 有病率、感度、特異度の事前分布はどうするかな? 有病率:一様分布 感度:最頻値0.6標準偏差0.1のベータ分布 特異度:最頻値0.9標準偏差0.05のベータ分布 として検査陽性数は有病率*感度+(1−有病率)*(1−特異度)の確率に従う二項分布 でいいか? >>618 再生産数から収束に至るまでの期間の計算法って公開されてたっけ? >>622 それは簡単に計算できるだろ?逆は難しくて、それが西浦先生の功績だと思うけど。 >>622 いっぱんに再生産数Rの伝染病が収束を始める免疫の比率は1-1/R。 R=2.5なら免疫獲得者が1-1/2.5=60%になれば拡大はしない。 なので60%になれば平衡状態になり、最低でもそれ以上減らさないとダメなのは明らか。 集団免疫の話読んだ事ある人間なら常識。 オレみにわかではあるが。 感染率が70%減、80%減になったとき、どれくらいの速度で落ちていくかもsirでだいたいわかる。 え? >>617 って皮肉じゃないの?w 西浦さんの8割削減って言葉がメディアに乗るころには拡大再生産数は 1を切ってたんだよね(4/1頃)。拡大再生産数が2を越えて高止まりして たのは3/15-3/25頃で、そこから1週間ほどで急落して1を切ってんだよね。 東京、大阪で自粛が叫ばれ、危機感がつのってたころだ。 >>625 それは今でも良くわかってないだろ。 何せ4月の初め頃って検査リソースの限界で東京とか陽性率が8割に近かった、もうほとんど陽性確定の人しか検査してなかった頃。 この“認知されたケース(confirmed case”となったケースのピークと実際の感染者数のピークは当然ずれてる。 実際の東京の罹患者は8万人を超えてるらしいので、相当数が検査されずに見過ごされてる。 おそらく見過ごされた多くのケースは陽性者数がピークを迎えた(陽性率が8割ぐらいあった)4/15前後に発症したケースだろうけど、それは明らかにconfirmed casesのピークより遅い。 >>627 4月の初め頃の検査対象はその1週間前に感染してんだから、 まさに専門家委員会が示したconfirmed casesの感染者数の ピークの頃だよ。 >>629 confermed casesのピークはもっと遅い4/15前後だよ。 その前と後に二つの最大のピークがある。 いずれにせよ、検査拒否の陽性率の最大がそのあたりにあって、どのくらいの見逃しがどのくらいあったのか科学的に推し量る方法などない。 その当時どんなコールセンターへの相談があり、どの程度怪しいケースを弾いていたのかなどコールセンターの資料をツブサに調べていくしかない。 そんなもんここで議論するような話ではない永遠の水かけ論にしかならない。 >>624 収束する時期は再生産数だけじゃでは計算できないのではという疑問。 >>628 有病率を10万分の3として、何人の検体を混ぜるのが一番効率がいいんだろ? 10万を3でわった3万3千人くらいかな? >>632 n/k+n/k(1-(1-p)^k)kを最小にするkだから、k=2Productlog(-0.5√-log(1-p))/log(1-p)だね。 p=3/100000を代入してk=183。 ちなみにk=10となるのはp=0.0111の時。そんなに陽性率が高いと見積もったわけでは無いだろうから、10ぐらいが混ぜて検査の限界なのだろう。 もっと究極的に検査回数を少なくすると>>419 の方法に行き着くと思う。 >>630 だ、か、ら、検査結果の判明した日付では4/15頃がピークだけど、検査の 待ち時間や発症までの潜伏期間を含めると、感染したのはその2週間以上 前になるっていうのが、専門家委員会の見積もり。 https://www.mhlw.go.jp/content/10900000/000630718.pdf >>635 confirmed cases の罹患→発症の平均時間と all cases の罹患→発症の平均時間にはズレが出る。 confirmed cases のそれが2週間以上でそのピークが4月より前でも all cases のそれはもっと小さいのでそのピークはもっと後ろになる。 具体例で実験してみればいい。 山がある曲線かく。 山をy軸方向に0.7倍したもの(軽症段階で相談した人)と0.3倍したもの(重篤化して認知されたもの)を作り、後者はx軸反転する。 前者は右に1週間、後者は左手に一週間ずらす。 この二つの曲線の差がその日に検査要請がかかった件数。 しかし検査リソースには限界があるのでx軸より下の重症者は全部検査されたとして、そこから一定幅より上の部分は検査拒否されたと考える、つまりy軸より下の部分を検査リソースの限界分だけ上に平行移動したものより下が実際に検査された軽症者の人数。 もとまった検査された軽症者の曲線を改めて一週間前にずらし、重症者のカステラの曲線と足してみる、コレがconfirmed casesに限った推定発症日の曲線。 元の山とできた山のピーク比べてみる。 ズレてるから。 >>636 だ、か、ら、>>627 で検査が滞って隠れ感染者が相対的に増えてたと 君がいう4月上旬の検査対象は3月下旬の感染者なんだから、全感染 者数のピーク時期はconfirmed casesの4/1頃のピークの前にこそなれ、 後ろにはならんでしょ? >>637 なるよ。 >>635 の実験やってみた? 必ずconfirmed casesのピークはall casesのピークより前にズレる。 つまりall casesのピークはconfirmed casesのピークより後ろにズレる。 >>637 前提に問題がありすぎ。検査数は逐次増えてて、感染者のピーク前後では ほぼ検査数に比例して陽性者数も増えてる(陽性率は十数%前後で大きな 変化はない)。キャパ一定で推移してるわけではない。 >>640 前にも書いたけどその話をしだすとココでは絶対答えが出ない水掛け論にしかならない。 結局症状が出てても検査を断られる4月頭時点での検査体制では、実際どれくらいのズレが発生したのかはハッキリとは言えない。 それは後年の調査である程度はわかるかもしれないけど。 何せ1日数万件もあった相談と検査拒否の実態を掴むのは容易ではない。 いずれにせよ、実際の感染のピークが confirmed casesのピーク日 -conformed cases内の(確定日-推定罹患日)の平均 より後ろにズレるのは確実。 何故ならどの程度の検査拒否があったにせよ、この相談日-罹患日の値が小さいグループ(症状の軽いグループ)の方が検査拒否された割合が有意に高いので。 まぁ実際confirmed casesのピーク日-(xxx) のxxxに入れるべき数値は2週間よりはずっと小さい。 職場で最初に発症した人が感染源と決めつけてはならないことを 数学的には確認。 あるクリニックで院長が新コロナを発症。 翌日2人の看護師が発症、その次の日に3人の看護師が発症。 各人が最初に感染していた確率を計算せよ。 こういう計算ができると、誰が移したとか責めたりせずにすむ。 (春節ウェルカムした安倍が悪い、で職員の心が一つになるw 計算に必要な数値は #--- incubation period --- # http://www.nejm.org/doi/full/10.1056/NEJMoa2001316 # Li et al NEJM 2020 # 潜伏期間は対数正規分布でパラメータは ln_par1 = 1.434065 ln_par2 = 0.6612 RとStanでやってみた。 > data.frame(p=round(sapply(1:n,fi),2)) p 1 0.26 2 0.18 3 0.18 4 0.13 5 0.13 6 0.13 コードはここ https://egg.5ch.net/test/read.cgi/hosp/1583407025/199 >>632 ウイルスを検出できる最低限度の濃度とかあるんじゃないかな 制約付きで何人分を混ぜるかを最大化する問題になると思う >>641 しつこいな。どこからそんな式が出てくるんだよw 現状で分かってる事実から推測して全感染者数の ピークが後ろにくる理由はないよ。むしろ前にくるほうが自然。 >>642 発症前から感染力があり、潜伏期間にも幅があるんだから、 その情報だけでだれが感染源かなんて断定できるわけない のは当たり前だろ。それでも最初の発症者が一番怪しい ってことに変わりはない。 常識で考えろ。 >>646 院長と3日めの発症した看護師とどれくらい確率が違うかは 常識では出せないだろ。 こういう計算には数値化が必要。 ある開業医が新型コロナ肺炎に罹患したとする。 行動調査によって発症前にキャバクラに行っており接客したキャバ嬢が開業医発症の2日後に発症していたことがわかった。 キャバ嬢は開業医から移されたと主張して1億円の賠償を求めている。 潜伏期間には幅がありキャバ嬢から移された可能性もあると主張してその確率を計算して賠償金を値切りたい。 いくら値切れるか計算せよ。 >>642 Mathematicaだと cond = Distributed[{a, b, c, d, e, f},ProductDistribution[{LogNormalDistribution[1.434065, 0.6612], 6}]]; NProbability[a > Max[b - 1, c - 1, d - 2, e - 2, f - 2], cond] NProbability[b - 1 > Max[a, c - 1, d - 2, e - 2, f - 2], cond] NProbability[d - 2 > Max[a, b - 1, c - 1, e - 2, f - 2], cond] >0.257392 >0.178947 >0.128264 >>649 レスありがとうございます。 有効数字2桁で同じ結果なので自分のコードの検証になりました。 中国湖北省武漢市(人口約1100万人)が今月中旬から全市民を対象に実施していたPCR検査(遺伝子検査)がおおむね終了した。 地元紙「湖北日報」(電子版)によると、直近に検査を済ませた人や乳児を除く900万人以上から検体を採取。 武漢では1日に100万人以上の検査が可能で、24日までに約650万人の検査が終わり、無症状感染者が218人見つかった。 事前分布に有病率は一様分布、感度は40-60%、特異度は90−100%のベータ分布を設定してStanでMCMCしてみた。 mean se_mean sd 2.5% 50% 97.5% n_eff Rhat spc 0.999977 0.000000 0.000008 0.999965 0.999976 0.999994 12593 1.00001 sen 0.452612 0.002116 0.220226 0.073638 0.444200 0.871391 10834 1.00000 prev 0.000041 0.000001 0.000118 0.000001 0.000023 0.000183 8206 1.00004 p 0.000034 0.000000 0.000002 0.000030 0.000034 0.000039 36412 0.99998 spc:特異度、sen:感度 prev:有病率 p:陽性率 900万人での陽性者数の予想 最頻値 306.15 95%信頼区間 lower upper 266.46 346.70 1100万人のうちの有病者数予想 > summary(ms$prev*11e6) ; MODE(ms$prev*11e6)[1] ; hdi(ms$prev*11e6)[1:2] Min. 1st Qu. Median Mean 3rd Qu. Max. 0 108 251 449 483 137796 x 179.24 lower upper 0.020213 1360.523395 岩田健太郎 Kentaro Iwata @georgebest1969 日本が第一波をかなりうまく乗り切ろうとしているのだけど、最大の功労者の一人は西浦博先生だよ。それは絶対に間違いない。 午後5:44 2020年5月20日 >>652 俺は志村けんとフィリピンパブ親父だと思う。 政府:「接触8割減を目指して下さい。」 マスコミ:「今朝の○○駅の映像をご覧下さい。人出8割減にはほど遠いですね...」 政府広報: 「接触8割減には、(1-x)^2=1-0.8 → x=0.552786 だから 人出55%減でいいんだけど、 馬○なマスコミが、勝手に目標を上げてくれてる。あえて指摘する必要ないよね...」 >>653 藤浪も仲間に入れてやって。 岡江久美子も5月10日前後の感染急減には貢献したかも。 しかし、岩田健太郎に擁護されてもむしろ逆効果でお気の毒>8割おじさん >>648 賠償金ゼロまで値切れる。 俺じゃない可能性が十分あるのならそれでいい。 期待値で賠償金を計算されちゃかなわん。 >>656 期待値を求めよ、では食いつきが悪いと思って金の問題にアレンジして開業医スレに書いてみたんだが、 答られる開業医は0だったよ。 ■ このスレッドは過去ログ倉庫に格納されています
read.cgi ver 07.5.4 2024/05/19 Walang Kapalit ★ | Donguri System Team 5ちゃんねる