多变量线性回归算法

其实所谓的多变量的线性回归（Linear Regression with multiple variables ）本质上将与单变量的线性回归没啥差别。因此我们完全可以用上一节中的梯度下降算法来解决，只需要在每一次迭代的时候多考虑几个变量而已。所以这一节就稍微介绍一下了，不再用例子分析。不过毕竟多了一些变量，在对多变量跑梯度下降算法时，显然对参数的调节就更加重要了，因此我们首先得学会一些参数调节的技巧。这些技巧在实际的操作过程中尤为重要。

调节参数

特征值缩放（Feature Scaling）

我们知道在进行梯度下降算法时，经常会有这样一个问题，就是如果两个特征变量x_i和x_j的范围差别很大，那么每一次迭代时，他们下降的相对速率就不一样，这就直接导致了范围大的那个变量下降的过慢。为了解决这个问题，我们引入了特征值缩放这个技巧。其实也很简单，就是把每一个参数的变化范围大概控制在一起，这样他们的下降速率就不会受到相互的牵连。具体的操作就是将他所有参数的值，都缩放在大小差不多的一个区域内（通常是(-1,1)）。为了进行这样的缩放，我们就需要接下来的方法。

标准化（Mean Normalization）

为了缩放到一个区域（(-1,1)）内，我们要做两件事：第一，把平均值对应到0点；第二，把最大最小值对应到-1和1这两点。想想也很好理解，用公式表示就是：

这里的x_i表示需要标准化的数据，u表示在数据集中这种参数的平均值，\sigma表示数据集中这种参数的标准差。看着很眼瘦？没错，这个就是我们求正态分布的那个标准化转换函数。经过这样的转换，每一个变量都会将他的值缩放到（-1,1）中了，进而方便我们进行梯度下降。

学习因子$\alpha$的选取

我们知道学习因子既不能过大也不能过小。实际上，对于一个学习因子是好是坏我们可以这样判断。对于一个特定的学习因子\alpha，我们可以画出J函数与迭代次数的关系图。对于一个好的\alpha，他的函数图像应该是一个类似y=\frac1x的图像。如果J函数下降的速度过慢，显然我们需要增加学习因子使他下降的更快；如果J函数是增函数，或者在某些区间内有上升的趋势，那么显然我们需要减小学习因子防止学习过度。

多项式回归（Polynomial Regression ）

对于某些不能用线性回归的问题，我们有时候可以试着用多项式来进行回归拟合。其实多项式回归完全可以看成是多变量的线性回归问题，因为我们完全可以把其中的x^i看成是第i个独立的变量，只不过他的值是由x推出来的而已。原理很简单，但是如果想不到那就头大了0.0。

公式法（Normal equation）

介绍

对于多变量的线性回归，除了用我们之前学的GD算法，我们其实还有另外一个直接套公式的算法（卧槽早说）。这牵涉到线性代数的知识，我们需要做的只是将数据集组合成几个矩阵，然后运算一个公式即可，这个公式就叫 Normal equation （觉得叫成“正规方程”好难听）：

其实这个方法就是矩阵最小二乘法了。

其中X^T表示X的转置，X^{-1}表示X的逆。

这个公式，得到的结果就是对整个数据集最好的线性逼近。下面介绍下这个公式怎么用。

用法

假设我们的数据集有m个数据，我们的变量有n个（即每一个数据有n个特征值），那么我们可以把所有的输入数据写成一个m*（n+1）的矩阵。为什么是n+1列，因为我们知道，对n个参数的回归需要用到n+1个\theta。只是这当中的\theta_0一直取1罢了。因此这个矩阵的第一列全部都是1。

上面的这个矩阵，就是我们公式中的X了。

对于公式中的y，其实是由所有的输出数据组成的列向量，与输入一一对应。

最后，我们只需要计算上面这个公式的结果就能计算出\theta向量的值了，（很明显这个向量应该是一个n+1维的列向量）

至于这个公式的证明。。。貌似很麻烦，百度normal equation 可以找到需要的证明过程。

说明

这里有一个问题，就是如果X^TX没有逆元怎么办?

事实上正常情况下是不会没有逆元的，没有逆元的情况一般有两种：

1、参数中有线性相关的两组量。比如一个参数是用厘米表示的高度，另一个参数是用米表示的同样的那个高度。显然，这样选择参数是没有意义的，对于这种情况我们只要删掉其中一个参数就好了。

2、数据集的数量m小于参数的个数n。通常情况下当然是不会发生这种事的，毕竟都没几个数据还怎么好意思来拟合？不过万一出现这种状况，也是有对应的办法的，这种情况以后再说。

性能

我们注意到这里需要求矩阵的逆，而X^TX是一个n*n的矩阵，因此他的复杂度通常至少是O(n^3)（虽然可以优化，但是也差不多这么大）。这在n比较小的时候效果还挺不错的，而且比GD算法更加方便准确。但是当参数的个数变多了的时候，效率就会特别慢，相比之下GD算法的优势就显现出来了。

实现

接下来我就结合上一节的例子试着实现一下，虽然是单变量，但是道理是一样的。

import numpy as np
import theano.tensor as T
from pylab import *
import  matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

#read file

file=open("data","r")
brain=[]
body=[]
for line in file:
	num,_brain,_body=line.split()
	brain.append(_brain)
	body.append(_body)

#calculate
size=len(brain)
X=np.matrix(np.vstack((np.ones(size),brain))).reshape(2,size).T.astype(np.float)
Y=np.matrix(body).T.astype(np.float)
out=(X.T*X).I*X.T*Y
print out

#draw
out=np.array(out)
px=np.linspace(-100,600,100)
py=px*out[1]+out[0]

plt.plot(px,py)
plt.scatter(brain,body,s=100,alpha=.5)
plt.xlim(-100,600)
plt.ylim(-200,1400)
plt.show()

（这里写代码的时候要非常注意np.matrix的用法以及强制指定数据类型）

最后算得：\theta_0=55.95185719,\theta_1=1.24271276，从图上看还是拟合的挺不错的。