Рассмотрим множественную линейную регрессионую модель, которая предсказывает реальный валовой национальный продукт США (GNPR) с использованием линейной комбинации индекса промышленного производства (IPI), общая занятость (E), и реальная заработная плата (WR).

Для всех $$t$$, $${\epsilon }_t$$ - серия независимых гауссовских нарушений порядка со средним значением 0 и отклонением $${\sigma }^2$$.

Предположим, что эти предыдущие распределения:

Создайте предыдущую модель для байесовской линейной регрессии. Задайте количество предикторов 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 являются эквивалентными объектами модели.

Можно задать значения свойств writable для созданных моделей с помощью записи через точку. Установите имена коэффициентов регрессии в соответствующие имена переменных.

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 ® связывает имена переменных с коэффициентами регрессии на отображениях.

Рассмотрим линейную регрессионую модель в Create Private Model для регрессии Лассо Байесова.

Создайте предыдущую модель для выполнения регрессии Бейесова лассо. Задайте количество предикторов 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 для всех других коэффициентов.

Загрузите набор данных Нельсона-Плоссера. Создайте переменные для ряда отклика и предиктора. Поскольку lasso чувствителен к переменным шкалам, стандартизируйте все переменные.

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 в числовой вектор значений усадки.

Реализуйте регрессию байесовского лассо путем оценки маргинальных апостериорных распределений $$\beta$$ и $${\sigma }^2$$. Поскольку байесовская регрессия лассо использует марковскую цепь Monte Carlo (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 объект модели, который хранит черты из апостериорных распределений $$\beta$$ и $${\sigma }^2$$ учитывая данные. estimate отображает сводные данные маргинальных апостериорных распределений в командной строке. Строки сводных данных соответствуют коэффициентам регрессии и нарушения порядка отклонения, а столбцы соответствуют характеристикам апостериорного распределения. Характеристики включают:

По умолчанию 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 находится в области высокой плотности в их апостериорных распределениях.

Рассмотрим линейную регрессионую модель в Create Private Model для регрессии Бейесова Лассо и ее реализацию в Performable Selection Using Default Lasso Shrinkage.

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

Создайте предыдущую модель для выполнения регрессии Бейесова лассо. Задайте количество предикторов 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];

Реализуйте регрессию байесовского лассо путем оценки маргинальных апостериорных распределений $$\beta$$ и $${\sigma }^2$$. Поскольку байесовская регрессия лассо использует 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  
 

Рассмотрим линейную регрессионую модель в Create Private Model для регрессии Лассо Байесова.

Выполните байесовскую регрессию лассо:

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.