установка нижнего предела производной в R desolve - PullRequest
Новая
       49

установка нижнего предела производной в R desolve

0 голосов
/ 20 февраля 2020

Я создаю математическую модель, используя пакет R deSolve. У меня есть три производных: dx (доля зараженных хозяев), dY (доля зараженных векторов) и dm (это отношение векторов к хозяевам в популяции). Цель моей модели - показать влияние определенной обработки инсектицидами на популяцию (с эффектами, представленными как параметр "z". Чтобы включить этот зависимый от времени ковариат в модель, была использована функция приближения. Модель работает правильно , однако я хотел бы установить нижний предел для dm (при условии, что не все векторы в совокупности будут затронуты). Без установки нижнего предела мой код и график выглядят так: 

Initial vectors for days post treatment, % killed
x <- c(4, 30, 60, 90, 120, 210, 360)
z <- c(1.0, 0.99, 0.99, 0.79, 0.7, 0.02, 0) 

plot(z ~ x)

#=============================================
#  fit data with logistic curve
#       -extract fit values using equation y = Asym / (1 + exp((xmid - input) / scal))
#=============================================
fit2 <- nls(z ~ SSlogis(x, Asym, xmid, scal), data = data.frame(x, z))
summary(fit2)

lines(seq(0, 400, length.out = 400),
      predict(fit2, newdata = data.frame(x = seq(0.5, 400, length.out = 400))))

Asym<-summary(fit2)$parameters[1,1]
xmid<-summary(fit2)$parameters[2,1]
scal<-summary(fit2)$parameters[3,1]

times <- seq(0, 1000, by = 1)
signal <- data.frame(times = times, import = rep(0, length(times)))
signal$import=  Asym / (1 + exp((xmid - times) / scal))

#Force time dependent covariate into the model
input <- approxfun(signal, rule = 2)


RMTx2 <- function(times, stateTx2, parametersTx2)   
{
  with(
    as.list(c(stateTx2, parametersTx2)), 
    {
      z <- input(times)
      dX <- ((m*a*b*Y)+(p*k*(a*m*z*Y)))*(1-X)-r*X 
      dY <- a*c*X*(exp(-g*n)-Y)-((g*(1-m/K)*Y)+(m*a*z*Y)) 
      dm <- ((R*(1-m/K)*m )+(-m*a*z))
      return(list(c(dX, dY, dm)))
    }
  )
}


initTx2 <- c(X = 0.01, Y= 0, m=40) 
parametersTx2 <- c(a=1/14, b=0.00068, n=45, g= 0.005, c=0.28, k= 0.10, r= 1/(3*365), p=0, K=40, R= 0.09)
outTx2 <- as.data.frame(ode(y = initTx2, times = times, func = RMTx2, parms = parametersTx2))
RESULTS2<-data.frame(outTx2$X,outTx2$Y)
RESULTS2m <-data.frame(outTx2$m, outTx2$Y*outTx2$Y)
timesTx2 <- seq(0, 10000, by = 1)

И вот график. То, что я хотел бы сделать, это ограничить падение дм во времени, чтобы оно не упало ниже определенного значения во время лечения enter image description here

Я пытаюсь установить параметр dm, чтобы, например, значение не могло опуститься ниже 15. Я попытался сделать несколько кодов, включая: 

RMTx2 <- function(times, stateTx2, parametersTx2)   
{
  with(
    as.list(c(stateTx2, parametersTx2)), 
    {
      z <- input(times)
      dX <- ((m*a*b*Y)+(p*k*(a*m*z*Y)))*(1-X)-r*X 
      dY <- a*c*X*(exp(-g*n)-Y)-((g*(1-m/K)*Y)+(m*a*z*Y))
      dm <- if (isTRUE (((R*(1-m/K)*m )+(-m*a*z)) > MM)) ((R*(1-m/K)*m )+(-m*a*z)) else 15 
      return(list(c(dX, dY, dm)))
    }
  )
}

initTx2 <- c(X = 0.01, Y= 0, m=40) 
parametersTx2 <- c(a=1/14, b=0.00068, n=45, g= 0.005, c=0.28, k= 0.10, r= 1/(3*365), p=0, K=40, R= 0.09, MM= 15)
outTx2 <- as.data.frame(ode(y = initTx2, times = times, func = RMTx2, parms = parametersTx2))
RESULTS2<-data.frame(outTx2$X,outTx2$Y)
RESULTS2m <-data.frame(outTx2$m, outTx2$Y*outTx2$Y)

К сожалению, по какой-то причине это приводит к тому, что население просто увеличиваться до бесконечности: enter image description here

Есть ли что-то фундаментальное в этом подходе, которое не будет к? Или это скорее ошибка кодирования? Спасибо!

1 Ответ

0 голосов
/ 21 февраля 2020

Вы можете попробовать для этого защиту типа Монода. MM - ваш порог, а km - постоянная монода. Его значение довольно некритично, но оно должно быть не слишком маленьким, чтобы избежать численных трудностей для решателя. 

RMTx2 <- function(times, stateTx2, parametersTx2)   
{
  with(
    as.list(c(stateTx2, parametersTx2)), 
    {
      z <- input(times)
      dX <- ((m*a*b*Y)+(p*k*(a*m*z*Y)))*(1-X)-r*X 
      dY <- a*c*X*(exp(-g*n)-Y)-((g*(1-m/K)*Y)+(m*a*z*Y)) 

      #dm <- (R*(1-m/K)*m )+(-m*a*z)
      dm <- (R*(1-m/K)*m )+(-m*a*z) * (m - MM) / (km + (m - MM))  

      list(c(dX, dY, dm), dm=dm) ## return derivative dm as additional output
    }
  )
}


initTx2 <- c(X = 0.01, Y= 0, m=40) 

parameters_old <- c(a=1/14, b=0.00068, n=45, g= 0.005, c=0.28, k= 0.10, r= 1/(3*365), p=0, K=40, R= 0.09, MM = 0, km=0)
parameters_new <- c(a=1/14, b=0.00068, n=45, g= 0.005, c=0.28, k= 0.10, r= 1/(3*365), p=0, K=40, R= 0.09, MM = 15, km=0.1)

out_old <- ode(y = initTx2, times = times, func = RMTx2, parms = parameters_old)
out_new <- ode(y = initTx2, times = times, func = RMTx2, parms = parameters_new)

plot(out_old, out_new)

Simulation without (black) and with (red) safeguard

...