Этот пример показывает, как выполнить выбор переменной при помощи Байесовой регрессии лассо.
Регрессия лассо является методом линейной регрессии, который комбинирует регуляризацию и выбор переменной. Регуляризация помогает предотвратить сверхподбор кривой путем уменьшения значения коэффициентов регрессии. Частотное представление регрессии лассо отличается от того из других методов регуляризации, такой как, гребенчатая регрессия, потому что лассо приписывает значение точно 0 к коэффициентам регрессии, соответствующим предикторам, которые незначительны или избыточны.
Считайте это несколькими модель линейной регрессии:
вектор n ответов.
матрица p соответствующих наблюдаемых переменных прогноза.
вектор p коэффициентов регрессии.
прерывание.
длина n вектор-столбец из единиц.
случайный вектор iid Гауссовых воздействий со средним значением 0 и отклонением
.
Целевая функция частотного представления регрессии лассо
штраф (или уменьшение) термин и, в отличие от других параметров, не является подходящим к данным. Необходимо задать значение для перед оценкой, и лучшая практика состоит в том, чтобы настроить ее. После того, как вы зададите значение,
минимизирован относительно коэффициентов регрессии. Получившиеся коэффициенты являются оценками лассо. Для получения дополнительной информации на частотном представлении регрессии лассо, см. [108].
В данном примере рассмотрите создание прогнозирующей линейной модели для уровня безработицы США. Вы хотите модель, которая делает вывод хорошо. Другими словами, вы хотите минимизировать сложность модели путем удаления всех избыточных предикторов и всех предикторов, которые являются некоррелироваными с уровнем безработицы.
Загрузите США макроэкономический набор данных 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');
Переменные имеют довольно подобные шкалы.
Чтобы определить, как хорошо модели временных рядов делают вывод, выполните эту процедуру:
Разделите данные в оценку и предскажите выборки.
Соответствуйте моделям к выборке оценки.
Используйте подобранные модели, чтобы предсказать ответы в горизонт прогноза.
Оцените среднеквадратическую ошибку прогноза (FMSE) для каждой модели.
Выберите модель с самым низким 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 регрессии лассо при помощи последовательности значений для. Поэтому LassoBetaEstimates 13 100 матрица оценок коэффициента регрессии. Строки соответствуют переменным прогноза в X, и столбцы соответствуют значениям. FitInfo является структурой, которая включает значения (FitInfo.Lambda), среднеквадратической ошибки для каждой подгонки (FitInfo.MSE) и предполагаемые прерывания (FitInfo.Intercept).
Вычислите FMSE каждой модели, возвращенной lasso.
yFLasso = FitInfo.Intercept + XF*LassoBetaEstimates; fmseLasso = sqrt(mean((yF - yFLasso).^2,1));
Постройте значение коэффициентов регрессии относительно значения уменьшения.
figure; 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 и шкалой
, где фиксированный параметр уменьшения и
. Для получения дополнительной информации смотрите lassoblm.
Как с частотным представлением регрессии лассо, если вы увеличиваетесь, затем разреженность получившейся модели увеличивается монотонно. Однако различающееся частотное лассо, Байесово лассо имеет незначительные или избыточные содействующие режимы, которые не являются точно 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. Раздуйте значения уменьшения фактором
, куда
MSE лассо запущен
, = 1 через количество значений уменьшения.
ismissing = any(isnan(X),2);
n = sum(~ismissing); % Effective sample size
lambda = FitInfo.Lambda*n./sqrt(FitInfo.MSE);
Считайте путь к регуляризации возвращенным lasso. Цикл через значения уменьшения в каждой итерации:
Оцените апостериорные распределения коэффициентов регрессии и отклонения воздействия, учитывая уменьшение и данные. Поскольку шкалы предикторов близки, приписывают то же уменьшение каждому предиктору и приписывают уменьшение 0.01 к прерыванию. Сохраните следующие средние значения и 95%-е вероятные интервалы.
Для графика лассо, если 95%-й вероятный интервал включает 0, то обнуленный соответствующее следующее среднее значение.
Предскажите в горизонт прогноза.
Оцените 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
estimate | forecast | simulate