В этом примере показано, как улучшить сэмплер среза для следующей оценки и вывода Байесовой модели линейной регрессии.
При оценке следующего, состоявшего из вероятности данных и полусопряженных или пользовательских предшествующих моделей, estimate
использует сэмплер MCMC. Если график трассировки выборки показывает переходное поведение или очень высоко последовательную корреляцию, или вы хотите сохранить немного выборок от следующего, то можно задать выборку выжигания дефектов или утончение.
Рассмотрите модель многофакторной линейной регрессии, которая предсказывает США действительный валовой национальный продукт (GNPR
) использование линейной комбинации индекса промышленного производства (IPI
), общая занятость (E
), и действительная заработная плата (WR
).
Навсегда точки, серия независимых Гауссовых воздействий со средним значением 0 и отклонение. Примите, что предшествующие распределения:
4-мерное t распределение с 50 степенями свободы для каждого компонента и единичной матрицей для корреляционной матрицы. Кроме того, распределение сосредоточено в, и каждый компонент масштабируется соответствующими элементами вектора.
.
bayeslm
обработки эти предположения и вероятность данных, как будто следующее соответствие аналитически тяжело.
Объявите функцию MATLAB® что:
Принимает значения и вместе в вектор-столбце и значениях гиперпараметров.
Возвращает значение объединенного предшествующего распределения, учитывая значения и.
function logPDF = priorMVTIG(params,ct,st,dof,C,a,b) %priorMVTIG Log density of multivariate t times inverse gamma % priorMVTIG 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. priorMVTIG returns the log of the product of % the two evaluated densities. % % 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 (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)); sigma2 = params(end); tVal = (beta - ct)./st; mvtDensity = mvtpdf(tVal,C,dof); igDensity = sigma2^(-a-1)*exp(-1/(sigma2*b))/(gamma(a)*b^a); logPDF = log(mvtDensity*igDensity); end
Создайте анонимную функцию, которая действует как priorMVTIG
, но принимает только значения параметров и содержит гиперзначения параметров, зафиксированные в их значениях.
dof = 50; C = eye(4); ct = zeros(4,1); st = ones(4,1); a = 3; b = 1; logPDF = @(params)priorMVTIG(params,ct,st,dof,C,a,b);
Создайте пользовательскую объединенную предшествующую модель для параметров линейной регрессии. Задайте количество предикторов, p
. Кроме того, задайте указатель на функцию для priorMVTIG
и имена переменных.
p = 3; PriorMdl = bayeslm(p,'ModelType','custom','LogPDF',logPDF,... 'VarNames',["IPI" "E" "WR"]);
PriorMdl
customblm
Байесов объект модели линейной регрессии, представляющий предшествующее распределение коэффициентов регрессии и отклонения воздействия.
Загрузите набор данных Нельсона-Плоссера. Создайте переменные для ряда предиктора и ответа.
load Data_NelsonPlosser X = DataTable{:,PriorMdl.VarNames(2:end)}; y = DataTable{:,'GNPR'};
Оцените крайние апостериорные распределения и. Чтобы определить соответствующую электротермотренировку, не задайте электротермотренировку. Задайте ширину для сэмплера среза, который является близко к следующему стандартному отклонению параметров, принимающих рассеянную предшествующую модель. Чертите 15000
выборки.
width = [20,0.5,0.01,1,20]; numDraws = 15e3; rng(1) % For reproducibility PosteriorMdl1 = estimate(PriorMdl,X,y,'BurnIn',0,'Width',width,... 'NumDraws',numDraws);
Warning: High autocorrelation of MCMC draws:IPI,E,WR,Sigma2. Method: MCMC sampling with 15000 draws Number of observations: 62 Number of predictors: 4 | Mean Std CI95 Positive Distribution ----------------------------------------------------------------------- Intercept | -0.6669 1.5325 [-3.202, 1.621] 0.309 Empirical IPI | 4.6196 0.1021 [ 4.437, 4.839] 1.000 Empirical E | 0.0004 0.0002 [ 0.000, 0.001] 0.995 Empirical WR | 2.6135 0.3025 [ 1.982, 3.152] 1.000 Empirical Sigma2 | 48.9990 9.1926 [33.147, 67.993] 1.000 Empirical
PosteriorMdl1
empiricalblm
объект модели. BetaDraws
и Sigma2Draws
свойства хранят выборки от сэмплера среза.
Чертите трассировку и графики автокорреляции для Цепей Маркова.
figure; for j = 1:4 subplot(2,3,j) plot(PosteriorMdl1.BetaDraws(j,:)); axis tight xlabel('Sample index') title(sprintf(['Trace Plot ' char(8212) ' %s'],PriorMdl.VarNames{j})); end subplot(2,3,5) plot(PosteriorMdl1.Sigma2Draws); axis tight xlabel('Sample index') title(['Trace Plot ' char(8212) ' Sigma2']); figure; for j = 1:4 subplot(2,3,j) autocorr(PosteriorMdl1.BetaDraws(j,:)); axis tight title(sprintf(['Correlogram ' char(8212) ' %s'],PriorMdl.VarNames{j})); end subplot(2,3,5) autocorr(PosteriorMdl1.Sigma2Draws); axis tight title(['Correlogram ' char(8212) ' Sigma2']);
Все параметры, кроме прерывания, показывают значительные переходные эффекты и высокую автокорреляцию. estimate
проверки на высокую корреляцию, и если есть выдают предупреждения. Кроме того, параметры, кажется, не сходились к их стационарному распределению.
Оцените апостериорное распределение снова. Чтобы уменьшать переходные эффекты, используйте электротермотренировку по умолчанию (5000
, который является значением по умолчанию). Уменьшать влияние последовательной корреляции, тонкой на коэффициент 20 и задавать 15e4/thin
эффективное количество ничьих. Кроме того, численно стабилизируйтесь, оценка определением повторно параметрируют к. Постройте графики автокорреляции и трассировка.
thin = 20; numDraws = 15e4/thin; PosteriorMdl2 = estimate(PriorMdl,X,y,'Thin',thin,'Width',width,... 'Reparameterize',true,'NumDraws',numDraws); figure; for j = 1:4 subplot(2,3,j) plot(PosteriorMdl2.BetaDraws(j,:)); axis tight xlabel('Sample index') title(sprintf(['Trace Plot ' char(8212) '%s'],PriorMdl.VarNames{j})); end subplot(2,3,5) plot(PosteriorMdl2.Sigma2Draws); axis tight xlabel('Sample index') title(['Trace Plot ' char(8212) ' Sigma2']); figure; for j = 1:4 subplot(2,3,j) autocorr(PosteriorMdl2.BetaDraws(j,:)); title(sprintf(['Correlogram ' char(8212) ' %s'],PriorMdl.VarNames{j})); end subplot(2,3,5) autocorr(PosteriorMdl2.Sigma2Draws); title(['Correlogram ' char(8212) ' Sigma2']);
Method: MCMC sampling with 7500 draws Number of observations: 62 Number of predictors: 4 | Mean Std CI95 Positive Distribution ----------------------------------------------------------------------- Intercept | -0.4777 1.2225 [-2.911, 1.955] 0.338 Empirical IPI | 4.6506 0.1128 [ 4.434, 4.884] 1.000 Empirical E | 0.0005 0.0002 [ 0.000, 0.001] 0.991 Empirical WR | 2.5243 0.3468 [ 1.815, 3.208] 1.000 Empirical Sigma2 | 51.3444 9.3170 [36.319, 72.469] 1.000 Empirical
Графики трассировки предполагают что Цепи Маркова:
Достигли их стационарных распределений, потому что распределение выборок имеет относительно постоянное среднее значение и отклонение. А именно, распределение точек, кажется, не изменяется.
Смешиваются справедливо быстро, потому что последовательная корреляция является умеренной, для некоторых параметров, к очень низко.
Не показывайте переходное поведение.
Эти результаты предполагают, что можно возобновить анализ чувствительности, или следующий вывод.