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

Скрипт Open Live Script

В этом примере показано, как реализовать стохастический поисковый выбор переменной (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');

Figure contains an axes object. The axes object with title Variable Box Plots contains 98 objects of type line.

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

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

Разделите данные в оценку и предскажите выборки.
Подбирайте модели к выборке оценки и задайте для всех $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

Figure contains an axes object. The axes object with title Log GDP, V1 = 10, V2 = 0.5 contains an object of type line.

Figure contains an axes object. The axes object with title Log GDP, V1 = 50, V2 = 0.5 contains an object of type line.

Figure contains an axes object. The axes object with title Log GDP, V1 = 100, V2 = 0.5 contains an object of type line.

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

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