Возможная проблема с индексированием в ODE?

Question

Следуя предложению @tpetzoldt, я открываю это как вопрос после предыдущего обсуждения (значения параметров как функция другого вектора. DeSolve ).

Я пытаюсь достичь того, чтобы иметь возможность интегрировать модель на каждом временном шаге по вектору DailyTemperature а затем соответствующие значения параметров для каждого дня являются функцией значений из фрейма данных других выходных значений температуры.

      library(deSolve)
set.seed(1)

deriv <- function(t, state, pars) {
  
  pars <- parameters[match(DailyTemperature[floor(t + 1)],parameters$TraitTemperature),2:5]
  #print(pars)
  
  with(as.list(c(state, pars)), {
    d_x <- alpha * x - beta * x * y
    d_y <- delta * beta * x * y - gamma * y
    list(c(x = d_x, y = d_y), alpha=alpha, beta=beta, gamma=gamma, delta=delta)
  })
}


state <- c(x = 1000, y = 10)
times = seq(0, 50, by = 1)

# Parameter datasets 
parameters <- data.frame(
  TraitTemperature = seq(0.1,40,0.1),
  alpha = rtruncnorm(40,a=0,b=1,mean = 1,sd=2),
  beta =  rtruncnorm(40,a=0,b=1,mean = 1,sd=2),
  delta =  rtruncnorm(40,a=0,b=1,mean = 1,sd=2),
  gamma = seq(0.025,1,0.025)
)

# random daily temperature dataset 
DailyTemperature <- round(runif(51, 0, 40),1) # one more because start zero
DailyTemperature

out <- ode(y = state, times = times, func = deriv, parms = pars)
plot(out)
out

На самом деле я начинаю думать, что это проблема со значениями параметров, а не с кодом. В любом случае, мне было бы интересно узнать, правильно ли мое индексирование?

Answer 1

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

• Линия с rtruncnormкажется полностью изуродованным. Здесь я предполагаю, что альфа, бета, гамма линейно зависят от температуры плюс случайная составляющая (которую я в данном случае не очень понимаю), но здесь это не имеет значения, если мы сосредоточимся на программном подходе.
• Состояния и значения параметров относительно велики. Поскольку члены в основном экспоненциальны (лучше всего видно с alpha * xто есть экспоненциальный рост), популяция, вероятно, взорвется или вымрет. Этого можно избежать путем тщательной балансировки параметров и переменных состояния.
• matchмогут возникнуть проблемы из-за ошибок округления, даже если roundа также truncиспользуются. Вместо проверки на равенство обычно лучше проверять минимальное расстояние. Это можно сделать, например, с помощью which.min

Собрав это вместе, мы могли бы сделать это следующим образом:

      library(deSolve)
set.seed(1)

deriv <- function(t, state, pars) {
  pars <- parameters[
    which.min(abs(DailyTemperature[floor(t + 1)] - parameters$TraitTemperature)), 
    1:5
  ]
  with(as.list(c(state, pars)), {
    d_x <- alpha * x - beta * x * y
    d_y <- delta * beta * x * y - gamma * y

    list(c(x = d_x, y = d_y), 
         alpha = alpha, beta = beta, gamma = gamma, delta = delta, 
         ## some extra output for debugging
         DailyTemperature = DailyTemperature[floor(t + 1)],
         TraitTemperature = TraitTemperature
    )
  })
}


state <- c(x = 10, y = 5)
times = seq(0, 150, by = 1)
TraitTemperature = seq(0.1, 40, 0.1)

## here we assume a linear increase of the parameters with temperature
parameters <- data.frame(
  TraitTemperature = TraitTemperature,
  alpha = 1.0 + 0.01 * TraitTemperature + rnorm(40, mean = 0, sd = 0.01),
  beta =  0.2 + 0.01 * TraitTemperature + rnorm(40, mean = 0, sd = 0.01),
  delta = 0.2 + 0.01 * TraitTemperature + rnorm(40, mean = 0, sd = 0.01),
  gamma = seq(0.025, 1, 0.025)
)

## which.min will also work without rounding the temperature

#DailyTemperature <- round(runif(length(times), 0, 40), 1)
DailyTemperature <- runif(length(times), 0, 40)
DailyTemperature

out <- ode(y = state, times = times, func = deriv, parms = pars)
plot(out)
head(out)

Я хочу добавить, что есть и другие (значительно более быстрые) способы поиска в таблице.