SIR in the R

RstudioでSIRモデルA

免疫を獲得しないで回復する率を0と置いた場合は以下のようになる。

 

線型代数,ベクトル解析,慣性モーメント,解析力学,微分方程式,フーリエ解析,物理学,数学

 

この第2方程式より以下が求まる。

 

さらに基本再生産数を次のように表す。

 

時刻時刻tにおける総人口を線型代数,ベクトル解析,慣性モーメント,解析力学,微分方程式,フーリエ解析,物理学,数学と置いた場合、どのようなモデルであっても

線型代数,ベクトル解析,慣性モーメント,解析力学,微分方程式,フーリエ解析,物理学,数学

のような等式が成り立つ。

 

上記関係式に関して右辺の総人口線型代数,ベクトル解析,慣性モーメント,解析力学,微分方程式,フーリエ解析,物理学,数学は、実際の総人口数ではなく全人口を1とおいて考える。

 

rstudioを起動させ次のように打ち込む(または以下のコードをコピペする)

 

sir202104 <- function(t, y, parameters){
with(
as.list(c(parameters, y)),
list(c
(S = -lambda * S * I,
I = lambda * S * I - gamma * I,
R = lambda * I)
)
)
}

線型代数,ベクトル解析,慣性モーメント,解析力学,微分方程式,フーリエ解析,物理学,数学

S、I、Rのそれぞれの変数に値を代入していく。
ちなみに代入している値はあくまで架空の数値なので本気にしないように。
y_init <- c(S = 0.998537432, I = 0.000213768, R = 0.001248782)

 

out202104にy_init、経過日数の範囲、生成データ、線型代数,ベクトル解析,慣性モーメント,解析力学,微分方程式,フーリエ解析,物理学,数学線型代数,ベクトル解析,慣性モーメント,解析力学,微分方程式,フーリエ解析,物理学,数学をdeSolveのパッケージ、“ode”を使って代入。
このode(常微分方程式の略?)が入っていないと動作しないのでライブラリ読み込みしていない場合はlibrary(deSolve)と打ち込んで読み込んでから出力をする。
out202104 <- ode(y = y_init, times = 1:100, func = sir202104, parms = c(lambda = 0.9, gamma = 0.1))

> out202104
time      S     I     R
1 1 0.9985374320 0.0002137680 0.001248782
2 2 0.9982412177 0.0004770163 0.001545476
3 3 0.9975856819 0.0010595619 0.002202387
4 4 0.9961291960 0.0023537033 0.003663488
5 5 0.9929072736 0.0052156544 0.006903230
6 6 0.9858225492 0.0115047026 0.014064316
7 7 0.9704548361 0.0251266669 0.029776055
8 8 0.9380492518 0.0537585492 0.063739373
9 9 0.8736448587 0.1102596172 0.134869298
10 10 0.7594138407 0.2089210159 0.274995872
11 11 0.5922624058 0.3484511232 0.523587821
12 12 0.4048325562 0.4936057964 0.904064409
13 13 0.2467751147 0.5966631001 1.399065649
14 14 0.1408441040 0.6402803569 1.959889434
15 15 0.0790277623 0.6378925161 2.537727076
16 16 0.0450236191 0.6093852100 3.100330120
17 17 0.0264805084 0.5689529912 3.631108086
18 18 0.0161876436 0.5245581509 4.123297431
19 19 0.0103008363 0.4802167218 4.575351559
20 20 0.0068163558 0.4378200691 4.988281758
21 21 0.0046803568 0.3981807049 5.364260026
22 22 0.0033258799 0.3615730656 5.705919073
23 23 0.0024393110 0.3280056426 6.016005000
24 24 0.0018414285 0.2973589587 6.297206097
25 25 0.0014271865 0.2694545577 6.552073884
26 26 0.0011329254 0.2440918293 6.782986791
27 27 0.0009191595 0.2210663166 6.992140296
28 28 0.0007606108 0.2001795920 7.181547756
29 29 0.0006407561 0.1812437970 7.353048604
30 30 0.0005486198 0.1640836679 7.508318993
31 31 0.0004766873 0.1485373620 7.648883138
32 32 0.0004197378 0.1344562589 7.776125612
33 33 0.0003740788 0.1217044258 7.891303040
34 34 0.0003370439 0.1101578867 7.995555207
35 35 0.0003066903 0.0997037918 8.089915243
36 36 0.0002815819 0.0902395553 8.175319348
37 37 0.0002606343 0.0816720100 8.252615783
38 38 0.0002430238 0.0739166188 8.322572798
39 39 0.0002281154 0.0668967004 8.385886240
40 40 0.0002154126 0.0605427163 8.443186422
41 41 0.0002045267 0.0547916623 8.495043880
42 42 0.0001951499 0.0495864596 8.541975097
43 43 0.0001870350 0.0448753981 8.584447684
44 44 0.0001799824 0.0406116438 8.622884946
45 45 0.0001738293 0.0367527845 8.657670058
46 46 0.0001684425 0.0332604159 8.689149856
47 47 0.0001637116 0.0300997667 8.717638278
48 48 0.0001595449 0.0272393567 8.743419468
49 49 0.0001558656 0.0246506874 8.766750605
50 50 0.0001526092 0.0223079605 8.787864455
51 51 0.0001497209 0.0201878236 8.806971681
52 52 0.0001471543 0.0182691386 8.824262946
53 53 0.0001448696 0.0165327726 8.839910802
54 54 0.0001428326 0.0149614082 8.854071415
55 55 0.0001410139 0.0135393713 8.866886115
56 56 0.0001393881 0.0122524758 8.878482807
57 57 0.0001379329 0.0110878822 8.888977245
58 58 0.0001366292 0.0100339705 8.898474185
59 59 0.0001354600 0.0090802237 8.907068428
60 60 0.0001344106 0.0082171241 8.914845770
61 61 0.0001334679 0.0074360578 8.921883851
62 62 0.0001326205 0.0067292291 8.928252935
63 63 0.0001318584 0.0060895831 8.934016609
64 64 0.0001311724 0.0055107350 8.939232416
65 65 0.0001305547 0.0049869066 8.943952430
66 66 0.0001299983 0.0045128688 8.948223778
67 67 0.0001294968 0.0040838895 8.952089106
68 68 0.0001290446 0.0036956860 8.955587006
69 69 0.0001286367 0.0033443828 8.958752406
70 70 0.0001282688 0.0030264727 8.961616909
71 71 0.0001279367 0.0027387815 8.964209118
72 72 0.0001276370 0.0024784370 8.966554916
73 73 0.0001273663 0.0022428400 8.968677726
74 74 0.0001271219 0.0020296380 8.970598743
75 75 0.0001269011 0.0018367024 8.972337151
76 76 0.0001267016 0.0016621067 8.973910307
77 77 0.0001265213 0.0015041078 8.975333920
78 78 0.0001263585 0.0013611279 8.976622205
79 79 0.0001262112 0.0012317395 8.977788026
80 80 0.0001260781 0.0011146505 8.978843024
81 81 0.0001259578 0.0010086917 8.979797736
82 82 0.0001258491 0.0009128052 8.980661694
83 83 0.0001257507 0.0008260336 8.981443523
84 84 0.0001256618 0.0007475105 8.982151032
85 85 0.0001255813 0.0006764517 8.982791284
86 86 0.0001255086 0.0006121478 8.983370674
87 87 0.0001254428 0.0005539567 8.983894987
88 88 0.0001253833 0.0005012972 8.984369458
89 89 0.0001253295 0.0004536435 8.984798825
90 90 0.0001252808 0.0004105198 8.985187377
91 91 0.0001252368 0.0003714954 8.985538992
92 92 0.0001251969 0.0003361808 8.985857183
93 93 0.0001251609 0.0003042231 8.986145126
94 94 0.0001251283 0.0002753034 8.986405697
95 95 0.0001250988 0.0002491328 8.986641498
96 96 0.0001250721 0.0002254500 8.986854883
97 97 0.0001250479 0.0002040185 8.987047984
98 98 0.0001250261 0.0001846243 8.987222729
99 99 0.0001250063 0.0001670737 8.987380862
100 100 0.0001249884 0.0001511917 8.987523961

得られたデータをCSVの格納

得られたデータをCSVファイルにする。
write.csv(out202104, "out202104.csv")
csvファイルの読み込み
csvファイルはカレントディレクトリに生成されても読み込みを行わないと認識されない。
次のように入力してEnter
out202104_csv <- read.csv("out202104.csv",header=TRUE,sep=",")
これで読み込みができる。
次のように入力してEnterを押せばデータが表示される。 out202104_csv X time S I R 1 1 1 0.9985374320 0.0002137680 0.001248782 2 2 2 0.9982412177 0.0004770163 0.001545476 3 3 3 0.9975856819 0.0010595619 0.002202387 99 99 99 0.0001250063 0.0001670737 8.987380862 100 100 100 0.0001249884 0.0001511917 8.987523961

 

生成されたCSVファイルのデータをout202104_dataに格納。
out202104_data <- read.csv(file = "out202104.csv", header = TRUE, sep = ",")

out2021004_data_1にデータフレームを格納。

out202104_data_1 <- as.data.frame(out202104_data)

線型代数,ベクトル解析,慣性モーメント,解析力学,微分方程式,フーリエ解析,物理学,数学

 

3つのグラフを同時に描画

withを使ってSIRを描画。次のように打ち込む。

with(out202104_data_1, {
plot(time, S, type = "l", col = "blue", xlab = "経過日数", ylab = "人口")
lines(time, I, col = "red")
lines(time, R, col = "green")
})
legend("right", c("感染人口", "感染者", "回復者"), col = c("blue", "red", "green"), lty = 1, bty = "n")

上記のように打ち込んでEnterを押せば以下のように描画ができる。

線型代数,ベクトル解析,慣性モーメント,解析力学,微分方程式,フーリエ解析,物理学,数学

 

線型代数,ベクトル解析,慣性モーメント,解析力学,微分方程式,フーリエ解析,物理学,数学

> flu_out202104_data <- read.csv(file = "flu_out202104.csv", header = TRUE, sep = ",")

 

 

> flu_out202104_data
X time S I cases R
1 1 1 0.998537432 0.000213768 21376.80 0.001248782
2 2 2 0.998241218 0.000477016 47701.63 0.001545476
3 3 3 0.997585682 0.001059562 105956.19 0.002202387
4 4 4 0.996129196 0.002353703 235370.33 0.003663488
5 5 5 0.992907274 0.005215654 521565.44 0.006903230
6 6 6 0.985822549 0.011504703 1150470.26 0.014064316
7 7 7 0.970454836 0.025126667 2512666.69 0.029776055
8 8 8 0.938049252 0.053758549 5375854.92 0.063739373
9 9 9 0.873644859 0.110259617 11025961.72 0.134869298
10 10 10 0.759413841 0.208921016 20892101.59 0.274995872
11 11 11 0.592262406 0.348451123 34845112.32 0.523587821
12 12 12 0.404832556 0.493605796 49360579.64 0.904064409
13 13 13 0.246775115 0.596663100 59666310.01 1.399065649
14 14 14 0.140844104 0.640280357 64028035.69 1.959889434
15 15 15 0.079027762 0.637892516 63789251.61 2.537727076
16 16 16 0.045023619 0.609385210 60938521.00 3.100330120
17 17 17 0.026480508 0.568952991 56895299.12 3.631108086
18 18 18 0.016187644 0.524558151 52455815.09 4.123297431
19 19 19 0.010300836 0.480216722 48021672.18 4.575351559
20 20 20 0.006816356 0.437820069 43782006.91 4.988281758
21 21 21 0.004680357 0.398180705 39818070.49 5.364260026
22 22 22 0.003325880 0.361573066 36157306.56 5.705919073
23 23 23 0.002439311 0.328005643 32800564.26 6.016005000
24 24 24 0.001841429 0.297358959 29735895.87 6.297206097
25 25 25 0.001427187 0.269454558 26945455.77 6.552073884
26 26 26 0.001132925 0.244091829 24409182.93 6.782986791
27 27 27 0.000919160 0.221066317 22106631.66 6.992140296
28 28 28 0.000760611 0.200179592 20017959.20 7.181547756
29 29 29 0.000640756 0.181243797 18124379.70 7.353048604
30 30 30 0.000548620 0.164083668 16408366.79 7.508318993
31 31 31 0.000476687 0.148537362 14853736.20 7.648883138
32 32 32 0.000419738 0.134456259 13445625.89 7.776125612
33 33 33 0.000374079 0.121704426 12170442.58 7.891303040
34 34 34 0.000337044 0.110157887 11015788.67 7.995555207
35 35 35 0.000306690 0.099703792 9970379.18 8.089915243
36 36 36 0.000281582 0.090239555 9023955.53 8.175319348
37 37 37 0.000260634 0.081672010 8167201.00 8.252615783
38 38 38 0.000243024 0.073916619 7391661.88 8.322572798
39 39 39 0.000228115 0.066896700 6689670.04 8.385886240
40 40 40 0.000215413 0.060542716 6054271.63 8.443186422
41 41 41 0.000204527 0.054791662 5479166.23 8.495043880
42 42 42 0.000195150 0.049586460 4958645.96 8.541975097
43 43 43 0.000187035 0.044875398 4487539.81 8.584447684
44 44 44 0.000179982 0.040611644 4061164.38 8.622884946
45 45 45 0.000173829 0.036752784 3675278.45 8.657670058
46 46 46 0.000168443 0.033260416 3326041.59 8.689149856
47 47 47 0.000163712 0.030099767 3009976.67 8.717638278
48 48 48 0.000159545 0.027239357 2723935.67 8.743419468
49 49 49 0.000155866 0.024650687 2465068.73 8.766750605
50 50 50 0.000152609 0.022307960 2230796.05 8.787864455
51 51 51 0.000149721 0.020187824 2018782.36 8.806971681
52 52 52 0.000147154 0.018269139 1826913.86 8.824262946
53 53 53 0.000144870 0.016532773 1653277.26 8.839910802
54 54 54 0.000142833 0.014961408 1496140.82 8.854071415
55 55 55 0.000141014 0.013539371 1353937.13 8.866886115
56 56 56 0.000139388 0.012252476 1225247.57 8.878482807
57 57 57 0.000137933 0.011087882 1108788.22 8.888977245
58 58 58 0.000136629 0.010033970 1003397.05 8.898474185
59 59 59 0.000135460 0.009080224 908022.37 8.907068428
60 60 60 0.000134411 0.008217124 821712.41 8.914845770
61 61 61 0.000133468 0.007436058 743605.78 8.921883851
62 62 62 0.000132621 0.006729229 672922.91 8.928252935
63 63 63 0.000131858 0.006089583 608958.31 8.934016609
64 64 64 0.000131172 0.005510735 551073.50 8.939232416
65 65 65 0.000130555 0.004986907 498690.66 8.943952430
66 66 66 0.000129998 0.004512869 451286.88 8.948223778
67 67 67 0.000129497 0.004083890 408388.95 8.952089106
68 68 68 0.000129045 0.003695686 369568.60 8.955587006
69 69 69 0.000128637 0.003344383 334438.28 8.958752406
70 70 70 0.000128269 0.003026473 302647.27 8.961616909
71 71 71 0.000127937 0.002738781 273878.15 8.964209118
72 72 72 0.000127637 0.002478437 247843.70 8.966554916
73 73 73 0.000127366 0.002242840 224284.00 8.968677726
74 74 74 0.000127122 0.002029638 202963.80 8.970598743
75 75 75 0.000126901 0.001836702 183670.24 8.972337151
76 76 76 0.000126702 0.001662107 166210.67 8.973910307
77 77 77 0.000126521 0.001504108 150410.78 8.975333920
78 78 78 0.000126358 0.001361128 136112.79 8.976622205
79 79 79 0.000126211 0.001231739 123173.95 8.977788026
80 80 80 0.000126078 0.001114650 111465.05 8.978843024
81 81 81 0.000125958 0.001008692 100869.17 8.979797736
82 82 82 0.000125849 0.000912805 91280.52 8.980661694
83 83 83 0.000125751 0.000826034 82603.36 8.981443523
84 84 84 0.000125662 0.000747510 74751.05 8.982151032
85 85 85 0.000125581 0.000676452 67645.17 8.982791284
86 86 86 0.000125509 0.000612148 61214.78 8.983370674
87 87 87 0.000125443 0.000553957 55395.67 8.983894987
88 88 88 0.000125383 0.000501297 50129.72 8.984369458
89 89 89 0.000125329 0.000453643 45364.35 8.984798825
90 90 90 0.000125281 0.000410520 41051.98 8.985187377
91 91 91 0.000125237 0.000371495 37149.54 8.985538992
92 92 92 0.000125197 0.000336181 33618.08 8.985857183
93 93 93 0.000125161 0.000304223 30422.31 8.986145126
94 94 94 0.000125128 0.000275303 27530.34 8.986405697
95 95 95 0.000125099 0.000249133 24913.28 8.986641498
96 96 96 0.000125072 0.000225450 22545.00 8.986854883
97 97 97 0.000125048 0.000204019 20401.85 8.987047984
98 98 98 0.000125026 0.000184624 18462.43 8.987222729
99 99 99 0.000125006 0.000167074 16707.37 8.987380862
100 100 100 0.000124988 0.000151192 15119.17 8.987523961

 

 

 

 

with(flu_out202104_data, plot(time, cases, pch = 19, col = "red", ylim = c(0, 600)))
predictions <- sir_1(beta = 0.004, gamma = 0.5, S = 0.998537432, I = 0.000213768, R = 0.001248782, times = flu_out202104_data$time)
with(predictions, lines(time, I, col = "red"))

 

 

ブログランキング・にほんブログ村へ

RstudioでSIRモデルA関連ページ

感染率による感染者数増加の比較
当サイトは主に物理に関する数学など、その他周辺も含めた少々ごった煮のウェブサイトです。 数学分野に関しての趣旨としては、通常のテキストでは割愛されてしまう内容などを詳しく記述し、さらには難しい説明をするのではなく、わかりにくい内容をいかにわかりやすく伝えるか━など、そういったウェブコンテンツならではの利便性と機動性を生かしたサイト作成を主眼としています。
Rでggplot2によるグラフ描画
当サイトは主に物理に関する数学など、その他周辺も含めた少々ごった煮のウェブサイトです。 数学分野に関しての趣旨としては、通常のテキストでは割愛されてしまう内容などを詳しく記述し、さらには難しい説明をするのではなく、わかりにくい内容をいかにわかりやすく伝えるか━など、そういったウェブコンテンツならではの利便性と機動性を生かしたサイト作成を主眼としています。
RstudioでSIRモデル@
当サイトは主に物理に関する数学など、その他周辺も含めた少々ごった煮のウェブサイトです。 数学分野に関しての趣旨としては、通常のテキストでは割愛されてしまう内容などを詳しく記述し、さらには難しい説明をするのではなく、わかりにくい内容をいかにわかりやすく伝えるか━など、そういったウェブコンテンツならではの利便性と機動性を生かしたサイト作成を主眼としています。

ホーム RSS購読 サイトマップ
TOP 線形代数 ベクトル解析 慣性モーメント 解析力学 微分方程式 NEへの道しるべ