动机

一直以来，使用机器学习的算法都是用他人写好的类库，总觉得云里雾里的，弄不清楚到底怎么回事。今天实现了一个最简单的线性回归分析，觉得收货很大。纸上得来终觉浅，绝知此事要躬行。

回归分析

数据

假设有一组数据,知道自变量和因变量的值，如下例：

3.0000   10.0000
3.1000   10.3000
3.2000   10.6000
3.3000   10.9000
3.4000   11.2000
3.5000   11.5000
3.6000   11.8000
3.7000   12.1000
3.8000   12.4000
3.9000   12.7000
4.0000   13.0000
4.1000   13.3000
4.2000   13.6000
4.3000   13.9000
4.4000   14.2000
4.5000   14.5000
4.6000   14.8000
4.7000   15.1000
4.8000   15.4000
4.9000   15.7000
5.0000   16.0000

第一列为自变量，第二列为因变量。 使用matlab绘制出一个图像。

plot(x,y,'o');

可以看到：

模型假设

通过观察我们发现，自变量和因变量大体位于一条线上，这样我们就假设他们遵循一个线性函数。

$f(x) = wx+b$

因为现在的自变量x是一个数字（或者看做1维的向量），所以w和b也是一个数字。

现在的情况是，我们知道自变量和变量，但不知道这个w和b。所以我们的目的就是要估计出一个最合理的w和b。

补充

* 入门读者 这个部分可以不看 *

在实际问题中，x可能是一个长度为n的向量，那么w的维度要和x一致。

$f(x) = w_1*x_1+w_2*x_2 + ... + w_n*x_n + b)$

这种写法非常的繁琐，可以用线性代数的写法，实际上这种写法也是一般所用的写法。

$f(x) = w\dot{} x +b$

这里w是一个横向量，x是一个列向量，b 是一个标量（就是一个数字）。

误差分析

为了刻画我们的模型的好坏，我们需要给出一个衡量标准。那么怎么衡量的。
通过一定的算法，我们得到了w和b的估计值 $\widehat{w}$和$\widehat{b}$
那么 y的估计值 ypred:

$ypred = \widehat{w} x + \widehat{b}$

这个y的估计值 ypred 和y 的差的绝对值 $| ypred - y |$越小，说明估计的参数越准确。我们将所有的样本(x,y) 的误差加起来就得到了总的误差 :

$$L = \sum_{(x,y)}{|ypred - y|}$$

在数学上，计算绝对值不是很方便，因为绝对值函数在0点不可导。我们希望把他给换成一个功能类似的可导函数。

为什么需要可导呢？ 后面会看到。

我们使用误差平方和函数来代替绝对值函数。这个又叫做二阶范式。

$$L = \frac{1}{2} (ypred - y)^2= \frac{1}{2}(wx+b - y)^2$$

note: 为什么要加二分之一？

这个函数是可导的。

$$y = \frac{x^2}{2}$$ , 其图像如下:

误差最小化

为了能得到最好的w和b 的估计，我们希望上面表示出的误差能最小。那么如何求一个函数的最小值呢？ 微积分告诉我们，当一个函数导数为0的时候，它去极大值或极小值。 那么函数的最小值就在导数为0的地方。

我们把误差函数看成是w的函数，对w进行求导，可以得到：

$$\frac{\partial{L}}{\partial{w}} = \sum_{(x,y)} (wx + b - y)*x$$

如果忘记了怎么求导，可以看看我这篇文章复习一下。
http://blog.csdn.net/taiji1985/article/details/72857554

误差函数对b求导，可以得到:

$$\frac{\partial{L}}{\partial{b}} = \sum_{(x,y)} (wx + b - y)$$

数学方法求误差最小

上面得到了L关于w和b的导数，另导数等于0，可以得到极值。

$$\frac{\partial{L}}{\partial{b}} = 0$$

$$\frac{\partial{L}}{\partial{w}} = 0$$

带入导数公式得到

$$\sum (wx + b - y) = 0$$
$$\sum (wx + b - y)*x = 0$$

求解可得最小二乘法公式，这里就不求了。

数值方法求解

上面的方法可以通过让导数为0得到结果，但对于一切比较复杂的模型，求解导数为0的公式是很困难的。这样就需要一些数值方法来学习。

数值方法一般使用迭代的方法来逼近最终答案，最常用的为梯度下降法。梯度下降法的原理是这样的。 随便选择一个x的初始点$x_0$ ,在该点沿着其导数的相反方向运动一小段，就可以靠近极小值点。（如果是沿着导数的方向移动就会得到极大值）。

具体做法是，在每次迭代时，计算w的新值

$$w' = w - a*\frac{\partial{L}}{\partial{w}} = w - a*\sum_{(x,y)} (wx + b - y)*x$$

$$b' = b - a*\frac{\partial{L}}{\partial{b}} = b - a*\sum_{(x,y)} (wx + b - y)$$

其中a称为学习因子,就是梯度每一步走多长。。如果学习因子太大，容易错过极值点，如果比较小，则收敛变慢。 学习因子一般通过经验给出。 可以先设置一个较小的数，逐渐调大，知道找到一个合适的值。

对于我们这个例子，直接让$a = 0.001$

算法

有了核心的思路就可以写算法了。先给出一个伪代码

    while(还没有收敛){计算误差计算误差和上一次误差的差如果这个差小于某个值，认为是收敛，退出。根据公式，更新w和b}

matlab代码如下：

function [w,b] =  predict(x,y)d = 10;n = length(x);w = rand(1,1);  % 随机生成 y = wx+b 中的w和b 做为初始值b = rand(1,1); olde = 0;       %上一次得到的误差i = 0;          %迭代计数器rate = 0.001;while d>0.001    % d表示两次迭代的误差的差i= i+1;%计算 yp = w*x+b;  %根据公式计算，当前w和b下，计算出的y的估计值e = sum((yp-y).*(yp-y)); %将估计值和真实值的差的平方和作为误差评估函数d = abs(olde  - e);     %计算上一次循环和这一次之间的误差的差fprintf('%d iter e = %f , d = %f \n',i,e,d);  if d<0.01break;endolde = e;%修改w和bw = w - rate*sum((yp-y).*x);b = b - rate*sum(yp-y);endhold off;plot(x,y,'o');hold on;xp = min(x):0.01:max(x);yp = w*xp+b;plot(xp,yp);end

写一个函数来调用这个学习函数

x = (3:0.1:5)'; 
w = 3; % 这是真实的参数值
b = 1;
y = w*x+b;
seed = 333;  % 随机数种子，为了让这个实验在重新运行时能重现所有现在的结果
rand('seed',seed)
y = y+rand(size(x))*0.7;
[wp,bp] = predict(x,y)

运行结果

1 iter e = 3335.957691 , d = 3335.957691 
2 iter e = 1348.893476 , d = 1987.064216 
3 iter e = 545.781076 , d = 803.112400 
4 iter e = 221.186873 , d = 324.594203 
5 iter e = 89.995526 , d = 131.191346 
6 iter e = 36.971877 , d = 53.023650 
7 iter e = 15.541291 , d = 21.430586 
8 iter e = 6.879684 , d = 8.661607 
9 iter e = 3.378919 , d = 3.500765 
10 iter e = 1.964014 , d = 1.414905 
11 iter e = 1.392151 , d = 0.571863 
12 iter e = 1.161021 , d = 0.231130 
13 iter e = 1.067605 , d = 0.093416 
14 iter e = 1.029850 , d = 0.037756 
15 iter e = 1.014590 , d = 0.015260 
16 iter e = 1.008422 , d = 0.006168 wp =2.9827bp =1.4161

图中直线为根据估计出的w和b画出的函数曲线。