Байесовская линейная регрессия

Байесовская статистика
Теория
Приемлемое правило решения Эффективность по Байесу Байесовская вероятность Интерпретация вероятности Теорема Байеса / Коэффициент Байеса Байесовский вывод / Байесовская сеть Априорная вероятность Апостериорная вероятность Правдоподобие Сопряжённое априорное распределение Апостериорное прогнозируемое распределение Гиперпараметр / Гиперприорное распределение Принцип безразличия Максимум энтропии Эмпирический метод Правило Кромвеля Теорема Бернштейна — фон Мизеса Критерий Шварца Доверительный интервал Оценка апостериорного максимума Радикальный пробабилизм
Техники
Байесовская линейная регрессия Байесовская оценка решения Приближенные байесовские вычисления Марковская цепь Монте-Карло
Портал:Статистика

Байесовская линейная регрессия — это подход в линейной регрессии, в котором статистический анализ проводится в контексте байесовского вывода: когда регрессионная модель имеет ошибки^[англ.], имеющие нормальное распределение, и, если принимается определённая форма априорного распределения, доступны явные результаты для апостериорных распределений вероятностей параметров модели.

Рассмотрим стандартную задачу линейной регрессии, в которой для ${\displaystyle i=1,...,n}$ мы указываем среднее условное распределение величины ${\displaystyle y_{i}}$ для заданного вектора ${\displaystyle k\times 1}$ предсказаний ${\displaystyle \mathbf {x} _{i}}$:

${\displaystyle y_{i}=\mathbf {x} _{i}^{\rm {T}}{\boldsymbol {\beta }}+\epsilon _{i},}$

где ${\displaystyle {\boldsymbol {\beta }}}$ является ${\displaystyle k\times 1}$ вектором, а ${\displaystyle \epsilon _{i}}$ являются независимыми и одинаково распределёнными нормально случайными величинами:

${\displaystyle \epsilon _{i}\sim N(0,\sigma ^{2}).}$

Это соответствует следующей функции правдоподобия:

${\displaystyle \rho (\mathbf {y} |\mathbf {X} ,{\boldsymbol {\beta }},\sigma ^{2})\propto (\sigma ^{2})^{-n/2}e^{-{\frac {1}{2{\sigma }^{2}}}(\mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})^{\rm {T}}(\mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})}.}$

Решение обычного метода наименьших квадратов является оценкой вектора коэффициентов с помощью псевдообратной матрицы Мура — Пенроуза:

${\displaystyle {\hat {\boldsymbol {\beta }}}=(\mathbf {X} ^{\rm {T}}\mathbf {X} )^{-1}\mathbf {X} ^{\rm {T}}\mathbf {y} }$

где ${\displaystyle \mathbf {X} }$ является ${\displaystyle n\times k}$ матрицей плана^[англ.], каждая строка которой является вектором предсказаний ${\displaystyle \mathbf {x} _{i}^{\rm {T}}}$, а ${\displaystyle \mathbf {y} }$ является вектор-столбцом r ${\displaystyle [y_{1}\;\cdots \;y_{n}]^{\rm {T}}}$.

Это является частотным^[англ.] подходом, и предполагается, что существует достаточно измерений для того, чтобы сказать что-то осмысленное о ${\displaystyle {\boldsymbol {\beta }}}$. В байесовском подходе данные сопровождаются дополнительной информацией в виде априорного распределения вероятности. Априорные убеждения о параметрах комбинируются с функцией правдоподобия данных согласно теореме Байеса для получения апостериорной уверенности о параметрах ${\displaystyle {\boldsymbol {\beta }}}$ и ${\displaystyle \sigma }$. Априорные данные могут принимать различные формы в зависимости от области применения и информации, которая доступна a priori.

Для любого априорного распределения, может не существовать аналитического решения для апостериорного распределения. В этом разделе мы рассмотрим так называемое сопряжённое априорное распределение, для которого апостериорное распределение можно вывести аналитически.

Априорное распределение ${\displaystyle \rho ({\boldsymbol {\beta }},\sigma ^{2})}$ является сопряжённым функции правдоподобия, если оно имеет ту же функциональную форму с учётом ${\displaystyle {\boldsymbol {\beta }}}$ и ${\displaystyle \sigma }$. Поскольку логарифмическое правдоподобие квадратично от ${\displaystyle {\boldsymbol {\beta }}}$, его перепишем так, что правдоподобие становится нормальным от ${\displaystyle ({\boldsymbol {\beta }}-{\hat {\boldsymbol {\beta }}})}$. Запишем

${\displaystyle {\begin{aligned}(\mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})^{\rm {T}}(\mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})&=(\mathbf {y} -\mathbf {X} {\hat {\boldsymbol {\beta }}})^{\rm {T}}(\mathbf {y} -\mathbf {X} {\hat {\boldsymbol {\beta }}})\\&+({\boldsymbol {\beta }}-{\hat {\boldsymbol {\beta }}})^{\rm {T}}(\mathbf {X} ^{\rm {T}}\mathbf {X} )({\boldsymbol {\beta }}-{\hat {\boldsymbol {\beta }}}).\end{aligned}}}$

Правдоподобие теперь переписывается как

${\displaystyle {\begin{aligned}\rho (\mathbf {y} |\mathbf {X} ,{\boldsymbol {\beta }},\sigma ^{2})&\propto (\sigma ^{2})^{-v/2}e^{-{\frac {vs^{2}}{2{\sigma }^{2}}}}(\sigma ^{2})^{-(n-v)/2}\\&\times e^{-{\frac {1}{2{\sigma }^{2}}}({\boldsymbol {\beta }}-{\hat {\boldsymbol {\beta }}})^{\rm {T}}(\mathbf {X} ^{\rm {T}}\mathbf {X} )({\boldsymbol {\beta }}-{\hat {\boldsymbol {\beta }}})},\end{aligned}}}$

где

${\displaystyle vs^{2}=(\mathbf {y} -\mathbf {X} {\hat {\boldsymbol {\beta }}})^{\rm {T}}(\mathbf {y} -\mathbf {X} {\hat {\boldsymbol {\beta }}})\quad }$ и ${\displaystyle \quad v=n-k}$,

где ${\displaystyle k}$ является числом коэффициентов регрессии.

Это указывает на вид априорного распределения:

${\displaystyle \rho ({\boldsymbol {\beta }},\sigma ^{2})=\rho (\sigma ^{2})\rho ({\boldsymbol {\beta }}|\sigma ^{2}),}$

где ${\displaystyle \rho (\sigma ^{2})}$ является обратным гамма-распределением^[англ.]

${\displaystyle \rho (\sigma ^{2})\propto (\sigma ^{2})^{-{\frac {v_{0}}{2}}-1}e^{-{\frac {v_{0}s_{0}^{2}}{2{\sigma }^{2}}}}.}$

В обозначениях, введённых в статье Обратное гамма-распределение^[англ.], это плотность распределения ${\displaystyle {\text{Inv-Gamma}}(a_{0},b_{0})}$ с ${\displaystyle a_{0}={\tfrac {v_{0}}{2}}}$ и ${\displaystyle b_{0}={\tfrac {1}{2}}v_{0}s_{0}^{2}}$, где ${\displaystyle v_{0}}$ и ${\displaystyle s_{0}^{2}}$ являются априорными значениями ${\displaystyle v}$ и ${\displaystyle s^{2}}$ соответственно. Эквивалентно, эту плотность можно описать как масштабированное обратное распределение хи-квадрат^[англ.] ${\displaystyle {\mbox{Scale-inv-}}\chi ^{2}(v_{0},s_{0}^{2}).}$

Далее, условная априорная плотность ${\displaystyle \rho ({\boldsymbol {\beta }}|\sigma ^{2})}$ является нормальным распределением,

${\displaystyle \rho ({\boldsymbol {\beta }}|\sigma ^{2})\propto (\sigma ^{2})^{-{\frac {k}{2}}}e^{-{\frac {1}{2{\sigma }^{2}}}({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{0})^{\rm {T}}\mathbf {\Lambda } _{0}({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{0})}.}$

В обозначениях нормального распределения условное априорное распределение равно ${\displaystyle {\mathcal {N}}\left({\boldsymbol {\mu }}_{0},\sigma ^{2}\mathbf {\Lambda } _{0}^{-1}\right).}$

При указанном априорным распределении апостериорное распределение можно выразить как

${\displaystyle \rho ({\boldsymbol {\beta }},\sigma ^{2}|\mathbf {y} ,\mathbf {X} )\propto \rho (\mathbf {y} |\mathbf {X} ,{\boldsymbol {\beta }},\sigma ^{2})\rho ({\boldsymbol {\beta }}|\sigma ^{2})\rho (\sigma ^{2})}$

${\displaystyle \propto (\sigma ^{2})^{-n/2}e^{-{\frac {1}{2{\sigma }^{2}}}(\mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})^{\rm {T}}(\mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})}}$

${\displaystyle \times (\sigma ^{2})^{-k/2}e^{-{\frac {1}{2{\sigma }^{2}}}({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{0})^{\rm {T}}{\boldsymbol {\Lambda }}_{0}({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{0})}}$

${\displaystyle \times (\sigma ^{2})^{-(a_{0}+1)}e^{-{\frac {b_{0}}{{\sigma }^{2}}}}.}$

После некоторых преобразований^[1] апостериорная вероятность может быть переписана так, что апостериорное среднее ${\displaystyle {\boldsymbol {\mu }}_{n}}$ вектора параметров ${\displaystyle {\boldsymbol {\beta }}}$ может быть выражено в терминах оценки по методу наименьших квадратов ${\displaystyle {\hat {\boldsymbol {\beta }}}}$ и априорного среднего ${\displaystyle {\boldsymbol {\mu }}_{0}}$, где поддержка априорной вероятности выражается матрицей априорной точности ${\displaystyle {\boldsymbol {\Lambda }}_{0}}$

${\displaystyle {\boldsymbol {\mu }}_{n}=(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0})^{-1}(\mathbf {X} ^{\rm {T}}\mathbf {X} {\hat {\boldsymbol {\beta }}}+{\boldsymbol {\Lambda }}_{0}{\boldsymbol {\mu }}_{0}).}$

Для подтверждения, что ${\displaystyle {\boldsymbol {\mu }}_{n}}$ в действительности является апостериорным средним, квадратичные члены в экспоненте можно преобразовать к квадратичной форме^[англ.] от ${\displaystyle {\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{n}}$^[2].

${\displaystyle (\mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})^{\rm {T}}(\mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})+({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{0})^{\rm {T}}{\boldsymbol {\Lambda }}_{0}({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{0})=}$

${\displaystyle ({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{n})^{\rm {T}}(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0})({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{n})+\mathbf {y} ^{\rm {T}}\mathbf {y} -{\boldsymbol {\mu }}_{n}^{\rm {T}}(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0}){\boldsymbol {\mu }}_{n}+{\boldsymbol {\mu }}_{0}^{\rm {T}}{\boldsymbol {\Lambda }}_{0}{\boldsymbol {\mu }}_{0}.}$

Теперь апостериорное распределение можно выразить как нормальное распределение, умноженное на обратное гамма-распределение^[англ.]:

${\displaystyle \rho ({\boldsymbol {\beta }},\sigma ^{2}|\mathbf {y} ,\mathbf {X} )\propto (\sigma ^{2})^{-{\frac {k}{2}}}e^{-{\frac {1}{2{\sigma }^{2}}}({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{n})^{\rm {T}}(\mathbf {X} ^{\rm {T}}\mathbf {X} +\mathbf {\Lambda } _{0})({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{n})}}$

${\displaystyle \times (\sigma ^{2})^{-{\frac {n+2a_{0}}{2}}-1}e^{-{\frac {2b_{0}+\mathbf {y} ^{\rm {T}}\mathbf {y} -{\boldsymbol {\mu }}_{n}^{\rm {T}}(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0}){\boldsymbol {\mu }}_{n}+{\boldsymbol {\mu }}_{0}^{\rm {T}}{\boldsymbol {\Lambda }}_{0}{\boldsymbol {\mu }}_{0}}{2{\sigma }^{2}}}}.}$

Поэтому апостериорное распределение можно параметризовать следующим образом.

${\displaystyle \rho ({\boldsymbol {\beta }},\sigma ^{2}|\mathbf {y} ,\mathbf {X} )\propto \rho ({\boldsymbol {\beta }}|\sigma ^{2},\mathbf {y} ,\mathbf {X} )\rho (\sigma ^{2}|\mathbf {y} ,\mathbf {X} ),}$

где два множителя соответствуют плотностям распределений ${\displaystyle {\mathcal {N}}\left({\boldsymbol {\mu }}_{n},\sigma ^{2}{\boldsymbol {\Lambda }}_{n}^{-1}\right)\,}$ и ${\displaystyle {\text{Inv-Gamma}}\left(a_{n},b_{n}\right)}$ с параметрами, задаваемыми выражениями

${\displaystyle {\boldsymbol {\Lambda }}_{n}=(\mathbf {X} ^{\rm {T}}\mathbf {X} +\mathbf {\Lambda } _{0}),\quad {\boldsymbol {\mu }}_{n}=({\boldsymbol {\Lambda }}_{n})^{-1}(\mathbf {X} ^{\rm {T}}\mathbf {X} {\hat {\boldsymbol {\beta }}}+{\boldsymbol {\Lambda }}_{0}{\boldsymbol {\mu }}_{0}),}$

${\displaystyle a_{n}=a_{0}+{\frac {n}{2}},\qquad b_{n}=b_{0}+{\frac {1}{2}}(\mathbf {y} ^{\rm {T}}\mathbf {y} +{\boldsymbol {\mu }}_{0}^{\rm {T}}{\boldsymbol {\Lambda }}_{0}{\boldsymbol {\mu }}_{0}-{\boldsymbol {\mu }}_{n}^{\rm {T}}{\boldsymbol {\Lambda }}_{n}{\boldsymbol {\mu }}_{n}).}$

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

${\displaystyle {\boldsymbol {\mu }}_{n}=(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0})^{-1}({\boldsymbol {\Lambda }}_{0}{\boldsymbol {\mu }}_{0}+\mathbf {X} ^{\rm {T}}\mathbf {X} {\hat {\boldsymbol {\beta }}})=(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0})^{-1}({\boldsymbol {\Lambda }}_{0}{\boldsymbol {\mu }}_{0}+\mathbf {X} ^{\rm {T}}\mathbf {y} ),}$

${\displaystyle {\boldsymbol {\Lambda }}_{n}=(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0}),}$

${\displaystyle a_{n}=a_{0}+{\frac {n}{2}},}$

${\displaystyle b_{n}=b_{0}+{\frac {1}{2}}(\mathbf {y} ^{\rm {T}}\mathbf {y} +{\boldsymbol {\mu }}_{0}^{\rm {T}}{\boldsymbol {\Lambda }}_{0}{\boldsymbol {\mu }}_{0}-{\boldsymbol {\mu }}_{n}^{\rm {T}}{\boldsymbol {\Lambda }}_{n}{\boldsymbol {\mu }}_{n}).}$

Обоснованность модели ${\displaystyle p(\mathbf {y} |m)}$ — это вероятность данных для данной модели ${\displaystyle m}$. Она известна также как предельное правдоподобие и как априорная предсказательная плотность. Здесь модель определяется функцией правдоподобия ${\displaystyle p(\mathbf {y} |\mathbf {X} ,{\boldsymbol {\beta }},\sigma )}$ и априорным распределением параметров, то есть, ${\displaystyle p({\boldsymbol {\beta }},\sigma )}$. Обоснованность модели фиксируется одним числом, показывающим, насколько хорошо такая модель объясняет наблюдения. Обоснованность модели байесовской линейной регрессии, представленная в этом разделе, может быть использована для сравнения конкурирующих линейных моделей путём байесовского сравнения моделей. Эти модели могут отличаться числом и значениями предсказывающих переменных, как и их априорными значениями в параметрах модели. Сложность модели принимается во внимание обоснованностью модели, поскольку она исключает параметры путём интегрирования ${\displaystyle p(\mathbf {y} ,{\boldsymbol {\beta }},\sigma |\mathbf {X} )}$ по всем возможным значениям ${\displaystyle {\boldsymbol {\beta }}}$ и ${\displaystyle \sigma }$.

${\displaystyle p(\mathbf {y} |m)=\int p(\mathbf {y} |\mathbf {X} ,{\boldsymbol {\beta }},\sigma )\,p({\boldsymbol {\beta }},\sigma )\,d{\boldsymbol {\beta }}\,d\sigma }$

Этот интеграл можно вычислить аналитически и решение задаётся следующим равенством^[3]

${\displaystyle p(\mathbf {y} |m)={\frac {1}{(2\pi )^{n/2}}}{\sqrt {\frac {\det({\boldsymbol {\Lambda }}_{0})}{\det({\boldsymbol {\Lambda }}_{n})}}}\cdot {\frac {b_{0}^{a_{0}}}{b_{n}^{a_{n}}}}\cdot {\frac {\Gamma (a_{n})}{\Gamma (a_{0})}}}$

Здесь ${\displaystyle \Gamma }$ означает гамма-функцию. Поскольку мы выбрали сопряжённое априорное распределение, предельное правдоподобие может быть легко вычислено путём решения следующего равенства для произвольных значений ${\displaystyle {\boldsymbol {\beta }}}$ и ${\displaystyle \sigma }$.

${\displaystyle p(\mathbf {y} |m)={\frac {p({\boldsymbol {\beta }},\sigma |m)\,p(\mathbf {y} |\mathbf {X} ,{\boldsymbol {\beta }},\sigma ,m)}{p({\boldsymbol {\beta }},\sigma |\mathbf {y} ,\mathbf {X} ,m)}}}$

Заметим, что это равенство является ни чем иным, как переформулировкой теоремы Байеса. Подстановка формулы для априорной вероятности, правдоподобия и апостериорной вероятности и упрощения получающегося выражения приводит к аналитическому выражению, приведённому выше.

В общем случае может оказаться невозможным или нецелесообразным получать апостериорное распределение аналитически. Однако можно аппроксимировать апостериорную вероятность методом приближенного байесовского вывода^[англ.], таким как выборка по методу Монте-Карло^[4] или вариационные байесовские методы^[англ.].

Частный случай ${\displaystyle {\boldsymbol {\mu }}_{0}=0,\mathbf {\Lambda } _{0}=c\mathbf {E} }$ называется гребневой регрессией.

Аналогичный анализ можно провести для общего случая множественной регрессии и частично для байесовской оценки ковариационной матрицы^[англ.] — см. Байесовская мультивариантная линейная регрессия^[англ.].

• Байесовская линейная статистика^[англ.]
• Регуляризованная линейной регрессия^[англ.]
• Метод регуляризации Тихонова
• Выбор спайковой переменной^[англ.]
• Байесовская интерпретация ядерной регуляризации^[англ.]

• George E. P. Box, Tiao G. C. Bayesian Inference in Statistical Analysis. — Wiley, 1973. — ISBN 0-471-57428-7.
• Bradley P. Carlin, Thomas A. Louis. Bayesian Methods for Data Analysis, Third Edition. — Boca Raton, FL: Chapman and Hall/CRC, 2008. — ISBN 1-58488-697-8.
• Fahrmeir L., Kneib T., Lang S. Regression. Modelle, Methoden und Anwendungen. — 2nd. — Heidelberg: Springer, 2009. — ISBN 978-3-642-01836-7. — doi:10.1007/978-3-642-01837-4.
• Fornalski K.W., Parzych G., Pylak M., Satuła D., Dobrzyński L. Application of Bayesian reasoning and the Maximum Entropy Method to some reconstruction problems // Acta Physica Polonica A. — 2010. — Т. 117, вып. 6. — С. 892—899. — doi:10.12693/APhysPolA.117.892.
• Krzysztof W. Fornalski. Applications of the robust Bayesian regression analysis // International Journal of Society Systems Science. — 2015. — Т. 7, вып. 4. — С. 314–333. — doi:10.1504/IJSSS.2015.073223.
• Andrew Gelman, John B. Carlin, Hal S. Stern, Donald B. Rubin. Bayesian Data Analysis, Second Edition. — Boca Raton, FL: Chapman and Hall/CRC, 2003. — ISBN 1-58488-388-X.
• Michael Goldstein, David Wooff. Bayes Linear Statistics, Theory & Methods. — Wiley, 2007. — ISBN 978-0-470-01562-9.
• Minka, Thomas P. (2001) Bayesian Linear Regression Архивная копия от 26 октября 2008 на Wayback Machine, Microsoft research web page
• Peter E. Rossi, Greg M. Allenby, Robert McCulloch. Bayesian Statistics and Marketing. — John Wiley & Sons, 2006. — ISBN 0470863676.
• Anthony O'Hagan. Bayesian Inference. — First. — Halsted, 1994. — Т. 2B. — (Kendall's Advanced Theory of Statistics). — ISBN 0-340-52922-9.
• Sivia, D.S., Skilling, J. Data Analysis - A Bayesian Tutorial. — Second. — Oxford University Press, 2006.
• Gero Walter, Thomas Augustin. Bayesian Linear Regression—Different Conjugate Models and Their (In)Sensitivity to Prior-Data Conflict // Technical Report Number 069, Department of Statistics, University of Munich. — 2009.

• Python
  • Bayesian Type-II Linear Regression code, tutorial Архивная копия от 18 декабря 2020 на Wayback Machine
  • ARD Linear Regression code Архивная копия от 1 марта 2017 на Wayback Machine
  • ARD Linear Regression with kernelized features code Архивная копия от 1 марта 2017 на Wayback Machine, tutorial Архивная копия от 18 декабря 2020 на Wayback Machine