1

我想使用 deSolve 创建蝴蝶生态的动态模型。模拟运行了几个模拟年，并且一些事件在一年中的某一天触发（所以我添加了一个状态变量days ）。为了触发这些事件，我想使用一个ifelse语句并且它工作正常，直到我尝试在ifelse语句中放入一个涉及另一个状态变量的操作：D.egg.sus=(ifelse(days<270,(400 * adult.sus),0)). 当我这样做时，模拟运行，但它似乎忽略了该ifelse语句。谁能帮帮我？这是我的完整代码：

days        = 1
egg.sus     = 0
larvae.sus  = 0 
pupae.sus   = 0 
adult.sus   = 1000

state = c(days = days, egg.sus=egg.sus, larvae.sus=larvae.sus,      
pupae.sus=pupae.sus, adult.sus=adult.sus)
model = function(t, state, parameters)
{ 
    with(as.list(c(state, parameters)), 
         {
    D.Days = 1
    D.egg.sus     =
        ( ifelse(days<270, (400*adult.sus) ,0))  ## This is the line causing trouble
        (- egg.sus/5) 
        (-  egg.sus * rbeta(1, 6.038892/5,1.4612593)*.95)                                                                                     
    D.larvae.sus  =
        (+ egg.sus/5) 
        (- larvae.sus * rbeta(1, 0.248531/14,0.2094379)*0.95)
        (- larvae.sus/14)                                                              
    D.pupae.sus   =  
        (+ larvae.sus/14)
        (- pupae.sus * rbeta(1, 0.022011/15, 1.43503))
        (- pupae.sus/15) 
    D.adult.sus   = 
        (+ pupae.sus/15) 
        (- adult.sus/30) 

    list(c( D.Days, D.egg.sus, D.larvae.sus,D.pupae.sus, D.adult.sus))
}
)}

events <- data.frame(var    = c('days'),
                 time   = seq(364,73000,by=365) ,
                 value  = 0,
                 method = "rep")

require(deSolve)

times = seq(1,900, by = 1) 
out = ode(y=state, times = times, func = model, parms = parameters,  events = list(data=events))

dev.cur()
plot(out, col = 2)
4

2 回答 2

1

问题中的模型有几个问题：

  1. 您可以直接使用仿真时间而不是状态变量days，因为函数中的仿真时间为t。然后只需使用模运算符%%，您就不再需要事件了。
  2. 参数都是硬编码的，所以parms=NULL在 ode 函数中使用。
  3. 换行是错误的。当（且仅当）它们在语法上尚未完成时，R会继续行。因此，删除过时的括号，例如，将-运算符放在行尾。
  4. 在 ODE 函数中使用随机数rgamma是一个非常糟糕的主意，特别是对于具有自动时间步长的求解器。ODE 根据定义是确定性的。可以考虑使用固定时间步长求解器，例如method="euler"使用非常小的时间步长或（更好）提供随机值作为外部输入（强制）。
  5. 如果您使用外部输入，则ifelse无论如何都可以避免。
于 2021-01-16T18:48:34.927 回答
1

我不知道大约五年前，但在撰写本文时ifelse使用deSolve. 您的问题似乎是您的条件的返回值没有按您的意愿返回。相反，您可能想要使用标志变量或将返回的值保存ifelse到您可以在模型中使用的变量中。

这是一个小示例，演示如何在模型参数中使用标志

library(deSolve)

# Our model function, first-order
# One parameter is a flag that is used by the ifelse to set Ka to zero if TRUE.

onecomp <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    
    Ka = ifelse(flag == TRUE, 0, Ka) # Use ifelse to check for negative values
    
    dX <- - X*Ka
    dY <- X*Ka - Y*Ke
    list(c(dX, dY))
  })
}

times <- seq(0, 24, by = 0.01)
parameters <- c(Ka = 0.8 , Ke = 0.2, flag = FALSE)
state <- c(X = 100 , Y = 0)

# Test for TRUE
out <- ode(y = state, times = times, func = onecomp, parms = parameters)
plot(out)


# Test for FALSE, where we expect no transfer.
parameters <- c(Ka = 0.8 , Ke = 0.2, flag = TRUE)
out <- ode(y = state, times = times, func = onecomp, parms = parameters)
plot(out)

由reprex 包（v0.3.0）于 2021-01-13 创建

于 2021-01-13T13:57:46.640 回答