deSolve R 中的时变参数矩阵

Question

我为此苦苦挣扎了这么久。我有一个逻辑增长函数，其中增长参数 r 是一个矩阵。该模型的构建方式是我将 N1 和 N2 作为输出两个 N。

我希望能够随着时间的推移更改 r 参数。当时间 < 50 我想要 r = r1 其中

r1=matrix(c(
  2,3),
  nrow=1, ncol=2

当时间 >= 50 我想要 r=r2 where

r2=matrix(c(
  1,2),
  nrow=1, ncol=2

这是我的功能。非常感谢任何帮助。

rm(list = ls())      
library(deSolve)

model <- function(time, y, params) {
  with(as.list(c(y,params)),{
    N = y[paste("N",1:2, sep = "")]
    
    dN <- r*N*(1-N/K)
    
    return(list(c(dN)))
  })
}

r=matrix(c(
  4,5),
  nrow=1, ncol=2)

K=100

params <- list(r,K)

y<- c(N1=0.1, N2=0.2)

times <- seq(0,100,1)

out <- ode(y, times, model, params)
plot(out)

我最想要这样的东西，但它不起作用

model <- function(time, y, params) {
  with(as.list(c(y,params)),{
    N = y[paste("N",1:2, sep = "")]
    
   r = ifelse(times < 10, matrix(c(1,3),nrow=1, ncol=2),
    ifelse(times > 10, matrix(c(1,4),nrow=1, ncol=2), matrix(c(1,2),nrow=1, ncol=2)))
     
    print(r)
    
    dN <- r*N*(1-N/K)
    
    return(list(c(dN)))
  })
}

感谢您的宝贵时间。

Answer 1

如果你想传递一个矩阵参数，你应该传递一个参数列表，你可以在超过时间限制时在模型中修改它（在下面的例子中你甚至不必传递 r矩阵到模型函数)

library(deSolve)

model <- function(time, y, params) {
  with(as.list(c(y,params)),{
    if(time < 3) r = matrix(c(2,3), nrow = 1, ncol = 2)
    else r = matrix(c(1,3), nrow = 1, ncol = 2)
    N = y[paste("N",1:2, sep = "")]
    
    dN <- r*N*(1-N/K)
    
    return(list(c(dN)))
  })
}

y <- c(N1=0.1, N2=0.2)

params <- list(r = matrix(c(0,0), nrow = 1, ncol = 2), K=100)
times <- seq(0,10,0.1)

out <- ode(y, times, model, params)
plot(out)

你可以看到这方面的例子，例如延迟微分方程 ?dede

Answer 2

这里是使用 approx 函数的扩展版本的通用方法。另请注意模型函数的一些进一步简化和参数值的附加图。

Edit根据Lewis Carter的建议修改了t=3的参数，这样就可以看到效果了

library(simecol) # contains approxTime, a vector version of approx

model <- function(time, N, params) {
    r <- approxTime(params$signal, time, rule = 2, f=0, method="constant")[-1]
    K <- params$K

    dN <- r*N*(1-N/K)
    
    return(list(c(dN), r))
}

signal <- matrix(
  # time, r[1, 2], 
  c(  0, 2, 3,
      3, 1, 2,
    100, 1, 2), ncol=3, byrow=TRUE
) 

## test of the interpolation
approxTime(signal, c(1, 2.9, 3, 100), rule = 2, f=0, method="constant")

params <- list(signal = signal, K = 100)

y <- c(N1=0.1, N2=0.2)

times <- seq(0, 10, 0.1)

out <- ode(y, times, model, params)
plot(out)

对于示例中的少量状态变量，使用包 stats 中带有 approxfun 的单独信号看起来不那么通用，但可能会稍微快一些。

作为进一步的改进，可以考虑用更平滑的过渡替换“硬”过渡。然后可以直接将其表示为一个函数，而不需要 approx、approxfun 或 approxTime.

编辑 2：

Package simecol 导入 deSolve，我们只需要它的一个小函数。因此，除了加载 simecol 之外，还可以在代码中显式包含 approxTime 函数。从数据帧到矩阵的转换提高了性能，但在这种情况下无论如何首选矩阵。

approxTime <- function(x, xout, ...) {
  if (is.data.frame(x)) {x <- as.matrix(x); wasdf <- TRUE} else wasdf <- FALSE
  if (!is.matrix(x)) stop("x must be a matrix or data frame")
  m <- ncol(x)
  y <- matrix(0, nrow=length(xout), ncol=m)
  y[,1] <- xout
  for (i in 2:m) {
    y[,i] <- as.vector(approx(x[,1], x[,i], xout, ...)$y)
  }
  if (wasdf) y <- as.data.frame(y)
  names(y) <- dimnames(x)[[2]]
  y
}

deSolve R 中的时变参数矩阵

Time varying parameter-matrix in deSolve R

r

if-statement

modeling

ode

differential-equations