第3课 多变量线性回归
多变量线性回归问题中，训练集一般包含n个特征$(x_1,x_2,\cdots,x_n)$，标签依然用$y$表示。
用$x^{(i)}$表示训练集中第$i$个样本的特征，用$x_j^{(i)}$表示第$i$个样本的第$j$个特征$(j=1,2,3,\cdots,n)$。
此时假设函数可以写为：
$h(x)=\theta_0+\sum_{i=1}^n\theta_ix_i$
写成向量的形式，令$x_0=1$，$\vec{x}\in\mathbb{R}^{n+1}$,$\vec{\theta}\in\mathbb{R}^{n+1}$
\(x=\begin{bmatrix}1\\ x_1\\ x_2\\ \vdots\\ x_n \end{bmatrix}\) \(\qquad \theta=\begin{bmatrix} \theta_0\\ \theta_1\\ \theta_2\\ \vdots\\ \theta_n \end{bmatrix}\)
假设函数可以写为：$h(x)=\theta^Tx$
利用梯度下降算法优化参数：
代价函数可以写为：
\[J(\theta_0,\theta_1,\cdots,\theta_n)=\frac{1}{2m}\sum_{i=1}^m(\theta^Tx^{(i)}-y^{(i)})^2\]迭代更新参数按照下面的法则：
\[\theta_j:=\theta_j-\frac{\alpha}{m}\sum_{i=1}^m(\theta^Tx^{(i)}-y^{(i)})x^{(i)}_j\qquad(j=0,1,\cdots,n)\]向量化操作的写法为：$\theta\in\mathbb{R}^n$,$X\in\mathbb{R}^{m\times n}$，$y\in\mathbb{R}^m$
\[\theta := \theta-\frac{\alpha}{m}(X\theta-y)^TX\]总体而言，处理多变量线性回归问题，可以以单变量线性回归的求解思路为基础，但是在实践过程中，经常会遇到不同特征之间数量级相差很大的情况，例如描述房屋的特征中，房间数量可能是个位数，而房屋总面积可能三位甚至四位数。这时候就要对特征进行缩放操作。
通过特征缩放，能够将所有特征变换到统一数量级，从而加快计算收敛的速度。
实施特征缩放可以利用以下两种思路：
- 对某一特征，找出所有样本中该特征的最大值，然后用所有样本的该特征值除以最大值；
-
或者利用均值归一化方式处理。如：
\(\frac{x^{(i)}_{j}-\mu}{std(x_j)}\) 或者 \(\frac{x^{(i)}_j-\mu}{x_{jmax}-x_{jmin}}\)
# 多变量线性回归的代码实现
# 载入数据
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
data = pd.read_csv('ex1data2.txt',header=None)
# 给数据框的列命名
data.columns = ['Area','RoomNums','Price']
# 提取特征矩阵和标签向量
X = data.loc[:,['Area','RoomNums']]
y = data.loc[:,['Price']]
m = len(y)
# 将数据框转换为数组
X = X.values
y = y.values
'''
================ 第一部分：特征缩放 ================
'''
# 两个特征之间的数量级相差较大，需要进行特征缩放。
X[1:10,:]
# 编写特征缩放函数
def FeatureNormalization(X):
'''
该函数用于对特征矩阵进行缩放
函数返回缩放后的特征矩阵，特征缩放前的各特征的平均值和标准差。
经过缩放后，每一个特征的均值变为0，标准差变为1.
'''
# 初始化变量
X_norm = X
mu = np.zeros(X.shape[1])
sigma = np.zeros(X.shape[1])
'''
首先，计算每一个特征的均值，将其存储在向量mu中，
然后，计算每一个特征的标准差，将其存储在向量sigma中。
需要注意，X的每一个特征均是按列存储，每一行代表一个样本。
'''
mu = np.mean(X, axis=0)
sigma = np.std(X,axis=0)
X_norm = (X - mu) / sigma
return X_norm, mu,sigma
# 对特征矩阵进行缩放
X_norm,mu,sigma = FeatureNormalization(X)
# 首先给X增加一列全1向量，作为x0
X = np.insert(X_norm, 0, values=1, axis=1)
'''
================ 第2部分：梯度下降 ================
在这一部分中，需要完成代价函数computeCost和梯度下降函数gradientDescent
然后，通过调整学习速率alpha来获得最好的模型训练结果。
最后，预测一下一幢1650平方英尺，3个房间房屋的价格。
'''
# 选择不同的学习速率
alpha = np.array([0.01,0.03,0.1,0.3])
num_iters = 400
# 代价函数
def computeCostMulti(X, y, theta):
'''
计算多变量线性回归模型的代价函数，并返回误差值J
J = computeCostMulti(X, y, theta) 利用特征矩阵X、标签向量y以及参数向量theta计算模型的误差
'''
#初始化变量
m = len(y) # 训练样本数量
J = 0 # 初始化误差值
J = np.sum((np.dot(X,theta) - y)**2)/2/m
return J
# 梯度下降函数
def gradientDescentMulti(X, y, theta, alpha, num_iters):
'''
本函数利用梯度下降算法优化参数向量theta，
theta,J_history = gradientDescentMulti(X, y, theta, alpha1, num_iters)通过num_iters次迭代
返回优化后的theta向量，并将每一次迭代的误差，记录在J_history中。
'''
# 初始化变量
m = len(y)
J_history = np.zeros((num_iters,))
# 迭代
for iter in range(num_iters):
"""
对参数数据框执行单步迭代。
"""
theta = theta - alpha / m * np.dot((np.dot(X,theta) - y[:,0]).T,X).T
J_history[iter] = computeCostMulti(X, y, theta)
return theta,J_history
# 初始化Theta，变换学习速率，然后利用梯度下降算法计算参数
theta = np.zeros((X.shape[1],len(alpha)))
J_history = np.zeros((num_iters,len(alpha)))
for i in range(len(alpha)):
theta[:,i],J_history[:,i] = gradientDescentMulti(X, y, theta[:,i], alpha[i], num_iters)
# 绘制J_history的变化图
from pylab import *
plt.rcParams['font.sans-serif'] = ['Arial Unicode MS'] # 使图形能够显示中文
plt.rcParams['axes.unicode_minus'] = False
labels = [0 for i in range(len(alpha))]
for i in range(len(alpha)):
labels[i] = 'alpha = ' + str(alpha[i])
plt.plot(np.arange(num_iters),J_history[:,i],label=labels[i])
plt.xlabel('迭代次数')
plt.ylabel('代价函数')
plt.title(u'代价函数迭代次数的变化趋势')
plt.legend()
# 假设一幢房屋的面积为1650，房间数量为3，预测其价格。
x_new = np.array([1650,3])
x_new = (x_new - mu) / sigma # 特征缩放
x_new = np.insert(x_new,0,values=1,axis=0) # 因为x0=1，所以可以先进行特征缩放，然后再插入x0列。
# 依据不同的学习速率得到的参数，能够得到不同的房价预测
price = np.dot(x_new,theta)
labels = [0 for i in range(len(alpha))]
# for i in range(len(alpha)):
# labels[i] = 'alpha = ' + str(alpha[i])
plt.bar(range(len(alpha)),price,tick_label=list(alpha))
plt.xlabel('学习速率')
plt.ylabel('预测房价')
plt.title('不同学习速率所预测房价')
plt.legend()
正规方程
除了采用梯度下降算法求解参数的方法之外，对于线性回归问题，还可以用正规方程的形式，直接利用解析式求解$\theta$。
其分析原理如下：
假设参数$\vec{\theta}\in\mathbb{R}^{n+1}$，线性回归问题的代价函数为：
$J(\vec{\theta})=\frac{1}{2m}\sum(X\theta-y)^2$
可以看到代价函数是关于$\theta$的一元二次方程。对于一般的一元二次方程
$J(\theta)=a\theta^2+b\theta+c$
并且已知函数$J(\theta)$为凸函数，则为了使$J(\theta)$取极值，$\frac{\mathrm{d} J(\theta)}{\mathrm{d} \theta} = 0$
$\therefore \theta=\frac{-b}{2a}$
以上是用解析法求解参数的基本原理，在具体实践中，以多变量线性回归问题为例，若$\theta\in\mathbb{R}^{n+1}$，$X\in\mathbb{R}^{m\times n+1}$，$y\in\mathbb{R}^m$
那么$\theta = (X^TX)^{-1}X^Ty$
方程中的$X^TX$项，在计算过程中可能会存在奇异矩阵的情况：
- 如果特征矩阵中两特征存在相关性
- 假如特征数量非常多，$n\geq m$
利用正规方程求解参数相对来说比较直观，但是因为要进行矩阵的乘法和求逆运算，因此在面对大量特征时，其计算速度会较为缓慢，在这种情况下，梯度下降算法仍然是优先考虑的对象。
'''
================ 正规方程的代码实现 ================
'''
# 载入数据
import pandas as pd
import numpy as np
data = pd.read_csv('ex1data2.txt',header=None)
X = data.iloc[:,0:-1]
y = data.iloc[:,-1]
X = X.values
y = y.values
m = len(y)
X = np.insert(X,0,values=1,axis=1)
# 正规方程计算函数
def normalEqn(X,y):
theta = np.zeros((X.shape[1],1))
theta = np.dot(np.dot(np.linalg.inv(np.dot(X.T,X)),X.T),y)
return theta
# 根据训练样本计算参数
theta = normalEqn(X, y)
# 预测1650平方英尺，3个房间的房屋价格。
price = np.dot([1,1650,3],theta)
print('面积为1650平方英尺，3个房间的房屋价格为 %.2f 美元'%(price))
最后，在线性回归模型的基础上，还可以进行更复杂的多项式回归假设。例如：
- 二次多项式：$\theta_0+\theta_1x+\theta_2x^2$
- 三次多项式：$\theta_0+\theta_1x+\theta_2x^2+\theta_3x^3$
- 平方根多项式：$\theta_0+\theta_1x+\theta_2\sqrt{x}$
在利用多项式进行回归问题处理时，更要注意对特征进行缩放。
