西电-机器学习-线性回归

线性回归

回归问题是非常常见的一类问题，目的是寻找变量之间的关系。比如要从数据中寻找房屋面积与价格的关系，年龄和身高的关系，气体压力和体积的关系等等。而机器学习要做的正是要让机器自己来学习这些关系，并为对未知的情况做出预测。

对于线性回归，假设变量之间的关系是线性的，也即：

$$h_{\theta }(x)=mx+b$$
其中 $\theta$ 就是学习算法需要学习的参数,在线性回归的问题上，就是$m$和$b$，而 $x$ 是我们对于问题锁选取的特征，也即输入。$h$表示算法得到的映射。

代价函数的表示

为了找到这个算法中合适的参数，我们需要制定一个标准。一般而言算法拟合出来的结果与真实的结果误差越小越好，试想一下如果算法拟合出来的结果与真实值的误差为零，那么就是说算法完美地拟合了数据。所以可以根据“真实值与算法拟合值的误差”来表示算法的“合适程度”。在线性回归中，我们经常使用最小二乘的思路构建平方差函数：

$$J(m,b)=\sum _{i=1}^n[y_{\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{e}}-y_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}{]}^2$$

这里 $y_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}$ 有算法预测得出，对线性回归，为 $mx+b$，这样得到误差函数为：
$$J(m,b)=\sum _{i=1}^n[y_{\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{e}}-(mx+b){]}^2$$

误差函数的值越小，则代表算法拟合结果与真实结果越接近。

梯度下降

梯度下降算法沿着误差函数的反向更新$\theta$的值，知道代价函数收敛到最小值。梯度下降算法更新${\theta }_i$的方法为：

$${\theta }_i={\theta }_i-\alpha \frac{\partial }{\partial {\theta }_i}J(\theta )$$

其中 $\alpha$表示学习率。对于线性回归的的参数，可以根据代价函数求出其参数更新公式：

$$\frac{\partial J}{\partial m}=-\frac{2}{n}\sum _{i=1}^n{x}_i(y_{\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{e}}-(m{x}_i+b))$$

$$\frac{\partial J}{\partial b}=-\frac{2}{n}\sum _{i=1}^n(y_{\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{e}}-(m{x}_i+b))$$

$$m=m+\alpha \frac{\partial \mathrm{l}\mathrm{o}\mathrm{s}\mathrm{s}}{\partial m}$$

$$b=b+\alpha \frac{\partial \mathrm{l}\mathrm{o}\mathrm{s}\mathrm{s}}{\partial b}$$

代码实现

现在让我们开始动手实现，首先让我们回顾一下numpy和matplotlib：

import numpy as np
import matplotlib.pyplot as plt
# matplotlib inline


def warm_up_exercise():
    """热身练习"""
    A = None
    A = np.eye(5,5)
    return A

# 当你的实现正确时，下面会输出一个单位矩阵：
print(warm_up_exercise())

[[1. 0. 0. 0. 0.]
 [0. 1. 0. 0. 0.]
 [0. 0. 1. 0. 0.]
 [0. 0. 0. 1. 0.]
 [0. 0. 0. 0. 1.]]

你需要实现绘制数据集中图像的函数，当你的实现|正确时，你应该会得到如下的图像：

def plot_data(x, y):
    """绘制给定数据x与y的图像"""
    plt.figure()
    # ====================== 你的代码 ==========================
    # 绘制x与y的图像
    # 使用 matplotlib.pyplt 的命令 plot, xlabel, ylabel 等。
    # 提示：可以使用 'rx' 选项使数据点显示为红色的 "x"，
    #       使用 "markersize=8, markeredgewidth=2" 使标记更大
    # 给制数据
    plt.plot(x,y,'rx',markersize = 8,markeredgewidth=2)
    # 设置y轴标题为 'Profit in $10,000s'
    plt.ylabel("Profit in $10,000s")
    # 设置x轴标题为 'Population of City in 10,000s'
    plt.xlabel("Population of City in 10,000s")
    # =========================================================
    # plt.show()
    
# 让我们测试一下你的实现是否正确
# 从txt中加载数据
print('Plotting Data ...\n')
data = np.loadtxt('./data/data5984/PRML_LR_data.txt', delimiter=',')
x, y = data[:, 0], data[:, 1]

# 绘图
plot_data(x, y)
plt.show()

Plotting Data ...

现在运用所学的知识，对上述数据利用线性回归进行拟合。首先我们对要学习的参数和数据做一个准备:

# Add a column of ones to x
m = len(y)
X = np.ones((m, 2))
X[:, 1] = data[:, 0]

# initialize fitting parameters
theta = np.zeros((2, 1))

# Some gradient descent settings
iterations = 1500
alpha = 0.01

计算初始误差函数的值，你需要实现误差函数的计算：

def compute_cost(X, y, theta):
    """计算线性回归的代价。"""
    m = len(y)
    J = 0.0
    # ====================== 你的代码 ==========================
    # 计算给定 theta 参数下线性回归的代价
    # 请将正确的代价赋值给 J
    # sum = 0.0
    # for i in range(m):
    #     sum += (y[i]-(theta[1]*X[i,1]+theta[0]))**2
    # J = sum/(2*m)
    
    J = 1/(2*m) * sum((y-(theta[1]*X[:,1]+theta[0]))**2)
    # =========================================================
    return J

# compute and display initial cost
# Expected value 32.07
J0 = compute_cost(X, y, theta)
print(J0)

32.072733877455654

现在你验证了代价计算的正确性，接下来就需要实现最核心的部分：梯度下降。在实现这一部分之前，确定你理解了上述各种变量及其表示。你需要完成梯度下降的核心代码部分：

def gradient_descent(X, y, theta, alpha, num_iters):
    """执行梯度下降算法来学习参数 theta。"""
    m = len(y)
    J_history = np.zeros((num_iters,))

    for iter in range(num_iters):
        # ====================== 你的代码 ==========================
        # 计算给定 theta 参数下线性回归的梯度，实现梯度下降算法
        temp_m = -2/m * sum( X[:,1] * (y-(theta[1]*X[:,1]+theta[0])) )
        temp_b = -2/m * sum(y-(theta[1]*X[:,1]+theta[0]))
        
        theta[1] -= alpha * temp_m
        theta[0] -= alpha * temp_b
        # =========================================================
        # 将各次迭代后的代价进行记录
        J_history[iter] = compute_cost(X, y, theta)
    return theta, J_history

# run gradient descent
# Expected value: theta = [-3.630291, 1.166362]
theta, J_history = gradient_descent(X, y, theta,
                                    alpha, iterations)

print(theta)
plt.figure()
plt.plot(range(iterations),J_history)
plt.xlabel("iterations")
plt.ylabel(r'$J(\theta)$')
plt.show()

[[-3.87813769]
 [ 1.19126119]]

为了验证梯度下降方法实现的正确性，你需要把学习的到的直线绘制出来，确定你的实现是否正确。前面你已经绘制了数据集中的点，现在你需要在点的基础上绘制一条直线，如果你的实现正确，那么得到的图像应该是如下这样：
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-QGBws3VX-1584200136136)(./fit.png)]
现在你已经正确实现了线性回归，你可能会对误差函数的优化过程比较好奇。为了更好地理解这个过程，你可以将损失函数的图像绘制出来。为此你需要将需要优化的参数的各个取值时误差函数的取值在图像上绘制出来，以下代码需要你进行填写。

plot_data(x, y)
y_predict = theta[1]*x+theta[0]
plt.plot(x,y_predict)
plt.show()

def plot_visualize_cost(X, y, theta_best):
    """可视化代价函数"""

    # 生成参数网格
    theta0_vals = np.linspace(-10, 10, 101)
    theta1_vals = np.linspace(-1, 4, 101)
    t = np.zeros((2, 1))
    J_vals = np.zeros((101, 101))
    for i in range(101):
        for j in range(101):
            # =============== 你的代码 ===================
            # 加入代码，计算 J_vals 的值
            J_vals[i][j] = compute_cost(X,y,[theta0_vals[i],theta1_vals[j]])
            
            # ===========================================

    plt.figure()
    plt.contour(theta0_vals, theta1_vals, J_vals,
                levels=np.logspace(-2, 3, 21))
    plt.plot(theta_best[0], theta_best[1], 'rx',
             markersize=8, markeredgewidth=2)
    plt.xlabel(r'$\theta_0$')
    plt.ylabel(r'$\theta_1$')
    plt.title(r'$J(\theta)$')
    plt.show()
    

plot_visualize_cost(X, y, theta)
plt.show()

进阶

在实现中，你可能采取了像上面公式中给出的结果一样逐个样本计算代价函数，或者在梯度下降的更新时也采用了逐个样本计算的方式。但事实上，你可以采用numpy的矩阵函数一次性计算所有样本的代价函数。可以采用矩阵乘法(np.matmul())求和等方式（np.sum（））。利用你学到的线性代数知识，将其实现更改一下吧。

在梯度更新时，我们保留了代价的历史信息。在参数的学习过程中，代价函数的变化过程你也可以作一个图来查看。观察最后得到的$J(\theta )$的图像以及代价的变化过程，可以加深你的理解。在梯度下降的迭代中，我们设置终止条件为完成了固定的迭代次数，但是在迭代次数完成时，由于学习率等参数的设置，可能得到的参数并不是使得代价最低的值。你可以通过观察代价函数的变化过程，想办法调整学习率等参数或者改进程序，使得参数的取值为搜索到的最优结果。