Использование пакета desolve - могут ли состояния быть определены матрицей?

Question

Я пытаюсь создать модель с использованием deSolve с довольно большим количеством состояний. Одно из состояний - «foo» - фактически состоит из 15 различных состояний, состоящих из foo [1,1: 5], foo [2,1: 5] и foo [3,1: 5], поэтому я думал, что это будет Проще всего передать функции матрицу состояний, а не вводить их по отдельности, а затем я мог бы обратиться к ним с индексацией:

par <- rep(NA,3) 
    par_names <- c('alpha','prog','death_rate') 
    names(par) <-par_names 
    par['alpha'] <- 0.7 
    par['prog'] <- 0.8 
    par['death_rate'] <- 0.3 

    foo <- matrix(0,nrow = 3,ncol = 5) 

    states <- foo 

    my_func <- function(t,states,par){ 
    with(as.list(c(states,par)),{ 

     for (j in 1:5){ 
     dfoo[1,j] <- par['alpha']*par['prog']*foo[1,j] - par['death_rate']*foo[1,j] 
     dfoo[2,j] <- par['prog']*foo[1,j] - par['prog']*foo[2,j] - par['death_rate']*foo[2,j] 
     dfoo[3,j] <- par['prog']*foo[2,j] - par['prog']*foo[3,j] - par['death_rate']*foo[3,j] 
     } 


     list(c(
     dfoo[] 
    )) 
    }) 
    } 

    times <- seq(1,365,by=1) 

    library(deSolve) 

    alldata <- as.data.frame(ode(y=states,times=times,func=my_func,parms=par)) 

Я попытался исправить это, но я просто продолжаю получать ту же ошибку:

Error in dfoo[1, j] <- par["alpha"] * par["prog"] * foo[1, j] - par["death_rate"] * : 
    object 'dfoo' not found 

Так кто-нибудь знает, как это может быть сделано для работы или более простой способ сделать это?

Answer 1

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

Вы используете ненужные искажения для создания своих данных. Кроме того, как указано в комментариях, ваша функция, похоже, не инициализирует dfoo. Наконец, ваш цикл for в функции можно было бы более тщательно обрабатывать несколькими векторизованными операциями. Вот пример:

my_func <- function(t,states,par){ 
    foo <- matrix(states, nrow = 3, ncol = 5) 
    dfoo <- with(as.list(par), rbind( 
     (prog * alpha * foo[1,]) - (death_rate * foo[1,]), 
     (prog * foo[-nrow(foo),]) - (prog * foo[-1,]) - (death_rate * foo[-1,]) 
    )) 
    list(dfoo) 
} 

library(deSolve) 
par <- c(alpha = 0.7, prog = 0.8, death_rate = 0.3) 
states <- matrix(0,nrow = 3,ncol = 5) 
ode <- ode(y=states, times=1:365, func=my_func, parms=par) 
alldata <- as.data.frame(ode)