Байесова модель линейной регрессии с пользовательским объединенным предшествующим распределением
Байесов объект модели линейной регрессии customblm содержит журнал PDF объединенного предшествующего распределения (β, σ 2). Логарифмический PDF является пользовательской функцией, которую вы объявляете.
Вероятность данных где ϕ (yt; xtβ, σ 2) является Гауссовой плотностью вероятности, оцененной в yt со средним xtβ и отклонением σ 2. MATLAB® обрабатывает предшествующую функцию распределения, как будто это неизвестно. Поэтому получившиеся апостериорные распределения не аналитически послушны. Чтобы оценить или моделировать от апостериорных распределений, MATLAB реализует сэмплер среза.
В целом, когда вы создаете Байесов объект модели линейной регрессии, он задает объединенное предшествующее распределение и характеристики модели линейной регрессии только. Таким образом, объект модели является шаблоном, предназначенным для дальнейшего использования. А именно, чтобы включить данные в модель для анализа апостериорного распределения, передайте объект модели и данные к соответствующей объектной функции.
PriorMdl = customblm(NumPredictors,'LogPDF',LogPDF)PriorMdl = customblm(NumPredictors,'LogPDF',LogPDF,Name,Value)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 коэффициентами регрессии, но никакое прерывание.
estimate | Подходящие параметры Байесовой модели линейной регрессии к данным |
simulate | Моделируйте коэффициенты регрессии и отклонение воздействия Байесовой модели линейной регрессии |
forecast | Предскажите ответы Байесовой модели линейной регрессии |
plot | Визуализируйте предшествующую и следующую плотность Байесовых параметров модели линейной регрессии |
summarize | Статистика сводных данных распределения стандартной Байесовой модели линейной регрессии |
Считайте несколько моделью линейной регрессии, которая предсказывает США действительный валовой национальный продукт (GNPR) с помощью линейной комбинации индекса промышленного производства (IPI), общая занятость (E) и действительная заработная плата (WR).
Навсегда точки,
серия независимых Гауссовых воздействий со средним значением 0 и отклонение
. Примите эти предшествующие дистрибутивы:
4-D 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 = [-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'};
Оцените крайние апостериорные распределения и . Задайте ширину для сэмплера среза, который является близко к следующему стандартному отклонению параметров, принимающих рассеянную предшествующую модель. Уменьшайте последовательную корреляцию путем определения утончающегося фактора 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, чертит от апостериорных распределений и учитывая данные. 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'};
Оцените условные апостериорные распределения учитывая данные и . Задайте ширину для сэмплера среза, который является близко к следующему стандартному отклонению параметров, принимающих рассеянную предшествующую модель. Уменьшайте последовательную корреляцию путем определения утончающегося фактора 10, и, сократите эффективное количество по умолчанию ничьих фактором 10.
width = [20,0.5,0.01,1]; thin = 10; numDraws = 1e5/thin; rng(1) % For reproducibility [Mdl,condPostMeanBeta,CondPostCovBeta] = 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
Warning: Current syntax supports 6 output arguments, and will be removed in a future release. For supported output arguments, see <a href="matlab:helpview(fullfile(docroot,'econ','econ.map'),'blmestimate')">estimate</a>.
estimate возвращается 4 1 вектор средних значений и ковариационная матрица 4 на 4 условного апостериорного распределения в condPostMeanBeta и CondPostCovBeta, соответственно. Кроме того, estimate отображает сводные данные условного апостериорного распределения .
Предупреждение указывает, что в будущем релизе синтаксисы estimate изменятся. В это время не обновляйте свой код. Для получения дополнительной информации смотрите Заменяющий Нежелательные Синтаксисы оценки.
Отобразите 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 может создать любой поддерживаемый предшествующий объект модели для Байесовой линейной регрессии.