Значения параметров как функция другого вектора. deSolve

Question

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

Простая модель

       library(deSolve)
set.seed(1)

pars <- c(alpha = 1, beta = 0.2, delta = 0.5, gamma = 0.2)



lv_model <- function(pars, times = seq(0, 50, by = 1)) {
  # initial state 
  state <- c(x = 1, y = 2)
  # derivative
  deriv <- function(t, state, pars) {
    with(as.list(c(state, pars)), {
      d_x <- alpha * x - beta * x * y
      d_y <- delta * beta * x * y - gamma * y
      return(list(c(x = d_x, y = d_y)))
    })
  }
  # solve
  ode(y = state, times = times, func = deriv, parms = pars)
}
lv_results <- lv_model(pars = pars, times = seq(0, 50, by = 1))

Теперь я хочу использовать последовательность дневных температур. DailyTemperature<-floor(runif(50,0,40))

и сделаем значения параметров функциями температур

      TraitTemperature<-seq(1,40,1)

#trait responses to temperature
alpha<- abs(rnorm(40,mean = 0.5,sd=1))
beta<- abs(rnorm(40,mean = 0.2,sd=0.5))
delta<-abs(rnorm(40,mean=1,sd=2))
gamma<- seq(0.025,1,0.025)
parameters<-as.data.frame(cbind(TraitTemperature,alpha,beta,delta,gamma))

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

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

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

Answer 1

Здесь возможный подход с использованием в качестве форсирования, а затем фрейма данных в качестве таблицы поиска. Это не нужно matchздесь до тех пор, пока температуры являются целыми числами, а температуры во фрейме данных соответствуют индексу строки фрейма данных.

Я хочу добавить, что в принципе разрывные воздействия делают модель медленной, потому что ОДУ по определению является непрерывным. К счастью, поскольку решатели довольно надежны, для практических приложений это должно:

      library(deSolve)
set.seed(1)

deriv <- function(t, state, pars) {

  pars <- parameters[DailyTemperature[floor(t + 1)], 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 = 1, y = 2)
times = seq(0, 50, by = 1)

# pars <- c(alpha = 1, beta = 0.2, delta = 0.5, gamma = 0.2)
parameters <- data.frame(
  TraitTemperature = seq(1,40,1),
  alpha = abs(rnorm(40,mean = 0.5,sd=1)),
  beta = abs(rnorm(40,mean = 0.2,sd=0.5)),
  delta = abs(rnorm(40,mean=1,sd=2)),
  gamma = seq(0.025,1,0.025)
)


DailyTemperature <- floor(runif(51, 0, 40)) # one more because start zero

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

Параметры в виде списка

В приведенном выше примере переменная передается из odeпросто перезаписывается с parsпроизводные от глобальных переменных parametersа также DailyTemperature. Это работает, но можно также рассмотреть возможность передачи обоих как списка в derivфункция.

      deriv <- function(t, state, p) {
  
  parameters <- p[[1]]
  DailyTemperature <- p[[2]]

  parms <- parameters[DailyTemperature[floor(t + 1)], 2:5]
  # ...
}

а потом:

      out <- ode(y = state, times = times, func = deriv,
  parms = list(parameters, DailyTemperature))

Answer 2

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

Существует несколько методов интерполяции или индексации в 2D (время и температура). Я бы предпочел создать 2D-модель, а затем использовать метод 2D-интерполяции. Лучше всего это работает, если поверхность параметров будет гладкой, а не просто случайной, как в данном примере. Однако для упрощения можно также использовать округление и поиск по таблице. Если значения не являются целыми числами, точное сравнение часто не работает из-за эффектов округления и ограниченной точности (числовые форматы IEEE), поэтому вместо совпадения можно использовать :

      DailyTemperature <- round(runif(51, 0, 40), 1)

TraitTemperature <- seq(0, 40, by=0.1)
N <- length(TraitTemperature)

parameters <- data.frame(
  TraitTemperature = TraitTemperature,
  alpha = abs(rnorm(N, mean = 0.5, sd=1)),
  beta = abs(rnorm(N, mean = 0.2, sd=0.5)),
  delta = abs(rnorm(N, mean=1,sd=2)),
  gamma = seq(0.025, 1, length.out=N)
)


t <- 17
actualTemp <- DailyTemperature[floor(t+1)]
actualTemp
pars <- parameters[which.min(abs(actualTemp - parameters$TraitTemperature)), 1:5]

head(pars)

Answer 3

Кажется, это работает для индексации с использованием текущего воспроизводимого кода:

      set.seed(1)
deriv <- function(t, state, pars) {
  pars<- parameters[match(parameters$TraitTemperature[parameters[2:5]],DailyTemperature),]
  #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 = 1, y = 2)
times = seq(0, 50, by = 1)

# pars <- c(alpha = 1, beta = 0.2, delta = 0.5, gamma = 0.2)
parameters <- data.frame(
  TraitTemperature = seq(1,40,1),
  alpha = abs(rnorm(40,mean = 0.5,sd=1)),
  beta = abs(rnorm(40,mean = 0.2,sd=0.5)),
  delta = abs(rnorm(40,mean=1,sd=2)),
  gamma = seq(0.025,1,0.025)
)

DailyTemperature <- floor(runif(51, 0, 40)) # one more because start zero

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

но если я увеличу разрешение значений, например

      
# Parameter datasets 
parameters <- data.frame(
  TraitTemperature = seq(0.1,40,0.1),
  alpha = abs(rnorm(400,mean = 0.5,sd=1)),
  beta = abs(rnorm(400,mean = 0.2,sd=0.5)),
  delta = abs(rnorm(400,mean=1,sd=2)),
  gamma = seq(0.0025,1,0.0025)
)

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

Затем я получаю NA через определенные промежутки времени.