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

\forall $$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 эквивалентные объекты модели.

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

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

Реализуйте Байесовую регрессию лассо путем оценки крайних апостериорных распределений $$\beta$$ и $${\sigma }^2$$. Поскольку Байесова регрессия лассо использует Цепь Маркова Монте-Карло (MCMC) для оценки, установите seed случайных чисел воспроизводить результаты.

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

Полагайте, что модель линейной регрессии в Создает Предшествующую Модель для Байесовой Регрессии Лассо, и ее реализация в Выполняют Выбор переменной Используя Уменьшение Лассо По умолчанию.

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

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

PriorMdl.Lambda = [1e-5 0.1 1e4 0.1];

Реализуйте Байесовую регрессию лассо путем оценки крайних апостериорных распределений $$\beta$$ и $${\sigma }^2$$. Поскольку Байесова регрессия лассо использует MCMC для оценки, установите seed случайных чисел воспроизводить результаты.

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  
 

Полагайте, что модель линейной регрессии в Создает Предшествующую Модель для Байесовой Регрессии Лассо.

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

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.