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