章节4：线性回归
本章节，我们将介绍线性回归问题，机器学习中最基础的问题。
线性回归
线性回归是指在一组数据中拟合一条最合适的直线，使得数据大致落在这条直线上，考虑以下数据集：
import numpy as np
import matplotlib.pyplot as plt
# x, y
data = np.array([
[2.4, 1.7],
[2.8, 1.85],
[3.2, 1.79],
[3.6, 1.95],
[4.0, 2.1],
[4.2, 2.0],
[5.0, 2.7]
])
有了这些{x，y}集合对，就可以将它们绘制在散点图中，
x, y = data[:,0], data[:,1]
plt.figure(figsize=(4, 3))
plt.scatter(x, y)
线性回归的目标是找到一条直线：y=mx+b，使得它能很好的拟合数据。那么，用什么来判断数据点拟合的好坏程度呢？我们将比较以下三条随机拟合的直线，他们的函数方程如下：
def f1(x):
return 0.92 * x - 1.0
def f2(x):
return -0.21 * x + 3.4
def f3(x):
return 0.52 * x + 0.1
# try some examples
print("f1(-1.0) = %0.2f" % f1(-1))
print("f2( 0.0) = %0.2f" % f2(0))
print("f3( 2.0) = %0.2f" % f3(2))
f1(-1.0) = -1.92
f2( 0.0) = 3.40
f3( 2.0) = 1.14
绘制这三条直线方程，观察每条直线对数据的拟合程度，
min_x, max_x = min(x), max(x)
fig = plt.figure(figsize=(10,3))
fig.add_subplot(131)
plt.scatter(x, y)
plt.plot([min_x, max_x], [f1(min_x), f1(max_x)], 'k-')
plt.title("f1")
fig.add_subplot(132)
plt.scatter(x, y)
plt.plot([min_x, max_x], [f2(min_x), f2(max_x)], 'k-')
plt.title("f2")
fig.add_subplot(133)
plt.scatter(x, y)
plt.plot([min_x, max_x], [f3(min_x), f3(max_x)], 'k-')
plt.title("f3")
下面，我们绘制误差线，用红色的虚线表示，
min_x, max_x = min(x), max(x)
fig = plt.figure(figsize=(4,3))
plt.scatter(x, y) # original data points
plt.plot([min_x, max_x], [f1(min_x), f1(max_x)], 'k-') # line of f1
plt.scatter(x, f1(x), color='black') # points predicted by f1
for x_, y_ in zip(x, y):
plt.plot([x_, x_], [y_, f1(x_)], '--', c='red') # error bars
plt.title("error bars: $y_i-f_1(x_i)$")
# sum squared error
def cost(y_pred, y_actual):
return 0.5 * np.sum((y_actual-y_pred)**2)
x, y = data[:,0], data[:,1]
J1 = cost(f1(x), y)
J2 = cost(f2(x), y)
J3 = cost(f3(x), y)
print("J1=%0.2f, J2=%0.2f, J3=%0.2f" % (J1, J2, J3))
J1=1.18, J2=2.16, J3=0.15
正如我们所料，第三个函数在该数据集上的代价最小，第二个函数的代价最大。
有什么好的方法能帮助我们找到最佳的m和b，使得代价函数最小？最简单的方法就是暴力解决，让计算机进行数以万计的猜测，保留拥有最小代价的m和b。但是在实际的问题中，我们经常又数十，上百乃至数百万的参数需要同时优化，没有计算机能做到在合理的时间内找到最优的解决方案。
所以我们需要一个更正规的方案。或许我们能从代价函数的损失平面找到方法，我们将绘制出一定范围内m和b两个参数组合后的MSE（均方误差）函数，
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.gca(projection='3d')
# check all combinations of m between [-2, 4] and b between [-6, 8], to precision of 0.1
M = np.arange(-2, 4, 0.1)
B = np.arange(-6, 8, 0.1)
# get MSE at every combination
J = np.zeros((len(M), len(B)))
for i, m_ in enumerate(M):
for j, b_ in enumerate(B):
J[i][j] = cost(m_*x+b_, y)
# plot loss surface
B, M = np.meshgrid(B, M)
ax.plot_surface(B, M, J, rstride=1, cstride=1, cmap=plt.cm.coolwarm, linewidth=0, antialiased=False)
plt.title("cost for different m, b")
plt.xlabel("b")
plt.ylabel("m")
由此可以看到，我们的损失平面有点像一个加长型的碗，因为我们的目的是找到误差最小的参数m和b，所以，大致可以看到最低点在m=0，b=2处。
梯度下降
在这，我们介绍梯度下降来寻找最合适的参数，注意，有比梯度下降更好的方法来解决线性回归的问题。但我们引入梯度下降，是为了将其运用于后面要学习的神经网络中。所以，先介绍下梯度下降在简单问题中的运用。
梯度下降的基本概念是，先随机猜测一组参数，接着计算它们的损失函数的梯度。注意，损失函数(Loss function)是定义在单个训练样本上的，而代价函数(Cost function)是定义在整个训练集上面的，也就是所有样本的误差的总和的平均，也就是损失函数的总和的平均，所以在这是损失函数。章节2介绍过，一个函数的梯度就是包含每个变量偏导数的一个向量，如下：
怎样计算对各个参数的偏导数呢？已知，
import random
# get our data
x, y = data[:,0], data[:,1]
# it is a good idea to normalize the data
x = x / np.amax(x, axis=0)
y = y / np.amax(y, axis=0)
# choose a random initial m, b
m, b = random.random(), random.random()#0.8, -0.5
def F(x, m, b):
return m * x + b
# what is our error?
y_pred = F(x, m, b)
init_cost = cost(y_pred, y)
print("initial parameters: m=%0.3f, b=%0.3f"%(m, b))
print("initial cost = %0.3f" % init_cost)
# implement partial derivatives of our parameters
def dJdm(x, y, m, b):
return -np.dot(x, y - F(x, m, b))
def dJdb(x, y, m, b):
return -np.sum(y - F(x, m, b))
# choose the alpha parameter and number of iterations
alpha = 0.01
n_iters = 2000
# keep track of error
errors = []
for i in range(n_iters):
m = m - alpha * dJdm(x, y, m, b)
b = b - alpha * dJdb(x, y, m, b)
y_pred = F(x, m, b)
j = cost(y_pred, y)
errors.append(j)
# plot it
plt.figure(figsize=(16, 3))
plt.plot(range(n_iters), errors, linewidth=2)
plt.title("Cost by iteration")
plt.ylabel("Cost")
plt.xlabel("iterations")
# what is our final error rate
y_pred = F(x, m, b)
final_cost = cost(y_pred, y)
print("final parameters: m=%0.3f, b=%0.3f"%(m, b))
print("final cost = %0.3f" % final_cost)
initial parameters: m=0.181, b=0.682
initial cost = 0.043
final parameters: m=0.584, b=0.325
final cost = 0.009
画出找到的最合适的直线：
min_x, max_x = min(x), max(x)
fig = plt.figure(figsize=(3,3))
plt.scatter(x, y)
plt.plot([min_x, max_x], [m * min_x + b, m * max_x + b], 'k-')
plt.title("line of best fit")
上面所展示的程序就是梯度下降方法，然而对于线性回归问题可以用分析法更好的解决，我们引入梯度下降，是为了在下一章节将其运用在神经网络中。
