Байесова регрессия лассо

Открыть скрипт

В этом примере показано, как выполнить выбор переменной при помощи Байесовой регрессии лассо.

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

Рассмотрите эту модель многофакторной линейной регрессии:

$$y = \beta_0 1_n + X\beta + \varepsilon.$$

  • $y$ вектор из n ответов.

  • $X$ матрица p соответствующих наблюдаемых переменных предикторов.

  • $\beta$ вектор из p коэффициентов регрессии.

  • $\beta_0$ прерывание.

  • $1_n$ длина n вектор-столбец из единиц.

  • $\varepsilon$ случайный вектор из iid Гауссовых воздействий со средним значением 0 и отклонением$\sigma^2$.

Целевая функция частотного представления регрессии лассо

$$f\left( {\beta_0,\beta;y,X} \right) = 0.5\left\| {y - \beta_0 1_n + X\beta } \right\|_2^2 + \lambda \left\| \beta \right\|_1^2.$$

$\lambda$ штраф (или уменьшение) термин и, в отличие от других параметров, не является подходящим к данным. Необходимо задать значение для$\lambda$ перед оценкой, и лучшая практика состоит в том, чтобы настроить ее. После того, как вы зададите значение,$f$ минимизирован относительно коэффициентов регрессии. Получившиеся коэффициенты являются оценками лассо. Для получения дополнительной информации о частотном представлении регрессии лассо см. [186].

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

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

Загрузите США макроэкономический набор данных 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 = gca;
h.XTickLabelRotation = 45;
title('Variable Box Plots');

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

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

  1. Разделите данные в оценку и предскажите выборки.

  2. Подбирайте модели к выборке оценки.

  3. Используйте подобранные модели, чтобы предсказать ответы в горизонт прогноза.

  4. Оцените среднеквадратическую ошибку прогноза (FMSE) для каждой модели.

  5. Выберите модель с самым низким FMSE.

Создайте оценку и предскажите демонстрационные переменные для данных о предикторе и ответа. Задайте горизонт прогноза 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);

Подбирайте простую модель линейной регрессии

Подбирайте простую модель линейной регрессии к выборке оценки. Задайте имена переменных (имя переменной отклика должно быть последним элементом). Извлеките p-значения. Идентифицируйте незначительные коэффициенты при помощи 5%-го уровня значения.

SLRMdlFull = fitlm(X,y,'VarNames',DataTableLog.Properties.VariableNames)
slrPValue = SLRMdlFull.Coefficients.pValue;
isInSig = SLRMdlFull.CoefficientNames(slrPValue > 0.05)

SLRMdlFull = 


Linear regression model:
    UNRATE ~ [Linear formula with 14 terms in 13 predictors]

Estimated Coefficients:
                    Estimate       SE        tStat       pValue  
                    ________    ________    _______    __________

    (Intercept)       88.014      5.3229     16.535    2.6568e-37
    log_COE           7.1111      2.3426     3.0356     0.0027764
    log_CPIAUCSL     -3.4032      2.4611    -1.3828       0.16853
    log_GCE            -5.63      1.1581    -4.8615    2.6245e-06
    log_GDP           14.522      3.9406     3.6851    0.00030659
    log_GDPDEF        11.926      2.9298     4.0704    7.1598e-05
    log_GPDI         -2.5939     0.54603    -4.7504    4.2792e-06
    log_GS10         0.57467     0.29313     1.9605      0.051565
    log_HOANBS       -32.861      1.4507    -22.652    2.3116e-53
    log_M1SL         -3.3023     0.53185     -6.209     3.914e-09
    log_M2SL         -1.3845      1.0438    -1.3264       0.18647
    log_PCEC          -7.143      3.4302    -2.0824      0.038799
    FEDFUNDS        0.044542    0.047747    0.93287        0.3522
    TB3MS            -0.1577    0.064421     -2.448      0.015375


Number of observations: 185, Error degrees of freedom: 171
Root Mean Squared Error: 0.312
R-squared: 0.957,  Adjusted R-Squared: 0.954
F-statistic vs. constant model: 292, p-value = 2.26e-109

isInSig =

  1x4 cell array

    {'log_CPIAUCSL'}    {'log_GS10'}    {'log_M2SL'}    {'FEDFUNDS'}

SLRMdlFull LinearModel объект модели. P-значения предполагают что, на 5%-м уровне значения, CPIAUCSL, GS10, M2SL, и FEDFUNDS незначительные переменные.

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

idxreduced = ~ismember(predictornames,isInSig);
XReduced = X(:,idxreduced);
SLRMdlReduced = fitlm(XReduced,y,'VarNames',[predictornames(idxreduced) {'UNRATE'}]);

yFSLRF = predict(SLRMdlFull,XF);
yFSLRR = predict(SLRMdlReduced,XF(:,idxreduced));

fmseSLRF = sqrt(mean((yF - yFSLRF).^2))
fmseSLRR = sqrt(mean((yF - yFSLRR).^2))

fmseSLRF =

    0.6674


fmseSLRR =

    0.6105

FMSE упрощенной модели меньше FMSE полной модели, которая предполагает, что упрощенная модель делает вывод лучше.

Подбирайте частотную модель регрессии лассо

Подбирайте модель регрессии лассо к данным. Используйте путь к регуляризации по умолчанию. lasso стандартизирует данные по умолчанию. Поскольку переменные находятся по подобным шкалам, не стандартизируйте их. Возвратите информацию о подгонках вдоль пути к регуляризации.

[LassoBetaEstimates,FitInfo] = lasso(X,y,'Standardize',false,...
    'PredictorNames',predictornames);

По умолчанию, lasso соответствует временам модели 100 регрессии лассо при помощи последовательности значений для$\lambda$. Поэтому LassoBetaEstimates 13 100 матрица оценок коэффициента регрессии. Строки соответствуют переменным предикторам в X, и столбцы соответствуют значениям$\lambda$. FitInfo структура, которая включает значения$\lambda$ (FitInfo.Lambda), среднеквадратическая ошибка для каждой подгонки (FitInfo.MSE), и предполагаемые прерывания (FitInfo.Intercept).

Вычислите FMSE каждой модели, возвращенной lasso.

yFLasso = FitInfo.Intercept + XF*LassoBetaEstimates;
fmseLasso = sqrt(mean((yF - yFLasso).^2,1));

Постройте величину коэффициентов регрессии относительно значения уменьшения.

hax = lassoPlot(LassoBetaEstimates,FitInfo);
L1Vals = hax.Children.XData;
yyaxis right
h = plot(L1Vals,fmseLasso,'LineWidth',2,'LineStyle','--');
legend(h,'FMSE','Location','SW');
ylabel('FMSE');
title('Frequentist Lasso')

Модель с 7 предикторами (df = 7), кажется, балансирует минимальный FMSE и сложность модели хорошо. Получите коэффициенты, соответствующие модели, содержащей 7 предикторов и дающей к минимальному FMSE.

fmsebestlasso = min(fmseLasso(FitInfo.DF == 7));
idx = fmseLasso == fmsebestlasso;
bestLasso = [FitInfo.Intercept(idx); LassoBetaEstimates(:,idx)];
table(bestLasso,'RowNames',["Intercept" predictornames])

ans =

  14x1 table

                    bestLasso
                    _________

    Intercept         61.154 
    log_COE         0.042061 
    log_CPIAUCSL           0 
    log_GCE                0 
    log_GDP                0 
    log_GDPDEF         8.524 
    log_GPDI               0 
    log_GS10          1.6957 
    log_HOANBS       -18.937 
    log_M1SL         -1.3365 
    log_M2SL         0.17948 
    log_PCEC               0 
    FEDFUNDS               0 
    TB3MS           -0.14459

Частотный анализ лассо предполагает что переменные CPIAUCSL, GCE, GDP, GPDI, PCEC, и FEDFUNDS или незначительны или избыточны.

Подбирайте байесовую модель регрессии лассо

В представлении Bayesian регрессии лассо предшествующим распределением коэффициентов регрессии является Лаплас (удвойте экспоненциал), со средним значением 0 и шкалой$\sigma/\lambda$, где$\lambda$ фиксированный параметр уменьшения и$\sigma^2\sim IG(A,B)$. Для получения дополнительной информации смотрите lassoblm.

Как с частотным представлением регрессии лассо, если вы увеличиваетесь$\lambda$, затем разреженность получившейся модели увеличивается монотонно. Однако различающееся частотное лассо, Байесово лассо имеет незначительные или избыточные содействующие режимы, которые не являются точно 0. Вместо этого оценки имеют апостериорное распределение и незначительны или избыточны, если область приблизительно 0 имеют высокую плотность.

Когда вы реализуете Байесовую регрессию лассо в MATLAB®, знать о нескольких различиях между Statistics and Machine Learning Toolbox™ функционируют lasso и объект Econometrics Toolbox™ lassoblm и его присоединенные функции.

  • lassoblm часть объектной среды, тогда как lasso функция. Объектная среда оптимизировала эконометрические рабочие процессы.

  • В отличие от lasso, lassoblm не стандартизирует данные о предикторе. Однако можно предоставить различные значения уменьшения для каждого коэффициента. Эта функция помогает обеспечить interpretability содействующих оценок.

  • lassoblm применяет одно значение уменьшения к каждому коэффициенту; это не принимает путь к регуляризации как lasso.

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

Создайте Байесовую регрессию лассо предшествующая модель при помощи bayeslm. Задайте количество переменных предикторов и имен переменных. Отобразите значение уменьшения по умолчанию для каждого коэффициента, сохраненного в Lambda свойство модели.

PriorMdl = bayeslm(p,'ModelType','lasso','VarNames',predictornames);
table(PriorMdl.Lambda,'RowNames',PriorMdl.VarNames)

ans =

  14x1 table

                    Var1
                    ____

    Intercept       0.01
    log_COE            1
    log_CPIAUCSL       1
    log_GCE            1
    log_GDP            1
    log_GDPDEF         1
    log_GPDI           1
    log_GS10           1
    log_HOANBS         1
    log_M1SL           1
    log_M2SL           1
    log_PCEC           1
    FEDFUNDS           1
    TB3MS              1

PriorMdl lassoblm объект модели. lassoblm приписывает уменьшение 1 для каждого коэффициента кроме прерывания, которое имеет уменьшение 0.01. Эти значения по умолчанию не могут быть полезными в большинстве приложений; лучшая практика состоит в том, чтобы экспериментировать со многими значениями.

Считайте путь к регуляризации возвращенным lasso. Раздуйте значения уменьшения на коэффициент$n/\sqrt(MSE_j)$, куда$MSE_j$ MSE лассо запущен$j$,$j$ = 1 через количество значений уменьшения.

ismissing = any(isnan(X),2);
n = sum(~ismissing); % Effective sample size
lambda = FitInfo.Lambda*n./sqrt(FitInfo.MSE);

Считайте путь к регуляризации возвращенным lasso. Цикл через значения уменьшения в каждой итерации:

  1. Оцените апостериорные распределения коэффициентов регрессии и отклонения воздействия, учитывая уменьшение и данные. Поскольку шкалы предикторов близки, приписывают то же уменьшение каждому предиктору и приписывают уменьшение 0.01 к прерыванию. Сохраните следующие средние значения и 95%-е вероятные интервалы.

  2. Для графика лассо, если 95%-й вероятный интервал включает 0, то обнуленный соответствующее следующее среднее значение.

  3. Предскажите в горизонт прогноза.

  4. Оцените FMSE.

numlambda = numel(lambda);

% Preallocate
BayesLassoCoefficients = zeros(p+1,numlambda);
BayesLassoCI95 = zeros(p+1,2,numlambda);
fmseBayesLasso = zeros(numlambda,1);
BLCPlot = zeros(p+1,numlambda);

% Estimate and forecast
rng(10); % For reproducibility
for j = 1:numlambda
    PriorMdl.Lambda = lambda(j);
    [EstMdl,Summary] = estimate(PriorMdl,X,y,'Display',false);
    BayesLassoCoefficients(:,j) = Summary.Mean(1:(end - 1));
    BLCPlot(:,j) = Summary.Mean(1:(end - 1));
    BayesLassoCI95(:,:,j) = Summary.CI95(1:(end - 1),:);
    idx = BayesLassoCI95(:,2,j) > 0 & BayesLassoCI95(:,1,j) <= 0;
    BLCPlot(idx,j) = 0;
    yFBayesLasso = forecast(EstMdl,XF);
    fmseBayesLasso(j) = sqrt(mean((yF - yFBayesLasso).^2));
end

В каждой итерации, estimate запускает сэмплер Гиббса к выборке от полных условных выражений (см. Аналитически Послушное Последующее поколение), и возвращает empiricalblm модель EstMdl, который содержит ничьи и сводную таблицу оценки Summary. Можно задать параметры для сэмплера Гиббса, такие как количество ничьих или утончающейся схемы.

Хорошая практика должна диагностировать следующую выборку путем осмотра графиков трассировки. Однако, потому что 100 распределений чертились, этот пример доходы, не выполняя тот шаг.

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

L1Vals = sum(abs(BLCPlot(2:end,:)),1)/max(sum(abs(BLCPlot(2:end,:)),1));

figure;
plot(L1Vals,BLCPlot(2:end,:))
xlabel('L1');
ylabel('Coefficient Estimates');
yyaxis right
h = plot(L1Vals,fmseBayesLasso,'LineWidth',2,'LineStyle','--');
legend(h,'FMSE','Location','SW');
ylabel('FMSE');
title('Bayesian Lasso')

Модель имеет тенденцию делать вывод лучше как нормированные увеличения штрафа L1 прошлые 0.1. Чтобы сбалансировать минимальный FMSE и сложность модели, выберите самую простую модель с FMSE, самым близким к 0,68.

idx = find(fmseBayesLasso > 0.68,1);
fmsebestbayeslasso = fmseBayesLasso(idx);
bestBayesLasso = BayesLassoCoefficients(:,idx);
bestCI95 = BayesLassoCI95(:,:,idx);
table(bestBayesLasso,bestCI95,'RowNames',["Intercept" predictornames])

ans =

  14x2 table

                    bestBayesLasso           bestCI95       
                    ______________    ______________________

    Intercept            90.587          79.819       101.62
    log_COE              6.5591          1.8093       11.239
    log_CPIAUCSL        -2.2971         -6.9019       1.7057
    log_GCE             -5.1707         -7.4902      -2.8583
    log_GDP              9.8839          2.3848       17.672
    log_GDPDEF           10.281          5.1677       15.936
    log_GPDI            -2.0973         -3.2168     -0.96728
    log_GS10            0.82485         0.22748       1.4237
    log_HOANBS          -32.538         -35.589      -29.526
    log_M1SL            -3.0656         -4.1499      -2.0066
    log_M2SL            -1.1007         -3.1243      0.75844
    log_PCEC            -3.3342         -9.6496       1.8482
    FEDFUNDS           0.043672       -0.056192      0.14254
    TB3MS              -0.15682        -0.29385    -0.022145

Один способ определить, незначительна ли переменная или избыточная, путем проверки, является ли 0 в его соответствующем коэффициенте 95%-м вероятным интервалом. Используя этот критерий, CPIAUCSL, M2SL, PCEC, и FEDFUNDS незначительны.

Отобразите все оценки в таблице для сравнения.

SLRR = zeros(p + 1,1);
SLRR([true idxreduced]) = SLRMdlReduced.Coefficients.Estimate;
table([SLRMdlFull.Coefficients.Estimate; fmseSLRR],...
    [SLRR; fmseSLRR],...
    [bestLasso; fmsebestlasso],...
    [bestBayesLasso; fmsebestbayeslasso],'VariableNames',...
    {'SLRFull' 'SLRReduced' 'Lasso' 'BayesLasso'},...
    'RowNames',[PriorMdl.VarNames; 'MSE'])

ans =

  15x4 table

                    SLRFull     SLRReduced     Lasso      BayesLasso
                    ________    __________    ________    __________

    Intercept         88.014       88.327       61.154       90.587 
    log_COE           7.1111       10.854     0.042061       6.5591 
    log_CPIAUCSL     -3.4032            0            0      -2.2971 
    log_GCE            -5.63      -6.1958            0      -5.1707 
    log_GDP           14.522       16.701            0       9.8839 
    log_GDPDEF        11.926       9.1106        8.524       10.281 
    log_GPDI         -2.5939      -2.6963            0      -2.0973 
    log_GS10         0.57467            0       1.6957      0.82485 
    log_HOANBS       -32.861      -33.782      -18.937      -32.538 
    log_M1SL         -3.3023      -2.8099      -1.3365      -3.0656 
    log_M2SL         -1.3845            0      0.17948      -1.1007 
    log_PCEC          -7.143      -14.207            0      -3.3342 
    FEDFUNDS        0.044542            0            0     0.043672 
    TB3MS            -0.1577    -0.078449     -0.14459     -0.15682 
    MSE              0.61048      0.61048      0.79425      0.69639

После того, как вы выбрали модель, повторно оцените его с помощью целого набора данных. Например, чтобы создать прогнозирующую Байесовую модель регрессии лассо, создайте предшествующую модель и задайте уменьшение, дающее к самой простой модели с минимальным FMSE, затем оцените его с помощью целого набора данных.

bestLambda = lambda(idx);
PriorMdl = bayeslm(p,'ModelType','lasso','Lambda',bestLambda,...
    'VarNames',predictornames);
PosteriorMdl = estimate(PriorMdl,[X; XF],[y; yF]);

Method: lasso MCMC sampling with 10000 draws
Number of observations: 201
Number of predictors:   14
 
              |   Mean      Std          CI95         Positive  Distribution 
-----------------------------------------------------------------------------
 Intercept    |  85.9643  5.2710   [75.544, 96.407]     1.000     Empirical  
 log_COE      |  12.3574  2.1817   [ 8.001, 16.651]     1.000     Empirical  
 log_CPIAUCSL |   0.1498  1.9150   [-3.733,  4.005]     0.525     Empirical  
 log_GCE      |  -7.1850  1.0556   [-9.208, -5.101]     0.000     Empirical  
 log_GDP      |  11.8648  3.9955   [ 4.111, 19.640]     0.999     Empirical  
 log_GDPDEF   |   8.8777  2.4221   [ 4.033, 13.745]     1.000     Empirical  
 log_GPDI     |  -2.5884  0.5294   [-3.628, -1.553]     0.000     Empirical  
 log_GS10     |   0.6910  0.2577   [ 0.194,  1.197]     0.997     Empirical  
 log_HOANBS   | -32.3356  1.4941  [-35.276, -29.433]    0.000     Empirical  
 log_M1SL     |  -2.2279  0.5043   [-3.202, -1.238]     0.000     Empirical  
 log_M2SL     |   0.3617  0.9123   [-1.438,  2.179]     0.652     Empirical  
 log_PCEC     | -11.3018  3.0353  [-17.318, -5.252]     0.000     Empirical  
 FEDFUNDS     |  -0.0132  0.0490   [-0.109,  0.082]     0.392     Empirical  
 TB3MS        |  -0.1128  0.0660   [-0.244,  0.016]     0.042     Empirical  
 Sigma2       |   0.1186  0.0122   [ 0.097,  0.145]     1.000     Empirical

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

| |

Похожие темы

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

Поддержка

Для просмотра документации необходимо авторизоваться на сайте