Задайте градиент для сэмплера HMC

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

Этот пример показывает, как настроить Байесовую модель линейной регрессии для эффективной следующей выборки с помощью сэмплера Гамильтонова Монте-Карло (HMC). Предшествующее распределение коэффициентов является многомерным t, и отклонение воздействия имеет обратную предшествующую гамму.

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

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

  • 4-D t распределение с 30 степенями свободы для каждого компонента, корреляционная матрица C, местоположение ct и шкала st.

  • , с формой и шкалой.

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

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

  • Принимает значения и вместе в вектор-столбце и принимает значения гиперпараметров.

  • Возвращает журнал объединенного предшествующего распределения, учитывая значения и, и возвращает градиент журнала предшествующая плотность относительно и.

function [logPDF,gradient] = priorMVTIGHMC(params,ct,st,dof,C,a,b)
%priorMVTIGHMC Log density of multivariate t times inverse gamma and
%gradient
%   priorMVTIGHMC 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, and passes params(end) to the inverse
%   gamma distribution with shape a and scale b.  priorMVTIG returns the
%   log of the product of the two evaluated densities and its gradient.
%   After applying the log, terms involving constants only are dropped from
%   logPDF.
%   
%   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 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));
numCoeffs = numel(beta);
sigma2 = params(end);

tVal = (beta - ct)./st;
expo = -(dof + numCoeffs)/2;
logmvtDensity = expo*log(1 + (1/dof)*((tVal'/C)*tVal));
logigDensity = (-a - 1)*log(sigma2) - 1/(sigma2*b);
logPDF = logmvtDensity + logigDensity;

grad_beta = (expo/(1 + (1/dof)*(tVal'/C*tVal)))*((2/dof)*((tVal'/C)'./st));
grad_sigma = (-a - 1)/sigma2 + 1/(sigma2^2*b);
gradient = [grad_beta; grad_sigma];
end

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

rng(1); % For reproducibility
dof = 30;
V = rand(4,4);
Sigma = V'*V;
st = sqrt(diag(Sigma));
C = Sigma./(st*st');
ct = -10*rand(4,1);
a = 10*rand;
b = 10*rand;
logPDF = @(params)priorMVTIGHMC(params,ct,st,dof,C,a,b);

Создайте пользовательскую объединенную предшествующую модель для параметров линейной регрессии. Задайте количество предикторов, p. Кроме того, задайте указатель на функцию для priorMVTIGHMC и имен переменных.

p = 3;
PriorMdl = bayeslm(p,'ModelType','custom','LogPDF',logPDF,...
    'VarNames',["IPI" "E" "WR"]);

PriorMdl является customblm Байесов объект модели линейной регрессии, представляющий предшествующее распределение отклонения воздействия и коэффициентов регрессии.

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

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

Оцените крайние апостериорные распределения и использование сэмплера HMC. Чтобы контролировать прогресс сэмплера, задайте уровень многословия 2. Для числовой устойчивости задайте репараметризацию отклонения воздействия.

options = sampleroptions('Sampler','hmc','VerbosityLevel',2)
PosteriorMdl = estimate(PriorMdl,X,y,'Options',options,'Reparameterize',true);

options = 

  struct with fields:

                        Sampler: 'HMC'
           StepSizeTuningMethod: 'dual-averaging'
         MassVectorTuningMethod: 'iterative-sampling'
    NumStepSizeTuningIterations: 100
          TargetAcceptanceRatio: 0.6500
                  NumStepsLimit: 2000
                 VerbosityLevel: 2
                       NumPrint: 100

o Tuning mass vector using method: iterative-sampling 

o Initial step size for dual-averaging = 0.0625

|==================================================================================|
|   ITER   |    LOG PDF    |  STEP SIZE  |  NUM STEPS  |  ACC RATIO  |  DIVERGENT  |
|==================================================================================|
|      100 |  4.321661e+01 |   8.954e-02 |          56 |   6.500e-01 |           2 |
Finished mass vector tuning iteration 1 of 5. 

o Initial step size for dual-averaging = 0.5

|==================================================================================|
|   ITER   |    LOG PDF    |  STEP SIZE  |  NUM STEPS  |  ACC RATIO  |  DIVERGENT  |
|==================================================================================|
|      100 |  4.137069e+01 |   5.039e-01 |          10 |   6.000e-01 |           0 |

|==================================================================================|
|   ITER   |    LOG PDF    |  STEP SIZE  |  NUM STEPS  |  ACC RATIO  |  DIVERGENT  |
|==================================================================================|
|      100 |  4.074939e+01 |   3.641e-02 |           6 |   7.750e-01 |           0 |
Finished mass vector tuning iteration 2 of 5. 

o Initial step size for dual-averaging = 0.5

|==================================================================================|
|   ITER   |    LOG PDF    |  STEP SIZE  |  NUM STEPS  |  ACC RATIO  |  DIVERGENT  |
|==================================================================================|
|      100 |  4.115031e+01 |   2.924e-01 |          17 |   6.300e-01 |           0 |

|==================================================================================|
|   ITER   |    LOG PDF    |  STEP SIZE  |  NUM STEPS  |  ACC RATIO  |  DIVERGENT  |
|==================================================================================|
|      100 |  3.770035e+01 |   2.632e-01 |           4 |   7.950e-01 |           0 |
|      200 |  4.151461e+01 |   1.080e-01 |          15 |   8.467e-01 |           0 |
Finished mass vector tuning iteration 3 of 5. 

o Initial step size for dual-averaging = 0.25

|==================================================================================|
|   ITER   |    LOG PDF    |  STEP SIZE  |  NUM STEPS  |  ACC RATIO  |  DIVERGENT  |
|==================================================================================|
|      100 |  4.229619e+01 |   3.252e-01 |          15 |   6.200e-01 |           0 |

|==================================================================================|
|   ITER   |    LOG PDF    |  STEP SIZE  |  NUM STEPS  |  ACC RATIO  |  DIVERGENT  |
|==================================================================================|
|      100 |  4.349496e+01 |   1.639e-02 |           2 |   7.750e-01 |           0 |
|      200 |  3.973301e+01 |   2.927e-01 |          15 |   8.400e-01 |           0 |
Finished mass vector tuning iteration 4 of 5. 

o Initial step size for dual-averaging = 0.25

|==================================================================================|
|   ITER   |    LOG PDF    |  STEP SIZE  |  NUM STEPS  |  ACC RATIO  |  DIVERGENT  |
|==================================================================================|
|      100 |  4.341187e+01 |   3.193e-01 |          16 |   6.400e-01 |           0 |

|==================================================================================|
|   ITER   |    LOG PDF    |  STEP SIZE  |  NUM STEPS  |  ACC RATIO  |  DIVERGENT  |
|==================================================================================|
|      100 |  4.026322e+01 |   9.213e-02 |           4 |   7.950e-01 |           0 |
|      200 |  4.269868e+01 |   2.873e-01 |           5 |   8.567e-01 |           0 |
Finished mass vector tuning iteration 5 of 5. 

o Tuning step size using method: dual-averaging. Target acceptance ratio = 0.65

o Initial step size for dual-averaging = 0.5

|==================================================================================|
|   ITER   |    LOG PDF    |  STEP SIZE  |  NUM STEPS  |  ACC RATIO  |  DIVERGENT  |
|==================================================================================|
|      100 |  4.252154e+01 |   3.176e-01 |          16 |   6.200e-01 |           0 |

Method: MCMC sampling with 10000 draws
Number of observations: 61
Number of predictors:   4
 
           |   Mean     Std         CI95        Positive  Distribution 
-----------------------------------------------------------------------
 Intercept |  0.0077  0.0087  [-0.009,  0.024]    0.812     Empirical  
 IPI       | -0.1637  0.1457  [-0.477,  0.099]    0.122     Empirical  
 E         |  2.3703  0.5434  [ 1.399,  3.528]    1.000     Empirical  
 WR        | -0.2309  0.3010  [-0.828,  0.347]    0.219     Empirical  
 Sigma2    |  0.0036  0.0007  [ 0.002,  0.005]    1.000     Empirical

PosteriorMdl является объектом модели empiricalblm, хранящим ничьи от апостериорных распределений.

Просмотрите графики трассировки и автокорреляционную функцию (ACF) графики ничьих от последующего поколения.

for j = 1:4
    figure;
    subplot(2,1,1)
    plot(PosteriorMdl.BetaDraws(j,:));
    title(sprintf('Trace plot - Beta %d',j));
    xlabel('MCMC draw')
    ylabel('Simulation index')
    subplot(2,1,2)
    autocorr(PosteriorMdl.BetaDraws(j,:))
end

figure;
subplot(2,1,1)
plot(PosteriorMdl.Sigma2Draws);
title('Trace plot - Disturbance Variance');
xlabel('MCMC draw')
ylabel('Simulation index')
subplot(2,1,2)
autocorr(PosteriorMdl.Sigma2Draws)

Графики трассировки показывают соответствующее смешивание и никакое переходное поведение. Графики автокорреляции показывают низкую корреляцию среди последующих выборок.

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

Объекты

Функции

Похожие темы

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