Реплицируйте процедуру подгонки с помощью функции ‘optim’ в R с пакетом ‘fitode’, добавляя известные времязависимые переменные.

Несколько лет назад я опубликовал следующий вопрос. Я пытаюсь воспроизвести ответ, предоставленный через пакет R ‘fitode’, добавляя к формулировке модели две временно зависимые переменные, представляющие собой ставки вылова.

Модель является простой модифицированной моделью Лотка-Вольтерры для двух видов следующего вида:

# Формулировка модели
inv <- 1e-5 # Числовая погрешность, чтобы избежать отрицательных биомасс.
HS2LLV <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    H1 <- H1(time)
    H2 <- H2(time)
    db1.dt = b1*(r1+a11*b1+a12*(b2^2))-H1*b1+inv
    db2.dt = b2*(r2+a22*b2+a21*(b1^2))-H2*b2+inv
    list(c(db1.dt, db2.dt))
  })
}

Здесь ‘H1’ и ‘H2’ представляют собой ставки вылова, которые являются функциями времени, определенными в R через ‘approxfun’, для которых у меня есть следующие данные:

Fdata <- data.frame(Herring = c(0.1603, 0.17, 0.1579, 0.1469,0.1263, 0.1533,
                                0.1593, 0.1839, 0.1648, 0.2255, 0.2369, 0.2524,
                                0.2265, 0.2626, 0.2547, 0.3426, 0.3324, 0.3416,
                                0.299, 0.3367, 0.4076, 0.3872, 0.4094, 0.4564,
                                0.4848, 0.4127, 0.49, 0.4336, 0.3995, 0.3179,
                                0.2713, 0.2464, 0.2712, 0.2729, 0.2744, 0.2485,
                                0.2943, 0.2299, 0.1658, 0.1492, 0.2099, 0.2886,
                                0.3496, 0.352, 0.4474, 0.4951, 0.46),
                    Sprat = c(0.3696, 0.396, 0.3777, 0.336, 0.3408, 0.2627,
                              0.3005, 0.1828, 0.3033, 0.1358, 0.1796, 0.1659,
                              0.2091, 0.2689, 0.2391, 0.2155, 0.1369, 0.1734,
                              0.2004, 0.1625, 0.2638, 0.3485, 0.2984, 0.416,
                              0.3981, 0.37, 0.3193, 0.2866, 0.4008, 0.3958,
                              0.4933, 0.4536, 0.3848, 0.3793, 0.4115, 0.5091,
                              0.4403, 0.3968, 0.3462, 0.4192, 0.4341, 0.4538,
                              0.3455, 0.3961, 0.395, 0.3918, 0.368))
rownames(Fdata) <- 1974:2020

Здесь ‘H1’ и ‘H2’ действуют как классические временно зависимые события в моделях ‘deSolve’.

Я подгоняю эту модель через ‘optim’ к следующим наблюдаемым данным SSB:

# Данные оценки из Балтийского моря
ObsSSB <- data.frame(Herring = c(1932041, 1864349, 1672273, 1944689, 1905001,
                                 1814679, 1626604, 1438139, 1520902, 1421057,
                                 1266502, 1177136, 1090921, 1011718, 1013695,
                                 856321, 714642, 647277, 675487, 649195, 651524,
                                 539582, 483859, 452706, 417827, 363049, 354724,
                                 339899, 332746, 367432, 377935, 424888, 461102,
                                 478633, 477803, 536040, 565200, 557749, 598115,
                                 629339, 695343, 643874, 579087, 597533, 581015,
                                 460378, 364981),
                     Sprat = c(940000, 726000, 625000, 1044000, 695000, 377000,
                               227000, 199000, 254000, 394000, 616000, 605000,
                               570000, 461000, 403000, 423000, 556000, 775000,
                               1045000, 1360000, 1375000, 1429000, 1811000,
                               1777000, 1354000, 1353000, 1319000, 1196000,
                               942000, 829000, 1040000, 1324000, 1055000,
                               898000, 921000, 827000, 948000, 752000, 694000,
                               706000, 620000, 680000, 1077000, 1095000, 982000,
                               855000, 817000))
rownames(ObsSSB) <- 1974:2020

Процедура подгонки начинается с обоснованного предположения начальных параметров:

### Подгонка модифицированной модели Лотка-Вольтерры для 2 видов к наблюдаемым данным SSB ###
library(deSolve)
spp <- ObsSSB/10^6
# Первое предположение параметров
nospp <- ncol(spp)
data_lagged <- rbind(NA, spp) # Задержанные данные
db.bdt <- list()
for (i in 1:ncol(spp)) {
  db.bdt[[i]] <- log(append(spp[,i], NA)/data_lagged[,i])
}
modlst <- NULL # Линейные модели
for(i in 1:nospp) {
  moddat <- data.frame(db.bdt = db.bdt[[i]],
                       data_lagged)[-nrow(data_lagged),]
  modlst[[i]] <- eval(parse(text = paste("lm(db.bdt~",
                                         paste(colnames(data_lagged),
                                               collapse = "+"), ", data=moddat)",
                                         sep = "")))
}
r_start <- unname(sapply(modlst, function(dat){coef(dat)["(Intercept)"]}))
A_start <- unname(sapply(modlst, function(dat){coef(dat)[-1]}))
# Входные данные модели
params <- c(r1 = r_start[1], r2 = r_start[2],
            a11 = A_start[1,1], a12 = A_start[2,1],
            a21 = A_start[1,2], a22 = A_start[2,2])
ini <- c(b1 = spp[1,1], b2 = spp[1,2])
tmax <- nrow(spp)
t <- seq(1,tmax,1)

Передавая модели правильные значения ставок вылова на каждом временном шаге, строя/запуская оптимизатор и, наконец, решая модель:

# Данные Fdata
H1 <- approxfun(Fdata[,1])
H2 <- approxfun(Fdata[,2])
# Оптимизатор
LVmse = function(parms) {
  out = as.matrix(deSolve::ode(ini, 1:tmax, HS2LLV, parms, method="rk4")[,-1])
  RSS = sum(log(out / spp) ^ 2,na.rm = TRUE)
  return(RSS)
}
# Запуск модели и оптимизация параметров с оценками Fs
optimres1 <- optim(par = params, fn = LVmse, method = "Nelder-Mead",
                   control = list(maxit = 1000))
fit1 <- deSolve::ode(y = ini, time = 1:tmax, func = HS2LLV, parms = optimres1$par)
# График
extrafont::loadfonts(device = "all")
matplot(1974:2020, ObsSSB/10^6, type = "n", xlab = "Время (годы)",
        ylab = "SSB (1000 x тонн)", ylim = c(0,2.025))
grid()
matplot(1974:2020, ObsSSB/10^6, pch = 16, col = c("gray", "steelblue"),
        add = TRUE)
matplot(1974:2020, fit1[,2:3], type = "l", add = TRUE,
        col = c("gray", "steelblue"), lwd = 2)
legend("topright", legend = c("Сельдь","Шпрот"),
       lty = c(1,2), col = c("gray", "steelblue"), lwd = 2, bty = "n")

Как указано в обширном ответе, предоставленном в оригинальном посте; “хорошая подгонка к данным в данном случае никоим образом не гарантирует разумного ответа”, учитывая, что проекция модели на 5 временных шагов вперед привела к тому, что 1 из 100 подгонок была разумным набором параметров, обеспечивающим сосуществование всех видов в модели. Моя цель – воспроизвести оригинальный пост для этого конкретного случая, чтобы выяснить, существует ли разумный набор параметров, при котором оба вида сосуществуют в будущем, учитывая данные.