Настройте сэмплер среза для следующей оценки

Открыть скрипт

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

При оценке следующего, состоявшего из вероятности данных и полусопряженных или пользовательских предшествующих моделей, 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

Графики трассировки предполагают что Цепи Маркова:

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

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

  • Не показывайте переходное поведение.

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

Смотрите также

Функции

Объекты

Похожие темы

Документация Econometrics Toolbox
Поддержка
Для просмотра документации необходимо авторизоваться на сайте