Задайте градиент для сэмплера 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.289332e+01 |   8.916e-02 |          56 |   6.700e-01 |           3 |
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.290130e+01 |   4.978e-01 |          10 |   6.400e-01 |           0 |

|==================================================================================|
|   ITER   |    LOG PDF    |  STEP SIZE  |  NUM STEPS  |  ACC RATIO  |  DIVERGENT  |
|==================================================================================|
|      100 |  4.085894e+01 |   3.597e-02 |           6 |   8.000e-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 |  3.849768e+01 |   2.804e-01 |          18 |   6.400e-01 |           0 |

|==================================================================================|
|   ITER   |    LOG PDF    |  STEP SIZE  |  NUM STEPS  |  ACC RATIO  |  DIVERGENT  |
|==================================================================================|
|      100 |  3.808618e+01 |   2.523e-01 |           4 |   7.950e-01 |           0 |
|      200 |  4.087904e+01 |   1.035e-01 |          16 |   8.567e-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.237073e+01 |   3.545e-01 |          14 |   6.400e-01 |           0 |

|==================================================================================|
|   ITER   |    LOG PDF    |  STEP SIZE  |  NUM STEPS  |  ACC RATIO  |  DIVERGENT  |
|==================================================================================|
|      100 |  4.320936e+01 |   1.787e-02 |           2 |   7.750e-01 |           0 |
|      200 |  3.929757e+01 |   3.191e-01 |          14 |   8.433e-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.252996e+01 |   3.253e-01 |          15 |   6.800e-01 |           0 |

|==================================================================================|
|   ITER   |    LOG PDF    |  STEP SIZE  |  NUM STEPS  |  ACC RATIO  |  DIVERGENT  |
|==================================================================================|
|      100 |  4.198629e+01 |   9.387e-02 |           4 |   8.050e-01 |           0 |
|      200 |  4.339092e+01 |   2.928e-01 |           4 |   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.016755e+01 |   3.491e-01 |          14 |   6.600e-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.0086  [-0.009,  0.024]    0.812     Empirical  
 IPI       | -0.1632  0.1446  [-0.473,  0.097]    0.123     Empirical  
 E         |  2.3719  0.5413  [ 1.414,  3.525]    1.000     Empirical  
 WR        | -0.2311  0.2998  [-0.836,  0.343]    0.220     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
Поддержка
Для просмотра документации необходимо авторизоваться на сайте