Machine Learning

线性模型-(一)-最小二乘

  • 线性模型形式简单，具有很好的可解释性。虽然简单却蕴藏着重要的机器学习的基本思想。一般实践中，都会先使用线性回归方法尝试。
  • 把分类问题的样本点放在在坐标轴上，坐标轴每一个轴代表一个特征。而回归问题的样本点放到坐标轴上，x轴是连续的，可以是其中一个特征的连续值，其他轴代表不同的特征。
  • 线性回归的目标是度量出模型没有拟合住样本点的部分，此时的目标函数称为 loss fucntion。
  • 而有的算法中度量的是拟合的程度，此时成目标函数为效用函数 utility function 。
  • 通过分析问题，确定问题的目标函数，通过最优化目标(最小化loss 或最大化utility)，获得机器学习模型。
  • 为参数学习，找到一组参数，这组参数可以最小化loss 或最大化utility。涉及到最优化原理或凸优化理论，如常见的牛顿法，梯度下降，模拟退火等算法。许多计算机问题如最短路径，背包问题都属于最优化问题。
  • 对于(一元或多元)线性回归，直接用最小二乘法求解。
  • 使用线性模型的前提，是假设数据与目标间由线性关系。

最小二乘法解一元线性回归

基于最小化均方误差来进行模型求解的方法为最小二乘法，即试图找一条直线，使得所有样本点到该直线的欧氏距离之和最小。具体的结果是经过数学推导，而非迭代参数学习。

数据为一维向量：

1
2
3
# 小写X表示这是一个向量
x = np.array([1,2,3,4,5,6,7,8,9,10,])
y = np.array([3,4,6,5,5,6,8,10,9,11,])

根据最小二乘法结果公式的得：

1
2
3
4
5
6
7
8
9
10
11
12
# 最小二乘法
x_bar = np.mean(x)
y_bar = np.mean(y)

fenzi = 0.0
fenmu = 0.0
for i, j in zip(x, y):
    fenzi += (i - x_bar) * (j - y_bar)
    fenmu += (i - x_bar)**2

a = fenzi / fenmu
b = y_bar - a * x_bar

绘制拟合直线：

1
2
3
4
5
6
7
# 绘制直线
y_hat = a*x + b

plt.scatter(x, y)
plt.plot(x, y_hat, color='r')
plt.axis([0, 12, 0, 12])
plt.show()

结果：

使用模型，分别预测单个值，预测一组值：

1
2
3
4
5
6
7
# 预测单个值
def predict(x):
    return a * x + b

# 预测多个值
def predict_X(X):
    return np.asarray([a * x + b for x in X])

上述计算过程是遍历每一个x和y，部分相乘后相加，这样计算的效率并不高。另一种方式是使用矩阵的乘法，即内积(Dot Product 点乘)

向量化

使用矩阵的内积，便可以用python中的向量相乘法则快速计算：

1
2
3
4
5
x_mean = np.mean(x_train)
y_mean = np.mean(y_train)

a = (x_train - x_mean).dot(y_train - y_mean) / (x_train - x_mean).dot(x_train - x_mean)
b = y_mean - a * x_mean

用内积代替了循环遍历，计算性能上会有很大提升。对于1亿的数据量，两者使用时间：

循环遍历： 55.904422
内积：    2.4315210000000036

很多时候，算法原理的数学推导，最终要能变化成向量内积的形式，很重要。

评价现行模型的性能

评价线性模型主要用MSE, RMSE, MAE:

1
2
3
4
5
6
7
8
# 公式：
mse = np.sum((y_predict - y_test)**2 / len(y_test))
rmse = math.sqrt(mse)
mae = np.sum(np.absolute(y_predict - y_test)) / len(y_test)

# sklearn：
from sklearn.metrics import mean_squared_error
from sklearn.metrics import mean_absolute_error

相对来说，RMSE比MAE好，最小化RMSE，它表示错误样本中最大的错误值相应的比较小。而在线性回归中最好的指标要数R Squared。

R Squared

这个是线性回归最常用的性能指标。其定义为：

1
2
3
r_squared = 1 - (y_predict - y_test).dot(y_predict - y_test)
               /
               (y_means - y_test).dot(y_means - y_test)

其中分子可以表示，使用我们的模型预测产生的错误。
而分母表示使用baseline 模型y=y_means 预测产生的错误。

分母的预测是不论x是多少，我都把他预测成y_means，错误率自然多。而分子是我们训练模型的实际预测，即我的模型实际拟合住样本的地方。

所以

  • r_squared 是我的模型与baseline的比较。
  • 两者相除小于等于1，
  • r_squared 越大，表示我们的模型越好
  • 当r_squared = 1，表示，我的模型没有犯任何错。
  • 当r_squared = 0，表示，我的模型性能等同于baseline。
  • 当r_squared < 0，表示，我的模型不如baseline，很可能，原数据不存在现线性关系。

当r_squared 定义中分子分母同时除以测试样本个数，r_squared还可以写成：

1
r_squared = 1 - mse(y_predict, y_test)/var(y_test)

其中var(y_test)，表示y_test 数据方差。所以r_squared 是由统计意义的。

在sklearn中使用

1
2
3
from sklearn.metrics import r2_score

r2_score(y_predict, y_test)

在sklearn中的LinearRegression模型里的 score成员方法返回的就是r_squared：

1
2
3
4
# MODEL 表示任何模型的实例
def score(x_test, y_test):
    y_predict = MODEL.predict(x_test)
    return r2_score(y_test, y_predict)

多元线性回归与Normal Equation

依旧可以使用最小二乘法得到最终解，即normal equation。但是计算机求解normal equation 的时间复杂度为O(n^3)。优点是根据数学解方程得到，不需要归一化处理。

sklearn中多元线性回归套路：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
from sklearn import datasets

boston = datasets.load_boston()
# print(boston.DESCR)  # (503, 13)
# print(boston.feature_names)

X = boston.data
y = boston.target

# 除去异常点
X = X[y < 50]
y = y[y < 50]

from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=123)

lin_reg = LinearRegression()
lin_reg.fit(X_train, y_train)

print(lin_reg.coef_)      #  得到最终模型所有参数
print(lin_reg.intercept_)  # 包括截距
print(lin_reg.score(X_test, y_test))   # 应用在测试数据集上的得分

实现normal equation：

首先在原本X_train中加入一列全为1：

1
X_b = np.hstack([np.ones((len(X_train), 1)), X_train])

根据normal equation 的公式，其中intercept_表示直线上的截距，截距与样本特征无关，或者说，截距对应的特征是常数1。coef_表示coefficient，每个特征的系数：

1
2
3
4
_theta = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y_train)  # 完全根据公式

intercept_ = _theta[0]  # 参数的第一列为截距
coef_ = _theta[1:]

有了参数，就可以用于测试样本。写成函数：

1
2
3
4
5
6
7
8
def predict(X_sample):
    assert intercept_ is not None and coef_ in not None, \
        "Must fit before predict"
    assert X_sample.shape[1] == len(coef_), \
        "the feature number of X_sample must be equal to the lenght of coef_"
    
    X_b = np.hstack([np.ones((len(X_sample), 1)), X_sample])
    return X_b.dot(_theta)

创建样本，使用模型：

1
2
3
X_sample = np.array([[0.5, 20., 4., 0., 0.57, 7., 52., 2.8, 5., 264., 13., 390., 3.16]])

lr.predict(X_sample)  # 返回模型认为的预测值

线性回归的可解释性

对所有数据解normal equation 方程得所有系数：

1
2
3
4
[ -1.05574295e-01   3.52748549e-02  -4.35179251e-02   4.55405227e-01
  -1.24268073e+01   3.75411229e+00  -2.36116881e-02  -1.21088069e+00
   2.50740082e-01  -1.37702943e-02  -8.38888137e-01   7.93577159e-03
  -3.50952134e-01]

系数的正负表示，这个特征与预测指标(如房价)是正相关还是负相关。
正值表示正相关，越大表示这个特征越大房价越高。
负值表示负相关，越大表示这个特征越大，房价越低。
系数的绝对值的大小表示了该特征对房价的影响程度。可以把特征的值从小到大排列，分析什么特征的影响大，什么特和那个的影响小：

1
2
3
4
# 按系数排序，返回排序后的索引
sort_index = np.argsort(lin_reg2.coef_)
# 使用排序后的索引，给特征名称排序 
boston.feature_names[sort_index]

比如，最大值系数表示房间数目。房间数目与房价正相关，且越大房价越高。合理。
最小系数对应一氧化氮浓度。浓度与房价成负相关，且值越大，房价越低。合理。

这就是“现行模型的可解释性”，有针对的采集更多数据，帮助决策。

注意：

1
2
np.array([1, 2, 3, 4, 5, 6, 7, 8])  
np.array([[1, 2, 3, 4, 5, 6, 7, 8]])

前者是向量，大小为： (8,)

后者是矩阵，大小为： (2, 8)

敲黑板机器学习算法实现不是目的，一个算法的优劣是将它放在特定任务中，与其他算法比较得出的。也就是说，单单训练好一个模型，还没结束，而是与其他训练好的模型一起评价，比如使用由confusion matrix 得到的指标：查准率，查全率，F-alpha，P-R曲线， ROC曲线等。