Ifelse заявление внутри deSolve не работает

Question

Я хочу создать динамическую модель экологии бабочек с использованием 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)

Answer 1

Не знаю как лет пять назад, но на момент написания работало просто отлично с 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)

^{Создано 13 января 2021 г. пакетом reprex (v0.3.0)}

Answer 2

Модель в вопросе имеет несколько проблем:

1. Вы можете использовать время моделирования напрямую вместо переменной состояния days, так как время симуляции в функции задается как t. Затем просто используйте оператор по модулю %%и вам больше не нужны события.
2. все параметры жестко закодированы, поэтому используйте parms=NULLв функции оды.
3. разрывы строк неверны. R продолжает строки, если (и только если) они еще не завершены синтаксически. Поэтому удалите устаревшие скобки и, например, поставьте -оператор в конце строки.
4. Использование случайного числа, например rgammaвнутри функции ОДУ - очень плохая идея, особенно для решателей с автоматическими временными шагами. ОДУ детерминированы по определению. Вместо этого можно рассмотреть решатель с фиксированным временным шагом, например method="euler"с очень маленьким временным шагом или (гораздо лучше) предоставить случайные значения в качестве внешнего ввода (форсирование).
5. Если вы используете внешний вход, вы можете избежать ifelseтем не мение.