ODE 中的潜在索引问题？

Posted 2023-02-14

【中文标题】ODE 中的潜在索引问题？

【英文标题】：Potential indexing issue in ODE?

【发布时间】：2022-01-22 22:03:22

【问题描述】：

根据@tpetzoldt 的建议，我将在之前的讨论 (Parameter values as a function of another vector. 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

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

【问题讨论】：

感谢您使用更新的代码示例提出新问题，即使标题“索引问题”有些误导。这不是一个“问题”。正确的主题是，如何组织索引（或表查找）。另外，参数的定义确实有问题：rtruncnorm是什么？我猜你的意思是trunc(rnorm())。还有，a=0, b=1 是什么？

【参考方案1】：

上面的代码有几个问题：

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)

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