作者anovachen (囧)
看板R_Language
标题[问题] 以R软体解ODE (传染病数学模型)
时间Sat Apr 20 12:50:12 2013
[问题类型]:
程式谘询(我想用R 做某件事情，但是我不知道要怎麽用R 写出来)
[软体熟悉度]:
入门(写过其他程式，只是对语法不熟悉)
[问题叙述]:
某传染病流行病学课本的公式：
n: size of the population
S: proportion of n that is susceptible
I: proportion of n that is currently infected and infectious
R: proportion of n that is immune.
D: 从感染到痊癒所需的时间
β：接触患者後感染的机率
κ：单位时间内平均一个人会接触几个人
S_t+1 = S_t - βκS_tI_t
I_t+1 = I_t + βκS_tIt - I_t/D
R_t+1 = R_t + I_t/D
(↑这些叫做差分方程，对吧？)
dS/dt = - βκS_tI_t
dI/dt = βκS_tI_t - I_t/D
dR/dt = I_t/D
(↑要解出给定时间的S,I,R，应该就是解微分方程？)
请问如果要把dS/dt, dI/dt, dR/dt在同一张图上作图，该怎麽做？
另外还有S, I, R对时间作图。
最主要是想了解dI/dt是否会在某段时间之後趋近於0。
以下程式执行结果跟预期的差很多...
不知道是不是function没写对？
[程式范例]:
课本给的参数是beta=0.15(/人), kapp=12(人/周), D=1(周), 时间单位是周
初始条件：族群总人数为1000人，有1人感染。
http://pastie.org/7682058
library(deSolve)
pars <- c(beta=0.15, kappa=12, D=1, pop_size=1000)
yini <- c(I=1, S=999, R=0)
times <- seq(0, 25, by=1)
infection <- function(Time, State, Pars) {
with(as.list(c(State, Pars)), {
dS <- -beta*kappa*(S/pop_size)*(I/pop_size)
dI <- beta*kappa*(S/pop_size)*(I/pop_size)-((I/pop_size)/D)
dR <- (I/pop_size)/D
return(list(c(dS, dI, dR)))
})
}
out<- ode(func=infection,y=yini,parms=pars,times=times)
matplot(out[,"time"], out[,2:4], type = "l", xlab = "week", ylab = "Rate",
main = "Infection", lwd = 2)
legend("topright", c("dS/dt", "dI/dt", "dR/dt"), col = 1:3, lty = 1:2)
[关键字]:
R, ODE, infection epidemiology, mathematical model
PS:
我自己用excel做的简易计算结果在这边(主要是要看dI/dt对时间作图的粗略结果)
，可以参考看看。
http://www.filethief.com/download/1983/epi_math.xlsx.html
感谢指教！
--
※ 发信站: 批踢踢实业坊(ptt.cc)
◆ From: 111.255.5.34
1F：→ obarisk:感觉像chaos 04/20 13:47
2F：→ lin15:R好像有做SIR model的package 04/20 13:51
3F：推 clickhere:你的input和output放错了, 所以解得不是你给的ode. 04/21 11:28
4F：→ clickhere:虽然ode()可以给变数名,但不代表fortran就懂名字. 04/21 11:29
5F：→ clickhere:fortran只管变数的顺序的.你input yini 是 I, S, R. 04/21 11:29
6F：→ clickhere:但在ode()里却回传 c(dI, dS, dR). 04/21 11:30
7F：→ clickhere:改 yini <- c(S = ..., I= ..., R = ...) 即可. 04/21 11:32
8F：→ clickhere: 回传 c(dS, dI, dR). 04/21 11:33
感谢各位回答！
刚开始先改yini的顺序，但还是执行不出预期结果。
後来参考其他人用R语言实作SIR model的程式，
http://archives.aidanfindlater.com/blog/2010/04/20/the-basic-sir-model-in-r/
发现似乎直接把S,I,R定义成proportion就好...
(也就是要把S,I定成999/1000,1/1000)
(如我上述的S,I,R定义，这三个变数其实是proportion)
我一开始是把S,I,R定义成人数，再让它除以pop_size(族群总人口数)，
dS <- -beta*kappa*(S/pop_size)*(I/pop_size)
dI <- beta*kappa*(S/pop_size)*(I/pop_size)-((I/pop_size)/D)
dR <- (I/pop_size)/D
这样是不符合SIR model公式的。
以下ODE解出来的SIR都是proportion，
但是最後会把out矩阵中的SIR值再乘上pop_size，
以求出人数。
http://pastie.org/7682055
library(deSolve)
pop_size=1000
pars <- c(beta=0.15, kappa=12, D=1)
yini <- c(S=0.999, I=0.001, R=0)
times <- seq(0, 25, by=1)
infection <- function(Time, State, Pars) {
with(as.list(c(State, Pars)), {
dS <- -beta*kappa*S*I
dI <- beta*kappa*S*I-(I/D)
dR <- I/D
return(list(c(dS, dI, dR)))
})
}
out<- ode(func=infection,y=yini,parms=pars,times=times)
out[,2:4] <- pop_size*out[,2:4]
matplot(out[,"time"], out[,2:4], type = "l", xlab = "week", ylab = "numbers",
main = "Infection", lwd = 2)
legend("topright", c("Susceptibles", "Infecteds", "Recovereds"), col = 1:3,
lty = 1:2)
不过，现在还有个问题是...
怎麽作出dS/dt, dI/dt, dR/dt对时间作图的结果呢？
因为在这个课本的范例中，
大约在第二十周之後感染者人数几乎是零人，
而且从SIR的图也可看出，感染者增加的速率是先上升後下降的。(由切线斜率判断)
也就是说...我似乎需要先写出d^2X/dt^2 (X=S,I,R)的function才行？
这下变成纯粹的数学问题了= ="
I''= βκ(S'I+SI') –I'/D
这样写对吗？
有办法用R解这种型态的ODE吗？
9F：→ clickhere:有没有除pop_size只影响解的范围.基本上,是顺序解错了. 04/22 00:37
呃...会有影响。
因为这公式对S,I,R的定义就是proportion，
所以要写成d(S/pop_size)/dt才算得出来？
10F：→ clickhere:你都已有out[,c(2,3,4)]了.何需在解ODE? 04/22 00:38
11F：→ clickhere:out[-1,2] - out[-nrow(out),2] 不是 dS/dt 吗? 04/22 00:39
感谢指教...
可是out[-1,2] - out[-nrow(out),2]无法把times这个column留下来？
out<- ode(func=infection,y=yini,parms=pars,times=times)
out[,2:4] <- pop_size*out[,2:4]
dout<- out[-1,1:4] - out[-nrow(out),1:4]
week <- times[2:26]
dim(week) <-c(25,1)
#因为times都变成1，只好设定一个叫做week的矩阵代表时间
matplot(week, dout[,2:4], type = "l", xlab = "week", ylab = "numbers", main =
"Infection", lwd = 2)
legend("topright", c("dS/dt", "dI/dt", "dR/dt"), col = 1:3, lty = 1:2)
12F：推 clickhere:proportion!! scale initial kappa & D即可. 04/22 18:47
感谢回答！
※ 编辑: anovachen 来自: 111.255.244.92 (04/23 23:15)