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

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

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  
 

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

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