使用 optim 从具有不同常量的多个时间序列中估计参数

Question

我正在尝试使用多个时间序列来估计模型参数，其中序列之间的常数值不同。为了便于解释，我将使用逻辑增长模型作为示例。

我已经能够从具有相同常量 K 的多个时间序列（N1 和 N2）中估计参数 (r)。

#Data; needs to be matrix instead of data frame
 dat <- as.matrix(cbind(time=seq(0,166,by=16),
         N1=c(0.020,0.030,0.060,0.100,0.160,0.26,0.360,0.50,0.70,0.800,0.90),
         N2=c(0.015,0.033,0.062,0.106,0.162,0.26,0.306,0.51,0.76,0.821,0.91)))

 #dynamical model to estimate r
 dNdt.model=function(t,x,params){
   N <- x
   with(as.list(c(params)), {
     dN <- r*N*(1-(N/K))
     list(c(dN))
   })
 }

 #sse objective function
 sse.dNdt=function(dNparams, data, Kfix=1){
   t <- data[,1]
   N <- data[,-1]
   N0 <- data[1,-1] 
   K <- Kfix
   r0 <- dNparams
   out <- as.data.frame(ode(y=N0, times=t, func=dNdt.model, parms=c(r=r0, K=K)))
   sse <- sum((out[2]-N)^2) + sum((out[3]-N)^2) #SSE needs to be sum of all trajectories
 }

 #run optim
 dNparams <- c(.3) #initial value of r
 optim(dNparams, sse.dNdt, data=dat) #estimate r based on N1 and N2 trajectories

哪个returns预期输出：

$平价 [1] 0.0637207

$值 [1] 4.279062

$ 计数函数梯度 32 不适用

$收敛 [1] 0

$消息空

如何转换此代码以对每个时间序列采用不同的 K 值？将 Kfix 转换为包含 2 个值的向量 returns 出错：未找到对象 'K'。

Answer 1

这是一个解决方案，它使用 purrr::map2_dfc 到运行每个时间序列的模型和 K 值以及 returns 数据框，其中每一列对应一个时间系列。

# Load package
library(deSolve)

# Data for fitting
dat <- data.frame(time = seq(0, 166, by=16),
                  N1 = c(0.020, 0.030, 0.060, 0.100, 0.160, 
                         0.26, 0.360, 0.50, 0.70, 0.800, 0.90),
                  N2 = c(0.015, 0.033, 0.062, 0.106, 0.162, 
                         0.26, 0.306, 0.51, 0.76, 0.821, 0.91))

# Model
dNdt_model <- function(t, x, params){
  # State variable
  N <- x

  with(as.list(c(params)), {
    # Rate of change
    dN <- r*N*(1-(N/K))

    #Return list
    list(c(dN))
  })
}

下面的代码是唯一真正改变的部分（除了我强迫性的整理）。 purrr::map2_dfc 遍历每个初始条件和 K 值以及运行模型。结果被绑定到一个数据框中。最后，从模型结果中减去观测值并计算平方和。

#sse objective function
sse_dNdt <- function(dNparams, data, K){
  # Output times
  t <- data[,1]

  # Time series
  N <- data[,-1]

  # Run for each time series (different initial condition, different K value)
  out <- purrr::map2_dfc(.x = data[1,-1], .y = K, 
                         ~ ode(y = .x, times = t, func = dNdt_model, parms = c(r = dNparams, K = .y))[, 2])

  #SSE needs to be sum of all trajectories
  sum((out - N)^2) 
}

但是，purrr::map2_dfc 是做什么的？正如其名称所示，它映射两个变量 (map2)，即初始条件和 K 值，以及 returns 一个数据框 (df)，即通过绑定列生成（c，而不是行的 r）。

# Initial value of r
dNparams <- c(.3)

下面，我用 K 个值的数组调用 optim。

# Run with different K values for time series N1 and N1
optim(dNparams, sse_dNdt, data = dat, K = c(1, 2))

使用 optim 从具有不同常量的多个时间序列中估计参数

Using optim to estimate parameters from multiple time series with different constants

r

ode