exponenta event banner

conjugateblm

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

Описание

Объект модели байесовской линейной регрессии conjugateblm указывает, что совместное предварительное распределение коэффициентов регрессии и дисперсии возмущений, то есть (β, start2), является зависимой, нормально-обратной-гамма-сопряженной моделью. Условное предшествующее распределение β 'start2 является многомерным гауссовым со средним λ и дисперсионным σ2V. Предыдущее распределение start2 - обратная гамма с формой A и шкалой B.

Правдоподобие данных равно ∏t=1Tϕ (yt; xtβ, start2), где (yt; xtβ, start2) - гауссова плотность вероятностей, оцениваемая на yt со средним xtβ и дисперсией start2. Указанные приоры являются сопряженными для правдоподобия, и результирующие краевое и условное апостериорное распределения являются аналитически прослеживаемыми. Для получения подробной информации о заднем распределении см. Аналитически отслеживаемые задние части.

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

Создание

Описание

пример

PriorMdl = conjugateblm(NumPredictors) создает объект модели байесовской линейной регрессии (PriorMdl) состоит из NumPredictors предикторы и перехват, и устанавливает NumPredictors собственность. Совместное предварительное распределение (β, start2) - зависимая нормально-обратно-гамма-сопряжённая модель. PriorMdl является шаблоном, который определяет предыдущие распределения и размерность β.

пример

PriorMdl = conjugateblm(NumPredictors,Name,Value) задает свойства (кроме NumPredictors) с использованием аргументов пары имя-значение. Заключите каждое имя свойства в кавычки. Например, conjugateblm(2,'VarNames',["UnemploymentRate"; "CPI"]) задает имена двух переменных предиктора в модели.

Свойства

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

Значения свойств, доступные для записи, можно задать при создании объекта модели с помощью синтаксиса аргумента пара имя-значение или после создания объекта модели с помощью точечной нотации. Например, чтобы установить более диффузную предшествующую ковариационную матрицу для 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

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

Имена переменных предиктора для дисплеев, заданные как строковый вектор или вектор ячейки векторов символов. VarNames должен содержать NumPredictors элементы. VarNames(j) - имя переменной в столбце j набора данных предиктора, который задается во время оценки, моделирования или прогнозирования.

Значение по умолчанию: {'Beta(1)','Beta(2),...,Beta(p)}, где p - значение NumPredictors.

Пример: 'VarNames',["UnemploymentRate"; "CPI"]

Типы данных: string | cell | char

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

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

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

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

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

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

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

  • Если c является NumPredictors, то conjugateblm устанавливает предыдущую ковариационную матрицу как

    [1e50⋯00⋮V0].

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

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

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

V является предшествующей ковариацией β вплоть до множителя, равного

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

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

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

A должно быть не менее –(Intercept + NumPredictors)/2.

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

Функциональную форму обратного гамма-распределения см. в разделе Аналитически отслеживаемые апостериоры.

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

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

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

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

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

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

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

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

Примеры

свернуть все

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

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

Для всех t-временных точек αt - это ряд независимых гауссовых возмущений со средним значением 0 и дисперсией start2.

Предположим, что предыдущие распределения:

  • β'σ2∼N4 (М, λ 2В). M является вектором средства 4 на 1, а V является масштабированной матрицей положительной определенной ковариации 4 на 4.

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

Эти допущения и правдоподобие данных подразумевают нормально-обратно-гамма-сопряженную модель.

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

p = 3;
Mdl = bayeslm(p,'ModelType','conjugate')
Mdl = 
  conjugateblm 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      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 Beta(1)   |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 Beta(2)   |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 Beta(3)   |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 Sigma2    | 0.5000   0.5000    [ 0.138,  1.616]     1.000   IG(3.00,    1)        
 

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

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

Mdl.VarNames = ["IPI" "E" "WR"]
Mdl = 
  conjugateblm 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      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 IPI       |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 E         |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 WR        |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 Sigma2    | 0.5000   0.5000    [ 0.138,  1.616]     1.000   IG(3.00,    1)        
 

Рассмотрим модель линейной регрессии в «Создать нормальную-обратную-гамма-сопряженную предшествующую модель».

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

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

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

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

Оцените краевые апостериорные распределения β и start2.

PosteriorMdl = estimate(PriorMdl,X,y);
Method: Analytic posterior distributions
Number of observations: 62
Number of predictors:   4
Log marginal likelihood: -259.348
 
           |   Mean      Std          CI95        Positive       Distribution      
-----------------------------------------------------------------------------------
 Intercept | -24.2494  8.7821  [-41.514, -6.985]    0.003   t (-24.25, 8.65^2, 68) 
 IPI       |   4.3913  0.1414   [ 4.113,  4.669]    1.000   t (4.39, 0.14^2, 68)   
 E         |   0.0011  0.0003   [ 0.000,  0.002]    1.000   t (0.00, 0.00^2, 68)   
 WR        |   2.4683  0.3490   [ 1.782,  3.154]    1.000   t (2.47, 0.34^2, 68)   
 Sigma2    |  44.1347  7.8020   [31.427, 61.855]    1.000   IG(34.00, 0.00069)     
 

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

  • CI95, который содержит 95% байесовских равных достоверных интервалов для параметров. Например, апостериорная вероятность того, что коэффициент регрессии WR в [1,782, 3,154] равно 0,95.

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

  • Distribution, который содержит описания задних распределений параметров. Например, краевое заднее распределение IPI t со средним значением 4,39, стандартным отклонением 0,14 и 68 степенями свободы.

Свойства доступа к заднему распределению с помощью точечной нотации. Например, отображение краевого заднего средства путем доступа к Mu собственность.

PosteriorMdl.Mu
ans = 4×1

  -24.2494
    4.3913
    0.0011
    2.4683

Рассмотрим модель линейной регрессии в «Создать нормальную-обратную-гамма-сопряженную предшествующую модель».

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

p = 3;
PriorMdl = bayeslm(p,'ModelType','conjugate','VarNames',["IPI" "E" "WR"])
PriorMdl = 
  conjugateblm 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      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 IPI       |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 E         |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 WR        |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 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'};

Оцените условное апостериорное распределение β, учитывая данные, и start2 = 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.2494  1.8695  [-27.914, -20.585]    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.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 фиксируется, выводы о ней тривиальны.

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

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

  -24.2494
    4.3913
    0.0011
    2.4683

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

    3.4950    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 = 
  conjugateblm 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      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 IPI       |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 E         |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 WR        |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 Sigma2    | 0.5000   0.5000    [ 0.138,  1.616]     1.000   IG(3.00,    1)        
 

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

Рассмотрим модель линейной регрессии в разделе Оценка предельных задних распределений.

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

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

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

PosteriorMdl = estimate(PriorMdl,X,y);
Method: Analytic posterior distributions
Number of observations: 62
Number of predictors:   4
Log marginal likelihood: -259.348
 
           |   Mean      Std          CI95        Positive       Distribution      
-----------------------------------------------------------------------------------
 Intercept | -24.2494  8.7821  [-41.514, -6.985]    0.003   t (-24.25, 8.65^2, 68) 
 IPI       |   4.3913  0.1414   [ 4.113,  4.669]    1.000   t (4.39, 0.14^2, 68)   
 E         |   0.0011  0.0003   [ 0.000,  0.002]    1.000   t (0.00, 0.00^2, 68)   
 WR        |   2.4683  0.3490   [ 1.782,  3.154]    1.000   t (2.47, 0.34^2, 68)   
 Sigma2    |  44.1347  7.8020   [31.427, 61.855]    1.000   IG(34.00, 0.00069)     
 

Извлеките заднее среднее β из задней модели и заднюю ковариацию β из сводки оценки, возвращенной summarize.

estBeta = PosteriorMdl.Mu;
Summary = summarize(PosteriorMdl);
estBetaCov = Summary.Covariances{1:(end - 1),1:(end - 1)};

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

Нарисовать 1e6 образцы из краевого заднего распределения β.

NumDraws = 1e6;
rng(1);
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.5362

Апостериорная вероятность того, что коэффициент реальной заработной платы меньше 2,5, составляет около 0.54.

Краевое заднее распределение коэффициента WR является t68, но центрирован в 2,47 и масштабирован на 0,34. Непосредственно вычислить апостериорную вероятность того, что коэффициент WR менее 2,5.

center = estBeta(isWR);
stdBeta = sqrt(diag(estBetaCov));
scale = stdBeta(isWR);
t = (2.5 - center)/scale;
dof = 68;
directProb = tcdf(t,dof)
directProb = 0.5361

Задние вероятности почти идентичны.

Copyright 2018 The MathWorks, Inc.

Рассмотрим модель линейной регрессии в разделе Оценка предельных задних распределений.

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

p = 3;
PriorMdl = bayeslm(p,'ModelType','conjugate','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

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');
ylabel('rGNP');
xlabel('Year');
hold off

Figure contains an axes. The axes with title Real Gross National Product contains 3 objects of type patch, line. These objects represent Forecast Horizon, True GNPR, Forecasted GNPR.

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

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

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

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

Copyright 2018 The MathWorks, Inc.

Подробнее

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

Алгоритмы

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

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

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

Представлен в R2017a