Использование цикла в ODE для графического сравнения различных параметров R

Я использую пакет deSolve для построения пары дифференциальных уравнений (если интересно, прочтите http://www.maa.org/press/periodicals/loci/joma/the-sir-model-for-spread-модели-болезни-дифференциального-уравнения).

Моя конечная цель - создать итеративную функцию или процесс (для цикла), чтобы построить график того, как изменения определенных параметров (бета и гамма) повлияют на решение. Предпочтительным выходом будет список, содержащий все ode решения для каждого указанного значения бета в цикле. У меня возникают проблемы с интеграцией цикла в настройку, которая требуется пакету deSolve для функции ode.

В приведенном ниже коде я пытаюсь изобразить, как диапазон значений (от 1 до 2 с шагом 0,1) в параметре beta повлияет на график дифференциальных уравнений.

for(k in seq(1,2,by=0.1)){ #range of values for beta

    init <- c(S=1-1e-6, I=1e-6, R=0) #initial conditions for odes
    time <- seq(0,80,by=1) #time period
    parameters <- c(beta=k, gamma=0.15) #parameters in ode

SIR <- function(time,state,parameters){ #function containing equaations
    with(as.list(c(state,parameters)),{
        dS <- -beta*S*I
        dI <- beta*S*I-gamma*I
        dR <- gamma*I

        return(list(c(dS,dI,dR)))
    })
}

ode(y=init,times=time,func=SIR()[beta],parms=parameters[k])}

}

Первая ошибка, которую я получаю, гласит, что параметры аргумента в функции SIR отсутствуют.

Ошибка в as.list (c (init, parameters)): аргумент "параметры" отсутствует, значение по умолчанию отсутствует.

Я не понимаю, почему сообщается об этой ошибке, когда я назначил parameters в предыдущих строках.


r ode
person EJJ    schedule 25.12.2016    source источник

Ответы (1)


arrow_upward
7
arrow_downward

Вы также можете определить свою функцию градиента (и другие неизменяемые элементы) вне цикла:

SIR <- function(time,state,parameters) {
  with(as.list(c(state,parameters)),{
    dS <- -beta*S*I
    dI <- beta*S*I-gamma*I
    dR <- gamma*I
    return(list(c(dS,dI,dR)))
  })
}
init <- c(S=1-1e-6, I=1e-6, R=0) #initial conditions for odes
time <- seq(0,80,by=1) #time period

Теперь определите вектор значений, который нужно попробовать (не обязательно, но удобно):

betavec <- seq(1,2,by=0.1)

и определите список для хранения результатов:

res <- vector(length(betavec),mode="list")

library(deSolve)
for (k in seq_along(betavec)){ #range of values for beta
    res[[k]] <- ode(y=init,times=time,func=SIR,
           parms=c(beta=betavec[k], gamma=0.15))
}

Теперь у вас есть список, каждый элемент которого содержит результаты одного запуска. Вы можете sapply или lapply над этим списком, например чтобы получить матрицу последних состояний каждого прогона:

t(sapply(res,tail,1))

Или, если вы хотите получить результаты в виде одного длинного фрейма данных ...

names(res) <- betavec  ## to get beta value incorporated in results
dd <- dplyr::bind_rows(lapply(res,as.data.frame),.id="beta")
dd$beta <- as.numeric(dd$beta)

do.call(rbind,...) будет работать почти так же хорошо, как bind_rows(), но аргумент bind_rows .id удобен для добавления значений beta в каждый фрейм данных. Вы также можете оставить результаты в виде списка и перебирать их при построении графика с отдельными вызовами lines() или (например) связать вместе только инфекционные столбцы и использовать matplot() для построения их всех одновременно. Это просто вопрос стиля и идиомы.

library(ggplot2); theme_set(theme_bw())
library(viridis)
ggplot(dd,aes(x=time,y=I,colour=beta))+
    geom_line(aes(group=beta))+
    scale_color_viridis()+
    scale_y_log10()

введите здесь описание изображения

person Ben Bolker    schedule 25.12.2016
comment
Спасибо за ответ! Я не знаком с пакетом dplyr, похоже, он делает то, что сделал бы rbind. Я предполагаю, что с использованием фрейма данных и списка было бы легче построить различные бета-версии? - person EJJ; 26.12.2016