Считайте несколько моделью линейной регрессии, которая предсказывает США действительный валовой национальный продукт (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$$ и $${\sigma }^2$$. Задайте ширину для сэмплера среза, который является близко к следующему стандартному отклонению параметров, принимающих рассеянную предшествующую модель. Уменьшайте последовательную корреляцию путем определения утончающегося фактора 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$$ и $${\sigma }^2$$ учитывая данные. 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$$ учитывая данные и $${\sigma }^2=2$$. Задайте ширину для сэмплера среза, который является близко к следующему стандартному отклонению параметров, принимающих рассеянную предшествующую модель. Уменьшайте последовательную корреляцию путем определения утончающегося фактора 10, и, сократите эффективное количество по умолчанию ничьих фактором 10.

width = [20,0.5,0.01,1];
thin = 10;
numDraws = 1e5/thin;
rng(1) % For reproducibility
[Mdl,condPostMeanBeta,CondPostCovBeta] = 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  
 

Warning: Current syntax supports 6 output arguments, and will be removed in a future release. For supported output arguments, see <a href="matlab:helpview(fullfile(docroot,'econ','econ.map'),'blmestimate')">estimate</a>.

estimate возвращается 4 1 вектор средних значений и ковариационная матрица 4 на 4 условного апостериорного распределения $$\beta$$ в condPostMeanBeta и CondPostCovBeta, соответственно. Кроме того, estimate отображает сводные данные условного апостериорного распределения $$\beta$$.

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

Отобразите 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, вместо этого задающий тот же seed случайных чисел, прежде чем вызов моделирует.

Считайте модель линейной регрессии в Оценке Крайними Апостериорными распределениями.

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

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 периодов данных из оценки, таким образом, можно использовать их, чтобы предсказать действительный GNP. Выключите отображение оценки.

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