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

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

数学統計に詳しい人が語るコロナウイルス - 暇つぶし2ch■コピペモード
□スレを通常表示
□オプションモード
□このｽﾚｯﾄﾞのURL
■項目テキスト

800:１３２人目の素数さん
 20/07/08 11:05:25  XD7Ql8W/.net
>>752 
>751の論文に95%CIが書いてあるのでそれを使って階層モデルを作ろうかと思って原著を読んだら95%CIはブートストラップ法で算出したと書いてあった。 
 
Rtの方は期待値がCI境界の平均なので正規分布を期待したのだが。 
dispersion parameterは非対称の分布みたい。

801:１３２人目の素数さん
 20/07/08 20:15:37  I3BoIViR.net
>>752 
レスありがとうございます。 
 
パラメータの分布を勝手に設定して階層モデルもどきでシミレーションしてみた。数字のお遊びでしょうけど。 
 
From the empirical offspring distribution and 
fitted negative binomial distribution shown in Figure 2B, 
we estimated an observed reproductive number (R) of 0.58 (95% CI: 0.45 - 0.71) 
and dispersion parameter (k) of 0.45 (95% CI: 0.31 - 0.76). 
 
Rは正規分布、kはガンマ分布として、期待値と95%CIが上記に適合するようなパラメータを 
ニュートン・ラフソン法で算出すると 
正規分布の標準偏差sd＝0.0663277494 
> pnorm(0.71,0.58,0.0663277494)-pnorm(0.45,0.58,0.0663277494) 
[1] 0.95 
ガンマ分布の形状パラメータsh＝24.38466,尺度パラメータ=0.01845 
> pgamma(0.76,sh,scale=sc)-pgamma(0.31,sh,scale=sc) 
[1] 0.95 
が得られた。 
 
各パラメータが上記の分布に従うとして 
R ~ normal(mean=0.58,sd=sd) 
k ~ gamma(shape=sh,scale=sc) 
Rt ~ nbinom(dispersion=k,mean=R) 
のモデルでシミレーション 
 
既述のコードに 
RR=rnorm(1e5,R,sd) # R ~ dnorm(R,sd) 
kk=rgamma(1e5,sh,scale=sc) # k ~ dgamma(shape=sh,scale=sc) 
RRt <- function(x){ # Rt ~ dnbinom(k,mu=R) 
  k=sample(kk,1) 
  R=sample(RR,1) 
  rnbinom(x,k,mu=R) # pick x random Rt numbers 
} 
を追加して 
while loop内のinfecteeの数を 
infectee=sum(RRt(infector)) # number of infectee  
に変更。 
 
最初の10人に何人が感染していたかの95％推定 
> HDInterval::hdi(spreader)[1:2] # 95% credibility interval 
lower upper  
    3     9  
> BEST::plotPost(spreader,xlim=c(1

次ページ