0

遵循@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

我实际上开始认为这是参数值而不是代码的问题。无论如何，我很想知道我的索引是否正确？

4

1 回答 1

0

上面的代码有几个问题：

  • 这条线rtruncnorm似乎完全被破坏了。在这里，我假设 alpha、beta、gamma 线性依赖于温度加上一个随机分量（在这种情况下我并不真正理解），但如果我们专注于编程方法，这并不重要。
  • 状态和参数值都比较大。由于这些术语本质上是指数的（最好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)

我想补充一点，还有其他（相当快的）查表方法。

于 2021-12-21T21:20:12.520 回答