run deSolve несколько раз, изменяя изменяющийся во времени параметр

Question

Я хотел бы, чтобы этот код выполнялся многократно, создавая единый выходной набор данных с другой переменной столбца для каждого запуска. Прямо сейчас, код работает и позволяет мне вставлять разные события в разное время. Тем не менее, я хотел бы быть в состоянии изменить величину события,

IPT <- ifelse (t<210, IPT, 0.35*exp(-(t-209)/21)) 

варьируя от 0,35 до 0,4, 0,5, 0,6 и т.д. Я попытался Для петель, но не мог заставить его работать. Мой код ниже:

library(deSolve) 
##Simple parameter list 
params <- c(b = 0.477, bs = .4, bsv = 0.1, nets = 0.4767, betah = 0.2, 
     rhos = 179, Bthetas = 0.2, psi = 14,phis = 0.5, gamma =14, 
     thetas = 0.79,piv = 1/19, betav = 0.09122, nu = 0.2085, sigma = 12, 
     muv = 1/19, IPT = 0, IPT2 = 0, IPT3 = 0) 
dt  <- seq(0, 5000, 7) 
inits  <- c(Ss = 30000, Is = 0, As = 0, Rs = 0, 
      Sv = 29999, Ev = 0, Iv = 1) 
Nh <- 30000 
Nv <- 30000 


## Create an SIR function 
sir1 <- function(t, x, params) { 

with(as.list(c(params, x)), { 

IPT <- ifelse (t<210, IPT, 0.35*exp(-(t-209)/21)) 

dSs <- -((b*bs*(1-nets))+(b*bs*nets*0.78))*betah*Iv*Ss   /Nh         + As*(1/rhos)*(1-Bthetas)                   + Rs*(1/psi) 
dIs <- ((b*bs*(1-nets))+(b*bs*nets*0.78))*betah*Iv*Ss*(1-phis)/Nh - 1/gamma * Is               - Is*(IPT + IPT2 + IPT3) 
dAs <- ((b*bs*(1-nets))+(b*bs*nets*0.78))*betah*Iv*Ss*( phis)/Nh + 1/gamma * Is*(1-thetas) - As*(1/rhos)*(1-Bthetas) - As*(2/rhos)*Bthetas       - As*(IPT + IPT2 + IPT3) 
dRs <-                1/gamma * Is*( thetas)       + As*(2/rhos)*Bthetas + Is*(IPT2 + IPT3+ IPT) + As*(IPT + IPT2 + IPT3) - Rs*(1/psi) 

dSv <- piv*Nv - Sv*betav*b*(nu*(
    ((bsv*(1-nets))+(bsv*nets*0.78))*As)+ 
    ((bsv*(1-nets))+(bsv*nets*0.78))*Is/Nh) - Sv*muv 
dEv <-   Sv*betav*b*(nu*(
    ((bsv*(1-nets))+(bsv*nets*0.78))*As)+ 
    ((bsv*(1-nets))+(bsv*nets*0.78))*Is/Nh) - Ev*(1/sigma + muv) 
dIv <-  Ev*(1/sigma)- Iv*muv 


der <- c(dSs, dIs, dAs, dRs, 
     dSv, dEv, dIv) 
list(der) 
}) 
} 

out <- as.data.frame(lsoda(inits, dt, sir1, parms = params)) 

out$prev <- with(out, Is+As/Nh) 

Я хотел бы окончательный набор данных, чтобы иметь несколько столбцов ПРЕД, по одному для каждого прогона с различными значениями события.

Любая помощь будет оценена, спасибо!

Answer 1

Потенциальным решением было бы иметь величину вместо параметра константа (здесь я называю это mag).

library(deSolve) 
##Simple parameter list 
params <- c(b = 0.477, bs = .4, bsv = 0.1, nets = 0.4767, betah = 0.2, 
     rhos = 179, Bthetas = 0.2, psi = 14,phis = 0.5, gamma =14, 
     thetas = 0.79,piv = 1/19, betav = 0.09122, nu = 0.2085, sigma = 12, 
     muv = 1/19, IPT = 0, IPT2 = 0, IPT3 = 0, mag=0.35) 
dt  <- seq(0, 5000, 7) 
inits  <- c(Ss = 30000, Is = 0, As = 0, Rs = 0, 
      Sv = 29999, Ev = 0, Iv = 1) 
Nh <- 30000 
Nv <- 30000 

Тогда мы можем настроить функцию sir1 принять параметр mag ...

## Create an SIR function 
sir1 <- function(t, x, params) { 

with(as.list(c(params, x)), { 

IPT <- ifelse (t<210, IPT, mag*exp(-(t-209)/21)) 

dSs <- -((b*bs*(1-nets))+(b*bs*nets*0.78))*betah*Iv*Ss   /Nh         + As*(1/rhos)*(1-Bthetas)                   + Rs*(1/psi) 
dIs <- ((b*bs*(1-nets))+(b*bs*nets*0.78))*betah*Iv*Ss*(1-phis)/Nh - 1/gamma * Is               - Is*(IPT + IPT2 + IPT3) 
dAs <- ((b*bs*(1-nets))+(b*bs*nets*0.78))*betah*Iv*Ss*( phis)/Nh + 1/gamma * Is*(1-thetas) - As*(1/rhos)*(1-Bthetas) - As*(2/rhos)*Bthetas       - As*(IPT + IPT2 + IPT3) 
dRs <-                1/gamma * Is*( thetas)       + As*(2/rhos)*Bthetas + Is*(IPT2 + IPT3+ IPT) + As*(IPT + IPT2 + IPT3) - Rs*(1/psi) 

dSv <- piv*Nv - Sv*betav*b*(nu*(
    ((bsv*(1-nets))+(bsv*nets*0.78))*As)+ 
    ((bsv*(1-nets))+(bsv*nets*0.78))*Is/Nh) - Sv*muv 
dEv <-   Sv*betav*b*(nu*(
    ((bsv*(1-nets))+(bsv*nets*0.78))*As)+ 
    ((bsv*(1-nets))+(bsv*nets*0.78))*Is/Nh) - Ev*(1/sigma + muv) 
dIv <-  Ev*(1/sigma)- Iv*muv 


der <- c(dSs, dIs, dAs, dRs, 
     dSv, dEv, dIv) 
list(der) 
}) 
} 

... и мы можем изменить params вектор в цикле, который также запускает модель, получает выводит prev и сохраняет его в файле данных out.

out <- as.data.frame(lsoda(inits, dt, sir1, parms = params)) 
magz <- seq(0.2, 0.5, length.out=10) 

for(i in 1:length(magz)){ 
    params['mag'] <- magz[i] 
    tmp <- as.data.frame(lsoda(inits, dt, sir1, parms = params)) 
    nm <- paste('prev', round(params['mag'],2), sep='') 
    out[,nm] <- with(tmp, Is+As/Nh) 
} 

Есть, вероятно, лучшие способы сделать то, что вы хотите сделать, но это потенциальное решение.