semiconjugateblm

Байесова модель линейной регрессии с полусопряженным, предшествующим для вероятности данных

Описание

Байесов объект модели линейной регрессии semiconjugateblm указывает, что условное предшествующее распределение β |σ2 многомерно Гауссов со средним μ и отклонением V, и предшествующее распределение σ 2 является обратной гаммой с формой A и шкала B. А именно, Байесовой моделью линейной регрессии является independent, normal-inverse-gamma semiconjugate model.

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

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

Создание

Описание

пример

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

пример

PriorMdl = semiconjugateblm(NumPredictors,Name,Value) дополнительные опции использования заданы одним или несколькими Name,Value парные аргументы. Name имя свойства, кроме NumPredictors, и Value соответствующее значение. Name должен появиться в одинарных кавычках (''). Можно задать несколько Name,Value парные аргументы в любом порядке как Name1,Value1,...,NameN,ValueN.

Свойства

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

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

PriorMdl.V = 100*eye(3);

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

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

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

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

  • Имена переменных (VarNames)

  • Предшествующее среднее значение β (Mu)

  • Предшествующая ковариационная матрица β (V)

Типы данных: 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

Средний параметр Гауссова предшествующего на β в виде числового скаляра или вектора.

Если Mu вектор, затем он должен иметь NumPredictors или NumPredictors + 1 элементы.

  • Для NumPredictors элементы, semiconjugateblm устанавливает предшествующее среднее значение NumPredictors предикторы только. Предикторы соответствуют столбцам в данных о предикторе (заданный во время оценки, симуляции, или предсказывающий). semiconjugateblm игнорирует прерывание в модели, то есть, semiconjugateblm задает предшествующее среднее значение по умолчанию к любому прерыванию.

  • Для NumPredictors + 1 элементы, первый элемент соответствует предшествующему среднему значению прерывания, и все другие элементы соответствуют предикторам.

Пример: 'Mu',[1; 0.08; 2]

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

Условная ковариационная матрица Гауссовых, предшествующих на β в виде c- c симметричная положительная определенная матрица. c может быть NumPredictors или NumPredictors + 1.

  • Если c NumPredictorsто semiconjugateblm устанавливает предшествующую ковариационную матрицу на

    [1e5000V0].

    semiconjugateblm приписывает предшествующие ковариации по умолчанию прерыванию и приписывает V к коэффициентам переменных предикторов в данных. Строки и столбцы V соответствуйте столбцам (переменные) в данных о предикторе.

  • Если c NumPredictors + 1то semiconjugateblm устанавливает целую предшествующую ковариацию на V. Первая строка и столбец соответствует прерыванию. Все другие строки и столбцы соответствуют столбцам в данных о предикторе.

Значением по умолчанию является flat prior. Для adaptive prior задайте diag(Inf(Intercept + NumPredictors,1)). Адаптивное уголовное прошлое указывает на нулевую точность для предшествующего распределения, чтобы иметь как можно меньше влияния на апостериорное распределение.

Пример: 'V',diag(Inf(3,1))

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

Сформируйте гиперпараметр обратной гаммы, предшествующей на σ 2 в виде числового скаляра.

A должен быть, по крайней мере, –(Intercept + NumPredictors)/2.

С B сохраненный зафиксированный, обратное гамма распределение становится более высоким и более сконцентрированным как A увеличения. Эта характеристика взвешивает предшествующую модель σ 2 в большей степени, чем вероятность во время следующей оценки.

Для функциональной формы обратного гамма распределения смотрите Аналитически Послушное Последующее поколение.

Пример: 'A',0.1

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

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

С A сохраненный зафиксированный, обратное гамма распределение становится более высоким и более сконцентрированным как B увеличения. Эта характеристика взвешивает предшествующую модель σ 2 в большей степени, чем вероятность во время следующей оценки.

Пример: 'B',5

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

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

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

Примеры

свернуть все

Рассмотрите модель многофакторной линейной регрессии, которая предсказывает американский действительный валовой национальный продукт (GNPR) использование линейной комбинации индекса промышленного производства (IPI), общая занятость (E), и действительная заработная плата (WR).

GNPRt=β0+β1IPIt+β2Et+β3WRt+εt.

\forall t моменты времени, εt серия независимых Гауссовых воздействий со средним значением 0 и отклонение σ2.

Примите, что предшествующие распределения:

  • β|σ2N4(M,V). M 4 1 вектор из средних значений, и V масштабированная положительная определенная ковариационная матрица 4 на 4.

  • σ2IG(A,B). A и B форма и шкала, соответственно, обратного гамма распределения.

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

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

p = 3;
Mdl = bayeslm(p,'ModelType','semiconjugate')
Mdl = 
  semiconjugateblm with properties:

    NumPredictors: 3
        Intercept: 1
         VarNames: {4x1 cell}
               Mu: [4x1 double]
                V: [4x4 double]
                A: 3
                B: 1

 
           |  Mean     Std           CI95         Positive     Distribution    
-------------------------------------------------------------------------------
 Intercept |  0       100    [-195.996, 195.996]    0.500   N (0.00, 100.00^2) 
 Beta(1)   |  0       100    [-195.996, 195.996]    0.500   N (0.00, 100.00^2) 
 Beta(2)   |  0       100    [-195.996, 195.996]    0.500   N (0.00, 100.00^2) 
 Beta(3)   |  0       100    [-195.996, 195.996]    0.500   N (0.00, 100.00^2) 
 Sigma2    | 0.5000  0.5000    [ 0.138,  1.616]     1.000   IG(3.00,    1)     
 

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

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

Mdl.VarNames = ["IPI" "E" "WR"]
Mdl = 
  semiconjugateblm with properties:

    NumPredictors: 3
        Intercept: 1
         VarNames: {4x1 cell}
               Mu: [4x1 double]
                V: [4x4 double]
                A: 3
                B: 1

 
           |  Mean     Std           CI95         Positive     Distribution    
-------------------------------------------------------------------------------
 Intercept |  0       100    [-195.996, 195.996]    0.500   N (0.00, 100.00^2) 
 IPI       |  0       100    [-195.996, 195.996]    0.500   N (0.00, 100.00^2) 
 E         |  0       100    [-195.996, 195.996]    0.500   N (0.00, 100.00^2) 
 WR        |  0       100    [-195.996, 195.996]    0.500   N (0.00, 100.00^2) 
 Sigma2    | 0.5000  0.5000    [ 0.138,  1.616]     1.000   IG(3.00,    1)     
 

Полагайте, что модель линейной регрессии в Создает Нормальную Обратную Гамму Полусопряженная Предшествующая Модель.

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

p = 3;
PriorMdl = bayeslm(p,'ModelType','semiconjugate','VarNames',["IPI" "E" "WR"]);

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

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

Оцените крайние апостериорные распределения β и σ2.

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  
 

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

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

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

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

По умолчанию, estimate чертит и отбрасывает выборку выжигания дефектов размера 5,000. Однако это - хорошая практика, чтобы смотреть график трассировки ничьих для соответствующего смешивания и отсутствия быстротечности. Постройте график трассировки ничьих для каждого параметра. Можно получить доступ к ничьим, которые составляют распределение, то есть, свойства BetaDraws и Sigma2Draws, использование записи через точку.

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

figure;
plot(PosteriorMdl.Sigma2Draws);
title('Sigma2');

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

Полагайте, что модель линейной регрессии в Создает Нормальную Обратную Гамму Полусопряженная Предшествующая Модель.

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

p = 3;
PriorMdl = bayeslm(p,'ModelType','semiconjugate','VarNames',["IPI" "E" "WR"])
PriorMdl = 
  semiconjugateblm with properties:

    NumPredictors: 3
        Intercept: 1
         VarNames: {4x1 cell}
               Mu: [4x1 double]
                V: [4x4 double]
                A: 3
                B: 1

 
           |  Mean     Std           CI95         Positive     Distribution    
-------------------------------------------------------------------------------
 Intercept |  0       100    [-195.996, 195.996]    0.500   N (0.00, 100.00^2) 
 IPI       |  0       100    [-195.996, 195.996]    0.500   N (0.00, 100.00^2) 
 E         |  0       100    [-195.996, 195.996]    0.500   N (0.00, 100.00^2) 
 WR        |  0       100    [-195.996, 195.996]    0.500   N (0.00, 100.00^2) 
 Sigma2    | 0.5000  0.5000    [ 0.138,  1.616]     1.000   IG(3.00,    1)     
 

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

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

Оцените условное апостериорное распределение β учитывая данные и σ2=2, и возвратите сводную таблицу оценки, чтобы получить доступ к оценкам.

[Mdl,Summary] = 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.2452  1.8693  [-27.909, -20.581]    0.000   N (-24.25, 1.87^2) 
 IPI       |   4.3914  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.4683  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 фиксируется в 2 во время оценки, выводы на ней тривиальны.

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

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

  -24.2452
    4.3914
    0.0011
    2.4683

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

    3.4944    0.0349   -0.0001    0.0241
    0.0349    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 = 
  semiconjugateblm with properties:

    NumPredictors: 3
        Intercept: 1
         VarNames: {4x1 cell}
               Mu: [4x1 double]
                V: [4x4 double]
                A: 3
                B: 1

 
           |  Mean     Std           CI95         Positive     Distribution    
-------------------------------------------------------------------------------
 Intercept |  0       100    [-195.996, 195.996]    0.500   N (0.00, 100.00^2) 
 IPI       |  0       100    [-195.996, 195.996]    0.500   N (0.00, 100.00^2) 
 E         |  0       100    [-195.996, 195.996]    0.500   N (0.00, 100.00^2) 
 WR        |  0       100    [-195.996, 195.996]    0.500   N (0.00, 100.00^2) 
 Sigma2    | 0.5000  0.5000    [ 0.138,  1.616]     1.000   IG(3.00,    1)     
 

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

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

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

p = 3;
PriorMdl = bayeslm(p,'ModelType','semiconjugate','VarNames',["IPI" "E" "WR"]);

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

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  
 

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

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

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

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

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

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

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

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

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

probWRLT2p5 = mean(isWRLT2p5)
probWRLT2p5 = 0.5283

Апостериорная вероятность, что коэффициент действительной заработной платы меньше 2.5, о 0.53.

Copyright 2018 The MathWorks, Inc.

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

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

p = 3;
PriorMdl = bayeslm(p,'ModelType','semiconjugate','VarNames',["IPI" "E" "WR"]);

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

rng(1); % For reproducibility
PosteriorMdl = estimate(PriorMdl,X,y,'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: 1909 - 1970');
ylabel('rGNP');
xlabel('Year');
hold off

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

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

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

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

Больше о

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

Алгоритмы

Можно сбросить все свойства модели с помощью записи через точку, например, PriorMdl.V = diag(Inf(3,1)). Для сброса свойства, semiconjugateblm делает минимальную проверку ошибок значений. Минимизация проверки ошибок имеет преимущество сокращения накладных расходов на симуляции Монте-Карло Цепи Маркова, который приводит к эффективному осуществлению алгоритма.

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

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

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