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

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

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

Создайте предыдущую модель регрессии лассо Байеса. Задайте количество предикторов, тип предыдущей модели и имена переменных. Задайте следующие усадки:

Порядок усадок соответствует порядку заданных имен переменных, но первый элемент является усадкой точки пересечения.

p = 3;
PriorMdl = bayeslm(p,'ModelType','lasso','Lambda',[0.01; 10; 1e5; 10],...
    'VarNames',["IPI" "E" "WR"]);

PriorMdl является lassoblm Байесовский объект линейной регрессионной модели, представляющий предшествующее распределение коэффициентов регрессии и отклонение нарушения порядка.

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

load Data_NelsonPlosser
X = DataTable{:,PriorMdl.VarNames(2:end)};
y = DataTable{:,"GNPR"};

Выполните регрессию байесовского лассо, передав предыдущую модель и данные в estimate, то есть путем оценки апостериорного распределения $\beta$ и ${\sigma }^2$. Байесовская регрессия лассо использует марковскую цепь Monte Carlo (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 | -1.3472   6.8160  [-15.169, 11.590]    0.427     Empirical  
 IPI       |  4.4755   0.1646   [ 4.157,  4.799]    1.000     Empirical  
 E         |  0.0001   0.0002   [-0.000,  0.000]    0.796     Empirical  
 WR        |  3.1610   0.3136   [ 2.538,  3.760]    1.000     Empirical  
 Sigma2    | 60.1452  11.1180   [42.319, 85.085]    1.000     Empirical  
 

PosteriorMdl является empiricalblm объект модели, который хранит черты из апостериорных распределений $$\beta$$ и $${\sigma }^2$$ учитывая данные. estimate отображает сводные данные маргинальных апостериорных распределений в командной строке MATLAB ®. Строки сводных данных соответствуют коэффициентам регрессии и нарушения порядка отклонения, а столбцы соответствуют характеристикам апостериорного распределения. Характеристики включают:

Постройте график апостериорных распределений.

plot(PosteriorMdl)

Учитывая усадку, распределение E довольно плотный около 0. Поэтому E возможно, не является важным предиктором.

По умолчанию 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');

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

Рассмотрим регрессионую модель в Выборе Переменных с Использованием Регрессии Бейесова Лассо.

Создайте предыдущую модель для выполнения выбора переменной стохастического поиска (SSVS). Предположим, что $\beta$ и ${\sigma }^2$являются зависимыми (модель сопряженной смеси). Задайте количество предикторов p и имена коэффициентов регрессии.

p = 3;
PriorMdl = mixconjugateblm(p,'VarNames',["IPI" "E" "WR"]);

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

load Data_NelsonPlosser
X = DataTable{:,PriorMdl.VarNames(2:end)};
y = DataTable{:,'GNPR'};

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

rng(1);
PosteriorMdl = estimate(PriorMdl,X,y);

Method: MCMC sampling with 10000 draws
Number of observations: 62
Number of predictors:   4
 
           |   Mean      Std           CI95        Positive  Distribution  Regime 
----------------------------------------------------------------------------------
 Intercept | -18.8333  10.1851  [-36.965,  0.716]    0.037     Empirical   0.8806 
 IPI       |   4.4554   0.1543   [ 4.165,  4.764]    1.000     Empirical   0.4545 
 E         |   0.0010   0.0004   [ 0.000,  0.002]    0.997     Empirical   0.0925 
 WR        |   2.4686   0.3615   [ 1.766,  3.197]    1.000     Empirical   0.1734 
 Sigma2    |  47.7557   8.6551   [33.858, 66.875]    1.000     Empirical    NaN   
 

PosteriorMdl является empiricalblm объект модели, который хранит черты из апостериорных распределений $$\beta$$ и $${\sigma }^2$$ учитывая данные. estimate отображает сводные данные маргинальных апостериорных распределений в командной строке. Строки сводных данных соответствуют коэффициентам регрессии и нарушения порядка отклонения, а столбцы соответствуют характеристикам апостериорного распределения. Характеристики включают:

Принимая, что переменные с Regime < 0,1 следует удалить из модели, результаты показывают, что можно исключить уровень безработицы из модели.

По умолчанию 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');

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

Рассмотрим регрессионую модель и предшествующее распределение в Select Variables Using Bayesian Lasso Regression.

Создайте байесову регрессию лассо предыдущую модель для 3 предикторов и задайте имена переменных. Задайте значения усадки 0.01, 10, 1e5, и 10 для точки пересечения и коэффициентов IPI, E, и WR.

p = 3;
PriorMdl = bayeslm(p,'ModelType','lasso','VarNames',["IPI" "E" "WR"],...
    'Lambda',[0.01; 10; 1e5; 10]);

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

load Data_NelsonPlosser
X = DataTable{:,PriorMdl.VarNames(2:end)};
y = DataTable{:,"GNPR"};

Оцените условное апостериорное распределение $$\beta$$ учитывая данные и что $${\sigma }^2=10$$и возвращает сводную таблицу оценок для доступа к оценкам.

rng(1); % For reproducibility
[Mdl,SummaryBeta] = estimate(PriorMdl,X,y,'Sigma2',10);

Method: lasso MCMC sampling with 10000 draws
Conditional variable: Sigma2 fixed at  10
Number of observations: 62
Number of predictors:   4
 
           |   Mean     Std          CI95        Positive  Distribution 
------------------------------------------------------------------------
 Intercept | -8.0643  4.1992  [-16.384,  0.018]    0.025     Empirical  
 IPI       |  4.4454  0.0679   [ 4.312,  4.578]    1.000     Empirical  
 E         |  0.0004  0.0002   [ 0.000,  0.001]    0.999     Empirical  
 WR        |  2.9792  0.1672   [ 2.651,  3.305]    1.000     Empirical  
 Sigma2    |   10      0       [10.000, 10.000]    1.000     Empirical  
 

estimate отображает сводные данные условного апостериорного распределения $$\beta$$. Поскольку ${\sigma }^2$ фиксируется на 10 во время оценки, выводы на ней тривиальны.

Отобразите Mdl.

Mdl

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      0.1000    [-0.200,  0.200]     0.500   Scale mixture  
 E         |  0      0.0000    [-0.000,  0.000]     0.500   Scale mixture  
 WR        |  0      0.1000    [-0.200,  0.200]     0.500   Scale mixture  
 Sigma2    | 0.5000  0.5000    [ 0.138,  1.616]     1.000   IG(3.00,    1) 
 

Потому что estimate вычисляет условное апостериорное распределение, оно возвращает входной параметр модели PriorMdl, не условный апостериор, в первой позиции списка выходных аргументов.

Отобразите сводную таблицу оценок.

SummaryBeta

SummaryBeta=5×6 table
                    Mean          Std                  CI95              Positive    Distribution                                   Covariances                              
                 __________    __________    ________________________    ________    _____________    _______________________________________________________________________

    Intercept       -8.0643        4.1992       -16.384       0.01837     0.0254     {'Empirical'}         17.633        0.17621    -0.00053724        0.11705              0
    IPI              4.4454      0.067949         4.312        4.5783          1     {'Empirical'}        0.17621      0.0046171    -1.4103e-06     -0.0068855              0
    E            0.00039896    0.00015673    9.4925e-05    0.00070697     0.9987     {'Empirical'}    -0.00053724    -1.4103e-06     2.4564e-08    -1.8168e-05              0
    WR               2.9792       0.16716        2.6506        3.3046          1     {'Empirical'}        0.11705     -0.0068855    -1.8168e-05       0.027943              0
    Sigma2               10             0            10            10          1     {'Empirical'}              0              0              0              0              0

SummaryBeta содержит условные апостериорные оценки.

Оцените условные апостериорные распределения $${\sigma }^2$$ учитывая, что $$\beta$$ - условное апостериорное среднее значение $$\beta |{\sigma }^2,X,y$$ (хранится в SummaryBeta.Mean(1:(end – 1))). Верните сводную таблицу оценок.

condPostMeanBeta = SummaryBeta.Mean(1:(end - 1));
[~,SummarySigma2] = estimate(PriorMdl,X,y,'Beta',condPostMeanBeta);

Method: lasso MCMC sampling with 10000 draws
Conditional variable: Beta fixed at -8.0643      4.4454  0.00039896      2.9792
Number of observations: 62
Number of predictors:   4
 
           |   Mean     Std          CI95        Positive  Distribution 
------------------------------------------------------------------------
 Intercept | -8.0643   0.0000  [-8.064, -8.064]    0.000     Empirical  
 IPI       |  4.4454   0.0000  [ 4.445,  4.445]    1.000     Empirical  
 E         |  0.0004   0.0000  [ 0.000,  0.000]    1.000     Empirical  
 WR        |  2.9792   0.0000  [ 2.979,  2.979]    1.000     Empirical  
 Sigma2    | 56.8314  10.2921  [39.947, 79.731]    1.000     Empirical  
 

estimate отображает сводные данные оценок условного апостериорного распределения $${\sigma }^2$$ учитывая данные и что $$\beta$$ является condPostMeanBeta. На отображении выводы $$\beta$$ тривиальны.

Рассмотрим регрессионую модель в Выборе Переменных с Использованием Регрессии Бейесова Лассо.

Создайте предыдущую модель для выполнения SSVS. Предположим, что $\beta$ и ${\sigma }^2$являются зависимыми (модель сопряженной смеси). Задайте количество предикторов p и имена коэффициентов регрессии.

p = 3;
PriorMdl = mixconjugateblm(p,'VarNames',["IPI" "E" "WR"]);

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

load Data_NelsonPlosser
X = DataTable{:,PriorMdl.VarNames(2:end)};
y = DataTable{:,'GNPR'};

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

rng(1);
[PosteriorMdl,Summary] = estimate(PriorMdl,X,y,'Display',false);

PosteriorMdl является empiricalblm объект модели, который хранит черты из апостериорных распределений $$\beta$$ и $${\sigma }^2$$ учитывая данные. Summary - таблица со столбцами, соответствующими апостериорным характеристикам, и строками, соответствующими коэффициентам (PosteriorMdl.VarNames) и отклонение нарушения порядка (Sigma2).

Отобразите оцененную ковариационную матрицу параметра (Covariances) и доля раз, в которую алгоритм включает каждый предиктор (Regime).

Covariances = Summary(:,"Covariances")

Covariances=5×1 table
                                              Covariances                              
                 ______________________________________________________________________

    Intercept        103.74         1.0486     -0.0031629         0.6791         7.3916
    IPI              1.0486       0.023815    -1.3637e-05      -0.030387        0.06611
    E            -0.0031629    -1.3637e-05     1.3481e-07    -8.8792e-05    -0.00025044
    WR               0.6791      -0.030387    -8.8792e-05        0.13066       0.089039
    Sigma2           7.3916        0.06611    -0.00025044       0.089039         74.911

Regime = Summary(:,"Regime")

Regime=5×1 table
                 Regime
                 ______

    Intercept    0.8806
    IPI          0.4545
    E            0.0925
    WR           0.1734
    Sigma2          NaN

Regime содержит предельную апостериорную вероятность включения переменной ($\gamma =1$ для переменной). Для примера, апостериорная вероятность того, что E должен быть включен в модель 0,0925.

Принимая, что переменные с Regime < 0,1 следует удалить из модели, результаты показывают, что можно исключить уровень безработицы из модели.