Задайте градиент для HMC Sampler

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

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

Рассмотрим множественную линейную регрессионую модель, которая предсказывает реальный валовой национальный продукт США (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-D t распределения с 30 степенями свободы для каждого компонента, корреляционная матрица C, местоположение ct, и масштабные st.

$\sigma^2 \sim IG(10r_1,10r_2)$ , с формой и шкалой.

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

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

Принимает значения $\beta$ и $\sigma^2$ вместе в векторе-столбце и принимает значения гиперпараметров.
Возвращает журнал предыдущего распределения соединений, $\pi\left(\beta,\sigma^2\right)$ учитывая значения $\beta$ и, $\sigma^2$ и возвращает градиент предыдущей плотности логарифмических соединений относительно и $\beta$ . $\gamma$

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

Оцените маргинальные апостериорные распределения $\beta$ и $\sigma^2$ использование HMC sampler. Чтобы контролировать прогресс дискретизатора, задайте уровень подробностей 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.340274e+01 |   9.440e-02 |          53 |   6.300e-01 |           4 |
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.323712e+01 |   5.075e-01 |          10 |   6.000e-01 |           0 |

|==================================================================================|
|   ITER   |    LOG PDF    |  STEP SIZE  |  NUM STEPS  |  ACC RATIO  |  DIVERGENT  |
|==================================================================================|
|      100 |  4.074007e+01 |   3.667e-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.125096e+01 |   3.086e-01 |          16 |   6.000e-01 |           0 |

|==================================================================================|
|   ITER   |    LOG PDF    |  STEP SIZE  |  NUM STEPS  |  ACC RATIO  |  DIVERGENT  |
|==================================================================================|
|      100 |  3.842133e+01 |   2.777e-01 |           4 |   7.700e-01 |           0 |
|      200 |  4.096947e+01 |   1.140e-01 |          15 |   8.400e-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.294008e+01 |   3.294e-01 |          15 |   6.300e-01 |           0 |

|==================================================================================|
|   ITER   |    LOG PDF    |  STEP SIZE  |  NUM STEPS  |  ACC RATIO  |  DIVERGENT  |
|==================================================================================|
|      100 |  4.294706e+01 |   1.660e-02 |           2 |   7.850e-01 |           0 |
|      200 |  3.957068e+01 |   2.965e-01 |          15 |   8.533e-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.241792e+01 |   3.173e-01 |          16 |   6.300e-01 |           0 |

|==================================================================================|
|   ITER   |    LOG PDF    |  STEP SIZE  |  NUM STEPS  |  ACC RATIO  |  DIVERGENT  |
|==================================================================================|
|      100 |  4.044149e+01 |   9.156e-02 |           4 |   7.800e-01 |           0 |
|      200 |  4.275379e+01 |   2.856e-01 |           5 |   8.467e-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.068391e+01 |   3.366e-01 |          15 |   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.0076  0.0086  [-0.009,  0.024]    0.808     Empirical  
 IPI       | -0.1635  0.1446  [-0.470,  0.102]    0.121     Empirical  
 E         |  2.3701  0.5429  [ 1.398,  3.526]    1.000     Empirical  
 WR        | -0.2304  0.3000  [-0.837,  0.342]    0.218     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)

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

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

Задайте градиент для HMC Sampler

См. также

Объекты

Функции

Похожие темы

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

Поддержка