所谓直线拟合,通常也叫做线性拟合、一元线性回归。指的是当我们有一批数据(xi,yi),这些数据在平面坐标系下落在一条直线上,或近似的落在一条直线上。我们就要求出这条直线的参数。如果这条直线可以写为:
y=kx+b
那么
k=∑(xi?xˉ)(yi?yˉ)∑(xi?xˉ)2
b=yˉ?kxˉ
这个关系式许多教科书上都有详细的推导,无需多说。
今天要说的是另一种情况,当我们的数据有可能落在一条竖直的直线上,也就是k 有可能为∞ 时,应该如何做拟合。这时我们肯定就不能用y=kx+b 了,但是可以将这个表达式变变形。我们知道
k=tanθ=sinθcosθ
那么原来的直线方程可以写为:
ycosθ=xsinθ+bcosθ
或者写为更一般的形式:
ax+by+c=0
同时满足附加条件:
a2+b2=1
下面就来说说这种形式的直线方程如何拟合。
点到直线的垂直距离
首先先要解决一个小问题,一个点 (xi,yi) 到这条直线的距离是多少。
直接计算有点麻烦,我们先考虑一种简单的形式,点 (0,0) 到这条直线的距离。
直线方程:
xsin(θ)?ycos(θ)+c=0
那么与它垂直的过原点的直线方程为:
xsin(θ+π/2)?ycos(θ+π/2)=0
化简后得到:
xcos(θ)+ysin(θ)=0
这两个直线的交点为两条直线方程组成的方程组的解:
{xsin(θ)?ycos(θ)+c=0xcos(θ)+ysin(θ)=0
简单计算后可以得到:
{x=?sin(θ)cy=cos(θ)c
容易看出交点到原点的距离为 c,而这就是原点到直线的距离。
下面就可以求点 (xi,yi) 到这条直线的距离了。只需要做个坐标变换。新坐标系下的坐标变量为 x′ 和 y′。 新旧坐标系的关系如下。
{x′=x?xiy′=y?yi
那么原始坐标系下的点(xi,yi) 在这个新的坐标系下称为了坐标原点(0,0)。
原始坐标系下的直线方程 ax+by+c=0 成为了:
a(x′+xi)+b(y′+yi)+c=0
整理一下成为:
ax′+by′+(axi+byi+c)=0
在这个坐标系下直线到原点的距离为: (axi+byi+c),这也就是我们要求的直线到点的垂直距离。
直线拟合
我们要求的直线方程的系数 a、b 和 c 就是在
a2+b2=1
的条件下,使得下面求和式:
f=∑(axi+byi+c)2
取得极小值的 a、b 和 c。
这个问题有标准的处理方法:拉格朗日乘子法。
f=∑(axi+byi+c)2?λ(a2+b2?1)
????????????????f?a=0?f?b=0?f?c=0?f?λ=0
这里先不急着把这几个式子都展开了。只来看看 ?f?c 这一项。
?f?c=2×∑(axi+byi+c)=0
变形之后可以得到:
axˉ+byˉ+c=0xˉ=∑xi/Nyˉ=∑yi/N
用这个式子, c 就确定了。可以将确定后的 c 带入 f 的表达式。
f=∑(a(xi?xˉ)+b(yi?yˉ))2?λ(x2+y2?1)
?f?a=2∑(a(xi?xˉ)+b(yi?yˉ))(xi?xˉ)?2aλ=2(aDxx+bDxy)?2aλ=0
?f?b=2∑(a(xi?xˉ)+b(yi?yˉ))(yi?xˉ)?2bλ=2(aDxy+bDyy)?2bλ=0
其中:
?????Dxx=∑(xi?xˉ)2Dxy=∑(xi?xˉ)(yi?yˉ)Dyy=∑(yi?yˉ)2
所以,我们的 a、b 和 λ 应该满足如下线性方程组。
{aDxx+bDxy=aλaDxy+bDyy=bλ
或者写成矩阵形式:
(DxxDxyDxyDyy)(ab)=λ(ab)
这样写就很清楚了。λ 是 (DxxDxyDxyDyy) 的特征值, a,b 是 (DxxDxyDxyDyy)
的特征向量。
但是这个矩阵有两个特征值,我们应该选用哪个特征值呢。
∑(a(xi?xˉ)+b(yi?yˉ))2=(ab)(DxxDxyDxyDyy)(ab)=λ(a2+b2)=λ
所以,我们应该选取两个特征值中较小的那个。
(Dxx?λDxyDxyDyy?λ)(ab)=0
有非零解的条件是:
∣∣∣Dxx?λDxyDxyDyy?λ∣∣∣=0
λ=(Dxx+Dyy)±(Dxx?Dyy)2+4D2xy?????????????????√2
较小的那个特征根是: λ=(Dxx+Dyy)?(Dxx?Dyy)2+4D2xy√2
(ab)=1D2xy+(λ?Dxx)2???????????????√(Dxyλ?Dxx)
至此,这个直线拟合问题的解决了。但是,关于这个直线拟合与传统的一元线性回归算法的区别我还想说几句。
传统的一元线性回归是认为数据点 (xi,yi) 中只有 yi 是有误差的,因此确定点到直线距离时用的是 y 方向的距离。 我本文中的算法认为数据点 (xi,yi) 都是有误差的,并且不确定度是相同的,因此,数据点到直线的距离用的是垂直距离。这两个异同可以参考下面的图片。当这条直线的斜率很小时,这两种方法求得直线方程很接近,当直线斜率很大时,两个结果可能有很大的区别。具体应该用哪种还是要根据数据点的性质来确定。
后记
写到这里关于直线拟合算法的实现方法就算是说清楚了。但其实还有一个重要的问题没有交代。我们做一元线性回归时会去求相关系数。相关系数反映的是数据点与拟合直线的吻合程度,对于我们的算法,也应该提出一个类似的指标。另外,求出的系数 a、b、c 也应给出不确定度或者置信区间。
关于这些问题,我准备再开一篇博客来详细的说说。(其实是因为有些问题还没想好如何解决(^.^))
版权声明:本文为博主原创文章,未经博主允许不得转载。