Байесов стохастический поисковый выбор переменной

Этот пример показывает, как реализовать стохастический поисковый выбор переменной (SSVS), Байесов метод выбора переменной для моделей линейной регрессии.

Введение

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

$y_{t} = \sum_{k} β_{k} x_{t k} + ε_{t} .$

Коэффициенты регрессии $β_{k} | σ^{2} \sim N (μ_{j}, σ^{2} V_{k})$ .
$k = 0, . . ., p$ .
Воздействия $ε_{t} \sim N (0, σ^{2})$ .
Отклонение воздействия $σ^{2} \sim I G (A, B)$ , где $I G (A, B)$ обратное гамма распределение с формой A и шкала B.

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

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

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

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

$γ_{k} = 1$ указывает на это $β_{k} \sim N (0, σ^{2} V_{1 k})$ и это $β_{k}$ включен в модель.
$γ_{k} = 0$ указывает на это $β_{k} \sim N (0, σ^{2} V_{2 k})$ и это $β_{k}$ исключен из модели.
$γ_{k} \sim B e r n o u l l i (g_{k})$ .
Выборочное пространство $γ_{k}$ имеет кардинальность $2^{p + 1}$ , и каждый элемент является a $p + 1$ - D вектор нулей или единиц.
$V_{2 k}$ маленькое, положительное число и $V_{1 k} > V_{2 k}$ .
Коэффициенты $β_{j}$ и $β_{k}$ , $j \neq k$ независимы, априорно.

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

$P (γ_{k} = 1 | β, σ^{2}, γ_{\neq k}) \propto g_{k} ϕ (β_{k}; 0, σ^{2} V_{1 k}),$

где $ϕ (μ, σ^{2})$ PDF Распределения Гаусса со скалярным средним значением $μ$ и отклонение $σ^{2}$ .

Econometrics Toolbox™ имеет две Байесовых модели линейной регрессии, которые задают предшествующие дистрибутивы для SSVS: mixconjugateblm и mixsemiconjugateblm. Среда представила, ранее описывает уголовное прошлое модели mixconjugateblm. Различие между моделями - это $β$ и $σ^{2}$ независимы, априорно, для моделей mixsemiconjugateblm. Поэтому предшествующее отклонение $β_{k}$ $V_{k 1}$ ( $γ_{k} = 1$ ) или $V_{k 2}$ ( $γ_{k} = 0$ ).

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

Intercept, логический скаляр, задающий, включать ли прерывание в модель.
Mu, (p + 1)-by-2 матрица, задающая предшествующие Гауссовы средние значения смеси $β$ . Первый столбец содержит средние значения для соответствия компонента $γ_{k} = 1$ , и второй столбец содержит среднее соответствие $γ_{k} = 0$ . По умолчанию все средние значения 0, который задает реализацию SSVS.
V, (p + 1)-by-2 матрица, задающая предшествующие Гауссовы факторы отклонения смеси (или отклонения) $β$ . Столбцы соответствуют столбцам Mu. По умолчанию отклонением первого компонента является 10, и отклонением второго компонента является 0.1.
Correlation, (p + 1) (p + 1) положительная определенная матрица, задающая предшествующую корреляционную матрицу $β$ для обоих компонентов. Значением по умолчанию является единичная матрица, которая подразумевает, что коэффициенты регрессии являются некоррелироваными, априорно.
Probability, (p + 1)-D вектор априорных вероятностей переменного включения ( $g_{k}$ , k = 0..., _p _) или указатель на функцию к пользовательской функции. $γ_{j}$ и $γ_{k}$ , $j \neq k$ , независимы, априорно. Однако с помощью указателя на функцию (@functionname), можно предоставить пользовательское предшествующее распределение, которое задает зависимости между $γ_{j}$ и $γ_{k}$ . Например, можно задать принуждение $x_{2}$ из модели, если $x_{4}$ включен.

После того, как вы создаете модель, передаете ее и данные к estimate. Функция estimate использует сэмплер Гиббса для выборки от полных условных выражений и оценочные характеристики апостериорных распределений $β$ и $σ^{2}$ . Кроме того, estimate возвращает следующие оценки $g_{k}$ .

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

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

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

load Data_USEconModel

Набор данных включает расписание MATLAB® DataTable, который содержит 14 переменных, измеренных от 1 квартала 1947 до 1 квартала 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');

Переменные имеют довольно подобные шкалы.

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

Разделите данные в оценку и предскажите выборки.
Соответствуйте моделям к выборке оценки и задайте для всех $k$ , $V_{1 k} = {10, 50, 100}$ и $V_{2 k} = {0 . 05, 0.1, 0.5}$ .
Используйте подобранные модели, чтобы предсказать ответы в горизонт прогноза.
Оцените прогноз MSE для каждой модели.
Выберите модель с самым низким прогнозом MSE.

Джордж и Маккалок предлагают другой способ настроить предшествующие отклонения $β$ в [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 и определения количества предикторов, типа модели, имен предиктора и факторов отклонения компонента. AssumeThat $β$ и $σ^{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

Предшествующие дистрибутивы $β$ имейте форму скачка-и-плиты. Когда 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 представляет апостериорную вероятность переменного включения ( $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 граничит с тем, чтобы быть незначительным.

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

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

Ссылки

[1] Джордж, E. I. и Р. Э. Маккалок. "Выбор переменной Через Гиббс, Выбирающий". Журнал американской Статистической Ассоциации. Издание 88, № 423, 1993, стр 881–889.

Смотрите также

estimate | sampleroptions

Документация