【人工智能基础】线性回归实验分析_线性回归详解实验

实验使用到的库：numpy、matplotlib、scikit-learn
实验使用的开发环境：anaconda、jupyter

一、线性回归

线性回归就是使用一个线性函数（多项式回归可以是曲线）去拟合给定的训练集，测试时，对输入的x值，返回这个线性函数的y值。最终目标是找到y=Θ₀ + Θ₁x₁ + Θ₂x₂ + …… + Θ_nx_n 函数式中的Θ₀、Θ₁、Θ₂、……、Θ_n值

由于上图中我们能看到存在一个Θ₀并不是x_i的系数，而Θ_i是从输入的标准化数据中获得的，所以在标准化数据中我们需要单独添加一列常数项1

实验步骤

1.产生数据

# 产生数据
import numpy as np
# 产生[0,2)的随机数据，作为x值
x = 2 * np.random.rand(100,1)
# 给定一个函数获得y值
y = 4 + 3 * x + np.random.randn(100,1)
1
2
3
4
5
6

2.展示产生的数据（这块只是一个数据展示，不影响整体流程）

# 绘制数据散点图图像
import matplotlib.pyplot as plt
plt.plot(x,y,'b.')
plt.xlabel('x1')
plt.ylabel('y')
plt.axis([0,2,0,15])
plt.show()
1
2
3
4
5
6
7

3.为数据添加常数项1，并求Θ

numpy.c_[矩阵1,矩阵2]可以讲矩阵1和矩阵2拼接成为一个矩阵
numpy.ones((100,1)) 可以生成一个100行1列的纯1矩阵

# 添加偏置项1
x_b= np.c_[np.ones((100,1)),x]
# 使用公式Θ = (xT ⋅ x)-1 ⋅ xT ⋅ y 求Θ 
# xT表示矩阵x的转置
# (xT ⋅ x)-1 表示 xT ⋅ x的逆
theta_best = np.linalg.inv(x_b.T.dot(x_b)).dot(x_b.T).dot(y)
1
2
3
4
5
6

4.绘制函数图像

# 选取两个点（函数的两个端点）
x_new = np.array([[0],[2]])

# 拼接偏置项1
x_new_b = np.c_[np.ones((2,1)),x_new]

# 获得预测值
y_predict = x_new_b.dot(theta_best)

# 数据绘制
plt.plot(x_new,y_predict,'r--')
plt.plot(x,y,'b.')
plt.xlabel('x1')
plt.ylabel('y')
plt.axis([0,2,0,15])
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

二、机器学习下的线性回归

在上面的实验中，我们直接使用公式Θ = (xT ⋅ x)-1 ⋅ xT ⋅ y得到了Θ值，但是在实际中我们的机器需要通过如梯度下降之类的算法，逐步减少误差，获取到Θ的值，这个获得Θ的过程就是学习。

scikit-learn给我们提供了线性回归的工具

from sklearn.linear_model import LinearRegression
lin_reg = LinearRegression()
lin_reg.fit(x,y)
print(lin_reg.coef_)
print(lin_reg.intercept_)
1
2
3
4
5

通过这段代码我们最终能得到：

• lin_reg.coef_,学习结束后得到的Θ值，即每个特征值
• lin_reg.intercept_,函数的截距值，即Θ₀

三、学习率的影响

我们知道。梯度下降就是沿着误差函数对Θ的偏导的负方向找到最小值Θ_imin

求导公式为 2/m * (Θ^T ⋅ x - y) ⋅ x_i
其中：Θ^T 表示Θ的转置
注：Θ 和 x 在代码中都是矩阵形式，使用Θ^T ⋅ x 就相当于对Θ_ix_i求和。

梯度下降实现

# 学习率
alpha = 0.1
# 迭代次数
n_iter = 1000
# 样本个数
count = 100
# 初始化theta位置
theta = np.random.randn(2,1)
for iter in range(n_iter):
    # 对thetaJ求导 = 2/m * 求和(thetaT * x - y) * xJ
    gradients = 2 / count * x_b.T.dot(x_b.dot(theta) - y)
    # 更新theta
    theta = theta - alpha * gradients
print(theta)
1
2
3
4
5
6
7
8
9
10
11
12
13
14

最后theta是一个2*1的矩阵存储着Θ₀和x对应的权重Θ

定义学习率影响因素的函数

# 这个代码也是全量梯度下降的代码（全部的Θ值都参与梯度下降计算）
# theta_path_bgd存储每次梯度下降得到的Θ值
theta_path_bgd = []
def plot_gradigent_descent(theta, alpha, theta_path = None):
    count = len(x_b)
    plt.plot(x,y,'b.')
    n_iter = 1000
	# 进行1000次迭代计算
    for iter in range(n_iter):
        y_predict = x_new_b.dot(theta)
        plt.plot(x_new,y_predict,'r--')
		# 计算梯度
        gradients = 2 / count * x_b.T.dot(x_b.dot(theta) - y)
		# 新的Θ值为 旧Θ值 - 学习率 * 梯度
        theta = theta - alpha * gradients
        if theta_path is not None:
            theta_path.append(theta)
	# 绘制图像
    plt.xlabel('x_1')
    plt.axis([0,2,0,15])
    plt.title('alpha={}'.format(alpha))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

调用函数绘制图像

theta = np.random.randn(2,1)
plt.figure(figsize=(10,4))
plt.subplot(131)
plot_gradigent_descent(theta,alpha = 0.02)
plt.subplot(132)
# 存储全量梯度下降的过程
plot_gradigent_descent(theta,0.1,theta_path_bgd)
plt.subplot(133)
plot_gradigent_descent(theta,alpha = 0.5)
plt.show()
1
2
3
4
5
6
7
8
9
10

从图中可以看出：

• 学习率过小（0.02）可以获得正确的Θ值，但是每次学习效率极低
• 学习率过大（0.5）虽然效率高，但是计算结果有误

四、不同梯度下降的对比

全量梯度下降已经在学习率中做过了

随机梯度下降

每次去一个Θ值进行梯度下降

theta_path_sgd = []
count = len(x_b)
# 迭代次数
n_epochs = 50

t0 = 5
t1 = 50

# 学习率修改
def learning_schedule(t) :
    return t0 / (t1 + t)

theta=np.random.randn(2,1)

for epoch in range(n_epochs):
    # 遍历样本
    for i in range(count):
        # 只画第一批前二十个
        if epoch == 0 and i < 20:
             y_predict = x_new_b.dot(theta)
             plt.plot(x_new,y_predict,'r--')
        random_index = np.random.randint(count)
		# 选中一个值
        xi = x_b[random_index:random_index + 1]
        yi = y[random_index:random_index + 1]
		# 对选中值进行梯度下降
        gradients = 2 * xi.T.dot(xi.dot(theta) - yi)
        alpha = learning_schedule(epoch * count + i)
        theta = theta - alpha * gradients
        theta_path_sgd.append(theta)
plt.plot(x,y,'b.')
plt.axis([0,2,0,15])
plt.show()
print(theta)
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
31
32
33
34

小批量梯度下降

即选中一批数据进行梯度下降

theta_path_mgd = []
count = len(x_b)
# 迭代次数
n_epochs = 50
# 一批数量
minibatch = 16

theta=np.random.randn(2,1)

t = 0
for epoch in range(n_epochs):
    # 洗牌
    shuffled_indices = np.random.permutation(count)
    x_b_shuffled = x_b[shuffled_indices]
    y_shuffled = y[shuffled_indices]
    for i in range(0,count,minibatch):
        t += 1
		# 获取一批数据
        xi = x_b_shuffled[i:i + minibatch]
        yi = y_shuffled[i:i + minibatch]
		# 对获取到的这批数据进行梯度下降
        gradients = 2 / minibatch * xi.T.dot(xi.dot(theta) - yi)
        alpha = learning_schedule(t)
        theta = theta - alpha * gradients
        theta_path_mgd.append(theta)

print(theta)
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

三种梯度下降方式的对比

数据处理

theta_path_bgd = np.array(theta_path_bgd)
theta_path_sgd = np.array(theta_path_sgd)
theta_path_mgd = np.array(theta_path_mgd)
1
2
3

图像绘制

plt.figure(figsize=(8,4))
plt.plot(theta_path_sgd[:,0],theta_path_sgd[:,1],'r-s',linewidth=1,label='SGD')
plt.plot(theta_path_mgd[:,0],theta_path_mgd[:,1],'b-+',linewidth=2,label='MGD')
plt.plot(theta_path_bgd[:,0],theta_path_bgd[:,1],'g-o',linewidth=3,label='BGD')
plt.legend(loc='upper left')
plt.axis([2,4,1,5.0])
plt.show()
1
2
3
4
5
6
7

可以看到：

• 全量梯度下降一直在向目标逼近
• 小批量梯度下降略有波动
• 随机梯度下降的波动很剧烈

五、多项式回归（非线性回归）

在实际情况中，我们得到的数据不一定能通过一次函数拟合，需要更高次的函数才能更好的拟合训练数据

创建数据

# 多项式回归
count = 100
# x(-3,+3)
x = 6 * np.random.rand(count,1) - 3
# y = 0.5 * x^2 + x
y = 0.5 * x**2 + x + np.random.randn(count,1)

plt.plot(x,y,'r.')
plt.show()
1
2
3
4
5
6
7
8
9

数据处理

PolynomialFeatures可以把我们的输入的数据处理为更高幂次的形式

from sklearn.preprocessing import PolynomialFeatures
poly_features = PolynomialFeatures(degree = 2, include_bias = False)

x_poly = poly_features.fit_transform(x)

print(x[0])
print(x_poly[0], x[0][0]**2)
1
2
3
4
5
6
7

这里我们把原本的x处理成x和x^2的形式

二次函数图像绘制

numpy.linspace(-3,3,100)生成[-3,3]的100个等距数值数组

from sklearn.linear_model import LinearRegression
lin_reg = LinearRegression()
lin_reg.fit(x_poly, y)

x_new = np.linspace(-3,3,100).reshape(100,1)
# 多项式转换
x_new_poly = poly_features.transform(x_new)
y_new = lin_reg.predict(x_new_poly)

plt.plot(x,y,'b.')
plt.plot(x_new,y_new,'r--',label='predcition')
plt.legend()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13

六、模型复杂度的影响

这里的复杂度就是指模型的最高幂次

Pipeline用于将多个数据处理步骤和机器学习模型串联起来，形成一个有序的工作流程，即“机器学习流水线”

StandardScaler用于数据的标准化

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

plt.figure(figsize=(12,6))

# 分别绘制50次，2次、1次的函数图像拟合训练数据
for style, width, degree in (('g-',1,50),('b--',2,2),('r-+',3,1)):
    poly_features = PolynomialFeatures(degree = degree, include_bias = False)
    std = StandardScaler()
    lin_reg = LinearRegression()
    pr= Pipeline([
        ('poly_features', poly_features),
        ('StandardScaler', std),
        ('regression', lin_reg)
    ])
    pr.fit(x,y)
    y_new_2 = pr.predict(x_new)
    plt.plot(x_new,y_new_2,style, label = str(degree), linewidth=width)
plt.plot(x,y,'b.')

plt.axis([-3,3,-5,10])

plt.legend()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

从下图中可以看到高复杂度的模型试图拟合更多的数据导致波动极大，在测试集中的表现可能还不如低复杂度的模型，这种情况称为过拟合

七、样本数量的影响（过拟合的避免）

样本量影响实验

def plot_learning_curves(model,x,y):
    x_train,x_val,y_train,y_val = train_test_split(x,y,test_size=20,random_state=0)
	# 存储训练方差和实际方差
    train_errors,val_errors = [],[]
    # 取1个 到 count个 样本量
    for m in range(1, len(x_train)):
        model.fit(x_train[:m],y_train[:m])
		# 截取样本
        y_train_predict = model.predict(x_train[:m])
        y_val_predict = model.predict(x_val)
		# 存储结果
		train_errors.append(mean_squared_error(y_train[:m],y_train_predict[:m]))
        val_errors.append(mean_squared_error(y_val,y_val_predict))

    plt.plot(np.sqrt(train_errors),'r-+',linewidth=2,label='train_error')
    plt.plot(np.sqrt(val_errors),'b-',linewidth=3,label='val_error')
    plt.legend()
	
lin_reg = LinearRegression()
plot_learning_curves(lin_reg,x,y)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

在低样本数据量的时候两者的差距很大

而大样本量的时候两个差距会缩小

多项式的过拟合风险

# 多项式回归的拟合风险
# 使用10次函数尝试拟合
pr= Pipeline([
        ('poly_features', PolynomialFeatures(degree = 10, include_bias = False)),
        ('regression', LinearRegression())
    ])
plot_learning_curves(pr,x,y)
plt.axis([0,80,0,45])
plt.show()
1
2
3
4
5
6
7
8
9

八、正则化

正则化的目的就是当两个模型对训练集给出的结果基本相同时，选取更稳定的模型
正则化通过对误差函数增加一个带权重的额外的值（岭回归中使用Θ^2求和，Lasso使用Θ的绝对值）

岭回归

from sklearn.linear_model import Ridge
np.random.seed(42)
count = 20
x = 3 * np.random.rand(count,1)
y = 0.5 * x + np.random.randn(count,1) / 1.5 + 1
x_new = np.linspace(0,3,100).reshape(100,1)

# alpha在这里表示正则的权重，越大则倾向于越稳定的模型
def plot_model(model_class, polynomial, alphas, ** model_kargs):
    # alphas[0]-’b-‘:alphas[0]-'g--':alphas[2]-’r.‘
    for alpha,style in zip(alphas,('b-','g--','r.')):
        model = model_class(alpha,**model_kargs)
        if polynomial:
            model = Pipeline([
                ('poly_features', PolynomialFeatures(degree = 10, include_bias = False)),
                ('StandardScaler', StandardScaler()),
                ('regression', model)
            ])
        model.fit(x,y)
        y_new_regul = model.predict(x_new)
        width = 2 if alpha > 0 else 1
        plt.plot(x_new,y_new_regul,style,linewidth = width, label='alpha={}'.format(alpha))
    plt.plot(x,y,'b.',linewidth=3)
    plt.legend()
plt.figure(figsize=(10,5))
plt.subplot(121)
plot_model(Ridge,polynomial=False,alphas = (0,10,100))
plt.subplot(122)
plot_model(Ridge,polynomial=True,alphas = (0,10**-5,1))
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