Exponenta Banner
exponenta event banner

Байесовская регрессия Лассо

Открыть сценарий

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

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

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

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

  • $y$ является вектором n ответов.

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

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

  • $\beta_0$ - перехват.

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

  • $\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 переменных, измеренных с 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 = 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 являются незначительными переменными.

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

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 являются либо незначительными, либо избыточными.

Подгонка байесовской регрессионной модели лассо

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

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

При внедрении байесовской регрессии лассо в MATLAB ® следует помнить о нескольких различиях между функцией статистики и функцией Toolbox™ машинного обученияlasso и объект Toolbox™ Econometrics lassoblm и связанные с этим функции.

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

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

  • 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$1$j$ от числа значений усадки.

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 постройте график FMSE.

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

Поддержка