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

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

Конфигурация модели

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

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

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

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

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

$\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 }})}.$

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

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

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

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

Регрессия с сопряжёнными распределениями

Сопряжённое априорное распределение

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

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

${\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}}$

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

${\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}}$

где

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

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

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

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

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

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

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

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

$\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})}.$

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

Апостериорное распределение

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

$\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})$

$\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 }})}$

$\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})}$

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

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

${\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}).$

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

$(\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})=$

$({\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}.$

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

$\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})}$

$\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}}}}.$

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

$\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} ),$

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

${\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}),$

$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}).$

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

${\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} ),$

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

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

$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}).$

Обоснованность модели

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

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

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

$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})}}$

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

$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)}}$

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

Другие случаи

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

Частный случай ${\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