常用算法分析——最小二乘法

作者: JamieCh | 来源:发表于2019-10-10 12:11 被阅读0次

常用算法分析——最小二乘法

1、引言

最小二乘法应该是我们最早接触的一种数值估计算法。它的特殊形式——一元线性回归，被广泛地应用于多种数值统计分析场合。例如，在验证欧姆定律（U = IR）时，通常的实验方法是分别测量出多个不同电压U_i下，通过电阻的电流值I_i，然后将这些(U_i, I_i)观测点，代入到一元最小二乘公式(1-1)中，便可计算出 $\hat{R}$ 。

$\begin{cases}a&=&\frac{\sum{xy}-\frac{1}{N}\sum{x}\sum{y}}{\sum{x^2}-\frac{1}{N}(\sum{x})^2}\\ b&=&\frac{1}{N}\sum{y}-\frac{a}{N}\sum{x}\end{cases}$ (1-1)

由此可得出线性拟合式(1-2)，

$\hat{y}=a\hat{x}+b$ (1-2)

其中， $\hat{y}=\hat{U},\ \hat{x}=\hat{I},\ a=\hat{R},\ b$ 是残差。通过此方法将观测点及拟合曲线绘制在同一个直角坐标系中，正常情况下可以直观地看到，观测点会均匀分布在直线附近，且每个点的残差平方和（即方差）最小。
“最小二乘法”由此得名。

2、普通最小二乘法（OLS）

最小二乘法显然不只是一元线性回归那么简单，它还可以应用于多元参数的拟合。本节将对普通最小二乘法（Ordinary Least Squares）的原理进行简单的推导和证明。

2.1、高斯—马尔可夫定理

高斯—马尔可夫定理（the Gauss–Markov theorem，简称G-M定理）　在给定经典线性回归的假定下，最小二乘估计量是具有最小方差的线性无偏估计量（即Best Linear Unbiased Estimator，简称BLUE)。

G-M定理共对OLS普通线性方程提出5个假设：
假设1（线性关系）：要求所有的母集团参数（population parameters）为常数，用来保证模型为线性关系。即，如果母集团方程为 $a,b_1,b_2,...,b_k$ 必须为常数，ε为无法检测的误差项，即实验过程中模型没有包含的因素。
假设2（随机采样）：假设我们有n个调查的样本，那么这n个样本必须是从母集团里面随机抽样得出的。
假设3（无完全共线关系）：在样本（母集团）中，没有独立变量是常数，且独立变量之间不能有完全共线性。
假设4（零条件均值）：母集团方程的误差项均值为0，且均值不受到独立变量的影响，即E(ε | X₁, X₂, …, X_k)＝0。
假设5（同方差性）：误差项ε的方差为一个固定不变的常数，且不受到独立变量的影响，即Var(ε | X₁, X₂, …, X_k)＝σ²或Var(ε)＝σ²I。

G-M定理的意义在于，当上述假定成立时，我们不需要再去寻找其它无偏估计量，因为没有一个会优于普通最小二乘估计量。换言之，如果存在一个好的线性无偏估计量，那么这个估计量的方差最多与普通最小二乘估计量的方差一样小，而不会小于普通最小二乘估计量的方差。

2.2、最小二乘估计量推导

假设使用多输入单输出（MISO）线性系统对未知系统进行估计，则估计系统的输入输出应满足

$y=\hat{\beta}_0+\hat{\beta}_1x_1+\hat{\beta}_2x_2+...+\hat{\beta}_kx_k=\sum_{j=0}^k\hat{\beta}_jx_j\ \ \ (x_0=1)$ (2-1)

此式称为系统的估计量。接下来，定义第 i 次观测值 y_i 应满足

$y_i=\beta_0+\beta_1x_{1i}+\beta_2x_{2i}+...+\beta_kx_{ki}+\varepsilon_i$ (2-2)

其中， $\small\beta_k$ 的估计量，ε_i 是第 i 次观测时无法预测的误差项。
由此可得

$\begin{array}{rcl}D&=& \sum_{i=1}^n \varepsilon_i^2\\ &=&\sum_{i=1}^n(y_i-\hat{y}_i)^2\\&=&\sum_{i=1}^n(y_i-\sum_{j=0}^k\hat{\beta}_jx_{ji})^2\end{array}$ (2-3)

若想求得使得式(2-3)取得极小值，那么需要对D分别求出 $\small\hat{\beta}_r$ 的偏导数，并令其等于零，即，

$\frac{\partial D}{\partial\beta_r} =2\sum_{i=1}^{n}(y_i-\sum_{j=1}^k\hat{\beta}_jx_{ji})x_{ri}=2\sum_{i=1}^{n}x_{ri}y_i-2\sum_{i=1}^{n}(x_{ri}\sum_{j=0}^k\hat{\beta}_jx_{ji})=0$
$\Rightarrow \sum_{i=1}^{n}x_{ri}y_i=\sum_{i=1}^{n}(x_{ri}\sum_{j=0}^k\hat{\beta}_jx_{ji})$ (2-4)

其中，r＝1, 2, 3, …, k，联立这k个关于 $\hat{\beta}_r$ 的方程，即为法方程组（或正规方程组），而法方程组的解即为式(2-1)的系数，通过此方法得到的估计量则称为系统的最小二乘估计量。可以证明，对于满足上述5个假设的OLS方程组，存在唯一实数解。
特殊地，当k＝1时，MISO退化成单输入单输出（SISO）的一元系统，则其法方程组可写作：

$\left\{ \begin{array}{rcl} \sum_{i=1}^n(\hat{\beta}_0+\hat{\beta}_1x_i)&=& \sum_{i=1}^ny_i \\ \sum_{i=1}^n(\hat{\beta}_0x_i+\hat{\beta}_1x_i^2)&=&\sum_{i=1}^nx_iy_i \end{array}\right. \\ \Rightarrow \left\{ \begin{array}{rcl} n\hat\beta_0+(\sum_{i=1}^nx_i)\hat{\beta}_1&=& \sum_{i=1}^ny_i \\ (\sum_{i=1}^nx_i)\hat{\beta}_0+(\sum_{i=1}^nx_i^2)\hat{\beta}_1&=&\sum_{i=1}^nx_iy_i \end{array}\right.$

求解后，可得出式(1-1)的结果。

2.3、矩阵形式

上节已推导出计算最小二乘估计量的一般方法，但显然，式(2-4)的形式过于繁琐。对于多元线性回归的场合，法方程组将会变得异常复杂，难以求解。同时，这种代数形式的公式也不便于得出更为一般的规律。因此，本节将用矩阵的形式，再次推导OLS的计算过程，并且尝试讨论法方程组解的存在问题。
首先，定义估计量

$\hat{y}={x}^\text{T}\hat{\boldsymbol{\beta}}$ (2-5)

其中， $x=(1,x_1,x_2,...,x_k)^\text{T},\ \hat{\boldsymbol{\beta}}=(\hat{\beta}_0,\hat{\beta}_1,\hat{\beta}_2,...,\hat{\beta}_k)^\text{T}$ 。
误差

$\boldsymbol{\varepsilon}=\boldsymbol{Y}-\hat{\boldsymbol{Y}}=\boldsymbol{Y}-\boldsymbol{X}\hat{\boldsymbol{\beta}}$ (2-6)

其中，
$\boldsymbol{\varepsilon}=(\varepsilon_1,\varepsilon_2,...,\varepsilon_n)^\text{T}$ ， $\boldsymbol{Y}=(y_1,y_2,...,y_n)^\text{T}$ ， $\boldsymbol{X}=(\boldsymbol{x}_1,\boldsymbol{x}_2,...,\boldsymbol{x}_n)^\text{T}=\left[\begin{matrix} 1&x_{11} &\cdots& x_{k1}\\1&x_{12}&&x_{k2}\\ \vdots&&\ddots&\vdots\\ 1&x_{1n} &\cdots& x_{kn}\end{matrix}\right]$ 。
取误差向量的2-范数^[1]的平方

$\begin{array}{rcl} \boldsymbol{D}&=&\left\|\boldsymbol{\varepsilon}\right\|_2^2=\boldsymbol{\varepsilon}^\text{T}\boldsymbol{\varepsilon} \\ &=& (\boldsymbol{Y}-\boldsymbol{X}\hat{\boldsymbol{\beta}})^\text{T}(\boldsymbol{Y}-\boldsymbol{X}\hat{\boldsymbol{\beta}}) \\ &=& (\boldsymbol{Y}^\text{T}-\hat{\boldsymbol{\beta}}^\text{T}\boldsymbol{X}^\text{T})(\boldsymbol{Y}-\boldsymbol{X}\hat{\boldsymbol{\beta}}) \\ &=& \boldsymbol{Y}^\text{T}\boldsymbol{Y}-\hat{\boldsymbol{\beta}}^\text{T}\boldsymbol{X}^\text{T}\boldsymbol{Y}-\boldsymbol{Y}^\text{T}\boldsymbol{X}\hat{\boldsymbol{\beta}}+\hat{\boldsymbol{\beta}}^\text{T}\boldsymbol{X}^\text{T}\boldsymbol{X}\hat{\boldsymbol{\beta}} \\ &=& \boldsymbol{Y}^\text{T}\boldsymbol{Y}-2\boldsymbol{Y}^\text{T}\boldsymbol{X}\hat{\boldsymbol{\beta}}+\hat{\boldsymbol{\beta}}^\text{T}\boldsymbol{X}^\text{T}\boldsymbol{X}\hat{\boldsymbol{\beta}} \end{array}$ (2-7)

然后对 D 求 $\hat{\boldsymbol{\beta}}$ 的偏导，并令其等于零，可得^[2]

$\begin{array}{rcl} \frac{\partial\boldsymbol{D}}{\partial\hat{\boldsymbol{\beta}}}=-2\boldsymbol{X}^\text{T}\boldsymbol{Y}+2\boldsymbol{X}^\text{T}\boldsymbol{X}\hat{\boldsymbol{\beta}}&=&0\\ \Rightarrow \boldsymbol{X}^\text{T}\boldsymbol{X}\hat{\boldsymbol{\beta}}&=&\boldsymbol{X}^\text{T}\boldsymbol{Y}\end{array}$ (2-8)

式(2-8)即为法方程组的矩阵形式。为了进一步证明极小值存在，可对D求 $\hat{\boldsymbol{\beta}}$ 的二阶导数

$\frac{\partial^2\boldsymbol{D}}{\partial\hat{\boldsymbol{\beta}}^2}=2\boldsymbol{X}^\text{T}\boldsymbol{X}$ (2-9)

因此，式(2-9)是正定阵，故该问题的最优解存在。

接下来，分情况讨论法方程组的解。
1）Rank( X )＝k＋1，即 X 列满秩。
令 C ＝ X^TX ，显然，Rank(C)＝k＋1，且C是k阶正定阵。因此，C^-1存在，同时，可知式(2-9)正定，故极值点存在，得

$\hat{\boldsymbol{\beta}}=[(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1}\boldsymbol{X}]^\text{T}\boldsymbol{Y}$ (2-10)

2）Rank( X )＝ r，0＜r＜k＋1。
此时，方程组并不符合假设3的条件，即存在完全共线的变量，但可以证明，方程组是相容的^[3]，因此，方程组有唯一的极小范数解^[4]，为

$\hat{\boldsymbol{\beta}}=[(\boldsymbol{X}^\text{T}\boldsymbol{X})^{+}\boldsymbol{X}]^\text{T}\boldsymbol{Y}$ (2-11)

这种情况下，估计量将无法保证是BLUE，此种情况已超出本文讨论的范围，从略。此时，可将共线的变量合并后，重新计算。
由此也可看出，若分析场景中无法保证所选变量完全共线时，建议先通过其他相关性分析算法，对数据进行一次前期的清洗，选择那些对因变量影响大、关联性强的自变量，剔除那些关联性弱或与其他自变量存在共线（或近似共线）的自变量。

综上所述，式(2-10)给出了计算最小二乘估计量的一般公式，根据此式，便可实现具体的最小二乘估计算法。

2.4、最小二乘估计量的统计特性

1）线性性
式(2-10)已表明， $\hat{\boldsymbol{\beta}}$ 是 Y 的线性组合。

2）无偏性
若证明 $\text{E}(\hat{\boldsymbol{\beta}})=\boldsymbol{\beta}$ 。
证明：

$\begin{array}{rcl}\text{E}(\hat{\boldsymbol{\beta}}) &=& \text{E}([(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1}\boldsymbol{X}]^\text{T}\boldsymbol{Y})\\&=&\text{E}([(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1}\boldsymbol{X}]^\text{T}(\boldsymbol{X}\hat{\boldsymbol{\beta}}+\boldsymbol{\varepsilon}))\\&=&\text{E}(\hat{\boldsymbol{\beta}}+[(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1}\boldsymbol{X}]^\text{T}\boldsymbol{\varepsilon})\\ &=& \text{E}(\boldsymbol{\beta}) \\ &=& \boldsymbol{\beta}\end{array}$ (2-12)

证毕。█

3）有效性
有效性是个定性的特性，即估计量的方差越小，则说明此估计量越有效。而G-M定理已指出，最小二乘估计量是最优的，即 $\text{Var}(\hat{\boldsymbol{\beta}})$ 是所有估计量中最小的。现在简单证明这一结论。
证明：
首先求 $\text{Var}(\hat{\boldsymbol{\beta}})$ ^[5]。

$\begin{array}{rcl}\text{Var}(\hat{\boldsymbol{\beta}}) &=& \text{E}([(\hat{\boldsymbol{\beta}}-\boldsymbol{\beta})(\hat{\boldsymbol{\beta}}-\boldsymbol{\beta})^\text{T}]\\&=&\text{E}([(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1}\boldsymbol{X}]^\text{T}\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}^\text{T}\boldsymbol{X}(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1})\\&=& [(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1}\boldsymbol{X}]^\text{T}\text{E}(\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}^\text{T})\boldsymbol{X}(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1}\\ &=& [(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1}\boldsymbol{X}]^\text{T}(\sigma^2\boldsymbol{I})\boldsymbol{X}(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1} \\ &=& \sigma^2(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1}\end{array}$ (2-13)

再证明式(2-13)是最小的。
令 $\tilde{\boldsymbol{\beta}}$ 是无偏的，则

$\text{E}(\boldsymbol{CY})=\text{E}(\boldsymbol{CX\beta}+\boldsymbol{C\varepsilon})=\boldsymbol{\beta}$

这表明，CX＝I。把 C 代入到式(2-10)、式(2-13)中可得，

$\text{Var}(\tilde{\boldsymbol{\beta}})=\sigma^2\boldsymbol{CC}^\text{T}$

令 $\boldsymbol{D}=\boldsymbol{C}-[(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1}\boldsymbol{X}]^\text{T}$ ，则

$\text{Var}(\tilde{\boldsymbol{\beta}})=\sigma^2\boldsymbol{DD}^\text{T}+\sigma^2[(\boldsymbol{D}+(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1}\boldsymbol{X})^\text{T})(\boldsymbol{D}+(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1}\boldsymbol{X})^\text{T})^\text{T}]$

由于 $\boldsymbol{DX}=\boldsymbol{\Theta}$ ，所以

$\text{Var}(\tilde{\boldsymbol{\beta}})=\sigma^2\boldsymbol{DD}^\text{T}+\sigma^2(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1}=\text{Var}(\hat{\boldsymbol{\beta}})+\sigma^2\boldsymbol{DD}^T$ (2-14)

因为DD^T是非负定矩阵^[6]，所以 $\text{Var}(\hat{\boldsymbol{\beta}})$ 是最小的。█

综上所述，最小二乘估计量是最小方差线性无偏估计量（BLUE）。

2.5、OLS推广

实际上，我们经常使用OLS对非线性系统进行回归分析。通常，有以下几种情形：多项式回归、对数回归、分段回归等。
1）多项式回归
令x＝(x⁰, x¹, x², …, x^k)^T，代入方程组进行求解。
但需注意，多项式的阶数不能过高，一般地，当k>3时，拟合误差会变得很大，拟合结果基本无法使用。因此，对于复杂曲线的拟合，可将曲线进行分段，然后采用低阶多项式进行拟合计算。
2）对数回归
令z＝lnx（或z＝e^x），代入方程组进行求解。
采用对数回归时，应注意x的定义域，若∃x≤0，则不能使用此方法。另外，对数有“缩小”效果，指数有“放大”效果。当x的范围跨度过大（lg(max(x)－min(x))＞3），且无明显线性特征时，可选用对数变换；当x的变化幅度过小（(max(x)－min(x))/mean(x)≪1），且无明显线性特征时，可选用指数变换。

3、OLS实现

对于简单的低阶（k≤2）估计而言，直接使用式(2-10)即可。但注意到，式(2-10)中需要计算矩阵 X^TX 的逆，而随着 k 的增长，这一步将变得十分低效，所以，应考虑使用某种方法避开矩阵求逆，于是我们想到了矩阵的QR分解。

3.1、QR分解

定理3.1 　设 $\boldsymbol{A}\in\mathbb{R}_n^{n\times n}$ ，^[7]则存在正交阵^[8] Q 及正线上三角阵 R ，使得 A ＝ QR 。
证明：因为 A 是满秩方阵，将 A 写成列分块形式： A ＝(a₁, a₂, …, a_n)，则a₁, a₂, …, a_n 均为 $\frak{B}=\lbrace \boldsymbol{u}_1, \boldsymbol{u}_2,..., \boldsymbol{u}_n\rbrace$ ，且 (a₁, a₂, …, a_n)＝(u₁, u₂, …, u_n)R 成立，其中 R 为正线上三角阵。
但由于 $\small a_i\in\mathbb{R}^n,u_i\in\mathbb{R}^n$ ，故有

$\boldsymbol{A} ＝ \boldsymbol{QR}$ (3-1)

注意 Q＝( u ₁, u ₂, …, u _n )的各列标准正交，故 Q^TQ＝I ，即 Q 为一个正交阵。█
Q有以下良好的性质：
1）Q 是非奇异的，且det(Q)＝±1；
2）Q ^-1＝ Q ^T；
3）正交阵的乘积仍是正交阵；
4）Q 的特征值的模为1。
定理3.2（Householder变换）　对于任意 $\boldsymbol{H}=\boldsymbol{I}_n-2\boldsymbol{uu}^\text{T}$ ，使得

$\boldsymbol{Hx}=a\boldsymbol{e}$ (3-2)

成立。其中 $\boldsymbol{e}$ 为n维自然基^[9]， $a$ 是实数，且 $|a|=‖\boldsymbol{u}‖_2$ 。
计算可得，

$\boldsymbol{u}=\frac{\boldsymbol{x}-a\boldsymbol{e}}{\left\|\boldsymbol{x}-a\boldsymbol{e}\right\|_2}\ \ (|a|=\left\|\boldsymbol{x}\right\|_2)$ (3-3)

定理3.2说明，任意列向量被Householder阵左乘后，可变换成与某自然基的共线向量。
推论1（Householder变换法） 对于满秩方阵 A ，将其按列分块得 A ＝(α₁, α₂, …, α_n)，取

$\boldsymbol{u}_1=\frac{\boldsymbol{\alpha}_1-a_1\boldsymbol{e}_1}{\left\|\boldsymbol{\alpha}_1-a_1\boldsymbol{e}_1\right\|_2}\ \ (|a_1|=\left\|\boldsymbol{\alpha}_1\right\|_2)$ (3-4)

则 $\boldsymbol{H}_1=\boldsymbol{I}-2\boldsymbol{u}_1\boldsymbol{u}_1^\text{T}$ ，有

$\boldsymbol{H}_1\boldsymbol{A}=(\boldsymbol{H}_1\boldsymbol{\alpha}_1+\boldsymbol{H}_1\boldsymbol{\alpha}_2+\ldots+\boldsymbol{H}_1\boldsymbol{\alpha}_n)=\left[\begin{matrix} a_1 & * & \ldots &* \\0& \\ \vdots & &\boldsymbol{B}_1 & \\0& \end{matrix} \right]$ (3-5)

将矩阵 $\boldsymbol{B}_1\in\mathbb{R}^{(n-1)\times (n-1)}$ 按列分块，得 B₁ ＝(β₂, …, β_n)，取

$\boldsymbol{u}_2=\frac{\boldsymbol{\beta}_2-a_2\boldsymbol{e}_2}{\left\|\boldsymbol{\beta}_2-a_2\boldsymbol{e}_2\right\|_2}\ \ (|a_2|=\left\|\boldsymbol{\beta}_2\right\|_2)$ (3-6)

$\tilde{\boldsymbol{H}}_2=\boldsymbol{I}-2\boldsymbol{u}_2\boldsymbol{u}_2^\text{T}\ \ \ \boldsymbol{H}_2=\left[\begin{matrix}1&\boldsymbol{\theta}^\text{T}\\\boldsymbol{\theta}&\tilde{\boldsymbol{H}}_2\end{matrix}\right]$ (3-7)

则

$\boldsymbol{H}_2\boldsymbol{H}_1\boldsymbol{A}=\left[\begin{matrix} a_1 & * & \ldots& \ldots &* \\0&a_2&*&\ldots&* \\ \vdots &0 \\ \vdots&\vdots&&\boldsymbol{C}_2 & \\0&0 \end{matrix} \right]$ (3-8)

依次进行下去，得到第n－1阶的Householder阵H_n－1，得

$\boldsymbol{H}_{n-1}\ldots\boldsymbol{H}_2\boldsymbol{H}_1\boldsymbol{A}=\left[\begin{matrix} a_1 & * & \ldots &* \\&a_2&\ldots&* \\ &&\ddots&\vdots\\&&&a_n \end{matrix} \right]=\boldsymbol{R}$ (3-9)

由于H^-1＝H，所以令Q＝H₁H₂…H_n-1，则A＝QR。█

由此可见，使用推论1提供的方法，对于任意满秩方阵，只需经过向量加减、矩阵相乘运算和基本四则运算就能实现QR分解。

现将X进行QR分解，即令X＝QR，则式(2-8)得

$\begin{aligned}\ \ \ \ \boldsymbol{X}^\text{T}\boldsymbol{X}\boldsymbol{\hat{\beta}}&=\boldsymbol{X}^\text{T}\boldsymbol{Y}\\ \Leftrightarrow (\boldsymbol{QR})^\text{T}\boldsymbol{QR}\boldsymbol{\hat{\beta}}&=(\boldsymbol{QR})^\text{T}\boldsymbol{Y}\\ \Leftrightarrow \boldsymbol{R}^\text{T}\boldsymbol{Q}^\text{T}\boldsymbol{QR}\boldsymbol{\hat{\beta}}&=\boldsymbol{R}^\text{T}\boldsymbol{Q}^\text{T}\boldsymbol{Y} \\ \Leftrightarrow \boldsymbol{R}^\text{T}\boldsymbol{R}\boldsymbol{\hat{\beta}}&=\boldsymbol{R}^\text{T}\boldsymbol{Q}^\text{T}\boldsymbol{Y} \\ \Leftrightarrow \boldsymbol{R}\boldsymbol{\hat{\beta}}&=\boldsymbol{Q}^\text{T}\boldsymbol{Y}\end{aligned}$ (3-10)

求解方程组(3-10)显然不需要求 R^-1 了，这是因为 R 是上三角阵，所以只需先计算出 β_n ，然后将其代入到第n－1个方程中，计算出β_n－1；再将β_n、β_n－1代入到第n－2个方程中，计算出β_n-2，以此类推，直至求出所有β为止。

综上所述，借助QR分解的方法，可以避免矩阵求逆，有效地降低了最小二乘算法的复杂度。

3.2、Apache Commons Math

Apache公共包中有一个专门用于数学计算的工具包Math（org.apache.commons.math3）^[10]，它提供了丰富且高效的数学分析算法，这其中就包含普通最小二乘法（org.apache.commons.math3.stat.regression.OLSMultipleLinearRegression类）。
OLSMultipleLinearRegression类的求解分两步：
1）按照推论1给出的方法，对X进行QR分解；
2）借助org.apache.commons.math3.linear.QRDecomposition类，求解方程组(3-10)。
此外，还提供评价OLS效果的一些系数计算方法，如SSR、SSTO、r²，调整的r²。

3.3、应用实例

基于Math包，我们可以轻易地开发出应用最小二乘法进行数据分析的功能。下面就以对一个时间序列，进行二次多项式拟合为例，说明OLSMultipleLinearRegression类的具体使用方法。
例：设某时间序列为 S＝(1.6, 4.4, 9.4, 16.6, 25.5)^[11]，猜测其符合s(k)＝b₀＋b₁k＋b₂k²，k＝1, 2, 3, 4, 5。现使用最小二乘法对 S 进行二次多项式拟合。
解：构造
$\boldsymbol{X}=\left[\begin{matrix}1&1&1\\1&2&4\\1&3&9\\1&4&16\\1&5&25\end{matrix}\right],\ \boldsymbol{Y}=\left[\begin{matrix}1.6\\4.4\\9.4\\16.6\\25.5\end{matrix}\right]$
然后传给如下代码：

 OLSMultipleLinearRegression regression = new OLSMultipleLinearRegression();
 regression.newSampleData(Y, X);
 double[] b = regression.estimateRegressionParameters();
 // s=b[0]+b[1]*k+b[2]*k*k
 double r2 = regression.calculateAdjustedRSquared();

所得数组b，即为拟合系数。█

4、广义最小二乘法（GLS）简介

在实际情况中，OLS成立的5个假设常常难以满足。当假设5（同方差性）不能满足时，即ε的方差随输入的变化而变化，那么之前对OLS的推导过程就不再成立了。对此，可以通过一些线性变换，使得模型满足同方差的假设，即又转换成OLS模型了。
具体地，设线性系统 Y ＝ Xβ ＋ ε ，若Var(ε)＝σ²Ω，其中，Ω是正定对称阵。
则令Ω^-1＝P^TP，则Ω＝P^-1(P^T)^-1，那么，

$P\Omega P^\text{T}=P(P^{-1}(P^\text{T})^{-1}) P^\text{T}=I_n$

在式(2-6)等号两端同乘P，有

$P\varepsilon=PY-PX\hat{\beta}$ (4-1)

进一步地，上式可表示为

$\varepsilon^{*}=Y^{*}-X^{*}\hat{\beta}$ (4-2)

其中，ε^*＝Pε，Y^*＝PY，X^*＝PX。变换后的随机干扰项具有均齐方差：

$Var(\varepsilon^{*})=Var(P\varepsilon)=\sigma^2P\Omega P^\text{T}=\sigma^2I$

可以看出，模型(4-2)已满足OLS的所有假设，因此可以采用OLS进行估计。模型(4-2)的OLS估计称为原模型的GLS估计量：

$\hat{\beta}=(X^{*T}X^*)^{-1}X^{*\text{T}}Y^*=(X^\text{T}\Omega ^{-1}X)^{-1}X^\text{T}\Omega ^{-1}Y$ (4-3)

式(4-3)的结论是在 Ω 已知的情况下给出的，但遗憾的是，在绝大多数情况下，Ω 是未知的。因此，我们需要采用以下两步法进行估计：
第一步：估计 Ω 中的未知参数，得到 $\hat{\Omega}$ ；
第二步：利用第一步得到的 $\hat{\Omega}$ 进行GLS估计：

$\hat{\beta}=(X^\text{T}\hat{\Omega}^{-1}X)^{-1}X^\text{T}\hat{\Omega}^{-1}Y$

在实际应用过程中，我们往往将上述＂两步法＂重复进行多次，直到相邻两次的估计结果差别很小为止（称为＂收敛＂）。这种做法有助于提高估计的有效性。此方法称作可行的广义最小二乘估计法（FGLS）。

注：

在矩阵理论中，向量的p-范数都是等价的，此处之所以取2-范数，很大程度上是出于接下来求偏导简便的考虑。 ↩
∂x/∂x＝I，∂x^T/∂x＝I，∂(x^TA)/ ∂x＝A，∂(x^TA^TAx)/ ∂x＝2A^TAx。 ↩
线性方程组Ax＝b相容的充要条件是AA⁺b＝b，其中A⁺是A的M-P广义逆。 ↩
即满足 ‖x‖=min_(Ax=b)‖x‖ 的解，其中‖‖是Cⁿ中由内积诱导的范数。 ↩
向量方差与标量方差的求法不同，向量方差对应的是协方差矩阵：Var(b)=E[(b-E(b)(b-E(b))^T)]。 ↩
关于DD^T的二次型是q^TDD^Tq＝z^Tz≥0。 ↩
R_n^(n×n)表示n阶的实数方阵，矩阵的秩为n，即满秩方阵。 ↩
即QQ^T＝Q^TQ＝I。 ↩
只有一个元素为1，其他元素均为0的列向量。 ↩
官网：http://commons.apache.org/proper/commons-math/。 ↩
按s(k)＝0.5＋k²设计，存在随机误差。 ↩

常用算法分析——最小二乘法

常用算法分析——最小二乘法

目录

1、引言

2、普通最小二乘法（OLS）

2.1、高斯—马尔可夫定理

2.2、最小二乘估计量推导

2.3、矩阵形式

2.4、最小二乘估计量的统计特性

2.5、OLS推广

3、OLS实现

3.1、QR分解

3.2、Apache Commons Math

3.3、应用实例

4、广义最小二乘法（GLS）简介

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读