Подходящие параметры Байесовой модели линейной регрессии к данным
Чтобы выполнить выбор переменной прогноза для Байесовой модели линейной регрессии, смотрите estimate.
PosteriorMdl = estimate(PriorMdl,X,y)PosteriorMdl = estimate(PriorMdl,X,y,Name,Value)[PosteriorMdl,estBeta,EstBetaCov,estSigma2,estSigma2Var]
= estimate(___)[PosteriorMdl,estBeta,EstBetaCov,estSigma2,estSigma2Var,Summary]
= estimate(___)PosteriorMdl = estimate(PriorMdl,X,y)PriorMdl задает объединенное предшествующее распределение параметров и структуру модели линейной регрессии. X является данными о предикторе, и y является данными об ответе. PriorMdl и PosteriorMdl не могут быть тем же типом объекта.
Чтобы произвести PosteriorMdl, функция estimate обновляет предшествующее распределение с информацией о параметрах, которые это получает из данных.
NaN s в данных указывает на отсутствующие значения, которые estimate удаляет при помощи мудрого списком удаления.
PosteriorMdl = estimate(PriorMdl,X,y,Name,Value)
Если вы задаете Beta или Sigma2, то PosteriorMdl и PriorMdl равны.
Следующие два синтаксиса будут удалены в будущем релизе. Для получения дополнительной информации см. Вопросы совместимости.
[ использование любая из комбинаций входных аргументов в предыдущих синтаксисах, чтобы возвратиться:PosteriorMdl,estBeta,EstBetaCov,estSigma2,estSigma2Var]
= estimate(___)
estBeta — Предполагаемые коэффициенты регрессии, то есть, средний вектор апостериорного распределения β
EstBetaCov — Предполагаемая ковариационная матрица содействующих оценок, то есть, ковариационная матрица апостериорного распределения β
EstSigma2 — Предполагаемое отклонение воздействия, то есть, среднее значение апостериорного распределения σ 2
EstSigma2Var — Предполагаемое отклонение отклонения воздействия, то есть, отклонения апостериорного распределения σ 2
Если вы задаете Beta или Sigma2, то estimate возвращает условные следующие оценки. В противном случае estimate возвращает объединенные следующие оценки.
[ также возвращает таблицу, которая содержит следующее для каждого параметра: следующие средние значения и стандартные отклонения, 95%-е вероятные интервалы, апостериорная вероятность, что параметр больше, чем 0, и описание апостериорного распределения (если вы существуете).PosteriorMdl,estBeta,EstBetaCov,estSigma2,estSigma2Var,Summary]
= estimate(___)
Рассмотрите модель, которая предсказывает экономию топлива (в MPG) автомобиля, учитывая его объем двигателя и вес.
Загрузите набор данных carsmall.
load carsmall
x = [Displacement Weight];
y = MPG;Экономия топлива регресса на объем двигателя и вес, включая прерывание, чтобы получить оценки обычных наименьших квадратов (OLS).
Mdl = fitlm(x,y)
Mdl =
Linear regression model:
y ~ 1 + x1 + x2
Estimated Coefficients:
Estimate SE tStat pValue
__________ _________ _______ __________
(Intercept) 46.925 2.0858 22.497 6.0509e-39
x1 -0.014593 0.0082695 -1.7647 0.080968
x2 -0.0068422 0.0011337 -6.0353 3.3838e-08
Number of observations: 94, Error degrees of freedom: 91
Root Mean Squared Error: 4.09
R-squared: 0.747, Adjusted R-Squared: 0.741
F-statistic vs. constant model: 134, p-value = 7.22e-28
Mdl.MSE
ans = 16.7100
Создайте значение по умолчанию, рассейте предшествующее распределение для одного предиктора.
p = 2; PriorMdl = bayeslm(p);
PriorMdl является объектом модели diffuseblm.
Используйте опции по умолчанию, чтобы оценить апостериорное распределение.
PosteriorMdl = estimate(PriorMdl,x,y);
Method: Analytic posterior distributions
Number of observations: 94
Number of predictors: 3
| Mean Std CI95 Positive Distribution
--------------------------------------------------------------------------------
Intercept | 46.9247 2.1091 [42.782, 51.068] 1.000 t (46.92, 2.09^2, 91)
Beta(1) | -0.0146 0.0084 [-0.031, 0.002] 0.040 t (-0.01, 0.01^2, 91)
Beta(2) | -0.0068 0.0011 [-0.009, -0.005] 0.000 t (-0.01, 0.00^2, 91)
Sigma2 | 17.0855 2.5905 [12.748, 22.866] 1.000 IG(45.50, 0.0013)
PosteriorMdl является объектом модели conjugateblm.
Следующие средние значения и содействующие оценки OLS почти идентичны. Кроме того, следующие стандартные отклонения и стандартные погрешности OLS почти идентичны. Следующее среднее значение Sigma2 близко к среднеквадратической ошибке (MSE) OLS.
Считайте несколько моделью линейной регрессии, которая предсказывает США действительный валовой национальный продукт (GNPR) с помощью линейной комбинации индекса промышленного производства (IPI), общая занятость (E) и действительная заработная плата (WR).
Для всех
,
серия независимых Гауссовых воздействий со средним значением 0 и отклонение
. Примите эти предшествующие дистрибутивы:
4-D t распределение с 30 степенями свободы для каждого компонента, корреляционная матрица C, местоположение ct и шкала st.
, с формой
и шкалой
.
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, но принимает значения параметров только и содержит гиперзначения параметров, зафиксированные к произвольно выбранным значениям.
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)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'};
Оцените крайние апостериорные распределения и использование сэмплера Гамильтонова Монте-Карло (HMC). Задайте выборки рисунка 10,000, и электротермотренировка 1 000 чертит.
PosteriorMdl = estimate(PriorMdl,X,y,'Sampler','hmc','NumDraws',1e4,... 'Burnin',1e3);
Method: MCMC sampling with 10000 draws
Number of observations: 62
Number of predictors: 4
| Mean Std CI95 Positive Distribution
--------------------------------------------------------------------------
Intercept | -23.8681 3.1357 [-30.155, -17.617] 0.000 Empirical
IPI | 4.3819 0.1089 [ 4.165, 4.596] 1.000 Empirical
E | 0.0011 0.0002 [ 0.001, 0.002] 1.000 Empirical
WR | 2.4998 0.3355 [ 1.849, 3.167] 1.000 Empirical
Sigma2 | 39.9245 6.7843 [28.684, 55.345] 1.000 Empirical
PosteriorMdl является объектом модели empiricalblm, хранящим ничьи от апостериорных распределений.
Просмотрите график трассировки и график ACF ничьих от следующего из
(например), и отклонение воздействия. Не стройте электротермотренировку.
figure; subplot(2,1,1) plot(PosteriorMdl.BetaDraws(2,1001:end)); title(['Trace Plot ' char(8212) ' \beta_1']); xlabel('MCMC Draw') ylabel('Simulation Index') subplot(2,1,2) autocorr(PosteriorMdl.BetaDraws(2,1001:end)) figure; subplot(2,1,1) plot(PosteriorMdl.Sigma2Draws(1001:end)); title(['Trace Plot ' char(8212) ' Disturbance Variance']); xlabel('MCMC Draw') ylabel('Simulation Index') subplot(2,1,2) autocorr(PosteriorMdl.Sigma2Draws(1001:end))
Выборка MCMC отклонения воздействия, кажется, смешивается хорошо.
Считайте модель регрессии в Оценке Следующей Используя гамильтонов Сэмплер Монте-Карло. Этот пример использует те же данные и контекст, но принимает рассеянную предшествующую модель вместо этого.
Создайте рассеянную предшествующую модель для параметров линейной регрессии. Задайте количество предикторов p и имена коэффициентов регрессии.
p = 3; PriorMdl = bayeslm(p,'ModelType','diffuse','VarNames',["IPI" "E" "WR"])
PriorMdl =
diffuseblm with properties:
NumPredictors: 3
Intercept: 1
VarNames: {4x1 cell}
| Mean Std CI95 Positive Distribution
-----------------------------------------------------------------------------
Intercept | 0 Inf [ NaN, NaN] 0.500 Proportional to one
IPI | 0 Inf [ NaN, NaN] 0.500 Proportional to one
E | 0 Inf [ NaN, NaN] 0.500 Proportional to one
WR | 0 Inf [ NaN, NaN] 0.500 Proportional to one
Sigma2 | Inf Inf [ NaN, NaN] 1.000 Proportional to 1/Sigma2
PriorMdl является объектом модели diffuseblm.
Загрузите набор данных Нельсона-Плоссера. Создайте переменные для ряда предиктора и ответа.
load Data_NelsonPlosser X = DataTable{:,PriorMdl.VarNames(2:end)}; y = DataTable{:,'GNPR'};
Оцените условные апостериорные распределения учитывая данные и это .
[Mdl,condPostMeanBeta,CondPostCovBeta] = estimate(PriorMdl,X,y,... 'Sigma2',2);
Method: Analytic posterior distributions
Conditional variable: Sigma2 fixed at 2
Number of observations: 62
Number of predictors: 4
| Mean Std CI95 Positive Distribution
--------------------------------------------------------------------------------
Intercept | -24.2536 1.8696 [-27.918, -20.589] 0.000 N (-24.25, 1.87^2)
IPI | 4.3913 0.0301 [ 4.332, 4.450] 1.000 N (4.39, 0.03^2)
E | 0.0011 0.0001 [ 0.001, 0.001] 1.000 N (0.00, 0.00^2)
WR | 2.4682 0.0743 [ 2.323, 2.614] 1.000 N (2.47, 0.07^2)
Sigma2 | 2 0 [ 2.000, 2.000] 1.000 Fixed value
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 вектор средних значений (condPostMeanBeta) и ковариационная матрица 4 на 4 (CondPostCovBeta) условного апостериорного распределения учитывая, что данные и это . Кроме того, estimate отображает сводные данные условного апостериорного распределения . Поскольку фиксируется во время оценки, выводы на ней тривиальны.
Предупреждение указывает, что в будущем релизе синтаксисы estimate изменятся. В это время не обновляйте свой код. Для получения дополнительной информации смотрите Заменяющий Нежелательные Синтаксисы оценки.
Отобразите Mdl.
Mdl
Mdl =
diffuseblm with properties:
NumPredictors: 3
Intercept: 1
VarNames: {4x1 cell}
| Mean Std CI95 Positive Distribution
-----------------------------------------------------------------------------
Intercept | 0 Inf [ NaN, NaN] 0.500 Proportional to one
IPI | 0 Inf [ NaN, NaN] 0.500 Proportional to one
E | 0 Inf [ NaN, NaN] 0.500 Proportional to one
WR | 0 Inf [ NaN, NaN] 0.500 Proportional to one
Sigma2 | Inf Inf [ NaN, NaN] 1.000 Proportional to 1/Sigma2
Поскольку estimate вычисляет условное апостериорное распределение, он возвращает исходную предшествующую модель, не следующее, в первом положении списка выходных аргументов.
Оцените условные апостериорные распределения учитывая, что condPostMeanBeta.
[~,~,~,condPostMeanSigma2,condPostVarSigma2] = estimate(PriorMdl,X,y,... 'Beta',condPostMeanBeta);
Method: Analytic posterior distributions
Conditional variable: Beta fixed at -24.2536 4.3913 0.00112035 2.46823
Number of observations: 62
Number of predictors: 4
| Mean Std CI95 Positive Distribution
--------------------------------------------------------------------------------
Intercept | -24.2536 0 [-24.254, -24.254] 0.000 Fixed value
IPI | 4.3913 0 [ 4.391, 4.391] 1.000 Fixed value
E | 0.0011 0 [ 0.001, 0.001] 1.000 Fixed value
WR | 2.4682 0 [ 2.468, 2.468] 1.000 Fixed value
Sigma2 | 48.5138 9.0088 [33.984, 69.098] 1.000 IG(31.00, 0.00069)
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 возвращает среднее значение (condPostMeanSigma2) и отклонение (CondPostVarSigma2) условного апостериорного распределения учитывая данные и это condPostMeanBeta. В отображении, выводах на тривиальны.
Считайте модель регрессии в Оценке Следующей Используя гамильтонов Сэмплер Монте-Карло. Этот пример использует те же данные и контекст, но принимает полусопряженную предшествующую модель вместо этого.
Создайте полусопряженную предшествующую модель для параметров линейной регрессии. Задайте количество предикторов p и имена коэффициентов регрессии.
p = 3; PriorMdl = bayeslm(p,'ModelType','semiconjugate',... 'VarNames',["IPI" "E" "WR"]);
PriorMdl является объектом модели semiconjugateblm.
Загрузите набор данных Нельсона-Плоссера. Создайте переменные для ряда предиктора и ответа.
load Data_NelsonPlosser X = DataTable{:,PriorMdl.VarNames(2:end)}; y = DataTable{:,'GNPR'};
Оцените крайние апостериорные распределения и . Затем получите сводную таблицу оценки:
Передача следующей модели к summarize
Извлечение поля MarginalDistributions из возвращенных сводных данных
rng(1); % For reproducibility
PosteriorMdl = estimate(PriorMdl,X,y);Method: Gibbs sampling with 10000 draws
Number of observations: 62
Number of predictors: 4
| Mean Std CI95 Positive Distribution
-------------------------------------------------------------------------
Intercept | -23.9922 9.0520 [-41.734, -6.198] 0.005 Empirical
IPI | 4.3929 0.1458 [ 4.101, 4.678] 1.000 Empirical
E | 0.0011 0.0003 [ 0.000, 0.002] 0.999 Empirical
WR | 2.4711 0.3576 [ 1.762, 3.178] 1.000 Empirical
Sigma2 | 46.7474 8.4550 [33.099, 66.126] 1.000 Empirical
Summary = summarize(PosteriorMdl); Summary = Summary.MarginalDistributions;
PosteriorMdl является объектом модели empiricalblm, потому что крайние апостериорные распределения полусопряженных моделей аналитически тяжелы, таким образом, estimate должен реализовать сэмплер Гиббса. Summary является таблицей, содержащей оценки, и заключает тот estimate отображения в командной строке.
Отобразите сводную таблицу.
Summary
Summary=5×5 table
Mean Std CI95 Positive Distribution
_________ __________ ________________________ ________ ____________
Intercept -23.992 9.052 -41.734 -6.1976 0.0053 'Empirical'
IPI 4.3929 0.14578 4.1011 4.6782 1 'Empirical'
E 0.0011124 0.00033976 0.00045128 0.0017883 0.9989 'Empirical'
WR 2.4711 0.3576 1.7622 3.1781 1 'Empirical'
Sigma2 46.747 8.455 33.099 66.126 1 'Empirical'
Доступ к 95% equitailed вероятный интервал коэффициента регрессии IPI.
Summary.CI95(2,:)
ans = 1×2
4.1011 4.6782
Если PriorMdl является объектом модели empiricalblm. Вы не можете задать Beta или Sigma2. Вы не можете оценить условные апостериорные распределения при помощи эмпирического предшествующего распределения.
Симуляция Монте-Карло подвергается изменению. Если estimate использует симуляцию Монте-Карло, то оценки и выводы могут отличаться, когда вы вызываете estimate многократно при на вид эквивалентных условиях. Чтобы воспроизвести результаты оценки, прежде, чем вызвать estimate, устанавливают seed случайных чисел при помощи rng.
Если estimate выдает ошибку при оценке, что апостериорное распределение с помощью пользовательской предшествующей модели, то пытается настроить начальные значения параметров при помощи BetaStart или Sigma2Start, или пытается настроить заявленный журнал предшествующая функция, и затем восстановить модель. Ошибка может указать, что журналом предшествующего распределения является –Inf в заданных начальных значениях.
Каждый раз, когда предшествующее распределение (PriorMdl) и вероятность данных приводит к аналитически послушному апостериорному распределению, estimate оценивает решения закрытой формы средств оценки Бейеса. В противном случае estimate обращается к симуляции Монте-Карло, чтобы оценить параметры и чертить выводы. Для получения дополнительной информации смотрите Следующую Оценку и Вывод.
Эта фигура иллюстрирует, как estimate уменьшает выборку Монте-Карло использование значений NumDraws, Thin и BurnIn.
Прямоугольники представляют последовательные ничьи от распределения. estimate удаляет белые прямоугольники из выборки Монте-Карло. Остающийся NumDraws черные прямоугольники составляет выборку Монте-Карло.