時変パラメータを変化させて複数回実行する。

Question

このコードを繰り返し実行して、実行ごとに異なる列変数を持つ単一の出力データセットを作成したいとします。今、コードが動作し、さまざまな時に異なるイベントを挿入することができます。しかし、私は、イベントの大きさを変更できるようになど0.4、0.5、0.6、0.35を変える

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

をしたいと思い、私はループのために試してみましたが、それはすべてで動作させることができませんでした。私のコードは以下の通りです：

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) 

最後のデータセットには、イベントの値が異なる実行ごとに1つずつ、複数の前の列があります。

ありがとうございます、ありがとうございます！

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 

その後、我々は... 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) 
}) 
} 

をsir1機能を調整することができます...そして、我々はまた、モデルを実行するループでparamsベクトルを変更することができ、取得しますprevを計算し、それをout data.frameに格納します。

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) 
} 

あなたがしたいことを行うためのより良い方法がありますが、これは潜在的な解決策です。