Байесова модель линейной регрессии с пользовательским объединенным предшествующим распределением
Байесов объект модели линейной регрессии customblm
содержит журнал PDF объединенного предшествующего распределения (β, σ 2). Логарифмическая PDF является пользовательской функцией, которую вы объявляете.
Вероятность данных где ϕ (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
, исключите любой термин прерывания из значения.
После создания модели, если вы изменяете значение NumPredictors
использование записи через точку, затем VarNames
возвращается к его значению по умолчанию.
Типы данных: double
Intercept
— Отметьте для включения прерывания модели регрессииtrue
(значение по умолчанию) | false
Отметьте для включения прерывания модели регрессии в виде значения в этой таблице.
Значение | Описание |
---|---|
false | Исключите прерывание из модели регрессии. Поэтому β является p - размерный вектор, где p значение NumPredictors . |
true | Включайте прерывание в модель регрессии. Поэтому β (p + 1) - размерный вектор. Эта спецификация заставляет T-by-1 вектор из единиц предварительно ожидаться к данным о предикторе во время оценки и симуляции. |
Если вы включаете столбец из единиц в данных о предикторе для термина прерывания, то установленный Intercept
к false
.
Пример: 'Intercept',false
Типы данных: логический
VarNames
— Имена переменного предиктораПеременный предиктор называет для отображений в виде вектора строки или вектора ячейки из векторов символов. VarNames
должен содержать NumPredictors
элементы. VarNames (
имя переменной в столбце j
)j
из набора данных предиктора, который вы задаете во время оценки, симуляции или прогнозирования.
Значением по умолчанию является {'Бета (1)', 'Бета (2)..., Бета (
, где p
)}p
значение NumPredictors
.
Пример: 'VarNames',["UnemploymentRate"; "CPI"]
Типы данных: string
| cell
| char
LogPDF
— Журнал объединенной функции плотности вероятности (β, σ 2)Журнал объединенной функции плотности вероятности (β, σ 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
).
Навсегда точки, серия независимых Гауссовых воздействий со средним значением 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,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 во время оценки, выводы на ней тривиальны.
Извлеките вектор средних значений и ковариационную матрицу условного выражения, следующего из из сводной таблицы оценки.
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 является лучше всего выполняющей моделью тех сравниваемых.
Bayesian linear regression model обрабатывает параметры β и σ 2 в модели yt многофакторной линейной регрессии (MLR) = xt β + εt как случайные переменные.
В течение многих времен t = 1..., T:
yt является наблюдаемым ответом.
xt является 1 на (p + 1) вектор-строка из наблюдаемых величин предикторов p. Вмещать прерывание модели, x 1t = 1 для всего t.
β (p + 1)-by-1 вектор-столбец коэффициентов регрессии, соответствующих переменным, которые составляют столбцы xt.
εt является случайным воздействием со средним значением нуля и Cov (ε) = σ 2IT×T, в то время как ε является T-by-1 вектор, содержащий все воздействия. Эти предположения подразумевают, что вероятность данных
ϕ (yt; xtβ, σ 2) является Гауссовой плотностью вероятности со средним xtβ и отклонением σ 2 оцененных в yt;.
Прежде, чем рассмотреть данные, вы налагаете предположение joint prior distribution на (β, σ 2). В Байесовом анализе вы обновляете распределение параметров при помощи информации о параметрах, полученных из вероятности данных. Результатом является joint posterior distribution (β, σ 2) или conditional posterior distributions параметров.
bayeslm
функция может создать любой поддерживаемый предшествующий объект модели для Байесовой линейной регрессии.
У вас есть модифицированная версия этого примера. Вы хотите открыть этот пример со своими редактированиями?
1. Если смысл перевода понятен, то лучше оставьте как есть и не придирайтесь к словам, синонимам и тому подобному. О вкусах не спорим.
2. Не дополняйте перевод комментариями “от себя”. В исправлении не должно появляться дополнительных смыслов и комментариев, отсутствующих в оригинале. Такие правки не получится интегрировать в алгоритме автоматического перевода.
3. Сохраняйте структуру оригинального текста - например, не разбивайте одно предложение на два.
4. Не имеет смысла однотипное исправление перевода какого-то термина во всех предложениях. Исправляйте только в одном месте. Когда Вашу правку одобрят, это исправление будет алгоритмически распространено и на другие части документации.
5. По иным вопросам, например если надо исправить заблокированное для перевода слово, обратитесь к редакторам через форму технической поддержки.