每天进步一点点《ML - 线性回归》_ml线性回归

本系列Machine Learning 的学习博文很大部分全是受教于 Andrew Ng 的Stanford录制的ML视频，我的男神。

致敬男神！
第一次写这样系列的博文，单纯是为了记录自己的所学，做好记录和总结，免得遗忘，水平也不会很高，适合我这种初学者哈，也希望能分享给其他初学者，哪里有写的不对或者错误的地方，请留言直接指正，后续随着学习的深入也会不断更新，谢谢！

一：前序
机器学习：根据一些给定实例X和预期Y，通过训练某个模式P，用于预测某个新的实例X_new的输出Y_new。
机器学习包括监督学习和非监督学习。
监督学习：给定每个实例的输出，训练某个模型，对新的实例做输出预测。这里还包括对连续值的预测和非连续值的预测，对连续值的预测称为回归（Regression），对非连续值的预测称为分类（Classification）。
非监督学习：不知道每个实例的具体输出，具体分布含义，通过算法自动做聚类，对某个新的实例自动去归类，自动发现实例之间的分布类别规律，帮助人们发现数据的类别规律。

为了方便一些概念的描述，我先提前做一些约定。

3. 参数向量用θ来表示. θ = (θ_0, θ_1, θ_2…… θ_n)T 是一个(n+1)*1的矩阵，列向量。
   其中θ_i代表了第i个参数。

4. 输出Y = (y_1, y_2…… y_m)T是一个m*1的矩阵。
   其中y_i代表了第i个实例的预期输出。

假设模型函数（Hypothesis）称为H(θ)

本文也会对一些关键名词标注英文单词，这样也能学习一些关键名词的英文词汇，我们还处在学习阶段。
线性回归是属于对连续数值进行假设预测的过程，属于监督学习。

二：线性回归（Linear Regression）
1：模式表示（Model Representation）
H(θ)(X) = Xθ （X是MN的矩阵，代表M个实例，θ是N1的列向量）

需要注意的是，每个实例的特征也许是自然采集的特征，有一些可能是经过认为计算后的特征，是有某几个特征组合而成（例如：房屋地皮的长度和宽度，还可以合并成一个特征面积，甚至含有更高次幂的项，比如立方，四次方，更高次方，混合项）
我们称这样的含有多项式的回归模型叫做，多项式回归（Polynomial Regression）。
例如下两多项式的例子：

所以，多项式的表达式的样子是很多很多的，当然构造的多项式越复杂，越能构造出更加复杂的模型曲线。越简单的多项式构造出的模型曲线越简单。

2：代价函数（Cost Function）
我们评估模型好坏，得要看模型计算的误差，我们最平常的方式是：

代价函数是关于θ的函数，将所有样本实例的预期输出和实际预测输出做差，求平方，再累加，既是对预测函数的功能上好坏的判断。这个也很直观：
1：当假设函数预测的非常准确，最适合预测这些样本，那么J(θ) 趋近于0值，即是趋近于无偏差， 越小，对样本的估计误差也越小，越符合要求。
2：当假设函数预测的非常不准，不适合预测这些样本，那么J(θ) 就会得到很大的值，这样的假设函数显然是不可取的，J(θ)越大，对样本的估计误差也越大，越不可取。
然而，我们说了一大堆的概念，到目前为止，我们预先写出来的回归函数，仅仅知道输入和输出，对于参数θ是多少一无所知，怎么求出来这个参数呢？这个才是重点。且看下一小章节的讲解。下面给两个直观的印象，用我们能画出来形式给出：
直观的印象一（含有一个参数）：

给出四个点（1，1.1），（2，2.3），（3，3.2），（4，4.6），画出J(θ)随着θ的函数图像

从数学上和图像都能分析出J(θ)是个抛物线，且当θ接近1的时候有最小的代价函数值，也就是有最小的误差.此时，这个模型简单的就出来了。这个很简单，能用眼睛看出来的，大部分问题是看不出来的。
计算过程如下：

% .* 逐个元素相乘法
% .^ 逐个元素取幂指数

clc;clear all;close all;
set(0,'defaultfigurecolor','w');
x = [1, 2, 3, 4];
y = [1.1, 2.3, 3.2, 4.6];

theta = -2:0.1:4;
cost = zeros(size(theta));
for m=1:size(theta)(2)
	cost(m) = sum((x .* theta(m) - y) .^ 2) / (2 * size(x)(2));
end;

figure();
subplot(1,1, 1)
scatter(theta, cost, 'r', 'linewidth', 3);
xlim([-2,7]);
ylim([-3,100]);
% axis equal;
grid on;
xlabel('theta');
ylabel('cost value');
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

直观的印象二（含有两个参数）：

给出四个点（1，1.1），（2，2.3），（3，3.2），（4，4.6），画出J(θ)随着θ的函数图像

随着θ1 和θ0 的变化，cost function value 也随着变化，终有取最小值的一系列点。
计算过程如下：

src_x = [1, 2, 3, 4];
src_y = [1.1, 2.3, 3.2, 4.6];
n = 50;
x = linspace (-5, 5, n)'; % 在x轴上取50点
y = linspace (-8, 10, n)'; % 在y轴上取50点
[xx, yy] = meshgrid (x, y); % xx和yy都是矩阵
cost = zeros(n, n);

for i=1:n
	for j=1:n
		cost(i,j) = sum((src_x .* xx(i,j) + yy(i,j) - src_y) .^ 2) / (2 * size(src_x)(2));
	end;
end;

%h= surf(xx,yy,cost); %画出立体曲面图
%colorbar;
%shading interp; %平滑网格形成平滑面

%contour(xx,yy,cost); %等高线
%colorbar;
%shading interp;

surfc(xx,yy,cost); %立体加等高线
colorbar;
shading interp; 

xlabel('theta1');
ylabel('theta0');
zlabel('cost Value');
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

直观的印象三（含有高次项）：

给出四个点（-3，9.1），（-2，4.1），（-1，1.2），（0，0.1），（1，1.1），（2，3.9），（3，9.2），画出J(θ)随着θ的函数图像

上图为三维面(碗状)和等高线，下图为等高线，可以发现在θ1= 1且θ0=0的点有最毒低点。
计算过程如下：

src_x = -3:1:3;
src_y = [9.1, 4.1, 1.2, 0.1, 1.1, 3.9, 9.2];

n = 50;
x = linspace (-5, 5, n)'; % 在x轴上取50点
y = linspace (-8, 10, n)'; % 在y轴上取50点
[xx, yy] = meshgrid (x, y); % xx和yy都是矩阵
cost = zeros(n, n);

for i=1:n
	for j=1:n
		cost(i,j) = sum(((src_x .^2) .* xx(i,j) + src_x .* yy(i,j) - src_y) .^ 2) / (2 * size(src_x)(2));
	end;
end;

%h= surf(xx,yy,cost); %画出立体曲面图
%colorbar;
%shading interp; %平滑网格形成平滑面

contour(xx,yy,cost); %等高线
colorbar;
%shading interp;

%surfc(xx,yy,cost); %立体加等高线
%colorbar;
%shading interp;

xlabel('theta1');
ylabel('theta0');
zlabel('cost Value');
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

为什么我要花这么多的精力讲解这个代价函数呢，就是为了让我们对其意义有形象的认识，这个也是后续的重中之重，因此我们可以得出结论，我们可以在大家函数的最低点，得到最优解的θ。
这个过程叫做最小化代价函数（Minimization of cost function）目的就是，通过最小化代价函数得到我们需要的最合适的参数，下一小章节我们讲解找到能够使代价函数最小化的参数的方法。

3：最小化过程（Minimization process）
我们知道了代价函数J(θ)，代价函数是一个关于θ的求值函数，也就是说每当我们给出一组θ的值，也就对应一个代价值J(θ)，我们希望代价最小，因为代价越小，才是越适合我们的输入数据，误差越小，因此，也就是要对这个函数求最小值，就能反得出对应的θ，问题转换为了纯数学期末考试题目，当θ是多少的时候，J(θ)有最小值。这样一来我们解题思路就多样化了。这样的解可能不止一种，所以涉及局部最优和全局最优的讨论。

1.偏微分求导（Derivative of partial differential）
如果你对这部分内容不熟悉，那么需要补一下高数的知识。当θ仅有一个参数的时候，在直角坐标系下就是一条抛物线，当θ有俩参数的时候，在三维坐标系下就是一条抛物面（碗状的），所以在最小值处，导数/偏导数是等于0的。因此我们需要求出各个偏导数，并让其等于0，即可得出。
For (i:n)
Let each Derivative of partial differential to θ(j) = 0
To solve this equation and get the value of θ(j).
End;

2.梯度下降（Gradient Descent Algorithm）
想当初第一次看到这个计算过程的时候，以及学习这门课程的时候，那时候还是大三下学期，我彻底醒悟了，真的彻底醒悟了，才真的发现数学是真的有天大的卵用的（我这个人笨，请大神尽情吐槽我），顺带吐槽一下大部分大学的授课方式，只管教授知识，不讲述知识的所用之处，我敢保证大部分学生都是期末考试结束了就立马下学期还给老师了，所以大部分的时候都会不清楚所学，后来我进入华为工作，现在的我重新自学一遍高数，发现现在学起来很清楚，拿到那个时候绝对能考很高分数，很有心得体会地去学，也能学的明白，不再是机械式的记忆。Better late than never，也不算太晚。

回到正题，大概描述下这个梯度下降的过程：
首先从给定某个具体值的，也就是给θ赋值初值。
一直去改变θ的值，去reduce J(θ)，直到我们达到了在期望上的最小值，也叫做收敛（convergence）。

How to reduce the J(θ)，我们可以想象一下，J(θ)函数是一座大山，山上每个点都是由坡度，对应于函数的导数/梯度，我们想要得到J(θ)的最小值，也就是去山脚谷底，策略就是每次走一小步，往坡度下降的方向走，再往坡度下降的地方走，总会走到一个谷底。对应于数学，我们每次都改变参数，沿着斜率/导数/梯度下降的方向变化，而且没变化一小点，就能到函数图像的小值点。下面亲看：

这些都是《数值计算方法》教程的内容，下面是具体操作过程。
令θ = (0,0,0…0)，假设给全部θ赋值初值0。

α 就是步长，每一移动一点点，
∂/∂x J(θ_j) 就是斜率下降的方向和倾斜程度。
给个图像简单说明一下吧：

假设J(θ)图像如上图所示，假设最小点是在θ=1的位置，我们赋初值是θ=0，在0~1的区间内，斜率是负值，因此在改变θ的值，θ会增大，也就是朝着θ=1的位置靠近，当我们赋初值在θ=1，在1~2的区间内，斜率是正值，因此在改变θ的值，θ减小，也就是朝着θ=1的位置靠近。
不管你出发位置在哪里，总是朝着最低点前进了。
1
2

更加确切地说呢是。每一步迭代更新的时候，仅仅只能用上一次迭代的θ值：
令θ = (0,0,0…0)，假设给全部θ赋值初值0。

那么 ∂/∂x J(θ_j)是多少呢？来，一起简单推倒下：

因此对θ_j求偏导可得。

所以，最小化过程为：

循环迭代，直到θ值不再改变，或者代价函数J(θ)的值不再改变，或到达某个小值区间既可认为不再改变。诶？说到这里啊，我们发现我们直接求偏导数，让∂/∂x J(θ_j) = 0，θ值不就再也不会改变了么？所以啊，和方法一是相通的。

3）正规方程（Normal Equation）

Y
X是样本输入矩阵MN，Y是样本的输出矩阵M1。这里我不打算给出证明和获得解过程，我们可以发现这个一个矩阵计算公式就可以搞出来，是不是特别简单？在octave中使用 inv 或者 pinv函数即可求解矩阵的逆。
接下来我会给出这个方法和梯度下降的区别：
梯度下降：需要选择一个α；需要很多次迭代；当特征数量很大的时候也能表现的很好。
正规方程：不需要选择一个α；不需要多次迭代；但是当特征数量极大的时候，求解矩阵的逆运算将变得很费力，因为 (XT * X)-1是一个N^2的矩阵。所以建议在特征值小的时候比如数量在1000以内的，可以选择该方法，但是当数量极大的时候，建议还是老实用梯度下降。

4）其他方法
Conjugate gradient
BFGS
L- BFGS
这几个算法，这篇文章暂时不展开，等候专门抽出一篇来介绍。暂时留下个印象，这几个计算过程不用人工选择α，通常运行速度比梯度下降要快，但是却更加复杂，篇幅有限，之后再详细描述了。

三：用例测试实践（TestCase）
为了能方便画图表述，我接下来会用一些简单维度的例子来描述，也会给出描述高维的公式来描述，但是就不能画图出来了。
举个例子：企业员工加班是很正常的现象，我们现在有一个统计数据是，员工的加班时长和工作进程贡献指数（别问我指数咋来的）。如下表的样本分布：
0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6
1．1 2.2 3.2 4.1 4.7 5.3 5.8 6.2 6.5 6.6 6.65 6.68

我们假设有好几种拟合模型：并且呢我想用更一般的梯度下降的方法来演示一遍，直接求导和nornal equation的方法我就不操作了。

模型一：先来看看单一的常数项的简单的线性回归
我们只有一个特征变量可得：Y = H(θ)(X) = Xθ = θ0

迭代次数和相应的代价函数值的曲线
计算过程如下：

clc;clear all;close all;
set(0,'defaultfigurecolor','w') ;
src_x = 0:0.5:6;
src_y = [0, 1.1, 2.2, 3.2, 4.1, 4.7, 5.3, 5.8, 6.2, 6.5, 6.6, 6.65, 6.68];

% model 1 : y = Xθ, θ = (θ0), X = (1)
X = ones(size(src_x)(2), 1);
X(:,1) = 1;
theta = zeros(size(X)(2), 1); % init value to 0

itr = 500; %iteration times
rate = 0.02; %step
cost = zeros(1, itr);
for it=1:itr
	% update each theta
	temp = theta;
	for j=1:size(theta)(1)
		temp(j) = theta(j) - rate * sum((X * theta .- src_y')' * diag(X(:,j))) / (size(X)(1));
	end;
	theta = temp;
	
	% compute the cost value
	cost(it) = sum((X * theta .- src_y') .^ 2) / (2 * size(X)(1));
end;
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

迭代了500次。每次都把代价函数计算了下来，可以发现其实在大概110次左右就已经开始收敛了。这次算法是成功的，诶？为什么这么说呢？这次运气好成功了，后面我将演示在不同的收缩速率步长和迭代次数下的表现。就可能会彻底影响结果。这种情况我们一看，误差挺大的，基本没啥用，这叫做欠拟合（under-fit）。
最后最小化的cost函数最小值是2.3836，此时theta = (4.5406)。h(x) = y = 4.5406;

模型二：先来看看带有一次项的简单的线性回归
我们只有一个特征变量可得：Y = H(θ)(X) = Xθ = θ0 +θ1X1

迭代次数和相应的代价函数值的曲线
核心计算过程如下：

% model 2 : y = Xθ, θ = (θ0, θ1), X = (1, x)
X = ones(size(src_x)(2), 2);
X(:,1) = 1;
X(:,1) = src_x';
theta = zeros(size(X)(2), 1); % init value to 0

itr = 500; %iteration times
rate = 0.02; %step
cost = zeros(1, itr);
for it=1:itr
	% update each theta
	temp = theta;
	for j=1:size(theta)(1)
		temp(j) = theta(j) - rate * sum((X * theta .- src_y')' * diag(X(:,j))) / (size(X)(1));
	end;
	theta = temp;
	
	% compute the cost value
	cost(it) = sum((X * theta .- src_y') .^ 2) / (2 * size(X)(1));
end;
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

这次发现，最后最小化的cost函数最小值是0.24417，发现这个模型可能更加适合，θ = [1.1599;1.1214]，h(x) = y = 1.1599 + x * 1.1214;

模型三：先来看看带有二次项的线性回归
我们只有一个特征变量可得：Y = H(θ)(X) = Xθ = θ0 +θ1X1+θ2X1^2

迭代次数和相应的代价函数值的曲线
核心计算过程如下：

% model 2 : y = Xθ, θ = (θ0, θ1, θ2), X = (1, x, x^2)
X = ones(size(src_x)(2), 3);
X(:,1) = 1;
X(:,2) = src_x';
X(:,3) = src_x' .^ 2;
theta = zeros(size(X)(2), 1); % init value to 0

itr = 3000; %iteration times
rate = 0.003; %step
cost = zeros(1, itr);
for it=1:itr
	% update each theta
	temp = theta;
	for j=1:size(theta)(1)
		temp(j) = theta(j) - rate * sum((X * theta .- src_y')' * diag(X(:,j))) / (size(X)(1));
	end;
	theta = temp;
	
	% compute the cost value
	cost(it) = sum((X * theta .- src_y') .^ 2) / (2 * size(X)(1));
end;
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

最小点是0.018891，误差又小了好多，看来这这个模型更加适合了，注意到我换掉了迭代次数和步长，因此按照之前的迭代次数和步长，会出现“跑飞”，不收敛的情况，这些后面都会讲解到。所以这俩值也是慢慢调节的过程，我增加了迭代次数，减少了步长。这里先提一下，我也是把这俩值调节了好几次。所以，收敛的过程也是受到迭代次数和步长的影响，怎么判断收敛呢？这个是一样的，画出每次迭代次数和相应的代价函数值的曲线。这种情况叫做合适（fit）.

模型四：先来看看带有高次项的线性回归
我们只有一个特征变量可得：Y = H(θ)(X) = Xθ = θ0 +θ1X1+θ2X1^2+θ3X13

核心过程如下：

clc;clear all;close all;
set(0,'defaultfigurecolor','w') ;
src_x = 0:0.5:6;
src_y = [0, 1.1, 2.2, 3.2, 4.1, 4.7, 5.3, 5.8, 6.2, 6.5, 6.6, 6.65, 6.68];

% model 4 : y = Xθ, θ = (θ0, θ1, θ2, θ3), X = (1, x, x^2, x^3)
X = ones(size(src_x)(2), 4);
X(:,1) = 1;
X(:,2) = src_x';
X(:,3) = src_x' .^ 2;
X(:,4) = src_x' .^ 3;
theta = zeros(size(X)(2), 1); % init value to 0

itr = 150000; %iteration times
rate = 0.00015; %step
cost = zeros(1, itr);
for it=1:itr
	% update each theta
	temp = theta;
	for j=1:size(theta)(1)
		temp(j) = theta(j) - rate * sum((X * theta .- src_y')' * diag(X(:,j))) / (size(X)(1));
	end;
	theta = temp;
	
	% compute the cost value
	cost(it) = sum((X * theta .- src_y') .^ 2) / (2 * size(X)(1));
end;
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

我们看到我们新增了好多迭代次数，步长也过大了。最小的值接近到了最小值 0.018714，θ = [0.462238, 1.661541, 0.064844, -0.029070] ，足够说明了一个问题就是，最小化过程受到迭代次数，步长的影响，存在高次项的情况下，很容易不收敛，我调整了好几次都是因为，形成了发散式的迭代，代价函数不收敛，每迭代一次代价函数数值增大一次，最后整个程序爆炸，不得不降低步长，但是一旦降低步长，就得调整增大迭代次数。最后就是这样，过程挺麻烦

以上，先学这么多，线性回归远不止这么多啊。
后面的文章会讲解到，其他最小化优化算法，正则化（预防过拟合），都是很实用的。