Tune Slice Sampler для задней оценки

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

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

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

$\texttt{GNPR}_t=\beta_0+\beta_1\texttt{IPI}_t+\beta_2\texttt{E}_t+\beta_3\texttt{WR}_t+\varepsilon_t.$

Для всех $t$ моментов времени $\varepsilon_t$ представляет собой ряд независимых гауссовых возмущений со средним значением 0 и дисперсией. $\sigma^2$ Предположим, что предыдущие распределения:

$\beta_j\vert\sigma^2$ является 4-мерным t-распределением с 50 степенями свободы для каждого компонента и единичной матрицей для корреляционной матрицы. Кроме того, распределение центрировано ${\left[ {\begin{array}{*{20}{c}} { - 25}&4&0&3 \end{array}} \right]^\prime }$ , и каждый компонент масштабируется соответствующими элементами вектора. ${\left[ {\begin{array}{*{20}{c}} {1}&1&1&1 \end{array}} \right]^\prime }$

$\sigma^2 \sim IG(3,1)$ .

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

Объявите функцию MATLAB ®, которая:

Принимает значения $\beta$ $\sigma^2$ вектора столбца и значения гиперпараметров.
Возвращает значение предыдущего распределения соединения,, $\pi\left(\beta,\sigma^2\right)$ учитывая значения $\beta$ и. $\sigma^2$

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

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

dof = 50;
C = eye(4);
ct = zeros(4,1);
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 Объект байесовской модели линейной регрессии, представляющий предварительное распределение коэффициентов регрессии и дисперсии возмущений.

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

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

Оцените краевые задние распределения $\beta$ и. $\sigma^2$ Для определения соответствующего периода записи укажите период записи. Укажите ширину образца среза, близкую к заднему стандартному отклонению параметров в предположении о диффузной предыдущей модели. Нарисовать 15000 образцы.

width = [20,0.5,0.01,1,20];
numDraws = 15e3;
rng(1) % For reproducibility
PosteriorMdl1 = estimate(PriorMdl,X,y,'BurnIn',0,'Width',width,...
    'NumDraws',numDraws);

Warning: High autocorrelation of MCMC draws:IPI,E,WR,Sigma2. 

Method: MCMC sampling with 15000 draws
Number of observations: 62
Number of predictors:   4
 
           |   Mean     Std         CI95        Positive  Distribution 
-----------------------------------------------------------------------
 Intercept | -0.6669  1.5325  [-3.202,  1.621]    0.309     Empirical  
 IPI       |  4.6196  0.1021  [ 4.437,  4.839]    1.000     Empirical  
 E         |  0.0004  0.0002  [ 0.000,  0.001]    0.995     Empirical  
 WR        |  2.6135  0.3025  [ 1.982,  3.152]    1.000     Empirical  
 Sigma2    | 48.9990  9.1926  [33.147, 67.993]    1.000     Empirical

PosteriorMdl1 является empiricalblm объект модели. BetaDraws и Sigma2Draws свойства хранят образцы из пробоотборника среза.

Нарисуйте контурные и автокорреляционные графики для цепей Маркова.

figure;
for j = 1:4
    subplot(2,3,j)
    plot(PosteriorMdl1.BetaDraws(j,:));
    axis tight
    xlabel('Sample index')
    title(sprintf(['Trace Plot ' char(8212) ' %s'],PriorMdl.VarNames{j}));
end
subplot(2,3,5)
plot(PosteriorMdl1.Sigma2Draws);
axis tight
xlabel('Sample index')
title(['Trace Plot ' char(8212) ' Sigma2']);

figure;
for j = 1:4
    subplot(2,3,j)
    autocorr(PosteriorMdl1.BetaDraws(j,:));
    axis tight
    title(sprintf(['Correlogram ' char(8212) ' %s'],PriorMdl.VarNames{j}));
end
subplot(2,3,5)
autocorr(PosteriorMdl1.Sigma2Draws);
axis tight
title(['Correlogram ' char(8212) ' Sigma2']);

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

Снова оцените апостериорное распределение. Для уменьшения переходных эффектов используйте период горения по умолчанию (5000, что является значением по умолчанию). Чтобы уменьшить влияние последовательной корреляции, тонкий в 20 раз и указать 15e4/thin эффективное количество розыгрышей. Кроме того, численно стабилизировать оценку $\sigma^2$ путем указания репарамеризации. $\log(\sigma^2)$ Графики трассировки и автокорреляции.

thin = 20;
numDraws = 15e4/thin;
PosteriorMdl2 = estimate(PriorMdl,X,y,'Thin',thin,'Width',width,...
    'Reparameterize',true,'NumDraws',numDraws);

figure;
for j = 1:4
    subplot(2,3,j)
    plot(PosteriorMdl2.BetaDraws(j,:));
    axis tight
    xlabel('Sample index')
    title(sprintf(['Trace Plot ' char(8212) '%s'],PriorMdl.VarNames{j}));
end
subplot(2,3,5)
plot(PosteriorMdl2.Sigma2Draws);
axis tight
xlabel('Sample index')
title(['Trace Plot ' char(8212) ' Sigma2']);

figure;
for j = 1:4
    subplot(2,3,j)
    autocorr(PosteriorMdl2.BetaDraws(j,:));
    title(sprintf(['Correlogram ' char(8212) ' %s'],PriorMdl.VarNames{j}));
end
subplot(2,3,5)
autocorr(PosteriorMdl2.Sigma2Draws);
title(['Correlogram  ' char(8212) ' Sigma2']);

Method: MCMC sampling with 7500 draws
Number of observations: 62
Number of predictors:   4
 
           |   Mean     Std         CI95        Positive  Distribution 
-----------------------------------------------------------------------
 Intercept | -0.4777  1.2225  [-2.911,  1.955]    0.338     Empirical  
 IPI       |  4.6506  0.1128  [ 4.434,  4.884]    1.000     Empirical  
 E         |  0.0005  0.0002  [ 0.000,  0.001]    0.991     Empirical  
 WR        |  2.5243  0.3468  [ 1.815,  3.208]    1.000     Empirical  
 Sigma2    | 51.3444  9.3170  [36.319, 72.469]    1.000     Empirical

Сюжеты следов позволяют предположить, что марковские цепи:

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

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

Документация

Tune Slice Sampler для задней оценки

См. также

Функции

Объекты

Связанные темы

Документация по инструментарию Econometrics

Поддержка