Функции плотности вероятности в R для прогнозирования следующего значения инцидентов

В широком смысле двумя основными подходами к моделированию являются так называемые «механистический» и «эмпирический» подходы. Оба имеют своих приверженцев (и своих недоброжелателей). Механистический подход утверждает, что моделирование должно исходить из понимания лежащих в основе явлений (механизма), которое затем преобразуется в некоторый тип математических уравнений, которые затем подгоняются к данным (для проверки механизма). Эмпирический подход собирает (обычно длинный) список моделей (уравнений) и пытается найти ту, которая «подходит лучше всего». Эмпирическое моделирование привлекательно, но опасно, потому что оценка того, когда вы «хорошо подходите», не является тривиальной задачей, хотя к ней часто относятся именно так.

Вы не дали нам достаточно информации, чтобы сформулировать механистическую модель, поэтому вот иллюстрация пары эмпирических моделей в качестве предостережения:

Модели сингулярности с конечным временем популярны для вашего типа данных. Среди прочего, эти модели используются для «предсказания» пузырей на фондовом рынке. (модель LPPL). Основная идея заключается в том, что грядет катастрофа (сингулярность), и мы хотим предсказать, когда. Поэтому мы используем функцию вида:

у = а × (с-х)^b

При b ‹ 0 y приближается к сингулярности при x -> c.

В коде R мы можем подогнать такую ​​модель следующим образом:

# Finite-Time Singularity Model
library(minpack.lm)
f <- function(par,x) {
  a <- par[1]
  b <- par[2]
  c <- par[3]
  a * (c - x)^b
}
resid   <- function(par,obs,xx) {obs-f(par,xx)}
start <- c(a=1, b=-1, c=2100)
nls.out <- nls.lm(par=start, fn=resid, obs =dat$incidents, xx=dat$year, 
                  control = nls.lm.control(maxiter=500))
coef(nls.out)
with(dat, plot(incidents~year, main="Finite-Time Singularity Model"))
lines(dat$year,f(coef(nls.out),year), col=2, lwd=2)

Это дает то, что кажется «довольно хорошим соответствием»:

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

with(dat,plot(year,resid(coef(nls.out),incidents,year),
              main="Residuals Plot", ylab="residuals"))

Другой подход отмечает, что ваши данные «подсчитываются» (например, количество инцидентов в год). Это предлагает обобщенную линейную модель в семействе Пуассона:

# generalized liner model, poisson family
fit.glm <- glm(incidents ~year,data=dat,family=poisson)
with(dat,plot(incidents~year))
lines(dat$year,predict(fit.glm,type="response"), col=2, lwd=2)
par(mfrow=c(2,2))
plot(fit.glm)

Это соответствие лучше, но все еще не очень хорошо, как показывают диагностические графики. Остатки следуют тренду, они не распределены нормально, а некоторые точки данных имеют неприемлемо высокое кредитное плечо.

dat <- data.frame(year,incidents)
with(dat, plot(incidents~year))

Итак, что-то изменилось... но что вызывает резкое увеличение числа инцидентов? Только у вас, ученого, есть ключ. Вероятно, вы можете предсказать, что в ближайшие год или два будет некоторое увеличение, но будет ли это увеличение следовать экспоненциальному или логистическому шаблону, определяется основной областью исследования. Логистическая модель не будет очень точной, если вы находитесь в так называемой «логарифмической фазе» роста, потому что верхний предел инцидентов в год неизвестен.