Подходящие параметры Байесовой модели линейной регрессии к данным
Чтобы выполнить выбор переменного предиктора для Байесовой модели линейной регрессии, смотрите estimate.
PosteriorMdl = estimate(PriorMdl,X,y)PriorMdl задает объединенное предшествующее распределение параметров и структуру модели линейной регрессии. X данные о предикторе и y данные об ответе. PriorMdl и PosteriorMdl не может быть тот же тип объекта.
Произвести PosteriorMdl, estimate функционируйте обновляет предшествующее распределение с информацией о параметрах, которые это получает из данных.
NaNs в данных указывают на отсутствующие значения, который estimate удаляет при помощи мудрого списком удаления.
PosteriorMdl = estimate(PriorMdl,X,y,Name,Value)
Если вы задаете Beta или Sigma2 аргумент пары "имя-значение", затем PosteriorMdl и PriorMdl равны.
[ использование любая из комбинаций входных аргументов в предыдущих синтаксисах, чтобы возвратить таблицу, которая содержит следующее для каждого параметра: следующее среднее и стандартное отклонение, 95%-й вероятный интервал, апостериорная вероятность, что параметр больше 0, и описание апостериорного распределения (если вы существуете). Кроме того, таблица содержит следующую ковариационную матрицу β и σ 2. Если вы задаете PosteriorMdl,Summary]
= estimate(___)Beta или Sigma2 аргумент пары "имя-значение", затем 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) использование линейной комбинации общей занятости (E) и действительная заработная плата (WR).
Для всех
,
серия независимых Гауссовых воздействий со средним значением 0 и отклонение
. Примите эти предшествующие распределения:
3-D t распределение с 10 степенями свободы для каждого компонента, корреляционная матрица 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, но принимает значения параметров только и содержит гиперзначения параметров, зафиксированные в произвольно выбранных значениях.
prednames = ["E" "WR"]; p = numel(prednames); numcoeff = p + 1; rng(1); % For reproducibility dof = 10; V = rand(numcoeff); Sigma = 0.5*(V + V') + numcoeff*eye(numcoeff); st = sqrt(diag(Sigma)); C = diag(1./st)*Sigma*diag(1./st); ct = rand(numcoeff,1); a = 10*rand; b = 10*rand; logPDF = @(params)priorMVTIG(params,ct,st,dof,C,a,b);
Создайте пользовательскую объединенную предшествующую модель для параметров линейной регрессии. Задайте количество предикторов p. Кроме того, задайте указатель на функцию для priorMVTIG и имена переменных.
PriorMdl = bayeslm(p,'ModelType','custom','LogPDF',logPDF,... 'VarNames',prednames);
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: 3
| Mean Std CI95 Positive Distribution
------------------------------------------------------------------------------
Intercept | -3.7453 5.6656 [-16.096, 6.164] 0.243 Empirical
E | -0.0056 0.0006 [-0.007, -0.004] 0.000 Empirical
WR | 15.2497 0.7697 [13.744, 16.736] 1.000 Empirical
Sigma2 | 1289.5491 242.0075 [907.338, 1832.866] 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,SummaryBeta] = 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
estimate отображает сводные данные условного апостериорного распределения . Поскольку фиксируется в 2 во время оценки, выводы на ней тривиальны.
Извлеките средний вектор и ковариационную матрицу условного выражения, следующего из из сводной таблицы оценки.
condPostMeanBeta = SummaryBeta.Mean(1:(end - 1))
condPostMeanBeta = 4×1
-24.2536
4.3913
0.0011
2.4682
CondPostCovBeta = SummaryBeta.Covariances(1:(end - 1),1:(end - 1))
CondPostCovBeta = 4×4
3.4956 0.0350 -0.0001 0.0241
0.0350 0.0009 -0.0000 -0.0013
-0.0001 -0.0000 0.0000 -0.0000
0.0241 -0.0013 -0.0000 0.0055
Отобразите 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.
[~,SummarySigma2] = 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)
estimate отображает сводные данные условного апостериорного распределения . Поскольку фиксируется к condPostMeanBeta во время оценки выводы на нем тривиальны.
Извлеките среднее значение и отклонение условного выражения, следующего из из сводной таблицы оценки.
condPostMeanSigma2 = SummarySigma2.Mean(end)
condPostMeanSigma2 = 48.5138
CondPostVarSigma2 = SummarySigma2.Covariances(end,end)
CondPostVarSigma2 = 81.1581
Считайте модель регрессии в Оценке Следующей Используя гамильтонов Сэмплер Монте-Карло. Этот пример использует те же данные и контекст, но принимает полусопряженную предшествующую модель вместо этого.
Создайте полусопряженную предшествующую модель для параметров линейной регрессии. Задайте количество предикторов 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'};
Оцените крайние апостериорные распределения и .
rng(1); % For reproducibility
[PosteriorMdl,Summary] = 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
PosteriorMdl empiricalblm объект модели, потому что крайние апостериорные распределения полусопряженных моделей аналитически тяжелы, таким образом, estimate должен реализовать сэмплер Гиббса. Summary таблица, содержащая оценки, и заключает тот estimate отображения в командной строке.
Отобразите сводную таблицу.
Summary
Summary=5×6 table
Mean Std CI95 Positive Distribution Covariances
_________ __________ ________________________ ________ _____________ ______________________________________________________________________
Intercept -23.992 9.052 -41.734 -6.1976 0.0053 {'Empirical'} 81.938 0.81622 -0.0025308 0.58928 0
IPI 4.3929 0.14578 4.1011 4.6782 1 {'Empirical'} 0.81622 0.021252 -7.1663e-06 -0.030939 0
E 0.0011124 0.00033976 0.00045128 0.0017883 0.9989 {'Empirical'} -0.0025308 -7.1663e-06 1.1544e-07 -8.4598e-05 0
WR 2.4711 0.3576 1.7622 3.1781 1 {'Empirical'} 0.58928 -0.030939 -8.4598e-05 0.12788 0
Sigma2 46.747 8.455 33.099 66.126 1 {'Empirical'} 0 0 0 0 71.487
Доступ к 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 черные прямоугольники составляют выборку Монте-Карло.