１．SIRモデルを統計言語[R]でシミュレーションする。
以下にRのプログラムを示す。なるべくシンプルなコーディングを目指したので、統計言語[R]の初歩的な知識があれば、コメント(＃)行をみれば、何をやっているのかがわかるであろう。
このプログラムは、自分のPCに[R]が既にインストールしてあれば、その上で動く。あるいは、手軽に試してみたい場合には、ネット上にオンラインで[R]を走らせるサイトがいくつかある。たとえば、
1) ここcompile R onlineをクリックして、以下のプログラムをコピペして、左下にある[run it]をクリックするか、自分のPCでF8を押せばよい。あるいは、
2) ここpaiza.ioでも同様なことができる。
試してみてはいかがだろうか？
# ==================================================
# [R]によるSIRモデルの計算 ＃#==================================================
##### シミュレーション期間の設定とS,I,Rベクトルの定義
No_of_Days=50; N1=No_of_Days+1; S=I=R=R_t=numeric(N1);
##### S,I,Rの初期値とパラメータ値の設定
S[1]=100; I[1]=1; R[1]=0; beta=0.01; gamma=0.1
##### SIRモデルの繰り返し計算
for (t in 1:No_of_Days){
S[t+1] = S[t] - beta*S[t]*I[t]
I[t+1] = I[t] + beta*S[t]*I[t] - gamma*I[t]
R[t+1] = R[t] + gamma*I[t]
R_t[t+1] = S[t+1]*beta/gamma
}
##### S,I,Rのグラフを描く（S,I,Rと横軸に日数をひとまとめのデータフレームとしてから）
output = data.frame( S, I, R, date=1:N1 )
library(ggplot2)
ggplot( output, aes(date)) +
geom_line(aes(y = S, colour = "S")) +
geom_line(aes(y = I, colour = "I")) +
geom_line(aes(y = R, colour = "R"))
##### 再生産数 R_tのグラフを描く
plot(R_t, type="l", lwd=2, xaxt="n")
# ==================================================
2. S,I,R の推移をグラフでしめす
上のプログラムを実行すれば、次のようSIRと実行再生産数R_tの推移をえることができる。
自分なりのシナリオを描いて、政策シミュレーションをやってみたい場合には、上のプログラムの
##### S,I,Rの初期値とパラメータ値の設定
S[1]=100; I[1]=1; R[1]=0; beta=0.01; gamma=0.1
で、それぞれの初期値S[1]=100; I[1]=1; R[1]=0; と 感染率と回復率のパラメータを変化させてることで可能である。
自分なりのシナリオを描いて、政策シミュレーションをやってみたい場合には、上のプログラムの
##### S,I,Rの初期値とパラメータ値の設定
S[1]=100; I[1]=1; R[1]=0; beta=0.01; gamma=0.1
で、それぞれの初期値S[1]=100; I[1]=1; R[1]=0; と 感染率と回復率のパラメータを変化させてることで可能である。
 |
| 統計言語[R]によるSIRモデルの計算結果 |
 |
| Rによる実行再生算数：R_t |
コメント
コメントを投稿
ここに感想をお願いいたします。