Введение

Рассмотрим модель байесовской линейной регрессии.

Цель выбора переменных состоит в том, чтобы включить в окончательную регрессионную модель только те предикторы, которые поддерживаются данными. Одним из способов сделать это является анализ $${\mathrm{}}^{\mathrm{2p}}$$ перестановок моделей, называемых режимами, где модели различаются по коэффициентам, которые включены. Если $$p$$ мал, то можно подогнать все перестановки моделей к данным, а затем сравнить модели с помощью показателей производительности, таких как goodness-of-fit (например, информационный критерий Акайке) или прогнозная среднеквадратичная ошибка (MSE). Однако для даже умеренных значений $$p$$оценка всех перестановок моделей неэффективна.

Байесовский вид выбора переменной - это коэффициент, исключаемый из модели, имеющий вырожденное заднее распределение. То есть исключенный коэффициент имеет дельта-распределение Дирака, которое имеет свою вероятностную массу, сконцентрированную на нуле. Чтобы обойти сложность, вызванную вырожденными вариациями, предшествующим для исключаемого коэффициента является гауссово распределение со средним значением 0 и малой дисперсией, например $$N\mathrm{(}\mathrm{}{\mathrm{}}^{\mathrm{}}0,0,12$$). Поскольку предшествующий уровень концентрируется около нуля, задний уровень также должен концентрироваться около нуля.

Предшествующим для включенного коэффициента может быть $$N(\lambda ,$$V), $$\mathrm{где}$$ V в достаточной степени отстоит от нуля $и$ λ обычно равно нулю. Эта структура подразумевает, что предшествующим для каждого коэффициента является модель гауссовой смеси.

Рассмотрим скрытые, двоичные случайные величины $${}_{\mathrm{\gamma \; k}}$$, $$k\mathrm{=}0,...$$, p, такие, что:

Одной из целей SSVS является оценка вероятности заднего режима $${}_{\mathrm{gk}}$$, оценки, которые определяют, должны ли соответствующие коэффициенты быть включены в модель. Учитывая $${}_{\mathrm{\beta k}}$$, $${}_{\mathrm{\gamma \; k}}$$ условно не зависит от данных. Поэтому для $$k\mathrm{=}0,...$$, p это уравнение представляет полное условное апостериорное распределение вероятности того, что переменная k включена в модель:

где $$\mathrm{start(}\lambda {,}^{\mathrm{}}$$start2) - pdf гауссова распределения со скалярным $$\mathrm{средним}$$ λ и $${\mathrm{дисперсией}}^{\mathrm{}}$$start2.

Эконометрика Toolbox™ имеет две байесовские модели линейной регрессии, которые определяют предыдущие распределения для SSVS: mixconjugateblm и mixsemiconjugateblm. Структура, представленная ранее, описывает предшествующие mixconjugateblm модель. Разница между моделями в том, что $$\beta$$ и $${}^{\mathrm{start2}}$$ независимы, априори, для mixsemiconjugateblm модели. Следовательно, предшествующая дисперсия $${}_{\mathrm{\beta k}}$$ равна $${}_{\mathrm{Vk1}}$$ (${}_{\mathit{\gamma \; k}}\mathrm{=}$1) $${\mathrm{или}}_{\mathrm{}}$$Vk2 ${}_{\mathit{(}}\mathrm{\gamma \; k}$= 0).

После того, как вы решите, какую предыдущую модель использовать, позвоните bayeslm для создания модели и задания значений гиперпараметров. Поддерживаемые гиперпараметры:

После создания модели передайте ее и данные в estimate. estimate функция использует пробоотборник Гиббса для выборки из полных условий и оценки характеристик задних распределений $$\beta$$ и $${}^{\mathrm{}}$$Также, estimate возвращает апостериорные оценки $${}_{\mathrm{gk}}$$.

Для этого примера рассмотрим создание прогностической линейной модели для уровня безработицы в США. Вы хотите модель, которая хорошо обобщается. Другими словами, необходимо минимизировать сложность модели, удалив все избыточные предикторы и все предикторы, которые не связаны с уровнем безработицы.

Загрузка и предварительная обработка данных

Загрузить набор макроэкономических данных США Data_USEconModel.mat.

load Data_USEconModel

Набор данных включает расписание MATLAB ®DataTable, который содержит 14 переменных, измеренных с Q1 1947 по Q1 2009; UNRATE - уровень безработицы в США. Для получения дополнительной информации введите Description в командной строке.

Постройте все серии на одном рисунке, но на отдельных вложенных графиках.

figure;
for j = 1:size(DataTable,2)
    subplot(4,4,j)
    plot(DataTable.Time,DataTable{:,j});
    title(DataTable.Properties.VariableNames(j));
end

Все серии, кроме FEDFUNDS, GS10, TB3MS, и UNRATE похоже, имеет экспоненциальный тренд.

Примените преобразование журнала к этим переменным с экспоненциальным трендом.

hasexpotrend = ~ismember(DataTable.Properties.VariableNames,...
    ["FEDFUNDS" "GD10" "TB3MS" "UNRATE"]);
DataTableLog = varfun(@log,DataTable,'InputVariables',...
    DataTable.Properties.VariableNames(hasexpotrend));
DataTableLog = [DataTableLog DataTable(:,DataTable.Properties.VariableNames(~hasexpotrend))];

DataTableLog является расписанием, как DataTable, но эти переменные с экспоненциальным трендом находятся в логарифмической шкале.

Коэффициенты, которые имеют относительно большие величины, имеют тенденцию доминировать над штрафом в целевой функции регрессии лассо. Поэтому важно, чтобы переменные имели одинаковый масштаб при реализации регрессии лассо. Сравнение масштабов переменных в DataTableLog путем построения графиков на одной и той же оси.

figure;
boxplot(DataTableLog.Variables,'Labels',DataTableLog.Properties.VariableNames);
h = gcf;
h.Position(3) = h.Position(3)*2.5;
title('Variable Box Plots');

Переменные имеют довольно схожие масштабы.

Чтобы настроить предыдущие коэффициенты дисперсии гауссовой смеси, выполните следующую процедуру:

Джордж и Маккаллох предлагают другой способ настроить предыдущие дисперсии $$\beta$$ в [1].

Создайте выборочные переменные оценки и прогноза для данных ответа и предиктора. Укажите горизонт прогноза 4 года (16 кварталов).

fh = 16;
y = DataTableLog.UNRATE(1:(end - fh));
yF = DataTableLog.UNRATE((end - fh + 1):end);
isresponse = DataTable.Properties.VariableNames == "UNRATE";
X = DataTableLog{1:(end - fh),~isresponse};
XF = DataTableLog{(end - fh + 1):end,~isresponse};
p = size(X,2); % Number of predictors
predictornames = DataTableLog.Properties.VariableNames(~isresponse);

Создание предыдущих байесовских моделей линейной регрессии

Создание предыдущих байесовских моделей линейной регрессии для SSVS путем вызова bayeslm и определение количества предикторов, типа модели, имен предикторов и коэффициентов дисперсии компонентов. Предположим, что $$\beta$$ и $${}^{\mathrm{start2}}$$ зависимы, априори (mixconjugateblm модель).

V1 = [10 50 100];
V2 = [0.05 0.1 0.5];
numv1 = numel(V1);
numv2 = numel(V2);

PriorMdl = cell(numv1,numv2); % Preallocate

for k = 1:numv2
    for j = 1:numv1
        V = [V1(j)*ones(p + 1,1) V2(k)*ones(p + 1,1)];
        PriorMdl{j,k} = bayeslm(p,'ModelType','mixconjugateblm',...
            'VarNames',predictornames,'V',V);
    end
end

PriorMdl является массивом ячеек 3 на 3, и каждая ячейка содержит mixconjugateblm объект модели.

Постройте график предварительного распределения log_GDP для моделей, в которых V2 составляет 0,5.

for j = 1:numv1
    [~,~,~,h] = plot(PriorMdl{j,3},'VarNames',"log_GDP");
    title(sprintf("Log GDP, V1 = %g, V2 = %g",V1(j),V2(3)));
    h.Tag = strcat("fig",num2str(V1(j)),num2str(V2(3)));
end

Предыдущие распределения $$\beta$$ имеют форму шипа и сляба. Когда V1 является низким, большая часть распределения сосредоточена вокруг 0, что затрудняет для алгоритма приписывание высокого значения для бета-версии. Однако переменные, которые алгоритм идентифицирует как важные, регуляризуются в том смысле, что алгоритм не относит высокую величину к соответствующим коэффициентам.

Когда V1 высокая, большая плотность происходит далеко от нуля, что облегчает для алгоритма приписывание ненулевых коэффициентов важным предикторам. Однако, если V1 слишком высока, тогда важные предикторы могут иметь завышенные коэффициенты.

Выполнение выбора переменных SSVS

Для выполнения SSVS оцените апостериорные распределения с помощью estimate. Используйте параметры по умолчанию для образца Гиббса.

PosteriorMdl = cell(numv1,numv2);
PosteriorSummary = cell(numv1,numv2);

rng(1); % For reproducibility
for k = 1:numv2
    for j = 1:numv1
        [PosteriorMdl{j,k},PosteriorSummary{j,k}] = estimate(PriorMdl{j,k},X,y,...
            'Display',false);
    end
end

Каждая ячейка в PosteriorMdl содержит empiricalblm объект модели, хранящий полные условные задние отрисовки из семплера Гиббса. Каждая ячейка в PosteriorSummary содержит таблицу задних оценок. Regime табличная переменная представляет заднюю вероятность включения переменной ($${}_{\mathrm{gk}}$$).

Отобразить таблицу задних оценок $${}_{\mathrm{gk}}$$.

RegimeTbl = table(zeros(p + 2,1),'RowNames',PosteriorSummary{1}.Properties.RowNames);
for k = 1:numv2
    for j = 1:numv1
        vname = strcat("V1_",num2str(V1(j)),"__","V2_",num2str(V2(k)));
        vname = replace(vname,".","p");
        tmp = table(PosteriorSummary{j,k}.Regime,'VariableNames',vname);
        RegimeTbl = [RegimeTbl tmp];
    end
end
RegimeTbl.Var1 = [];
RegimeTbl

RegimeTbl=15×9 table
                    V1_10__V2_0p05    V1_50__V2_0p05    V1_100__V2_0p05    V1_10__V2_0p1    V1_50__V2_0p1    V1_100__V2_0p1    V1_10__V2_0p5    V1_50__V2_0p5    V1_100__V2_0p5
                    ______________    ______________    _______________    _____________    _____________    ______________    _____________    _____________    ______________

    Intercept           0.9692                 1                 1            0.9501                1                 1           0.9487           0.9999                 1    
    log_COE             0.4686            0.4586            0.5102            0.4487           0.3919            0.4785           0.4575           0.4147            0.4284    
    log_CPIAUCSL        0.9713            0.3713            0.4088             0.971           0.3698            0.3856            0.962           0.3714            0.3456    
    log_GCE             0.9999                 1                 1            0.9978                1                 1           0.9959                1                 1    
    log_GDP             0.7895            0.9921            0.9982            0.7859           0.9959                 1           0.7908           0.9975            0.9999    
    log_GDPDEF          0.9977                 1                 1                 1                1                 1           0.9996                1                 1    
    log_GPDI                 1                 1                 1                 1                1                 1                1                1                 1    
    log_GS10                 1                 1            0.9991                 1                1            0.9992                1           0.9992             0.994    
    log_HOANBS          0.9996                 1                 1            0.9887                1                 1           0.9763                1                 1    
    log_M1SL                 1                 1                 1                 1                1                 1                1                1                 1    
    log_M2SL            0.9989            0.9993            0.9913            0.9996           0.9998            0.9754           0.9951           0.9983            0.9856    
    log_PCEC            0.4457            0.6366            0.8421            0.4435           0.6226            0.8342           0.4614            0.624              0.85    
    FEDFUNDS            0.0762            0.0386            0.0237            0.0951           0.0465            0.0343           0.1856           0.0953             0.068    
    TB3MS               0.2473            0.1788            0.1467            0.2014           0.1338            0.1095           0.2234           0.1185            0.0909    
    Sigma2                 NaN               NaN               NaN               NaN              NaN               NaN              NaN              NaN               NaN    

Используя произвольное пороговое значение 0,10, все модели соглашаются с тем, что FEDFUNDS является незначительным или избыточным предиктором. Когда V1 высокий, TB3MS граничит с незначительностью.

Прогнозирование ответов и вычисление прогнозных MSE с использованием расчетных моделей.

yhat = zeros(fh,numv1*numv2);
fmse = zeros(numv1,numv2);

for k = 1:numv2
    for j = 1:numv1
        idx = ((k - 1)*numv1 + j); 
        yhat(:,idx) = forecast(PosteriorMdl{j,k},XF);
        fmse(j,k) = sqrt(mean((yF - yhat(:,idx)).^2));
    end
end

Определите параметры коэффициента отклонения, которые дают минимальный прогноз MSE.

minfmse = min(fmse,[],'all');
[idxminr,idxminc] = find(abs(minfmse - fmse) < eps);
bestv1 = V1(idxminr)

bestv1 = 100

bestv2 = V2(idxminc)

bestv2 = 0.0500

Оцените модель SSVS, используя весь набор данных и настройки коэффициента дисперсии, которые дают минимальный прогноз MSE.

XFull = [X; XF];
yFull = [y; yF];
EstMdl = estimate(PriorMdl{idxminr,idxminc},XFull,yFull);

Method: MCMC sampling with 10000 draws
Number of observations: 201
Number of predictors:   14
 
              |   Mean      Std          CI95         Positive  Distribution  Regime 
-------------------------------------------------------------------------------------
 Intercept    |  29.4598  4.2723   [21.105, 37.839]     1.000     Empirical    1     
 log_COE      |   3.5380  3.0180   [-0.216,  9.426]     0.862     Empirical   0.7418 
 log_CPIAUCSL |  -0.6333  1.7689   [-5.468,  2.144]     0.405     Empirical   0.3711 
 log_GCE      |  -9.3924  1.4699  [-12.191, -6.494]     0.000     Empirical    1     
 log_GDP      |  16.5111  3.7131   [ 9.326, 23.707]     1.000     Empirical    1     
 log_GDPDEF   |  13.0146  2.3992   [ 9.171, 19.131]     1.000     Empirical    1     
 log_GPDI     |  -5.9537  0.6083   [-7.140, -4.756]     0.000     Empirical    1     
 log_GS10     |   1.4485  0.3852   [ 0.680,  2.169]     0.999     Empirical   0.9868 
 log_HOANBS   | -16.0240  1.5361  [-19.026, -13.048]    0.000     Empirical    1     
 log_M1SL     |  -4.6509  0.6815   [-5.996, -3.313]     0.000     Empirical    1     
 log_M2SL     |   5.3320  1.3003   [ 2.738,  7.770]     0.999     Empirical   0.9971 
 log_PCEC     |  -9.9025  3.3904  [-16.315, -2.648]     0.006     Empirical   0.9858 
 FEDFUNDS     |  -0.0176  0.0567   [-0.125,  0.098]     0.378     Empirical   0.0269 
 TB3MS        |  -0.1436  0.0762   [-0.299,  0.002]     0.026     Empirical   0.0745 
 Sigma2       |   0.2891  0.0289   [ 0.238,  0.352]     1.000     Empirical    NaN   
 

EstMdl является empiricalblm модель, представляющая результат выполнения SSVS. Вы можете использовать EstMdl прогнозировать уровень безработицы, например, с учетом будущих данных прогноза.