customblm

Байесова модель линейной регрессии с пользовательским объединенным предшествующим распределением

Описание

Байесов объект модели линейной регрессии customblm содержит журнал PDF объединенного предшествующего распределения (β, σ 2). Логарифмическая PDF является пользовательской функцией, которую вы объявляете.

Вероятность данных t=1Tϕ(yt;xtβ,σ2), где ϕ (yt; xtβ, σ 2) является Гауссовой плотностью вероятности, оцененной в yt со средним xtβ и отклонением σ 2. MATLAB® обрабатывает предшествующую функцию распределения, как будто это неизвестно. Поэтому получившиеся апостериорные распределения не аналитически послушны. Чтобы оценить или симулировать от апостериорных распределений, MATLAB реализует сэмплер среза.

В общем случае, когда вы создаете Байесов объект модели линейной регрессии, он задает объединенное предшествующее распределение и характеристики модели линейной регрессии только. Таким образом, объект модели является шаблоном, предназначенным для дальнейшего использования. А именно, чтобы включить данные в модель для анализа апостериорного распределения, передайте объект модели и данные к соответствующей объектной функции.

Создание

Описание

пример

PriorMdl = customblm(NumPredictors,'LogPDF',LogPDF) создает Байесов объект модели линейной регрессии (PriorMdl) состоявший из NumPredictors предикторы и прерывание и наборы NumPredictors свойство. LogPDF функция, представляющая журнал объединенного предшествующего распределения (β, σ 2). PriorMdl шаблон, который задает предшествующие распределения и размерность β.

пример

PriorMdl = customblm(NumPredictors,'LogPDF',LogPDF,Name,Value) свойства наборов (кроме NumPredictors) использование аргументов пары "имя-значение". Заключите каждое имя свойства в кавычки. Например, customblm(2,'LogPDF',@logprior,'Intercept',false) задает функцию, которая представляет журнал объединенной предшествующей плотности (β, σ 2), и задает модель регрессии с 2 коэффициентами регрессии, но никакое прерывание.

Свойства

развернуть все

Можно установить перезаписываемые значения свойств, когда вы создаете объект модели при помощи синтаксиса аргумента пары "имя-значение", или после того, как вы создаете объект модели при помощи записи через точку. Например, чтобы исключить прерывание из модели, войти

PriorMdl.Intercept = false;

Количество переменных предикторов в Байесовой модели многофакторной линейной регрессии в виде неотрицательного целого числа.

NumPredictors должен совпасть с количеством столбцов в ваших данных о предикторе, которые вы задаете во время оценки модели или симуляции.

При определении NumPredictors, исключите любой термин прерывания из значения.

После создания модели, если вы изменяете значение NumPredictors использование записи через точку, затем VarNames возвращается к его значению по умолчанию.

Типы данных: double

Отметьте для включения прерывания модели регрессии в виде значения в этой таблице.

ЗначениеОписание
falseИсключите прерывание из модели регрессии. Поэтому β является p- размерный вектор, где p значение NumPredictors.
trueВключайте прерывание в модель регрессии. Поэтому β (p + 1) - размерный вектор. Эта спецификация заставляет T-by-1 вектор из единиц предварительно ожидаться к данным о предикторе во время оценки и симуляции.

Если вы включаете столбец из единиц в данных о предикторе для термина прерывания, то установленный Intercept к false.

Пример: 'Intercept',false

Типы данных: логический

Переменный предиктор называет для отображений в виде вектора строки или вектора ячейки из векторов символов. VarNames должен содержать NumPredictors элементы. VarNames (j) имя переменной в столбце j из набора данных предиктора, который вы задаете во время оценки, симуляции или прогнозирования.

Значением по умолчанию является {'Бета (1)', 'Бета (2)..., Бета (p)}, где p значение NumPredictors.

Пример: 'VarNames',["UnemploymentRate"; "CPI"]

Типы данных: string | cell | char

Журнал объединенной функции плотности вероятности (β, σ 2) в виде указателя на функцию в форме @fcnName, где fcnName имя функции.

Предположим logprior имя функции MATLAB, задающей объединенное предшествующее распределение (β, σ 2). Затем logprior должен иметь эту форму.

function [logpdf,glpdf] = logprior(params)
	...
end
где:

  • logpdf числовой скаляр, представляющий журнал объединенной плотности вероятности (β, σ 2).

  • glpdf (Intercept + NumPredictors + 1)-by-1 числовой вектор, представляющий градиент logpdf. Элементы соответствуют элементам params.

    glpdf дополнительный выходной аргумент, и только гамильтонов сэмплер Монте-Карло (см. hmcSampler) применяет его. Если вы знаете аналитическую частную производную относительно некоторых параметров, но не других, то установленный элементы glpdf соответствие неизвестным частным производным к NaN. MATLAB вычисляет числовой градиент для пропавших без вести частных производных, который удобен, но замедляет выборку.

  • params (Intercept + NumPredictors + 1)-by-1 числовой вектор. Первый Intercept + NumPredictors элементы должны соответствовать значениям β, и последний элемент должен соответствовать значению σ 2. Первым элементом β является прерывание, если вы существуете. Все другие элементы соответствуют переменным предикторам в данных о предикторе, которые вы задаете во время оценки, симуляции или прогнозирования.

Пример: 'LogPDF',@logprior

Функции объекта

estimateОцените апостериорное распределение Байесовых параметров модели линейной регрессии
simulateСимулируйте коэффициенты регрессии и отклонение воздействия Байесовой модели линейной регрессии
forecastПредскажите ответы Байесовой модели линейной регрессии
plotВизуализируйте предшествующую и следующую плотность Байесовых параметров модели линейной регрессии
summarizeСтатистика сводных данных распределения стандартной Байесовой модели линейной регрессии

Примеры

свернуть все

Рассмотрите модель многофакторной линейной регрессии, которая предсказывает США действительный валовой национальный продукт (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 распределение с 50 степенями свободы для каждого компонента и единичной матрицей для корреляционной матрицы. Кроме того, распределение сосредоточено в${\left[ {\begin{array}{*{20}{c}} { - 25}&4&0&3 \end{array}} \right]^\prime }$, и каждый компонент масштабируется соответствующими элементами вектора${\left[ {\begin{array}{*{20}{c}} {1}&1&1&1 \end{array}} \right]^\prime }$.

  • $\sigma^2 \sim IG(3,1)$.

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

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

  • Принимает значения$\beta$ и$\sigma^2$ вместе в вектор-столбце и принимает значения гиперпараметров

  • Возвращает значение объединенного предшествующего распределения$\pi\left(\beta,\sigma^2\right)$, учитывая значения$\beta$ и $\sigma^2$

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 = [-25; 4; 0; 3];
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 with properties:

    NumPredictors: 3
        Intercept: 1
         VarNames: {4x1 cell}
           LogPDF: @(params)priorMVTIG(params,ct,st,dof,C,a,b)

The priors are defined by the function:
    @(params)priorMVTIG(params,ct,st,dof,C,a,b)

PriorMdl customblm Байесов объект модели линейной регрессии, представляющий предшествующее распределение коэффициентов регрессии и отклонения воздействия. В этом случае, bayeslm не отображает сводные данные предшествующих распределений в командной строке.

Полагайте, что модель линейной регрессии в Создает Пользовательскую Многомерную t Предшествующую Модель для Коэффициентов.

Создайте анонимную функцию, которая действует как priorMVTIG, но принимает значения параметров только и содержит гиперзначения параметров, зафиксированные в их значениях.

dof = 50;
C = eye(4);
ct = [-25; 4; 0; 3];
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 with properties:

    NumPredictors: 3
        Intercept: 1
         VarNames: {4x1 cell}
           LogPDF: @(params)priorMVTIG(params,ct,st,dof,C,a,b)

The priors are defined by the function:
    @(params)priorMVTIG(params,ct,st,dof,C,a,b)

Загрузите набор данных Нельсона-Плоссера. Создайте переменные для ряда предиктора и ответа.

load Data_NelsonPlosser
X = DataTable{:,PriorMdl.VarNames(2:end)};
y = DataTable{:,'GNPR'};

Оцените крайние апостериорные распределения β и σ2. Задайте ширину для сэмплера среза, который является близко к следующему стандартному отклонению параметров, принимающих рассеянную предшествующую модель. Уменьшайте последовательную корреляцию путем определения утончающегося фактора 10 и сократите эффективное количество по умолчанию ничьих на коэффициент 10.

width = [20,0.5,0.01,1,20];
thin = 10;
numDraws = 1e5/thin;
rng(1) % For reproducibility
PosteriorMdl = estimate(PriorMdl,X,y,'Width',width,'Thin',thin,...
    'NumDraws',numDraws);
Method: MCMC sampling with 10000 draws
Number of observations: 62
Number of predictors:   4
 
           |   Mean      Std          CI95         Positive  Distribution 
--------------------------------------------------------------------------
 Intercept | -25.0069  0.9919  [-26.990, -23.065]    0.000     Empirical  
 IPI       |   4.3544  0.1083   [ 4.143,  4.562]     1.000     Empirical  
 E         |   0.0011  0.0002   [ 0.001,  0.001]     1.000     Empirical  
 WR        |   2.5613  0.3293   [ 1.939,  3.222]     1.000     Empirical  
 Sigma2    |  47.0593  8.7570   [32.690, 67.115]     1.000     Empirical  
 

PosteriorMdl empiricalblm хранение объекта модели чертит от апостериорных распределений β и σ2 учитывая данные. estimate отображает сводные данные крайних апостериорных распределений к командному окну. Строки сводных данных соответствуют коэффициентам регрессии и отклонению воздействия и столбцам к характеристикам апостериорного распределения. Характеристики включают:

  • CI95, который содержит 95%-е Байесовы equitailed вероятные интервалы для параметров. Например, апостериорная вероятность, что коэффициент регрессии WR находится в [1.939, 3.222] 0.95.

  • Positive, который содержит апостериорную вероятность, что параметр больше 0. Например, вероятность, что прерывание больше 0, 0.

estimate выводит следующие характеристики из ничьих от апостериорных распределений, которые MATLAB® хранит как матрицы в свойствах BetaDraws и Sigma2Draws.

Чтобы контролировать смешивание и сходимость выборки MCMC, создайте графики трассировки. В BetaDraws свойство, чертит, соответствуют столбцам и параметрам к строкам.

figure;
for j = 1:4
    subplot(2,2,j)
    plot(PosteriorMdl.BetaDraws(j,:))
    title(sprintf('Trace Plot of %s',PosteriorMdl.VarNames{j}));
end

figure;
plot(PosteriorMdl.Sigma2Draws)
title('Trace Plot of Sigma2');

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

Полагайте, что модель линейной регрессии в Создает Пользовательскую Многомерную t Предшествующую Модель для Коэффициентов.

Создайте анонимную функцию, которая действует как priorMVTIG, но принимает значения параметров только и содержит зафиксированные гиперзначения параметров.

dof = 50;
C = eye(4);
ct = [-25; 4; 0; 3];
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 with properties:

    NumPredictors: 3
        Intercept: 1
         VarNames: {4x1 cell}
           LogPDF: @(params)priorMVTIG(params,ct,st,dof,C,a,b)

The priors are defined by the function:
    @(params)priorMVTIG(params,ct,st,dof,C,a,b)

Загрузите набор данных Нельсона-Плоссера. Создайте переменные для ряда предиктора и ответа.

load Data_NelsonPlosser
X = DataTable{:,PriorMdl.VarNames(2:end)};
y = DataTable{:,'GNPR'};

Оцените условное апостериорное распределение β учитывая данные и σ2=2, и возвратите сводную таблицу оценки, чтобы получить доступ к оценкам. Задайте ширину для сэмплера среза, который является близко к следующему стандартному отклонению параметров, принимающих рассеянную предшествующую модель. Уменьшайте последовательную корреляцию путем определения утончающегося фактора 10 и сократите эффективное количество по умолчанию ничьих на коэффициент 10.

width = [20,0.5,0.01,1];
thin = 10;
numDraws = 1e5/thin;
rng(1) % For reproducibility
[Mdl,Summary] = estimate(PriorMdl,X,y,'Sigma2',2,...
    'Width',width,'Thin',thin,'NumDraws',numDraws);
Method: MCMC sampling with 10000 draws
Conditional variable: Sigma2 fixed at   2
Number of observations: 62
Number of predictors:   4
 
           |   Mean      Std          CI95         Positive  Distribution 
--------------------------------------------------------------------------
 Intercept | -24.7820  0.8767  [-26.483, -23.054]    0.000     Empirical  
 IPI       |   4.3825  0.0254   [ 4.332,  4.431]     1.000     Empirical  
 E         |   0.0011  0.0000   [ 0.001,  0.001]     1.000     Empirical  
 WR        |   2.4752  0.0724   [ 2.337,  2.618]     1.000     Empirical  
 Sigma2    |    2       0       [ 2.000,  2.000]     1.000     Empirical  
 

estimate отображает сводные данные условного апостериорного распределения β. Поскольку σ2 фиксируется в 2 во время оценки, выводы на ней тривиальны.

Извлеките вектор средних значений и ковариационную матрицу условного выражения, следующего из β из сводной таблицы оценки.

condPostMeanBeta = Summary.Mean(1:(end - 1))
condPostMeanBeta = 4×1

  -24.7820
    4.3825
    0.0011
    2.4752

CondPostCovBeta = Summary.Covariances(1:(end - 1),1:(end - 1))
CondPostCovBeta = 4×4

    0.7686    0.0084   -0.0000    0.0019
    0.0084    0.0006    0.0000   -0.0015
   -0.0000    0.0000    0.0000   -0.0000
    0.0019   -0.0015   -0.0000    0.0052

Отобразите Mdl.

Mdl
Mdl = 
  customblm with properties:

    NumPredictors: 3
        Intercept: 1
         VarNames: {4x1 cell}
           LogPDF: @(params)priorMVTIG(params,ct,st,dof,C,a,b)

The priors are defined by the function:
    @(params)priorMVTIG(params,ct,st,dof,C,a,b)

Поскольку estimate вычисляет условное апостериорное распределение, оно возвращает исходную предшествующую модель, не следующее, в первом положении списка выходных аргументов. Кроме того, estimate не возвращает выборку MCMC. Поэтому, чтобы контролировать сходимость выборки MCMC, используйте simulate вместо этого и задайте тот же seed случайных чисел.

Считайте модель линейной регрессии в Оценке Крайними Апостериорными распределениями.

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

dof = 50;
C = eye(4);
ct = [-25; 4; 0; 3];
st = ones(4,1);
a = 3;
b = 1;
logPDF = @(params)priorMVTIG(params,ct,st,dof,C,a,b);

p = 3;
PriorMdl = bayeslm(p,'ModelType','custom','LogPDF',logPDF,...
    'VarNames',["IPI" "E" "WR"]);

load Data_NelsonPlosser
X = DataTable{:,PriorMdl.VarNames(2:end)};
y = DataTable{:,'GNPR'};  

width = [20,0.5,0.01,1,20];
thin = 10;
numDraws = 1e5/thin;

rng(1) % For reproducibility
PosteriorMdl = estimate(PriorMdl,X,y,'Width',width,'Thin',thin,...
    'NumDraws',numDraws,'Display',false);

Оцените статистику сводных данных апостериорного распределения для β при помощи ничьих от апостериорного распределения, сохраненного в следующей модели.

estBeta = mean(PosteriorMdl.BetaDraws,2);
EstBetaCov = cov(PosteriorMdl.BetaDraws');

Предположим что если коэффициент действительной заработной платы (WR) ниже 2.5, затем политика выполнена. Несмотря на то, что апостериорное распределение WR известен, и можно вычислить вероятности непосредственно, можно оценить вероятность с помощью симуляции Монте-Карло вместо этого.

Чертите 1e6 выборки от крайнего апостериорного распределения β.

NumDraws = 1e6;
BetaSim = simulate(PosteriorMdl,'NumDraws',NumDraws);

BetaSim 4-by-1e6 матрица, содержащая ничьи. Строки соответствуют коэффициенту регрессии и столбцам к последовательным ничьим.

Изолируйте ничьи, соответствующие коэффициенту WR, и затем идентифицируйте, который чертит, меньше 2.5.

isWR = PosteriorMdl.VarNames == "WR";
wrSim = BetaSim(isWR,:);
isWRLT2p5 = wrSim < 2.5;

Найдите крайнюю апостериорную вероятность что коэффициент регрессии WR ниже 2.5 путем вычисления пропорции ничьих, которые меньше 2.5.

probWRLT2p5 = mean(isWRLT2p5)
probWRLT2p5 = 0.4430

Апостериорная вероятность, что коэффициент WR меньше 2.5, о 0.4430.

Считайте модель линейной регрессии в Оценке Крайними Апостериорными распределениями.

Создайте предшествующую модель для коэффициентов регрессии и отклонения воздействия, затем оцените крайние апостериорные распределения. Протяните последние 10 периодов данных из оценки, таким образом, можно использовать их, чтобы предсказать действительный GNP. Выключите отображение оценки.

load Data_NelsonPlosser
VarNames = {'IPI'; 'E'; 'WR'};
fhs = 10; % Forecast horizon size
X = DataTable{1:(end - fhs),VarNames};
y = DataTable{1:(end - fhs),'GNPR'};
XF = DataTable{(end - fhs + 1):end,VarNames}; % Future predictor data
yFT = DataTable{(end - fhs + 1):end,'GNPR'};  % True future responses

dof = 50;
C = eye(4);
ct = [-25; 4; 0; 3];
st = ones(4,1);
a = 3;
b = 1;
logPDF = @(params)priorMVTIG(params,ct,st,dof,C,a,b);

p = 3;
PriorMdl = bayeslm(p,'ModelType','custom','LogPDF',logPDF,...
    'VarNames',VarNames);

width = [20,0.5,0.01,1,20];
thin = 10;
numDraws = 1e5/thin;
rng(1) % For reproducibility
PosteriorMdl = estimate(PriorMdl,X,y,'Width',width,'Thin',thin,...
    'NumDraws',numDraws,'Display',false);

Предскажите ответы с помощью следующего прогнозирующего распределения и будущих данных о предикторе XF. Постройте истинные значения ответа и предсказанных значений.

yF = forecast(PosteriorMdl,XF);

figure;
plot(dates,DataTable.GNPR);
hold on
plot(dates((end - fhs + 1):end),yF)
h = gca;
hp = patch([dates(end - fhs + 1) dates(end) dates(end) dates(end - fhs + 1)],...
    h.YLim([1,1,2,2]),[0.8 0.8 0.8]);
uistack(hp,'bottom');
legend('Forecast Horizon','True GNPR','Forecasted GNPR','Location','NW')
title('Real Gross National Product');
ylabel('rGNP');
xlabel('Year');
hold off

yF вектор 10 на 1 из будущих значений действительного GNP, соответствующего будущим данным о предикторе.

Оцените среднеквадратическую ошибку (RMSE) прогноза.

frmse = sqrt(mean((yF - yFT).^2))
frmse = 12.8148

Прогноз RMSE является относительной мерой точности прогноза. А именно, вы оцениваете несколько моделей с помощью различных предположений. Модель с самым низким прогнозом RMSE является лучше всего выполняющей моделью тех сравниваемых.

Больше о

развернуть все

Альтернативы

bayeslm функция может создать любой поддерживаемый предшествующий объект модели для Байесовой линейной регрессии.

Введенный в R2017a
Для просмотра документации необходимо авторизоваться на сайте