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 вектора-столбца таковых.

  • $\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$ минимизируется относительно коэффициентов регрессии. Получившиеся коэффициенты являются оценками lasso. Для получения дополнительной информации о частом взгляде на регрессию лассо смотрите [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 имеет высокую плотность.

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

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

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

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

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

Создайте байесову регрессию лассо предыдущую модель при помощи 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. Для lasso plot, если 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

Поддержка