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

Для всех моментов времени представляет собой ряд независимых гауссовых возмущений со средним значением 0 и дисперсией. Предположим, что эти предыдущие распределения:

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

function logPDF = priorMVTIG(params,ct,st,dof,C,a,b)
%priorMVTIG Log density of multivariate t times inverse gamma 
%   priorMVTIG passes params(1:end-1) to the multivariate t density
%   function with dof degrees of freedom for each component and positive
%   definite correlation matrix C. priorMVTIG returns the log of the product of
%   the two evaluated densities.
%   
%   params: Parameter values at which the densities are evaluated, an
%           m-by-1 numeric vector.
%
%       ct: Multivariate t distribution component centers, an (m-1)-by-1
%           numeric vector.  Elements correspond to the first m-1 elements
%           of params.
%
%       st: Multivariate t distribution component scales, an (m-1)-by-1
%           numeric (m-1)-by-1 numeric vector.  Elements correspond to the
%           first m-1 elements of params.
%
%      dof: Degrees of freedom for the multivariate t distribution, a
%           numeric scalar or (m-1)-by-1 numeric vector. priorMVTIG expands
%           scalars such that dof = dof*ones(m-1,1). Elements of dof
%           correspond to the elements of params(1:end-1).
%
%        C: Correlation matrix for the multivariate t distribution, an
%           (m-1)-by-(m-1) symmetric, positive definite matrix. Rows and
%           columns correspond to the elements of params(1:end-1).
% 
%        a: Inverse gamma shape parameter, a positive numeric scalar.
%
%        b: Inverse gamma scale parameter, a positive scalar.
%
beta = params(1:(end-1));
sigma2 = params(end);

tVal = (beta - ct)./st;
mvtDensity = mvtpdf(tVal,C,dof);
igDensity = sigma2^(-a-1)*exp(-1/(sigma2*b))/(gamma(a)*b^a);

logPDF = log(mvtDensity*igDensity);
end

dof = 50;
C = eye(4);
ct = [-25; 4; 0; 3];
st = ones(4,1);
a = 3;
b = 1;
logPDF = @(params)priorMVTIG(params,ct,st,dof,C,a,b);

p = 3;
PriorMdl = bayeslm(p,'ModelType','custom','LogPDF',logPDF,...
    'VarNames',["IPI" "E" "WR"])

PriorMdl = 

  customblm with properties:

    NumPredictors: 3
        Intercept: 1
         VarNames: {4x1 cell}
           LogPDF: @(params)priorMVTIG(params,ct,st,dof,C,a,b)

The priors are defined by the function:
    @(params)priorMVTIG(params,ct,st,dof,C,a,b)

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

Создать анонимную функцию, которая работает как priorMVTIG, но принимает только значения параметров и сохраняет значения гиперпараметров фиксированными по их значениям.

dof = 50;
C = eye(4);
ct = [-25; 4; 0; 3];
st = ones(4,1);
a = 3;
b = 1;
logPDF = @(params)priorMVTIG(params,ct,st,dof,C,a,b);

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

p = 3;
PriorMdl = bayeslm(p,'ModelType','custom','LogPDF',logPDF,...
    'VarNames',["IPI" "E" "WR"])

PriorMdl = 
  customblm with properties:

    NumPredictors: 3
        Intercept: 1
         VarNames: {4x1 cell}
           LogPDF: @(params)priorMVTIG(params,ct,st,dof,C,a,b)

The priors are defined by the function:
    @(params)priorMVTIG(params,ct,st,dof,C,a,b)

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

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

Оцените краевые апостериорные распределения $$\beta$$ и $${}^{\mathrm{start2}}$$. Укажите ширину образца среза, близкую к заднему стандартному отклонению параметров в предположении о диффузной предыдущей модели. Уменьшите последовательную корреляцию, указав коэффициент прореживания 10, и уменьшите фактическое количество розыгрышей по умолчанию в 10 раз.

width = [20,0.5,0.01,1,20];
thin = 10;
numDraws = 1e5/thin;
rng(1) % For reproducibility
PosteriorMdl = estimate(PriorMdl,X,y,'Width',width,'Thin',thin,...
    'NumDraws',numDraws);

Method: MCMC sampling with 10000 draws
Number of observations: 62
Number of predictors:   4
 
           |   Mean      Std          CI95         Positive  Distribution 
--------------------------------------------------------------------------
 Intercept | -25.0069  0.9919  [-26.990, -23.065]    0.000     Empirical  
 IPI       |   4.3544  0.1083   [ 4.143,  4.562]     1.000     Empirical  
 E         |   0.0011  0.0002   [ 0.001,  0.001]     1.000     Empirical  
 WR        |   2.5613  0.3293   [ 1.939,  3.222]     1.000     Empirical  
 Sigma2    |  47.0593  8.7570   [32.690, 67.115]     1.000     Empirical  
 

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

estimate извлекает задние характеристики из рисований из задних распределений, которые MATLAB ® сохраняет в качестве матриц в свойствахBetaDraws и Sigma2Draws.

Для мониторинга смешения и сходимости образца MCMC создайте графики трассировки. В BetaDraws свойства, рисунки соответствуют столбцам и параметры строкам.

figure;
for j = 1:4
    subplot(2,2,j)
    plot(PosteriorMdl.BetaDraws(j,:))
    title(sprintf('Trace Plot of %s',PosteriorMdl.VarNames{j}));
end

figure;
plot(PosteriorMdl.Sigma2Draws)
title('Trace Plot of Sigma2');

Графики трассировки указывают на адекватное смешение и сходимость, и нет временных эффектов для удаления.

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

Создать анонимную функцию, которая работает как priorMVTIG, но принимает только значения параметров и сохраняет фиксированные значения гиперпараметров.

dof = 50;
C = eye(4);
ct = [-25; 4; 0; 3];
st = ones(4,1);
a = 3;
b = 1;
logPDF = @(params)priorMVTIG(params,ct,st,dof,C,a,b);

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

p = 3;
PriorMdl = bayeslm(p,'ModelType','custom','LogPDF',logPDF,...
    'VarNames',["IPI" "E" "WR"])

PriorMdl = 
  customblm with properties:

    NumPredictors: 3
        Intercept: 1
         VarNames: {4x1 cell}
           LogPDF: @(params)priorMVTIG(params,ct,st,dof,C,a,b)

The priors are defined by the function:
    @(params)priorMVTIG(params,ct,st,dof,C,a,b)

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

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

Оцените условное апостериорное распределение $$\beta$$, учитывая данные, и $${}^{\mathrm{start2}}\mathrm{=}$$2, и верните сводную таблицу оценки для доступа к оценкам. Укажите ширину образца среза, близкую к заднему стандартному отклонению параметров в предположении о диффузной предыдущей модели. Уменьшите последовательную корреляцию, указав коэффициент прореживания 10, и уменьшите фактическое количество розыгрышей по умолчанию в 10 раз.

width = [20,0.5,0.01,1];
thin = 10;
numDraws = 1e5/thin;
rng(1) % For reproducibility
[Mdl,Summary] = estimate(PriorMdl,X,y,'Sigma2',2,...
    'Width',width,'Thin',thin,'NumDraws',numDraws);

Method: MCMC sampling with 10000 draws
Conditional variable: Sigma2 fixed at   2
Number of observations: 62
Number of predictors:   4
 
           |   Mean      Std          CI95         Positive  Distribution 
--------------------------------------------------------------------------
 Intercept | -24.7820  0.8767  [-26.483, -23.054]    0.000     Empirical  
 IPI       |   4.3825  0.0254   [ 4.332,  4.431]     1.000     Empirical  
 E         |   0.0011  0.0000   [ 0.001,  0.001]     1.000     Empirical  
 WR        |   2.4752  0.0724   [ 2.337,  2.618]     1.000     Empirical  
 Sigma2    |    2       0       [ 2.000,  2.000]     1.000     Empirical  
 

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

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

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

condPostMeanBeta = 4×1

  -24.7820
    4.3825
    0.0011
    2.4752

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

CondPostCovBeta = 4×4

    0.7686    0.0084   -0.0000    0.0019
    0.0084    0.0006    0.0000   -0.0015
   -0.0000    0.0000    0.0000   -0.0000
    0.0019   -0.0015   -0.0000    0.0052

Показ Mdl.

Mdl

Mdl = 
  customblm with properties:

    NumPredictors: 3
        Intercept: 1
         VarNames: {4x1 cell}
           LogPDF: @(params)priorMVTIG(params,ct,st,dof,C,a,b)

The priors are defined by the function:
    @(params)priorMVTIG(params,ct,st,dof,C,a,b)

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

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

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

dof = 50;
C = eye(4);
ct = [-25; 4; 0; 3];
st = ones(4,1);
a = 3;
b = 1;
logPDF = @(params)priorMVTIG(params,ct,st,dof,C,a,b);

p = 3;
PriorMdl = bayeslm(p,'ModelType','custom','LogPDF',logPDF,...
    'VarNames',["IPI" "E" "WR"]);

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

width = [20,0.5,0.01,1,20];
thin = 10;
numDraws = 1e5/thin;

rng(1) % For reproducibility
PosteriorMdl = estimate(PriorMdl,X,y,'Width',width,'Thin',thin,...
    'NumDraws',numDraws,'Display',false);

Оцените сводную статистику заднего распределения для $$\beta$$, используя розыгрыши из заднего распределения, сохраненного в задней модели.

estBeta = mean(PosteriorMdl.BetaDraws,2);
EstBetaCov = cov(PosteriorMdl.BetaDraws');

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

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

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

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

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

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

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

probWRLT2p5 = mean(isWRLT2p5)

probWRLT2p5 = 0.4430

Задняя вероятность того, что коэффициент WR меньше 2,5 составляет около 0.4430.

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

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

load Data_NelsonPlosser
VarNames = {'IPI'; 'E'; 'WR'};
fhs = 10; % Forecast horizon size
X = DataTable{1:(end - fhs),VarNames};
y = DataTable{1:(end - fhs),'GNPR'};
XF = DataTable{(end - fhs + 1):end,VarNames}; % Future predictor data
yFT = DataTable{(end - fhs + 1):end,'GNPR'};  % True future responses

dof = 50;
C = eye(4);
ct = [-25; 4; 0; 3];
st = ones(4,1);
a = 3;
b = 1;
logPDF = @(params)priorMVTIG(params,ct,st,dof,C,a,b);

p = 3;
PriorMdl = bayeslm(p,'ModelType','custom','LogPDF',logPDF,...
    'VarNames',VarNames);

width = [20,0.5,0.01,1,20];
thin = 10;
numDraws = 1e5/thin;
rng(1) % For reproducibility
PosteriorMdl = estimate(PriorMdl,X,y,'Width',width,'Thin',thin,...
    'NumDraws',numDraws,'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 = 12.8148

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