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

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

Для всех временных точек является рядом независимых Гауссовых нарушений порядка со средним значением 0 и отклонением. Примите, что предыдущие распределения:

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

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

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'};

Оцените маргинальные апостериорные распределения и. Чтобы определить подходящий период горения, не задайте период горения. Задайте ширину для дискретизатора среза, которая близка к апостериорному стандартному отклонению параметров, предполагающих диффузную предшествующую модель. Нарисуйте 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 объект модели. The 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 эффективное количество розыгрышей. Кроме того, численно стабилизируйте оценку путем определения reparameterize с. Постройте графики трассировки и автокорреляции.

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  
 

Следовые графики предполагают, что марковские цепи:

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