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