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

Рассмотрим стандартную задачу линейной регрессии, в которой для ${\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} }$ называется гребневой регрессией.

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

• 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