机器学习实战：回归

线性回归基本问题

线性回归模型函数与损失函数

• 问题
  m个样本数据，每个样本对应n个特征及一个输出，针对样本输入$x_i$，其输出预测$\widehat{y_i}$如何计算？
• 模型函数
  $$h_{\theta }(x_0,x_1,...,x_n)=\sum _{i=0}^n{\theta }_i\ast x_i$$
  其中$x_0=1$.写成矩阵形式，$h_{\theta }(X):(m,1)$、$X:(m,n)$、$\theta :(n,1)$形矩阵
  $$h_{\theta }(X)=\sum _{i=0}^{n}x\ast \theta$$
• 损失函数
  $$J({\theta }_0,{\theta }_1,...,{\theta }_n)=\sum _{i=1}^m(h_{\theta }(x_0,x_1,...,x_n)-y_i{)}^2$$$$J(\theta )=\frac{1}{2}(X\theta -Y{)}^T\ast (X\theta -Y)$$

线性回归算法

• 梯度下降
  $$\theta =\theta -\alpha X^T(X\theta -Y)$$

• 最小二乘法
  $$\theta =(X^{T}X{)}^{-1}{X}^{T}Y$$

• 判断回归效果
  计算相关系数R:预测序列$\hat{y}$与真实值序列$y$的匹配程度

标准线性回归

• 底层实现

数据集：ex0.txt

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.font_manager import FontProperties

font = FontProperties(fname=r"C:\Windows\Fonts\STFANGSO.TTF")
file_name = "D:\深度学习\机器学习\8回归\ex0.txt"


"""
函数名：load_dataset
函数说明：加载数据
input:
    file_name:数据路径
output:
    data_mat:特征集
    label_mat:标签集
"""
def load_dataset(file_name):
    numFeat = len(open(file_name).readline().split('\t'))-1  # 特征数
    dataMat = []; labelMat = []
    fr = open(file_name)
    for line in fr.readlines():
        linarr = []
        currentline = line.strip().split('\t')
        for i in range(numFeat):
            linarr.append(float(currentline[i]))  # float必须的，否则数据为str类型
        dataMat.append(linarr)
        labelMat.append(float(currentline[-1]))
    return dataMat,labelMat

def standRegression(dataMat,labelMat):
    data_x = np.mat(dataMat)
    data_y = np.mat(labelMat).T
    xTx = data_x.T * data_x
    if np.linalg.det(xTx) == 0.0:
        print('矩阵不可逆')
    ws = xTx.I*(data_x.T*data_y)
    return ws

if __name__ == '__main__':
    dataMat,labelMat = load_dataset(file_name)
    dataArr = np.array(dataMat)
    labelArr = np.array(labelMat)
    dataArr = dataArr[:,1]
    ws = standRegression(dataMat,labelMat)
    print(np.shape(ws))
    yHat = (np.mat(dataMat)*ws).flatten().A[0]
    fig = plt.figure()
    axs = fig.add_subplot(111)
    axs.scatter(dataArr,labelArr,label='真实')
    axs.plot(dataArr,yHat,label="预测")
    plt.xlabel("x")
    plt.ylabel('y')
    plt.legend(prop=font)
    plt.show()

• scikit-learn实现
  数据集：Folds5x2_pp.csv

from sklearn.linear_model import LinearRegression
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.model_selection import cross_val_predict
from sklearn.model_selection import cross_val_score
import pandas as pd
import numpy as np
from sklearn import metrics
data = pd.read_csv("D:\深度学习\机器学习\8回归\Folds5x2_pp.csv" )
# 观察数据是否有缺失
print(data.info())
print(data.head())
print(data.tail())

X = data[['AT','V','AP','RH']]
y = data['PE']
print(X.head())

# 训练集测试集划分
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=1)
print(np.shape(X_train))
print(np.shape(X_test))

# 定义模型
lin1 = LinearRegression()
# 训练集上训练数据集
lin1.fit(X_train, y_train)
w0 = lin1.intercept_
w = lin1.coef_
print(w0)
print(w)

# 模型我们需要评估我们的模型的好坏程度，对于线性回归来说，
# 我们一般用均方差（Mean Squared Error, MSE）或者均方根差(Root Mean Squared Error, RMSE)
# 在测试集上的表现来评价模型的好坏
# 在测试集上拟合
predicted = lin1.predict(X_test)
MSE = metrics.mean_squared_error(predicted, y_test)
RMSE = np.sqrt(metrics.mean_squared_error(predicted, y_test))
print('MSE:',MSE)
print('RMSE', RMSE)

# 如果我们用其他方法得到了不同的系数，需要选择模型时，就用MSE小的时候对应的参数。
# 比如这次我们用AT， V，AP这3个列作为样本特征。不要RH， 输出仍然是PE
# 观察到MSE值增大，模型拟合效果变差
lin2 = LinearRegression()
lin2.fit(X_train[['AT','V','AP']], y_train)
print(lin2.intercept_)
print(lin2.coef_)

y_pred = lin2.predict(X_test[['AT','V','AP']])
print('MSE:',metrics.mean_squared_error(y_pred, y_test))
print('RSE:',np.sqrt(metrics.mean_squared_error(y_pred, y_test)))

# 通过交叉验证来持续优化模型，代码如下，我们采用10折交叉验证
linreg = LinearRegression()
predicted = cross_val_predict(linreg, X, y, cv=10)
print('MSE:', metrics.mean_squared_error(predicted, y))
print('RMSE:',np.sqrt(metrics.mean_squared_error(predicted, y)))
# 评估交叉验证准确率
print(cross_val_score(linreg, X, y, cv=10))
# 绘图观察拟合情况，如果散点离直线越近，拟合效果越好，即y_predict越接近y
fig = plt.figure()
axs = fig.add_subplot(111)
axs.scatter(y, predicted)
axs.plot([y.min(), y.max()],[y.min(), y.max()], 'k--',lw=4)
plt.show()

局部加权线性回归

线性回归的一个问题是有可能出现欠拟合现象，因为它求的是具有最小均方误差的无偏估计。显而易见，如果模型欠拟合将不能取得最好的预测效果。所以有些方法允许在估计中引入一些偏差，从而降低预测的均方误差。 其中的一个方法是局部加权线性回归（Locally Weighted Linear Regression，LWLR）。在该算法中，我们给待预测点附近的每个点赋予一定的权重，算法解出回归系数w的形式如下：
$$\hat{\theta }=(X^T\omega X{)}^{-1}{X}^T\omega Y$$
LWLR使用“核”（与支持向量机中的核类似）来对附近的点赋予更高的权重①。核的类型可以自由选择，最常用的核就是高斯核，高斯核对应的权重如下：
$$\omega (i,i)=exp(\frac{|x_i-x|}{-2k^2})$$
构建了一个只含对角元素的权重矩阵w，并且点x与x(i)越近，w(i,i)将会越大，式包含一个需要用户指定的参数k，它决定了对附近的点赋予多大的权重。k越大，用于训练模型的数据点的数量越多；k越小，用于训练模型的数据点的数量越少，容易过拟合。
数据集：ex0.txt

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties

font = FontProperties(fname=r"C:\Windows\Fonts\STFANGSO.TTF")
"""
函数说明：加载数据
Input:
    file_path:文件路径
Outout:
    dataMat
    labelMat
"""
def load_dataset(file_path):
    numFeat = len(open(file_path).readline().split('\t'))-1
    dataMat = []
    labelMat = []
    fr = open(file_path)
    for line in fr.readlines():
        line_arr = []
        current_line = line.strip().split('\t')
        for i in range(numFeat):
            line_arr.append(float(current_line[i]))
        dataMat.append(line_arr)
        labelMat.append(float(current_line[-1]))
    return dataMat,labelMat

"""
函数说明：LWLR计算权重
Input：
    test_point: 测试点
    k
    dataMat
    labelMat
Returns:
    ws
"""
def LELR(test_point,k,dataMat,yMat):
    m = np.shape(dataMat)[0]
    w = np.mat(np.eye((m)))
    for j in range(m):
        diffMat = test_point - dataMat[j,:]
        w[j,j] = np.exp(diffMat*diffMat.T/(-2.0*k**2))
    xTx = dataMat.T*w*dataMat
    if np.linalg.det(xTx) == 0.0:
        print("矩阵不可逆")
        return
    ws = xTx.I*dataMat.T*w*yMat
    return ws
"""
函数说明：计算出yHat:
Input:
    dataMat
    labelMat
Output:
    yHat
"""
def caluc_yhat(test_arr,dataMat,yMat,k):
    n = np.shape(test_arr)[0]
    yHat = np.zeros(n)
    for i in range(np.shape(test_arr)[0]):
        test_point = test_arr[i]
        ws = LELR(test_point,k,dataMat,yMat)
        yHat[i] = test_point*ws
    return yHat
"""
函数说明：不同k值局部加权回归的效果图
Input:
    k_list
    filepath"""
def k_plot(k_list,filepath):
    dataMat,labelMat = load_dataset(filepath)
    dataMat = np.mat(dataMat)
    yMat = np.mat(labelMat).T
    strInd = dataMat[:,1].argsort(0)
    xSort = dataMat[strInd][:,0]
    print(np.shape(xSort[:,0]))
    fig,axs = plt.subplots(nrows=3,ncols=1)
    print(np.shape(dataMat[:,1]))
    print(np.shape(yMat.flatten().A[0]))
    for i in range(len(k_list)):
        yHat = caluc_yhat(test_arr=dataMat,dataMat=dataMat,yMat=yMat,k = k_list[i])
        axs[i].scatter(np.array(dataMat[:,1]),yMat.flatten().A[0],s=3)
        axs[i].plot(xSort[:,1],yHat[strInd],linewidth=1.5,color='orange')
        axs[i].set_title("k="+str(k_list[i]))
    plt.tight_layout()
    plt.show()

filepath = "D:\深度学习\机器学习\8回归\ex0.txt"
k_list = [1.0,0.01,0.003]
k_plot(k_list,filepath)

随着样本点与待预测点距离的递增，权重将以指数级衰减 。输入参数k控制衰减的速度。当k = 1.0时权重很大，如同将所有的数据视为等权重，得出的最佳拟合直线与标准的回归一致。使用k = 0.01得到了非常好的效果，抓住了数据的潜在模式。下图使用k = 0.003纳入了太多的噪声点，拟合的直线与数据点过于贴近，出现过拟合现象。
局部加权线性回归也存在一个问题，即增加了计算量，因为它对每个点做预测时都必须使用整个数据集。

预测鲍鱼年龄

数据集：abalone.txt

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties

font = FontProperties(fname=r"C:\Windows\Fonts\STFANGSO.TTF")
"""
函数说明：加载数据
Input:
    file_path:文件路径
Outout:
    dataMat
    labelMat
"""
def load_dataset(file_path):
    numFeat = len(open(file_path).readline().split('\t'))-1
    dataMat = []
    labelMat = []
    fr = open(file_path)
    for line in fr.readlines():
        line_arr = []
        current_line = line.strip().split('\t')
        for i in range(numFeat):
            line_arr.append(float(current_line[i]))
        dataMat.append(line_arr)
        labelMat.append(float(current_line[-1]))
    return dataMat,labelMat

"""
函数说明：LWLR计算权重
Input：
    test_point: 测试点
    k
    dataMat
    labelMat
Returns:
    ws
"""
def LELR(test_point,k,dataMat,yMat):
    m = np.shape(dataMat)[0]
    w = np.mat(np.eye((m)))
    for j in range(m):
        diffMat = test_point - dataMat[j,:]
        w[j,j] = np.exp(diffMat*diffMat.T/(-2.0*k**2))
    xTx = dataMat.T*w*dataMat
    if np.linalg.det(xTx) == 0.0:
        print("矩阵不可逆")
        return
    ws = xTx.I*dataMat.T*w*yMat
    return ws
"""
函数说明：计算出yHat:
Input:
    dataMat
    labelMat
Output:
    yHat
"""
def caluc_yhat(test_arr,dataMat,yMat,k):
    n = np.shape(test_arr)[0]
    yHat = np.zeros(n)
    for i in range(np.shape(test_arr)[0]):
        test_point = test_arr[i]
        ws = LELR(test_point,k,dataMat,yMat)
        yHat[i] = test_point*ws
    return yHat

file_path = "D:\深度学习\机器学习\8回归\\abalone.txt"
dataMat,labelMat = load_dataset(file_path)
dataMat = np.mat(dataMat)
labelMat = np.mat(labelMat).T
"""
函数说明：计算损失值
Input:
    yTrue - 真实输出
    yHat - 预测输出
Output:
    (bias**2).sum() - 损失值
"""
def loss(yTrue,yHat):
    bias = yTrue.flatten().A[0] - yHat
    return (bias**2).sum()
y01hat = caluc_yhat(test_arr=dataMat[0:99],dataMat=dataMat[0:99],yMat=labelMat[0:99],k=0.1)
y1hat = caluc_yhat(test_arr=dataMat[0:99],dataMat=dataMat[0:99],yMat=labelMat[0:99],k=1)
y10hat = caluc_yhat(test_arr=dataMat[0:99],dataMat=dataMat[0:99],yMat=labelMat[0:99],k=10)
loss01  = loss(yTrue=labelMat[0:99],yHat=y01hat)
loss1 = loss(yTrue=labelMat[0:99],yHat=y1hat)
loss10 = loss(yTrue=labelMat[0:99],yHat=y10hat)
print(loss01)
print(loss1)
print(loss10)
ay01hat = caluc_yhat(test_arr=dataMat[100:230],dataMat=dataMat[0:99],yMat=labelMat[0:99],k=0.1)
ay1hat = caluc_yhat(test_arr=dataMat[100:230],dataMat=dataMat[0:99],yMat=labelMat[0:99],k=1)
ay10hat = caluc_yhat(test_arr=dataMat[100:230],dataMat=dataMat[0:99],yMat=labelMat[0:99],k=10)
aloss01  = loss(yTrue=labelMat[100:230],yHat=ay01hat)
aloss1 = loss(yTrue=labelMat[100:230],yHat=ay1hat)
aloss10 = loss(yTrue=labelMat[100:230],yHat=ay10hat)
print(aloss01)
print(aloss1)
print(aloss10)

使用较小的核将在训练集上得到较低的误差。那么，为什么不在所有数据集上都使用最小的核呢？这是因为使用最小的核将造成过拟合，对新数据不一定能达到最好的预测效果。从上述结果可以看到，在上面的三个参数中，核大小等于10时的测试误差最小，但它在训练集上的误差却是最大的。

缩减系数来"理解"数据

当特征比样本点还多（n > m），也就是说输入数据的矩阵X不是满秩矩阵。非满秩矩阵在求逆时会出现问题。本部分主要介绍岭回归于lasso回归。
目前，也引入岭回归与lasso回归来防止过拟合。

岭回归

• 岭回归是在损失函数中引入L2正则化
  $$J(\theta )=\frac{1}{2}(X\theta -Y{)}^T(X\theta -Y)+\frac{1}{2}\alpha \Vert \theta {\Vert }_2^2$$

• 求得回归系数：
  $$\hat{\theta }=(X^{T}X+\lambda I{)}^{-1}{X}^T\omega Y$$
  其中I为单位矩阵，维度等于特征数n。

• 实现
  底层实现
  数据集：abalone.txt

import numpy as np
import matplotlib.pyplot as plt

"""
函数说明：加载数据集
Input：
    file_path - 文件路径
Output:
    dataMat,labelMat
"""
def load_dataset(file_path):
    numFeat = len(open(file_path).readline().split('\t'))-1  # 特征数
    dataMat = []; labelMat = []
    fr = open(file_path)
    for line in fr.readlines():
        lin_arr = []
        current_line = line.strip().split('\t')
        for i in range(numFeat):
            lin_arr.append(float(current_line[i]))
        dataMat.append(lin_arr)
        labelMat.append(float(current_line[-1]))
    return dataMat,labelMat


"""
函数说明：岭回归求解回归系数
Input: 
    xMat - 特征数据集
    yMat - 输出数据集
    lam - lamda,岭回归乘法因子
Output:
    ws - 回归系数
"""
def linregress(xMat,yMat,lam):
    xTx = xMat.T * xMat
    Imat = np.mat(np.eye(np.shape(xMat)[1]))
    demo = xTx + lam*Imat
    if np.linalg.det(demo) == 0.0:
        print("此行列式为0")
        return()
    ws = demo.I * xMat.T * yMat
    return ws


"""
函数说明：获取不同lamda值对应的ws
Input: 
    file_path - 文件路径
    lam_list - lamda值列表
Output:
    wMat - 回归系数矩阵
"""
def ridgetest(file_path):
    lam_list, xMat, yMat = creat_lam_list(file_path)
    wMat = np.mat(np.zeros((len(lam_list),np.shape(xMat)[1])))  # 回归系数矩阵初始化为0
    for i in range(len(lam_list)):
        ws = linregress(xMat,yMat,lam_list[i])
        wMat[i,:] = ws.T
    return wMat


"""
函数说明：给定lam_list
Input:
Output:
    lamb_list - 惩罚因子列表
"""
def creat_lam_list(file_path):
    dataMat, labelMat = load_dataset(file_path)
    yMat = np.mat(labelMat).T
    xMat = np.mat(dataMat)
    Xmean = np.mean(xMat, 0)
    Xvar = np.var(xMat, 0)
    xMat = (xMat - Xmean) / Xvar  # 特征标准化
    Ymean = np.mean(yMat, 0)
    yMat = yMat - Ymean
    lam_list = []
    numTests = 30
    for i in range(numTests):
        lam_list.append(np.exp(i-10))
    return lam_list, xMat, yMat


file = "D:\深度学习\机器学习\8回归\\abalone.txt"
weights = ridgetest(file)
lam_list, xMat, yMat = creat_lam_list(file)
plt.figure()
plt.plot(np.log(lam_list), weights)
plt.show()

在λ最小时，可以得到所有系数的原始值（与线性回归一致）；而λ非常大时，系数全部缩减成0；在中间部分的某值将可以取得最好的预测效果。

scikit-learn实现
数据集：Folds5x2_pp.csv

import numpy as np
from sklearn.linear_model import Ridge
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
import pandas as pd
from sklearn.linear_model import RidgeCV
from sklearn import metrics

data = pd.read_csv("D:\深度学习\机器学习\8回归\Folds5x2_pp.csv")
# 观察数据是否有缺失
print(data.info())
print(data.head())
print(data.tail())

X = data[['AT','V','AP','RH']]
y = data['PE']
print(X.head())

# 训练集、测试集划分
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state = 1)
print(np.shape(X_train))
print(np.shape(X_test))

# 定义模型
clf = Ridge(alpha=1)
clf.fit(X_train, y_train)
y_pred = clf.predict(X_test)
print(metrics.mean_squared_error(y_pred, y_test))
print(clf.intercept_)
print(clf.coef_)

"""选择最佳alpha"""
clf = RidgeCV(alphas = [0.5, 1, 1.5, 2, 3, 4, 5, 6, 7, 8])
clf.fit(X_train, y_train)
y_predict= clf.predict(X_test)
print(metrics.mean_squared_error(y_predict, y_test))
print(clf.intercept_)
print(clf.coef_)
# 输出最佳的alpha
print(clf.alpha_)

"""研究alpha和超参数的过程：
alpha越大，那么正则项惩罚的就越厉害，得到的回归系数θ就越小，最终趋近于0
而alpha越小，即正则化项越小，那么回归系数θ就越来越接近于普通的线性回归系数"""
# X is a 10x10 matrix
X = 1. / (np.arange(1, 11) + np.arange(0, 10)[:, np.newaxis])
# y is a 10 x 1 vector
y = np.ones(10)

n_alphas = 200
# alpha一共200个，都在10的-10次方和10的-2次方之间
alphas = np.logspace(-10, -2, n_alphas)

clf = Ridge(fit_intercept=False)
coefs = []
for i in alphas:
    clf.set_params(alpha=i)
    clf.fit(X,y)
    coefs.append(clf.coef_)

ax = plt.gca()
ax.plot(alphas, coefs)
ax.set_xscale('log')
# 翻转x轴的大小方向，让alpha从大到小显示
ax.set_xlim(ax.get_xlim()[::-1])
plt.xlabel('alpha')
plt.ylabel('weights')
plt.title('Ridge coefficients as a function of the regularization')
plt.axis('tight')
plt.show()

岭回归不抛弃任何一个特征,缩小了回归系数，使得模型相对而言比较的稳定，但和Lasso回归比，这会使得模型的特征留的特别多，模型解释性差。

lasso

• Lasso回归是在损失函数中引入L1正则化
  $$J(\theta )=\frac{1}{2}(X\theta -Y{)}^T(X\theta -Y)+\lambda \Vert \theta {\Vert }_1$$
  Lasso回归可以使得一些特征的系数变小，甚至还是一些绝对值较小的系数直接变为0，增强模型的泛化能力。Lasso求解比较困难，可以用前向逐步回归进行求解。

前向逐步回归

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
font = FontProperties(fname=r"C:\Windows\Fonts\STFANGSO.TTF")


"""
函数说明：加载数据
Input:
    file_path:文件路径
Outout:
    dataMat
    labelMat
"""
def load_dataset(file_path):
    numFeat = len(open(file_path).readline().split('\t'))-1
    dataMat = []
    labelMat = []
    fr = open(file_path)
    for line in fr.readlines():
        line_arr = []
        current_line = line.strip().split('\t')
        for i in range(numFeat):
            line_arr.append(float(current_line[i]))
        dataMat.append(line_arr)
        labelMat.append(float(current_line[-1]))
    return dataMat,labelMat


"""
函数说明：计算误差
Input:
    yMat - 真实输出
    yHat - 预测输出
Output:
    loss - 损失值
"""
def rssError(yMat, yHat):
    bias = yMat.flatten().A[0] - yHat
    loss = (bias**2).sum()
    return loss
"""
函数说明：前向逐步线性回归
Input:  
    dataMat - 特征数据
    labelMat - 输出数据
    eps - 改变的很小值
    numIt - 迭代次数
Output:
    RetrunMat - 权重矩阵
"""
def stageWise(dataMat, labelMat, eps, numIt):
    xMat = np.mat(dataMat); yMat = np.mat(labelMat).T
    Xmean = np.mean(xMat, 0)
    Xvar = np.var(xMat, 0)
    xMat = (xMat - Xmean) / Xvar  # 特征标准化
    Ymean = np.mean(yMat, 0)
    yMat = yMat - Ymean
    n = np.shape(xMat)[1]
    # 初始化
    ReturnMat = np.zeros((numIt, n))
    ws = np.zeros((n,1))
    wsMax = ws.copy()
    wsTest = ws.copy()
    for i in range(numIt):
        print(ws.T)
        lossMax = np.inf # 每次迭代设置误差为无穷大
        for j in range(n):
            # 增加或减小wsTest
            for sign in [-1,1]:
                wsTest = ws.copy()
                wsTest[j] += eps*sign
                yHat = xMat*wsTest
                loss = rssError(yMat, yHat)
                if loss < lossMax:
                    lossMax = loss
                    wsMax = wsTest.copy()
        ws = wsMax.copy()
        ReturnMat[i,:] = ws.T
    return ReturnMat


file = "D:\深度学习\机器学习\8回归\\abalone.txt"
dataMat, labelMat = load_dataset(file_path=file)
ReturnMat = stageWise(dataMat=dataMat, labelMat=labelMat, eps=0.01, numIt=200)

最后分享自己做的线性回归的思维导图
线性回归