처음에

대학의 세미나에서 감염증의 유행 모델인 SIR 모델 에 대해서 조사했으므로, 기사로서 남겨 둔다.
부호는 htps : // 기주 b. 이 m / cy-yon than / Shi R_p t에 놓여 있다.

실행 환경

windows10
Rstudio
R 4.0.1

주제

라이브러리 로드

SIR_plot.R

library(deSolve)
library(ggplot2)

미분 방정식을 나타내기 위해 deSolve, 시각화를 위해 ggplot 라이브러리를 읽습니다.

미분 방정식 및 초기 값 정의

SIR_plot.R

parameters <- c(b = 1/2, g = 1/3)# パラメータのセット
initial <- c(s = 1, i = 1.27*10^(-6), r = 0)# 初期条件
times <- seq(1, 180, 1)# 差分刻み

SIR <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    ds <- -b * s * i
    di <- (b * s - g)*i
    dr <- g * i
    list(c(ds, di, dr))
  })
}
out <- ode(y = initial, times = times, func = SIR, parms = parameters)

1960년대 NY에서 일어난 홍콩 감기 시뮬레이션을 실시한다. 각 파라미터나 초기치는 이 사이트 를 참고로 했다.
시각화를 위해 미분 방정식의 결과를 out에 대입해 둔다.

ggplot으로 시각화

SIR_plot.R

out_df <- as.data.frame(out)
SIRplot_base <- ggplot(data=out_df) +
  xlab("time")+ ylab("Y value") +
  ggtitle("SIR model")

SIRplot <- SIRplot_base +
  layer(
    mapping=aes(x=time, y=s, colour="s"), 
    geom="line", 
    stat="identity", 
    position="identity",
  ) +
  layer(
    mapping=aes(x=time, y=i, colour="i"), 
    geom="line", 
    stat="identity", 
    position="identity",
  ) +
  layer(
    mapping=aes(x=time, y=r,colour="r"), 
    geom="line", 
    stat="identity", 
    position="identity",
  )

미분 방정식의 결과를 data.frame으로 변환하여 S, I, R의 값을 플롯합니다.

결과

끝