贝叶斯线性回归（Bayesian linear regression）是使用统计学中贝叶斯推断（Bayesian inference）方法求解的线性回归（linear regression）模型。
贝叶斯线性回归将线性模型的参数视为随机变量（random variable），并通过模型参数（权重系数）的先验（prior）计算其后验（posterior）。贝叶斯线性回归可以使用数值方法求解，在一定条件下，也可得到解析型式的后验或其有关统计量。
贝叶斯线性回归具有贝叶斯统计模型的基本性质，可以求解权重系数的概率密度函数，进行在线学习以及基于贝叶斯因子（Bayes factor）的模型假设检验。
给定相互独立的 $N$ 组学习样本 $\boldsymbol{X}=\left\{\boldsymbol{X}_1, \boldsymbol{X}_2, \ldots, \boldsymbol{X}_N\right\} \in \mathbb{R}^N ， \boldsymbol{y}=\left\{\boldsymbol{y}_1, \boldsymbol{y}_2, \ldots, \boldsymbol{y}_N\right\}$ ，贝叶斯线性回归使用如下的多元线性回归模型 $f(\boldsymbol{X})$ : $$ f(\boldsymbol{X})=\boldsymbol{X}^{\top} \boldsymbol{w}, \quad \boldsymbol{y}=f(\boldsymbol{X})+\epsilon $$ 这里 $\boldsymbol{w}$ 为权重系数， $\epsilon$ 为残差。由于学习样本间相互独立，因此 $\epsilon$ 为独立同分布 (independent and identically distributed, iid）。贝叶斯线性回归假设残差服从正态分布，其方差服从逆Gamma分布 (Inverse-Gamma distribution) : $$ \begin{aligned} &p(\epsilon)=N\left(\epsilon \mid \mu_n, \sigma_n^2\right) \\ &\sigma_n^2=\operatorname{Inv}-\operatorname{Gamma}\left(\sigma_n^2 \mid a, b\right) \end{aligned} $$ 正态分布的均值 $\mu_n$ 和逆Gamma分布的系数 $(a, b)$ 要求预先指定。通常设定均值 $\mu_n=0$ ，对应白噪声 (white noise) 残差， 因此贝叶斯线性回归的模型本身至少包含2个超参数 (hyperparameter) 。以上的贝叶斯线性回归也可推广至广义线性模型 (Generalized Linear Model, GLM)，得到贝叶斯广义线性模型（Bayesian GLM) 。
根据线性模型的定义，权重系数 $\boldsymbol{w}$ 与观测数据 $\boldsymbol{X}$ 相互独立，也与残差的方差 $\sigma_n^2$ 相互独立，由贝叶斯定理 (Bayes' theorem）可推出，贝叶斯线性回归中权重系数的后验有如下表示: $$ p\left(\boldsymbol{w} \mid \boldsymbol{X}, \boldsymbol{y}, \sigma_n^2\right)=\frac{p\left(\boldsymbol{y} \mid \boldsymbol{X}, \boldsymbol{w}, \sigma_n^2\right) p(\boldsymbol{w})}{p\left(\boldsymbol{y} \mid \boldsymbol{X}, \sigma_n^2\right)} \quad(*) $$ 式中 $p\left(\boldsymbol{y} \mid \boldsymbol{X}, \boldsymbol{w}, \sigma_n^2\right)$ 称为似然（likelihood)，由线性回归模型完全决定，以模型残差服从iid的0均值正态分布为例，这里似然 也服从iid的正态分布 :
$$ p\left(\boldsymbol{y} \mid \boldsymbol{X}, \boldsymbol{w}, \sigma_n^2\right)=N\left(\boldsymbol{y} \mid \boldsymbol{X}^{\top} \boldsymbol{w}+0, \sigma_n^2\right)=\prod_{i=1}^N p\left(\boldsymbol{y}_i \mid \boldsymbol{X}_i, \boldsymbol{w}\right) $$$$ =\prod_{i=1}^N \frac{1}{\sqrt{2 \pi} \sigma_n} \exp \left[\frac{\left(\boldsymbol{y}_i-\boldsymbol{X}_i^{\top} \boldsymbol{w}\right)^2}{2 \sigma_n^2}\right]=\frac{1}{\sqrt{2 \pi} \sigma_n} \exp \left(-\frac{\left|\boldsymbol{y}-\boldsymbol{X}^{\top} \boldsymbol{w}\right|^2}{2 \sigma_n^2}\right) $$$(*)$ 式中的 $p\left(\boldsymbol{y} \mid \boldsymbol{X}, \sigma_n^2\right)$ 是 $\boldsymbol{y}$ 的边缘似然 (marginal likelihood)，在贝叶斯推断中也被称为模型证据 (model evidence)， 仅与观测数据有关，与权重系数相互独立 。 求解 $(*)$ 式要求预先给定权重系数的先验 $p(\boldsymbol{w})$ ，即一个连续概率分布（continuous probability distribution)，通常的选择 为 0 均值的正态分布 :
$$ p(\boldsymbol{w})=N\left(\boldsymbol{w} \mid 0, \sigma_w^2\right)=\frac{1}{\sqrt{2 \pi} \sigma_w} \exp \left(-\frac{\boldsymbol{w}^{\top} \boldsymbol{w}}{2 \sigma_w^2}\right) $$式中 $\sigma_w^2$ 为预先给定的超参数。和其它贝叶斯推断一样，根据似然和先验的类型，可用于 求解贝叶斯线性回归的方法有三类，即极大后验估计（Maximum A Posteriori estimation, MAP）、共轭先验 (conjugate prior) 求解和数值方法，前两者要求似然或先验满足特定 条件，数值方法没有特定要求，可以通过迭代逼近任意形式的后验分布。 在权重系数没有合理先验假设的问题中，贝叶斯线性回归可使用无信息先验 (uninformative prior)，即一个均匀分布（uniform distribution)，此时权重系数按均等 的机会取任意值 。
在贝叶斯线性回归中，MAP可以被视为一种特殊的贝叶斯估计（Bayesian estimator)，其求解步骤与极大似然估计 (maximum likelihood estimation) 类似。对给定的先验，MAP将 $(*)$ 式转化为求解 $\boldsymbol{w}$ 使后验概率最大的优化问题，并求得后验的众数 (mode)。由于正态分布的众数即是均值，因此MAP通常被应用于正态先验。
这里以的0均值正态先验为例介绍MAP的求解步骤，首先给定权重系数 $\boldsymbol{w}$ 的0均值正态分布先验: $p(\boldsymbol{w})=N\left(\boldsymbol{w} \mid 0, \sigma_w^2\right)$ 由于边缘似然与 $\boldsymbol{w}$ 相互独立，此时求解后验概率的极大值等价于求解似然和先验乘积的极大值 ${ }^{[3]}$ : $$ \arg \max _{\boldsymbol{w}} p\left(\boldsymbol{w} \mid \boldsymbol{X}, \boldsymbol{y}, \sigma_n^2\right) \Longleftrightarrow \arg \max _{\boldsymbol{w}} p\left(\boldsymbol{y} \mid \boldsymbol{X}, \boldsymbol{w}, \sigma_n^2\right) p(\boldsymbol{w}) $$ 对上述极值问题取自然对数并考虑正态分布的解析型式，可有如下推导:
$$ \begin{aligned} \arg \max _{\boldsymbol{w}} & \log \left[p\left(\boldsymbol{y} \mid \boldsymbol{X}, \boldsymbol{w}, \sigma_n^2\right)\right]+\log [p(\boldsymbol{w})] \\ \log \left[p\left(\boldsymbol{y} \mid \boldsymbol{X}, \boldsymbol{w}, \sigma_n^2\right)\right]+\log [p(\boldsymbol{w})] &=\log \left[\frac{1}{\sqrt{2 \pi} \sigma_n} \exp \left(-\frac{\left|\boldsymbol{y}-\boldsymbol{X}^{\top} \boldsymbol{w}\right|^2}{2 \sigma_n^2}\right)\right]+\log \left[\frac{1}{\sqrt{2 \pi} \sigma_w} \exp \left(-\frac{\boldsymbol{w}^{\top} \boldsymbol{w}}{2 \sigma_w^2}\right)\right] \\ &=-\frac{1}{2 \sigma_n^2}\left(\boldsymbol{y}-\boldsymbol{X}^{\top} \boldsymbol{w}\right)^{\top}\left(\boldsymbol{y}-\boldsymbol{X}^{\top} \boldsymbol{w}\right)-\frac{1}{2 \sigma_w^2} \boldsymbol{w}^{\top} \boldsymbol{w} \\ &=-\frac{1}{2 \sigma_n^2}\left\|\boldsymbol{y}-\boldsymbol{X}^{\top} \boldsymbol{w}\right\|^2-\frac{1}{2 \sigma_w^2}\|\boldsymbol{w}\|^2 \end{aligned} $$由于式中各项的系数均为负数，因此该极大值问题可转化为仅与 $\boldsymbol{w}$ 有关的极小值问题，并可通过线性代数得到 $\boldsymbol{w}$ 的解 $[14]$ : $$ \begin{aligned} &\arg \min _{\boldsymbol{w}}\left\|\boldsymbol{y}-\boldsymbol{X}^{\top} \boldsymbol{w}\right\|^2+\lambda\|\boldsymbol{w}\|^2, \quad \lambda=\frac{\sigma_n^2}{\sigma_w^2} \\ &\Rightarrow \boldsymbol{w}=\left(\boldsymbol{X}^{\top} \boldsymbol{X}+\lambda \boldsymbol{I}\right)^{-1} \boldsymbol{X}^{\top} \boldsymbol{y} \end{aligned} $$ 式中 $\lambda$ 为模型残差和权重系数方差的比值，由超参数直接计算， $\boldsymbol{I}$ 为单位矩阵。 除正态先验外，MAP也使用拉普拉斯分布（Laplace distribution）作为权重系数的先验: $$ p(\boldsymbol{w} \mid a, b)=\text { Laplace }(a, b)=\frac{1}{(2 b)^n} \prod_{i=1}^n \exp \left(-\frac{\left|w_i-a\right|}{b}\right) $$ 式中 $(a, b)$ 为超参数。取均值 $a=0$ ，在经过与正态先验类似的推导后，拉普拉斯先验的MAP可转化为如下优化问题 : $$ \arg \min _{\boldsymbol{w}}\left\|\boldsymbol{y}-\boldsymbol{X}^{\top} \boldsymbol{w}\right\|^2+\lambda\|\boldsymbol{w}\|, \quad \lambda=\frac{\sigma_n^2}{b} $$ 该优化问题对应一个泛定方程，在 $\lambda$ 足够大时，可以得到稀疏解（sparse solution）。 MAP是单点估计，不支持在线学习（online learning），也无法提供置信区间，但能够以很小的计算量求解贝叶斯线性回归 的权重系数，且可用于最常见的正态先验情形 $[3]$ 。使用正态先验和MAP求解的贝叶斯线性回归等价于岭回归 (ridge regression），优化目标中的 $\lambda\|\boldsymbol{w}\|^2$ 被称为 $L_2$ 正则化项；使用拉普拉斯分布先验的情形对应线性模型的LASSO（Least Absolute Shrinkage and Selection Operator)， $\lambda\|\boldsymbol{w}\|$ 被称为 $L_1$ 正则化项 。
一般地，贝叶斯推断的数值方法都适用于贝叶斯线性回归，其中最常见的是马尔可夫链蒙特卡罗 (Markov Chain Monte Carlo, MCMC)。这里以MCMC中的吉布斯采样 (Gibbs sampling) 算法为例介绍。 给定均值为 $\mu$ ，方差为 $\sigma_w^2$ 的正态先验 $p(\boldsymbol{w})=N\left(\boldsymbol{w} \mid \mu, \sigma_w^2\right)$ 和权重系数 $\boldsymbol{w}=\left\{w_1, w_2, \ldots, w_m\right\}$ ，吉布斯采样是一 个迭代算法，每个迭代都依次采样所有的权重系数 :
迭代步骤中的 $N\left(w_j \mid \boldsymbol{X}, \boldsymbol{y}, w_1, \ldots, w_{j-1}, w_{j+1}, \ldots, w_m, \sigma_n^2, \mu_1, \sigma_{w_1}^2\right)$ 和 $\sigma_n^2(i+1)=\operatorname{Inv}-\operatorname{Gamma}\left[\sigma_n^2(i+1) \mid \boldsymbol{X}, \boldsymbol{y}, w^{i+1}, a, b\right]$ 被称为条件采样分布 (conditional sampling distribution)。条 件采样的推导较为繁琐，这里给出一维情形下截距 $w_1$ 、斜率 $w_2$ 和残差方差 $\sigma_n^2$ 的条件采样分布: $$ \begin{aligned} &N\left(w_1 \mid X, y, w_2, \sigma_n^2, \mu_1, \sigma_{w_1}^2\right)=N\left[w_1 \mid \frac{\mu_1 \sigma_n^2+\sigma_{w_1}^2 \sum_{i=1}^N\left(y_i-w_2 X_i\right)}{\sigma_n^2+\sigma_{w_1}^2 N}, \frac{\sigma_{w_1}^2 \sigma_n^2}{\sigma_n^2+\sigma_{w_1}^2 N}\right] \\ &N\left(w_2 \mid X, y, w_1, \sigma_n^2, \mu_2, \sigma_{w_2}^2\right)=N\left[w_2 \mid \frac{\mu_2 \sigma_n^2+\sigma_{w_2}^2 \sum_{i=1}^N\left(y_i X_i-w_1 X_i\right)}{\sigma_n^2+\sigma_{w_2}^2 \sum_{i=1}^N X_i^2}, \frac{\sigma_{w_2}^2 \sigma_n^2}{\sigma_n^2+\sigma_{w_2}^2 \sum_{i=1}^N X_i^2}\right] \\ &\operatorname{Inv}-\operatorname{Gamma}\left(\sigma_n^2 \mid X, y, w, a, b\right)=\operatorname{Inv}-\operatorname{Gamma}\left(\sigma_n^2\left|a+\frac{N}{2}, b+\frac{1}{2} \sum_{i=1}^N\right| y_i-w_1-\left.w_2 X_i\right|^2\right) \end{aligned} $$
由MAP求解的贝叶斯线性回归可直接使用权重系数对预测数据进行计算。MAP是单点估计，因此预测结果不提供后验分布。 对共轭先验或数值方法求解的贝叶斯线性回归，可通过边缘化 (marginalized out) 模型权重，即按其后验积分得到预测结果 $[2]$ : $$ p\left(\boldsymbol{f}_* \mid \boldsymbol{X}, \boldsymbol{y}, \boldsymbol{X}_*, \sigma_n^2\right)=\int p\left(\boldsymbol{f}_* \mid \boldsymbol{X}_*, \boldsymbol{w}, \sigma_n^2\right) p(\boldsymbol{w} \mid \boldsymbol{X}, \boldsymbol{y}) d \boldsymbol{w}, \quad \boldsymbol{f}_* \triangleq f\left(\boldsymbol{X}_*\right) $$ 式中 $\boldsymbol{X}_*$ 为预测数据， $\boldsymbol{f}_*$ 为预测结果。对0均值的正态先验，由其共轭性可知，预测结果的后验也为正态分布: $$ p\left(\boldsymbol{f}_* \mid \boldsymbol{X}, \boldsymbol{y}, \boldsymbol{X}_*, \sigma_n^2\right)=N\left(\boldsymbol{f}_* \mid \frac{1}{\sigma_n^2} \boldsymbol{X}_*^{\top} \mathbf{\Lambda}^{-1} \boldsymbol{X} \boldsymbol{y}, \boldsymbol{X}_*^{\top} \boldsymbol{\Lambda}^{-1} \boldsymbol{X}_*\right) $$ 上式右侧的正态分布可以提供预测结果的置信区间，在本质上是线性模型使用所有潜在权重系数计算所得结果的概率密度函数。
边缘似然描述了似然和先验的组合在多大程度上解释了观测数据，因此边缘似然可以用于计算贝叶斯因子，与其它模型进行 比较并验证模型和先验的合理性 ${ }^{[1]}$ 。由全概率公式（law of total probability) 可知，边缘似然有如下积分形式: $$ p\left(\boldsymbol{y} \mid \boldsymbol{X}, \sigma_n^2\right)=\int p\left(\boldsymbol{y} \mid \boldsymbol{X}, \boldsymbol{w}, \sigma_n^2\right) p(\boldsymbol{w}) d \boldsymbol{w} $$ 由上式计算的边缘似然拥有描述模型复杂度的全部信息，包括先验假设的类型，权重系数的数量等，因此可以与任意复杂度的其 它模型计算贝叶斯因子进行比较。