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

Для всех $$\mathrm{t-временных}$$ точек $${}_{\mathrm{\alpha t}}$$ - это ряд независимых гауссовых возмущений со средним значением 0 и дисперсией $${}^{\mathrm{start2}}$$.

Предположим, что коэффициенты регрессии $$\beta ={}_{\mathrm{[}}\beta 0,..{.}_{\mathrm{,}}{}^{}$$β3] ′ и $${}^{\mathrm{}}$$дисперсию возмущений start2 являются случайными переменными, и у вас нет предварительных знаний об их значениях или распределении. То есть вы хотите использовать неинформативный Jeffreys previous: совместное предварительное распределение пропорционально $$\mathrm{}{}^{\mathrm{}}$$1/start2.

Эти допущения и вероятность данных подразумевают аналитически прослеживаемое апостериорное распределение.

Создайте диффузную предыдущую модель для параметров линейной регрессии, которая является типом предыдущей модели по умолчанию. Укажите количество предикторов p.

p = 3;
Mdl = bayeslm(p)

Mdl = 
  diffuseblm with properties:

    NumPredictors: 3
        Intercept: 1
         VarNames: {4x1 cell}

 
           | Mean  Std        CI95        Positive        Distribution       
-----------------------------------------------------------------------------
 Intercept |  0    Inf  [   NaN,    NaN]    0.500   Proportional to one      
 Beta(1)   |  0    Inf  [   NaN,    NaN]    0.500   Proportional to one      
 Beta(2)   |  0    Inf  [   NaN,    NaN]    0.500   Proportional to one      
 Beta(3)   |  0    Inf  [   NaN,    NaN]    0.500   Proportional to one      
 Sigma2    |  Inf  Inf  [   NaN,    NaN]    1.000   Proportional to 1/Sigma2 
 

Mdl является diffuseblm Объект байесовской модели линейной регрессии, представляющий предварительное распределение коэффициентов регрессии и дисперсии возмущений. В окне команд bayeslm отображает сводку предыдущих распределений. Поскольку предшествующий уровень является неинформативным, а данные еще не включены, резюме является тривиальным.

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

Mdl.VarNames = ["IPI" "E" "WR"]

Mdl = 
  diffuseblm with properties:

    NumPredictors: 3
        Intercept: 1
         VarNames: {4x1 cell}

 
           | Mean  Std        CI95        Positive        Distribution       
-----------------------------------------------------------------------------
 Intercept |  0    Inf  [   NaN,    NaN]    0.500   Proportional to one      
 IPI       |  0    Inf  [   NaN,    NaN]    0.500   Proportional to one      
 E         |  0    Inf  [   NaN,    NaN]    0.500   Proportional to one      
 WR        |  0    Inf  [   NaN,    NaN]    0.500   Proportional to one      
 Sigma2    |  Inf  Inf  [   NaN,    NaN]    1.000   Proportional to 1/Sigma2 
 

Рассмотрим модель линейной регрессии в разделе Создание диффузионной предыдущей модели.

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

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

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

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

Оцените краевые апостериорные распределения $$\beta$$ и $${}^{\mathrm{start2}}$$.

PosteriorMdl = estimate(PriorMdl,X,y);

Method: Analytic posterior distributions
Number of observations: 62
Number of predictors:   4
 
           |   Mean      Std           CI95        Positive       Distribution      
------------------------------------------------------------------------------------
 Intercept | -24.2536   9.5314  [-43.001, -5.506]    0.006   t (-24.25, 9.37^2, 58) 
 IPI       |   4.3913   0.1535   [ 4.089,  4.693]    1.000   t (4.39, 0.15^2, 58)   
 E         |   0.0011   0.0004   [ 0.000,  0.002]    0.999   t (0.00, 0.00^2, 58)   
 WR        |   2.4682   0.3787   [ 1.723,  3.213]    1.000   t (2.47, 0.37^2, 58)   
 Sigma2    |  51.9790  10.0034   [35.965, 74.937]    1.000   IG(29.00, 0.00069)     
 

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

Свойства доступа к заднему распределению с помощью точечной нотации. Например, отображение краевого заднего средства путем доступа к Mu собственность.

PosteriorMdl.Mu

ans = 4×1

  -24.2536
    4.3913
    0.0011
    2.4682

Рассмотрим модель линейной регрессии в разделе Создание диффузионной предыдущей модели.

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

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

PriorMdl = 
  diffuseblm with properties:

    NumPredictors: 3
        Intercept: 1
         VarNames: {4x1 cell}

 
           | Mean  Std        CI95        Positive        Distribution       
-----------------------------------------------------------------------------
 Intercept |  0    Inf  [   NaN,    NaN]    0.500   Proportional to one      
 IPI       |  0    Inf  [   NaN,    NaN]    0.500   Proportional to one      
 E         |  0    Inf  [   NaN,    NaN]    0.500   Proportional to one      
 WR        |  0    Inf  [   NaN,    NaN]    0.500   Proportional to one      
 Sigma2    |  Inf  Inf  [   NaN,    NaN]    1.000   Proportional to 1/Sigma2 
 

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

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

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

[Mdl,Summary] = estimate(PriorMdl,X,y,'Sigma2',2);

Method: Analytic posterior distributions
Conditional variable: Sigma2 fixed at   2
Number of observations: 62
Number of predictors:   4
 
           |   Mean      Std          CI95         Positive     Distribution    
--------------------------------------------------------------------------------
 Intercept | -24.2536  1.8696  [-27.918, -20.589]    0.000   N (-24.25, 1.87^2) 
 IPI       |   4.3913  0.0301   [ 4.332,  4.450]     1.000   N (4.39, 0.03^2)   
 E         |   0.0011  0.0001   [ 0.001,  0.001]     1.000   N (0.00, 0.00^2)   
 WR        |   2.4682  0.0743   [ 2.323,  2.614]     1.000   N (2.47, 0.07^2)   
 Sigma2    |    2       0       [ 2.000,  2.000]     1.000   Fixed value        
 

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

Извлеките средний вектор и ковариационную матрицу условного задника $\beta$ из сводной таблицы оценки.

condPostMeanBeta = Summary.Mean(1:(end - 1))

condPostMeanBeta = 4×1

  -24.2536
    4.3913
    0.0011
    2.4682

CondPostCovBeta = Summary.Covariances(1:(end - 1),1:(end - 1))

CondPostCovBeta = 4×4

    3.4956    0.0350   -0.0001    0.0241
    0.0350    0.0009   -0.0000   -0.0013
   -0.0001   -0.0000    0.0000   -0.0000
    0.0241   -0.0013   -0.0000    0.0055

Показ Mdl.

Mdl

Mdl = 
  diffuseblm with properties:

    NumPredictors: 3
        Intercept: 1
         VarNames: {4x1 cell}

 
           | Mean  Std        CI95        Positive        Distribution       
-----------------------------------------------------------------------------
 Intercept |  0    Inf  [   NaN,    NaN]    0.500   Proportional to one      
 IPI       |  0    Inf  [   NaN,    NaN]    0.500   Proportional to one      
 E         |  0    Inf  [   NaN,    NaN]    0.500   Proportional to one      
 WR        |  0    Inf  [   NaN,    NaN]    0.500   Proportional to one      
 Sigma2    |  Inf  Inf  [   NaN,    NaN]    1.000   Proportional to 1/Sigma2 
 

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

Рассмотрим модель линейной регрессии в разделе Оценка предельных задних распределений.

Создайте предыдущую модель для коэффициентов регрессии и дисперсии возмущений, а затем оцените предельные апостериорные распределения.

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

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

PosteriorMdl = estimate(PriorMdl,X,y);

Method: Analytic posterior distributions
Number of observations: 62
Number of predictors:   4
 
           |   Mean      Std           CI95        Positive       Distribution      
------------------------------------------------------------------------------------
 Intercept | -24.2536   9.5314  [-43.001, -5.506]    0.006   t (-24.25, 9.37^2, 58) 
 IPI       |   4.3913   0.1535   [ 4.089,  4.693]    1.000   t (4.39, 0.15^2, 58)   
 E         |   0.0011   0.0004   [ 0.000,  0.002]    0.999   t (0.00, 0.00^2, 58)   
 WR        |   2.4682   0.3787   [ 1.723,  3.213]    1.000   t (2.47, 0.37^2, 58)   
 Sigma2    |  51.9790  10.0034   [35.965, 74.937]    1.000   IG(29.00, 0.00069)     
 

Извлеките заднее среднее $$\beta$$ из задней модели и заднюю ковариацию β $$$$из сводки оценки, возвращенной summarize.

estBeta = PosteriorMdl.Mu;
Summary = summarize(PosteriorMdl);
estBetaCov = Summary.Covariances{1:(end - 1),1:(end - 1)};

Предположим, что если коэффициент реальной заработной платы ниже 2,5, то вводится политика. Хотя заднее распределение WR известен, и можно вычислить вероятности напрямую, можно оценить вероятность с помощью моделирования Монте-Карло.

Нарисовать 1e6 образцы из краевого заднего распределения $$\beta$$.

NumDraws = 1e6;
rng(1);
BetaSim = simulate(PosteriorMdl,'NumDraws',NumDraws);

BetaSim является 4-by- 1e6 матрица, содержащая черточки. Строки соответствуют коэффициенту регрессии, а столбцы - последовательным розыгрышам.

Выделить розыгрыши, соответствующие коэффициенту реальной заработной платы, а затем определить, какие розыгрыши меньше 2,5.

isWR = PosteriorMdl.VarNames == "WR";
wrSim = BetaSim(isWR,:);
isWRLT2p5 = wrSim < 2.5;

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

probWRLT2p5 = mean(isWRLT2p5)

probWRLT2p5 = 0.5341

Апостериорная вероятность того, что коэффициент реальной заработной платы меньше 2,5, составляет около 0.53.

Краевое заднее распределение коэффициента WR является $${}_{\mathrm{}\mathrm{t58}}$$, но центрирован в 2,47 и масштабирован на 0,37. Непосредственно вычислить апостериорную вероятность того, что коэффициент WR менее 2,5.

center = estBeta(isWR);
stdBeta = sqrt(diag(estBetaCov));
scale = stdBeta(isWR);
t = (2.5 - center)/scale;
dof = 68;
directProb = tcdf(t,dof)

directProb = 0.5333

Задние вероятности почти идентичны.

Рассмотрим модель линейной регрессии в разделе Оценка предельных задних распределений.

Создайте предыдущую модель для коэффициентов регрессии и дисперсии возмущений, а затем оцените предельные апостериорные распределения. Удерживайте последние 10 периодов данных из оценки, чтобы использовать их для прогнозирования реального ВНП. Выключите отображение оценки.

p = 3;
PriorMdl = bayeslm(p,'ModelType','diffuse','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

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');
ylabel('rGNP');
xlabel('Year');
hold off

yF является вектором 10 на 1 будущих значений реального GNP, соответствующего будущим данным предиктора.

Оценка прогнозируемой среднеквадратичной ошибки (RMSE).

frmse = sqrt(mean((yF - yFT).^2))

frmse = 25.5489

RMSE прогноза является относительным показателем точности прогноза. В частности, можно оценить несколько моделей с использованием различных допущений. Модель с самым низким прогнозом RMSE является лучшей из сравниваемых моделей.

Copyright 2018 The MathWorks, Inc.