Байесовская модель линейной регрессии с регуляризацией лассо
Объект модели байесовской линейной регрессии lassoblm задает совместное предварительное распределение коэффициентов регрессии и дисперсии возмущений (β, start2) для реализации байесовской регрессии лассо
[1]. Для j = 1,...,NumPredictors, условное предшествующее распределение βj 'start2 - это распределение Лапласа (двойное экспоненциальное) со средним значением 0 и шкалой λ 2/λ, где λ - параметр регуляризации лассо, или усадки. Предыдущее распределение start2 - обратная гамма с формой A и шкалой B.
Правдоподобие данных равно start2), где (yt; xtβ, start2) - гауссова плотность вероятностей, оцениваемая на yt со средним xtβ и дисперсией start2. Полученное апостериорное распределение не поддается аналитическому отслеживанию. Для получения подробной информации о заднем распределении см. Аналитически отслеживаемые задние части.
Как правило, при создании объекта модели байесовской линейной регрессии задается совместное предварительное распределение и характеристики только модели линейной регрессии. То есть объект модели является шаблоном, предназначенным для дальнейшего использования. В частности, для включения данных в модель для последующего анализа распределения и выбора признаков передайте объект модели и данные соответствующей функции объекта.
PriorMdl = lassoblm(NumPredictors)PriorMdl) состоит из NumPredictors предикторы и перехват, и устанавливает NumPredictors собственность. Совместное предварительное распределение (β, start2) подходит для реализации байесовской регрессии лассо [1]. PriorMdl является шаблоном, который определяет предыдущие распределения и задает значения параметра регуляризации лассо λ и размерность β.
estimate | Выполнение выбора переменных предиктора для байесовских моделей линейной регрессии |
simulate | Моделирование коэффициентов регрессии и дисперсии возмущений байесовской модели линейной регрессии |
forecast | Прогнозные отклики байесовской модели линейной регрессии |
plot | Визуализация предыдущих и задних плотностей параметров байесовской модели линейной регрессии |
summarize | Сводная статистика распределения байесовской модели линейной регрессии для выбора переменных предиктора |
Рассмотрим модель множественной линейной регрессии, которая предсказывает реальный валовой национальный продукт США (GNPR) с использованием линейной комбинации индекса промышленного производства (IPI), общая занятость (E) и реальная заработная плата (WR).
β3WRt + αt.
Для всех - это ряд независимых гауссовых возмущений со средним значением 0 и дисперсией .
Предположим, что эти предыдущие распределения:
Для j = 0..., 3, есть лапласовское распределение со средним из 0 и масштабом , где - параметр сжатия. Коэффициенты условно независимы.
B А В - форма и масштаб, соответственно, обратного гамма-распределения.
Создайте предыдущую модель для байесовской линейной регрессии. Укажите количество предикторов p.
p = 3; Mdl = lassoblm(p);
PriorMdl является lassoblm Объект байесовской модели линейной регрессии, представляющий предварительное распределение коэффициентов регрессии и дисперсии возмущений. В окне команд lassoblm отображает сводку предыдущих распределений.
Кроме того, можно создать предыдущую модель для байесовской регрессии лассо, передав число предикторов bayeslm и установка ModelType аргумент пары имя-значение для 'lasso'.
MdlBayesLM = bayeslm(p,'ModelType','lasso')
MdlBayesLM =
lassoblm with properties:
NumPredictors: 3
Intercept: 1
VarNames: {4x1 cell}
Lambda: [4x1 double]
A: 3
B: 1
| Mean Std CI95 Positive Distribution
---------------------------------------------------------------------------
Intercept | 0 100 [-200.000, 200.000] 0.500 Scale mixture
Beta(1) | 0 1 [-2.000, 2.000] 0.500 Scale mixture
Beta(2) | 0 1 [-2.000, 2.000] 0.500 Scale mixture
Beta(3) | 0 1 [-2.000, 2.000] 0.500 Scale mixture
Sigma2 | 0.5000 0.5000 [ 0.138, 1.616] 1.000 IG(3.00, 1)
Mdl и MdlBayesLM являются эквивалентными объектами модели.
Значения свойств создаваемых моделей, доступные для записи, можно задать с помощью точечной нотации. Установите имена коэффициентов регрессии в соответствующие имена переменных.
Mdl.VarNames = ["IPI" "E" "WR"]
Mdl =
lassoblm with properties:
NumPredictors: 3
Intercept: 1
VarNames: {4x1 cell}
Lambda: [4x1 double]
A: 3
B: 1
| Mean Std CI95 Positive Distribution
---------------------------------------------------------------------------
Intercept | 0 100 [-200.000, 200.000] 0.500 Scale mixture
IPI | 0 1 [-2.000, 2.000] 0.500 Scale mixture
E | 0 1 [-2.000, 2.000] 0.500 Scale mixture
WR | 0 1 [-2.000, 2.000] 0.500 Scale mixture
Sigma2 | 0.5000 0.5000 [ 0.138, 1.616] 1.000 IG(3.00, 1)
MATLAB ® связывает имена переменных с коэффициентами регрессии в дисплеях.
Рассмотрим модель линейной регрессии в разделе Создание предыдущей модели для байесовской регрессии лассо.
Создайте предыдущую модель для выполнения байесовской регрессии лассо. Укажите количество предикторов p и названия коэффициентов регрессии.
p = 3; PriorMdl = bayeslm(p,'ModelType','lasso','VarNames',["IPI" "E" "WR"]); shrinkage = PriorMdl.Lambda
shrinkage = 4×1
0.0100
1.0000
1.0000
1.0000
PriorMdl сохраняет значения усадки для всех предикторов в Lambda собственность. shrinkage(1) - усадка для перехвата и элементы shrinkage(2:end) соответствуют коэффициентам предикторов в Mdl.VarNames. Усадка по умолчанию для перехвата: 0.01, и по умолчанию 1 для всех остальных коэффициентов.
Загрузите набор данных Нельсона-Плоссера. Создайте переменные для последовательности ответа и предиктора. Поскольку лассо чувствителен к переменным масштабам, стандартизируйте все переменные.
load Data_NelsonPlosser X = DataTable{:,PriorMdl.VarNames(2:end)}; y = DataTable{:,'GNPR'}; X = (X - mean(X,'omitnan'))./std(X,'omitnan'); y = (y - mean(y,'omitnan'))/std(y,'omitnan');
Хотя в этом примере стандартизируются переменные, можно задать различные значения усадки для каждого коэффициента, установив Lambda имущество PriorMdl к числовому вектору значений усадки.
Реализуйте байесовскую регрессию лассо, оценивая краевые апостериорные распределения и . Поскольку байесовская регрессия лассо использует цепь Маркова Монте-Карло (MCMC) для оценки, установите начальное число случайных чисел, чтобы воспроизвести результаты.
rng(1); PosteriorMdl = estimate(PriorMdl,X,y);
Method: lasso MCMC sampling with 10000 draws
Number of observations: 62
Number of predictors: 4
| Mean Std CI95 Positive Distribution
-----------------------------------------------------------------------
Intercept | -0.4490 0.0527 [-0.548, -0.344] 0.000 Empirical
IPI | 0.6679 0.1063 [ 0.456, 0.878] 1.000 Empirical
E | 0.1114 0.1223 [-0.110, 0.365] 0.827 Empirical
WR | 0.2215 0.1367 [-0.024, 0.494] 0.956 Empirical
Sigma2 | 0.0343 0.0062 [ 0.024, 0.048] 1.000 Empirical
PosteriorMdl является empiricalblm объект модели, хранящий черпания из задних распределений и , учитывая данные. estimate отображает сводку по краевым задним распределениям в командной строке. Строки сводки соответствуют коэффициентам регрессии и дисперсии возмущений, а столбцы - характеристикам заднего распределения. Характеристики включают в себя:
CI95, который содержит 95% байесовских равных достоверных интервалов для параметров. Например, апостериорная вероятность того, что коэффициент регрессии E (стандартизировано) в [-0.110, 0,365] равно 0,95.
Positive, которая содержит апостериорную вероятность того, что параметр больше 0. Например, вероятность того, что перехват больше 0, равна 0.
По умолчанию estimate рисует и отбрасывает загоревшую выборку размером 5000. Тем не менее, хорошей практикой является проверка следового графика розыгрышей на предмет адекватного смешивания и отсутствия скороспелости. Постройте график трассировки черчений для каждого параметра. Доступ к чертежам, составляющим распределение (свойства BetaDraws и Sigma2Draws) с использованием точечной нотации.
figure; for j = 1:(p + 1) subplot(2,2,j); plot(PosteriorMdl.BetaDraws(j,:)); title(sprintf('%s',PosteriorMdl.VarNames{j})); end
figure;
plot(PosteriorMdl.Sigma2Draws);
title('Sigma2');
Графики следов показывают, что розыгрыши, по-видимому, хорошо смешиваются. График не показывает обнаруживаемой переходности или последовательной корреляции, и розыгрыши не прыгают между состояниями.
Постройте график апостериорных распределений коэффициентов и дисперсии возмущений.
figure; plot(PosteriorMdl)
E и WR может не быть важными предикторами, потому что 0 находится в области высокой плотности в их задних распределениях.
Рассмотрим модель линейной регрессии в разделе Создание предыдущей модели для байесовской регрессии лассо и ее реализацию в разделе Выполнение выбора переменных с использованием усадки лассо по умолчанию.
При реализации регрессии лассо обычной практикой является стандартизация переменных. Однако, если требуется сохранить интерпретацию коэффициентов, но переменные имеют разные масштабы, то можно выполнить дифференциальную усадку, указав различную усадку для каждого коэффициента.
Создайте предыдущую модель для выполнения байесовской регрессии лассо. Укажите количество предикторов p и названия коэффициентов регрессии.
p = 3; PriorMdl = bayeslm(p,'ModelType','lasso','VarNames',["IPI" "E" "WR"]);
Загрузите набор данных Нельсона-Плоссера. Создайте переменные для последовательности ответа и предиктора. Определите, имеют ли переменные экспоненциальные тренды, построив каждый график на отдельных рисунках.
load Data_NelsonPlosser X = DataTable{:,PriorMdl.VarNames(2:end)}; y = DataTable{:,'GNPR'}; figure; plot(dates,y) title('GNPR')
for j = 1:3 figure; plot(dates,X(:,j)); title(PriorMdl.VarNames(j + 1)); end
Переменные GNPR, IPI, и WR похоже, имеет экспоненциальный тренд.
Удаление экспоненциального тренда из переменных GNPR, IPI, и WR.
y = log(y); X(:,[1 3]) = log(X(:,[1 3]));
Все переменные предиктора имеют разные масштабы (для получения дополнительной информации введите Description в командной строке). Отображение среднего значения для каждого предиктора. Удаление наблюдений, содержащих ведущие отсутствующие значения, из всех предикторов.
predmeans = mean(X,'omitnan')predmeans = 1×3
104 ×
0.0002 4.7700 0.0004
Значения второго предиктора намного больше, чем значения других двух предикторов и ответа. Следовательно, коэффициент регрессии второго предсказателя может казаться близким к нулю.
Используя точечную нотацию, отнесите очень низкую усадку к перехвату, усадку 0,1 к первому и третьему предикторам и усадку 1000 ко второму предиктору.
PriorMdl.Lambda = [1e-5 0.1 1e4 0.1];
Реализуйте байесовскую регрессию лассо, оценивая краевые апостериорные распределения и . Поскольку байесовская регрессия лассо использует MCMC для оценки, задайте начальное число случайных чисел для воспроизведения результатов.
rng(1); PosteriorMdl = estimate(PriorMdl,X,y);
Method: lasso MCMC sampling with 10000 draws
Number of observations: 62
Number of predictors: 4
| Mean Std CI95 Positive Distribution
----------------------------------------------------------------------
Intercept | 2.0281 0.6839 [ 0.679, 3.323] 0.999 Empirical
IPI | 0.3534 0.2497 [-0.139, 0.839] 0.923 Empirical
E | 0.0000 0.0000 [-0.000, 0.000] 0.762 Empirical
WR | 0.5250 0.3482 [-0.126, 1.209] 0.937 Empirical
Sigma2 | 0.0315 0.0055 [ 0.023, 0.044] 1.000 Empirical
Рассмотрим модель линейной регрессии в разделе Создание предыдущей модели для байесовской регрессии лассо.
Выполнить байесовскую регрессию лассо:
Создайте предыдущую модель байесовского лассо для коэффициентов регрессии и дисперсии возмущений. Используйте усадку по умолчанию.
Удерживайте последние 10 периодов данных из оценки.
Оцените краевые задние распределения.
p = 3; PriorMdl = bayeslm(p,'ModelType','lasso','VarNames',["IPI" "E" "WR"]); load Data_NelsonPlosser fhs = 10; % Forecast horizon size X = DataTable{1:(end - fhs),PriorMdl.VarNames(2:end)}; y = DataTable{1:(end - fhs),'GNPR'}; XF = DataTable{(end - fhs + 1):end,PriorMdl.VarNames(2:end)}; % Future predictor data yFT = DataTable{(end - fhs + 1):end,'GNPR'}; % True future responses rng(1); % For reproducibility PosteriorMdl = estimate(PriorMdl,X,y,'Display',false);
Ответы прогноза с использованием апостериорного прогностического распределения и будущих данных предиктора XF. Постройте график истинных значений ответа и прогнозируемых значений.
yF = forecast(PosteriorMdl,XF); figure; plot(dates,DataTable.GNPR); hold on plot(dates((end - fhs + 1):end),yF) h = gca; hp = patch([dates(end - fhs + 1) dates(end) dates(end) dates(end - fhs + 1)],... h.YLim([1,1,2,2]),[0.8 0.8 0.8]); uistack(hp,'bottom'); legend('Forecast Horizon','True GNPR','Forecasted GNPR','Location','NW') title('Real Gross National Product: 1909 - 1970'); ylabel('rGNP'); xlabel('Year'); hold off
yF является вектором 10 на 1 будущих значений реального GNP, соответствующего будущим данным предиктора.
Оценка прогнозируемой среднеквадратичной ошибки (RMSE).
frmse = sqrt(mean((yF - yFT).^2))
frmse = 25.4831
RMSE прогноза является относительным показателем точности прогноза. В частности, можно оценить несколько моделей с использованием различных допущений. Модель с самым низким прогнозом RMSE является лучшей из сравниваемых моделей.
При выполнении байесовской регрессии лассо рекомендуется искать соответствующие значения усадки. Один из способов сделать это - оценить прогнозный RMSE по сетке значений усадки и выбрать усадку, которая минимизирует прогнозный RMSE.
Lambda является параметром настройки. Поэтому выполните байесовскую регрессию лассо, используя сетку значений усадки, и выберите модель, которая наилучшим образом уравновешивает критерий соответствия и сложность модели.
Для оценки, моделирования и прогнозирования MATLAB ® не стандартизирует данные предиктора. Если переменные в данных предиктора имеют различные масштабы, то укажите параметр усадки для каждого предиктора, предоставив числовой вектор дляLambda.
bayeslm может создавать любой поддерживаемый объект предыдущей модели для байесовской линейной регрессии.
[1] Парк, T. и Г. Казелла. «Байесовский Лассо.» Журнал Американской статистической ассоциации. т. 103, № 482, 2008, стр. 681-686.