8.3 Методы комплексации модельных прогнозов

Существует два возможных подхода к получению некоторого обобщенного предсказания на основе формального или неформального объединения частных прогнозов. В разделе 4.4 нами подробно рассматривалось агрегирование результатов моделирования, когда каждая модель имеет фактически одну и ту же структуру, но многократно обучается на слегка модифицированных исходных данных или при небольшой вариации управляющих параметров (bagging). Другой путь - построить на одних и тех же исходных данных некоторый “коллектив” (ensemble) из принципиально разнотипных моделей и выполнить комбинирование прогнозов (forecast combinations).

По аналогии с методами коллективного решения, столь эффективно используемым в человеческом обществе, считается, что суммарная эффективность системы коллективного прогнозирования теоретически будет в среднем значительно выше любого из его членов. Этот оптимизм мало чем подтверждается на практике, но нельзя не признать, что комбинированный прогноз значительно более устойчив к не вполне объяснимым “провалам”, которые свойственны отдельным индивидуальным методам, и поэтому может быть эффективным инструментом страхования на случай возможных стохастических рисков (в мире трейдинга это называется хеджированием, англ. hedging).

Пусть имеется \(m\) прогнозов \(\boldsymbol{y_1, y_2, ..., y_m}\) для переменной \(\mathbf{Y}\), полученных с помощью \(m\) различных моделей. Под коллективным прогнозом \(\boldsymbol{g}\) будем понимать прогноз той же переменной \(\mathbf{Y}\) как некоторую функции от индивидуальных прогнозов \(\boldsymbol{y_k}\): \(\boldsymbol{g = F(y_1, y_2, \dots, y_m; X)}\).

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

Коллективный прогноз \(\boldsymbol{g}\), построенный без адаптации, представляют в виде линейной комбинации из базового или суженного (наиболее информативного) множества частных прогнозов

\[g = \sum_{k=1}^m y_k \omega_k,\]

где \(\mathbf{\omega}\) - вектор неизвестных весовых коэффициентов, компоненты которого неизменны на всем диапазоне варьирования исходных данных \(\mathbf{X}\). Тогда задача комплексации эквивалентна определению совокупности значений вектора \(\omega_k\), удовлетворяющих заданным ограничениям и минимизирующих некоторый критерий качества. Очевидно, что качество (надежность) прогноза \(\boldsymbol{g}\) тем выше, чем “ближе” расчетные значения \(\mathbf{Y}\) к фактическим.

Обратимся вновь к примеру из предыдущих разделов, но будем теперь строить модели прогнозирования не для возрастных категорий, а непосредственно для числа колец в раковинах rings тасманских морских ушек. Заменим фактор пол на две метрические переменные полM и полI, после чего разделим исходную выборку на обучающую и проверочную в соотношении 2:1.

load(file = "data/abalone.RData")
# Формирование матрицы предикторов
X <- model.matrix(rings ~ ., data = abalone[, 1:9])[, -1]
Y <- abalone$rings
# Разделение на обучающую и проверочную выборки
train <- runif(nrow(X)) <= .66

На объектах обучающей выборки (2748 наблюдений) построим шесть моделей оценки rings с использованием функции train() из пакета caret:

library(caret)
library(party)

# Для параллельных вычислений:
require(doParallel)
cl <- makeCluster(detectCores() - 1)
registerDoParallel(cl)

myControl <- trainControl(method = 'cv', number = 10, 
                          savePredictions = TRUE, returnData = FALSE,
                          verboseIter = FALSE)
PP <- c('center', 'scale')

# Обучение избранных моделей
set.seed(202)
# Регрессия на k-ближайших соседей 
m.knn <- train(X[train, ], Y[train], method = 'knn',
               trControl = myControl, preProcess = PP)
# Линейная модель
m.lm <- train(X[train, ], Y[train], method = 'lm',
              trControl = myControl, preProcess = PP)
# Гребневая регрессия с регуляризацией
m.rlm <- train(X[train, ], Y[train], method = 'glmnet',
               trControl = myControl, preProcess = PP)
# Модель опорных векторов
m.svm <- train(X[train, ], Y[train], method = 'svmRadial',
               trControl = myControl, preProcess = PP)
#  Метод случайного леса
m.rf <- train(X[train, ], Y[train], method = 'rf',
              trControl = myControl)
# Бэггинг деревьев условного вывода
m.ctr <- train(X[train, ], Y[train], "bag", B = 10, 
               bagControl = bagControl(fit = ctreeBag$fit,
                                       predict = ctreeBag$pred, 
                                       aggregate = ctreeBag$aggregate))

stopCluster(cl)

Обратим внимание на модель m.ctr, построенную методом "bag". При построении этой модели используется специальная функция bag() из пакета caret, позволяющая выполнить бэггинг широкого класса различных моделей (daBag, plsBag, nbBag, ctreeBag, svmBag, nnetBag). Мы осуществили бэггинг деревьев условного вывода (см. раздел 4.3), агрегируя каждый раз по 10 деревьев.

Упакуем результаты моделирования в объект класса list и извлечем значения среднеквадратичной ошибки RMSE каждой индивидуальной модели на обучающей выборке:

# Создание списка всех моделей
all.models <- list(m.knn, m.lm, m.rlm, m.svm, m.rf, m.ctr)
names(all.models) <- sapply(all.models, function(x) x$method)
sort(sapply(all.models, function(x) min(x$results$RMSE)))
##        rf svmRadial        lm    glmnet       bag       knn 
##  2.142822  2.151066  2.189748  2.190679  2.213785  2.251248

В целом все модели показали близкие результаты, хотя можно отметить традиционное лидерство методов случайного леса (rf) и опорных векторов (svmRadial).

Сформируем таблицу частных прогнозов для объектов тестовой выборки (1429 наблюдений) и оценим ошибки индивидуальных моделей RMSE:

preds.all <- data.frame(sapply(all.models, function(x){predict(x, X)}))
head(preds.all) 
##         knn        lm    glmnet svmRadial        rf       bag
## 1  9.111111  9.271303  9.253349  8.479311 12.071000  9.063321
## 2  8.888889  7.829085  7.842020  8.365610  7.177400  8.979464
## 3 10.000000 11.169981 11.109249 10.796564  9.761967 12.347767
## 4  9.222222  9.613559  9.611096  9.568827  9.788767  9.714990
## 5  6.333333  6.690712  6.684828  6.237234  6.521167  6.646564
## 6  7.000000  7.808261  7.816015  7.287081  7.803167  7.670769
rmse <- function(x,y){sqrt(mean((x - y)^2))}
print("Среднеквадратичная ошибка на тестовой выборке:")
## [1] "Среднеквадратичная ошибка на тестовой выборке:"
print(sort(apply(preds.all[!train,], 2, rmse, y = Y[!train])), digits = 3)
##        rf svmRadial       bag        lm    glmnet       knn 
##      2.19      2.19      2.24      2.25      2.25      2.25

В целом незначительное (2-4%) увеличение ошибки RMSE на свежих наблюдениях свидетельствует о хорошем качестве моделей.

Рассмотрим теперь возможные алгоритмы комплексации частных прогнозов, реализованные в пакете ForecastCombinations. Два из них основаны только на самих предсказаниях и не требуют предварительного обучения:

  • простое среднее (simple) при \(k = const = 1/m\);
  • по лучшему частному прогнозу (best), т.е. \(k = 1\) для \(max(y_k)\) и \(k = 0\) для остальных \(y_k\).

Комплексация, основанная на дисперсиях (variance based), требует оценки на некоторой обучающей выборке значений среднеквадратичных отклонений \(MSE_k\) для каждой индивидуальной модели и тогда

\[\omega_k = (1/MSE_k)/\sum_{k=1}^m (1/MSE_k).\] Другим способом комплексации прогнозов разных моделей является использование линейной модели, в которой коэффициенты играют роль весов, определяющих вклад каждой из объединяемых моделей. С легкой руки букмекеров эту процедуру стали часто называть стэкингом (stacking).

Разумеется, обычный метод наименьших квадратов (ols) является одним из претендентов на создание такой “модели моделей”. Однако, поскольку прогнозы частных моделей сильно коррелируют между собой, то для получения комплексного результата предлагают использовать иные методы регрессии, которые, по мнению авторов, могут лучше справиться с этой напастью. В пакет ForecastCombinations включены метод наименьших абсолютных отклонений (robust) и наименьших квадратов с ограничениями (cls - Constrained Least Squares). Во втором случае на значения коэффициентов регрессии накладываются ограничения: все они неотрицательны и в сумме равны 1, т.е. по сути отражают вклад каждой частной модели в общий прогноз. Коэффициенты CLS-модели рассчитываются с использованием метода квадратичного программирования, представленного функцией solve.QP() из пакета quadprog.

На вход функции Forecast_comb() из пакета ForecastCombinations подаются объекты для обучения (матрица fhat со значениями частных прогнозов \(y_1, y_2,\dots, y_m\) и вектор obs с соответствующими наблюдаемыми величинами \(y\)) и прогноза (матрица fhat_new со значениями частных прогнозов, которые следует объединить с использованием метода Averaging_scheme). На выходе формируются векторы прогнозов pred и весов weights. С учетом выполненного нами разбиения исходной таблицы abalone, рассчитаем сравнительную оценку RMSE, полученную каждым методом комплексации на проверочной выборке:

library(ForecastCombinations)
scheme <- c("simple",  "variance based", "ols", "robust", "cls", "best")
combine_f <- list()
for (i in scheme) {
    combine_f[[i]] <- Forecast_comb(obs = Y[train],
                                    fhat = as.matrix(preds.all[train, ]),
                                    fhat_new = as.matrix(preds.all[!train, ]),
                                    Averaging_scheme = i)
    tmp <- round(sqrt(mean((combine_f[[i]]$pred - Y[!train])^2)), 3)
    cat(i, ":", tmp, "\n")
}
## simple : 2.155 
## variance based : 2.151 
## ols : 2.303 
## robust : 2.304 
## cls : 2.19 
## best : 2.19

Рассмотрим веса, которые получили каждые модели при реализации каждой схемы усреднения (приведем их предварительно к единой шкале):

w.list <- sapply(combine_f, function(sp.list) sp.list$weights)
w.list$ols <- w.list$ols[-1]
w.list$robust <- w.list$robust[-1]
weights <- as.data.frame.list(w.list)
rownames(weights) <- scheme  
wstan <- function(x) abs(x)/sum(abs(x)) 
weights[3, ] <- wstan(weights[3, ])
weights[4, ] <- wstan(weights[4, ])
print(round(weights, 3))
##                simple.1 simple.2 simple.3 simple.4 simple.5 simple.6
## simple            0.167    0.167    0.167    0.167    0.167    0.167
## variance based    0.167    0.167    0.167    0.167    0.167    0.167
## ols               0.123    0.123    0.123    0.123    0.123    0.123
## robust            0.104    0.104    0.104    0.104    0.104    0.104
## cls               0.167    0.167    0.167    0.167    0.167    0.167
## best              0.167    0.167    0.167    0.167    0.167    0.167
##                variance.based    ols robust cls best
## simple                  0.108 -0.142 -0.148   0    0
## variance based          0.094 -0.105 -0.193   0    0
## ols                     0.069  0.059  0.136   0    0
## robust                  0.065  0.150  0.164   0    0
## cls                     0.493  1.754  1.738   1    1
## best                    0.107 -0.390 -0.350   0    0

Комплексные прогнозы, созданные по схемам cls и best, эквивалентны мнению “непогрешимого” авторитета (а таковым оказался метод rf). Соответственно они дали ошибку, в точности равную погрешности наиболее качественной индивидуальной модели. Схемы ols и robust потерпели относительную неудачу, тогда как наиболее точным оказался прогноз, основанный на дисперсии (variance based).

П. Полищук (http://www.qsar4u.com) предлагает применять метод неотрицательных наименьших квадратов (non-negative least squares), используя пакет nnls:

require(nnls)
m.nnls <- nnls(as.matrix(preds.all[train, ]), Y[train])
coef(m.nnls)
## [1] 0.000000 0.000000 0.000000 0.000000 1.008573 0.000000

Этот метод дал те же значения весов, что и CLS-регрессия (впрочем, эти методы идентичны, по-видимому, и по своей сути).

Кроме линейного стэкинга для комплексации можно применять любые изощренные алгоритмы статистического моделирования. Например, набор частных прогнозов можно выразить через небольшое число главных компонент и на последующих шагах применить любой другой алгоритм взвешивания (Горелик, Френкель, 1983).9 П. Брусиловский и Г. Розенберг (1983) использовали для комплексации многорядный алгоритм МГУА, строили самообучающиеся иерархические модели и назвали этот метод “модельным штурмом”. Дело устойчиво идет к проблемам создания “коллектива коллективов”.

Пакет ForecastCombinations также позволяет реализовать эксперименты подобного рода. Из 6 векторов частных прогнозов, полученных нами по индивидуальным моделям, можно построить 53 модели регрессии с различными наборами переменных от одного до 5. Функция Forecast_comb_all() выполняет настройку всех этих моделей и сводит результаты в новую таблицу (т.е. вместо 6 прогнозов предлагается иметь дело с 53-мя прогнозами на основе прогнозов). Далее их можно просто усреднить:

combine_f_all <- Forecast_comb_all(Y[train],
                                   fhat = as.matrix(preds.all[train, ]),
                                   fhat_new = as.matrix(preds.all[!train, ]))
combine_f_all <- Forecast_comb_all(Y[train],
                                   fhat = as.matrix(preds.all[train, ]),
                                   fhat_new = as.matrix(preds.all[!train, ]))
# Усредняем комбинированные прогнозы по всем регрессиям
Combined_f_all_simple <- apply(combine_f_all$pred, 1, mean)
print(sqrt(mean((Combined_f_all_simple - Y[!train])^2)), digits = 3)
## [1] 2.17

Для моделей регрессии по каждому подмножеству отдельных моделей функция Forecast_comb_all() рассчитывает также набор информационных критериев (AIC, BIC и Mallow’s Cp), которые могут быть использованы для расчета взвешенного среднего:

# Комбинирование прогнозов с использованием  Mallow's Cp:
Combined_f_all_mal <- t(combine_f_all$mal %*% t(combine_f_all$pred))
print(sqrt(mean((Combined_f_all_mal - Y[!train])^2)), digits = 3)
## [1] 2.27

Итогом всех этих достаточно трудоемких вычислительных процедур оказался весьма скромный результат.

Альтернативой механическому усреднению прогнозов является учет структурных связей в ансамбле таким образом, чтобы положительные свойства той или иной модели (метода) доминировали при формировании коллективного отклика, а отрицательные - отфильтровывались (Растригин, Эренштейн, 1981). В частности, существенным влиянием должны обладать самые “непохожие” между собой модели, являющиеся носителями нетривиальных тенденций моделируемого процесса. Основным направлением развития в этой области является разработка методов адаптации, в основе которых лежат алгоритмы Бейтса-Гренджера (Bates, Grander, 1965), Ньюболда-Гренджера (Newbald, Grander, 1974), Ершова (1975) и др. К сожалению, нам не удалось обнаружить пакеты R, в которых были бы реализованы вменяемые методы комплексации прогнозов с адаптацией.

  1. Информация о всех литературных источниках этого раздела приведена в (Розенберг и др., 1994).↩