Введение

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

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

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

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

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

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

где $$\varphi (\mu ,{\sigma }^2)$$ - PDF Гауссова распределения со скалярным средним $$\mu$$ и отклонение $${\sigma }^2$$.

Econometrics Toolbox™ имеет две байесовские линейные регрессионые модели, которые задают предыдущие распределения для SSVS: mixconjugateblm и mixsemiconjugateblm. Представленные ранее среды описывают приоритеты mixconjugateblm модель. Различие между моделями заключается в том, что $$\beta$$ и $${\sigma }^2$$ являются независимыми, априори, для mixsemiconjugateblm модели. Поэтому предыдущее отклонение $${\beta }_k$$ является $$V_{k1}$$ (${\gamma }_{\mathit{k}}=1$) или $$V_{k2}$$ (${\gamma }_{\mathit{k}}=0$).

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

После того, как вы создаете модель, передайте ее и данные в estimate. The estimate функция использует семплер Гиббса для выборки из полных обусловленности и оценки характеристик апостериорных распределений $$\beta$$ и $${\sigma }^2$$. Кроме того, estimate возвращает апостериорные оценки $$g_k$$.

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

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

Загрузите набор макроэкономических данных США 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$$ и $${\sigma }^2$$ являются зависимыми, априори (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 содержит таблицу апостериорных оценок. The Regime табличная переменная представляет апостериорную вероятность включения переменной ($$g_k$$).

Отобразите таблицу апостериорных оценок $$g_k$$.

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 прогнозировать уровень безработицы, например, с учетом будущих прогнозирующих данных.