学习视频地址：第一章：线性回归原理推导 1-回归问题概述_哔哩哔哩_bilibili
一、线性回归原理推导
1、线性回归问题的数学表达可以写作如下方程式（1）
其中y(i) 和x(i)均是数据集已知的量，ε是需要推导的参数矩阵，ε是随机误差矩阵，独立且具有相同的愤怒，服从均值为0方差为θ2的高斯分布。
由于ε服从高斯分布，所以可以得到一个方程式（2）
将方程式（1）代入到方程式（2）中，可以消去误差ε，得到方程式（3）
我们可以找到一个关于θ的似然函数（4）
由于我们仅仅关注极值点的位置，并不关注极值的大小，所以可以给两边取对数，就可以将累乘变成累加，得到方程式（5）
提出常数项可以得到方程式（6）
方程式（6）可以看作是一个常数减去一个平方数的形式，所以要找方程式（6）的极大值，就可以看作找后面平方数的极小值，我们得到方程式（7）（最小二乘法的公式）
对J(θ）求偏导，并且令偏导为0，我们就可以求出了θ的值为(可以当作一个巧合，机器学习是需要迭代的，直接求出不符合机器学习的思想)
2、机器学习的套路：给机器一堆数据，然后告诉机器什么样的学习结果是对的（目标函数），让他朝着这个方向进行迭代。
3、梯度下降算法
首先由上面的推导，我们可以得出线性回归的目标函数为方程式（8）
其中m为样本的数量，我们对J求关于θ的偏导得到梯度方程式（9）
梯度下降算法原理就是每次迭代过程对参数向梯度反方向前进梯度个步数生成新的参数、直到找到最拟合目标函数的参数为止。
3.1批量梯度下降：每次前进总样本的平均梯度，容易得到最优解，但是速度很慢，在数据量大的时候一般不使用，参数迭代方程式（10）如下：
3.2随机梯度下降：每次找一个样本的梯度进行迭代，速度快，但是随机性强，每次迭代不一定向着收敛的方向，参数迭代方程式（11）如下：
3.3小批量梯度下降：每次使用一小部分的平均梯度进行迭代，速度快，迭代方向也比较好，经常被使用，参数迭代方程式（12）如下：
其中α为学习率
二、线性回归代码实现
数据集我们使用阿里天池的幸福度预测数据集
数据集地址：快来一起挖掘幸福感！_学习赛_天池大赛-阿里云天池 (aliyun.com)
下载好数据集后我们先不处理它，先来搭建模型
#!/usr/bin/env pytorch # -*- coding: UTF-8 -*- """ @Project :Follow-Tang-Yudi-machine-learning-algorithm-intensive-code-practice @File :linear_model.py @IDE :PyCharm @Author :张世航 @Date :2023/3/7 9:08 @Description :一个线性回归的模型 """ import numpy as np import prepare_train class LinearRegression: def __init__(self, data, labels, polynomial_degree=0, sinusoid_degree=0, normalize_data=True): """ 1、对数据进行预处理操作 2、得到特征个数 3、初始化参数矩阵 :param data:原始数据 :param labels:数据的标签 :param polynomial_degree:多项式维度 :param sinusoid_degree:正弦维度 :param normalize_data:是否进行归一化 """ (data_processed, features_mean, features_deviation) = prepare_train.prepare_for_train(data, polynomial_degree, sinusoid_degree, normalize_data) self.data = data_processed self.labels = labels self.features_mean = features_mean self.features_deviation = features_deviation self.polynomial_degree = polynomial_degree self.sinusoid_degree = sinusoid_degree self.normalize_data = normalize_data num_features = self.data.shape[1] self.theta = np.zeros((num_features, 1)) def train(self, alpha, num_epoch=500): """ 开始训练 :param alpha: 学习速率 :param num_epoch: 迭代次数 :return: 训练好的参数,损失值记录 """ cost_history = self.gradient_descent(alpha, num_epoch) return self.theta, cost_history def gradient_descent(self, alpha, num_epoch): """ 小批量梯度下降算法 :param alpha: 学习速率 :param num_epoch: 迭代次数 :return: 损失值的记录 """ cost_history = [] for i in range(num_epoch): self.gradient_step(alpha) cost_history.append(self.cost_function(self.data, self.labels)) return cost_history def gradient_step(self, alpha): """ 梯度下降参数更新方法 :param alpha: 学习率 :return: 无 """ num_examples = self.data.shape[0] prediction = LinearRegression.hypothesis(self.data, self.theta) delta = prediction - self.labels theta = self.theta theta = theta - alpha * (1 / num_examples) * (np.dot(delta.T, self.data)).T self.theta = theta def cost_function(self, data, labels): """ 计算损失值函数（最小二乘法） :param data:当前数据 :param labels:当前标签 :return:预测损失值 """ num_example = data.shape[0] delta = LinearRegression.hypothesis(data, self.theta) - labels cost = (1 / 2) * np.dot(delta.T, delta) return cost[0][0] @staticmethod def hypothesis(data, theta): """ 求预测值 :param data: 输入数据 :param theta: 当前参数 :return: 预测值 """ prediction = np.dot(data, theta) return prediction def get_cost(self, data, labels): """ 获取损失值 :param data: 传入的数据 :param labels: 数据的标签 :return: 当前模型预测数据的损失值 """ (data_processed, _, _) = prepare_train.prepare_for_train(data, self.polynomial_degree, self.sinusoid_degree, self.normalize_data) return self.cost_function(data_processed, labels) def predict(self, data): """ 预测函数 :param data: 输入数据 :return: 预测的值 """ (data_processed, _, _) = prepare_train.prepare_for_train(data, self.polynomial_degree, self.sinusoid_degree, self.normalize_data) predictions = LinearRegression.hypothesis(data_processed, self.theta) return predictions
其中在模型初始化时候对数据进行了预处理
#!/usr/bin/env pytorch # -*- coding: UTF-8 -*- """ @Project :Follow-Tang-Yudi-machine-learning-algorithm-intensive-code-practice @File :prepare_train.py @IDE :PyCharm @Author :张世航 @Date :2023/3/7 9:24 @Description :数据预处理 """ import numpy as np import normalize import generate_polynomials import generate_sinusoids def prepare_for_train(data, polynomial_degree=0, sinusoid_degree=0, normalize_data=True): """ 对数据进行预处理 :param data: 原始数据 :param polynomial_degree: 多项式维度 :param sinusoid_degree: 正弦维度 :param normalize_data: 是否进行归一化 :return: 处理后的数据,特征均值,特征方差 """ # 获取样本总数 num_examples = data.shape[0] data_processed = np.copy(data) # 预处理 features_mean = 0 features_deviation = 0 data_normalized = data_processed if normalize_data: data_normalized, features_mean, features_deviation = normalize.normalize(data_processed) data_processed = data_normalized # 特征变量正弦变换 if sinusoid_degree > 0: sinusoids = generate_sinusoids.generate_sinusoids(data_normalized, sinusoid_degree) data_processed = np.concatenate((data_processed, sinusoids), axis=1) # 特征变量多项式变换 if polynomial_degree > 0: polynomials = generate_polynomials.generate_polynomials(data_processed, polynomial_degree, normalize_data) data_processed = np.concatenate((data_processed, polynomials), axis=1) # 加一列1 data_processed = np.hstack((np.ones((num_examples, 1)), data_processed)) return data_processed, features_mean, features_deviation
预处理时候使用了正弦变换和多项式变换增加数据的维度，增加收敛速度，但是如果维度过高可能会过拟合
#!/usr/bin/env pytorch # -*- coding: UTF-8 -*- """ @Project :Follow-Tang-Yudi-machine-learning-algorithm-intensive-code-practice @File :generate_polynomials.py @IDE :PyCharm @Author :张世航 @Date :2023/3/8 8:34 @Description :进行多项式变换增加变量复杂度 """ import numpy as np import normalize def generate_polynomials(dataset, polynomials_degree, normalize_data=False): """ 变换方法：x1, x2, x1^2, x2^2, x1 * x2, x1 * x2^2, etc. :param dataset:原始数据 :param polynomials_degree:多项式的维度 :param normalize_data:是否归一化 :return:生成的多项式参数 """ features_split = np.array_split(dataset, 2, axis=1) dataset_1 = features_split[0] dataset_2 = features_split[1] (num_examples_1, num_features_1) = dataset_1.shape (num_examples_2, num_features_2) = dataset_2.shape if num_examples_1 != num_examples_2: raise ValueError("can not generate polynomials for two sets with different number") if num_features_1 == 0 and num_features_2 == 0: raise ValueError("can not generate polynomials for two sets with no colums") if num_features_1 == 0: dataset_1 = dataset_2 elif num_features_2 == 0: dataset_2 = dataset_1 num_features = num_features_1 if num_features_1 < num_examples_2 else num_features_2 dataset_1 = dataset_1[:, :num_features] dataset_2 = dataset_2[:, :num_features] polynomials = np.empty((num_examples_1, 0)) for i in range(1, polynomials_degree + 1): for j in range(i + 1): polynomial_feature = (dataset_1 ** (i - j)) * (dataset_2 ** j) polynomials = np.concatenate((polynomials, polynomial_feature), axis=1) if normalize_data: polynomials = normalize.normalize(polynomials)[0] return polynomials #!/usr/bin/env pytorch # -*- coding: UTF-8 -*- """ @Project :Follow-Tang-Yudi-machine-learning-algorithm-intensive-code-practice @File :generate_sinusoids.py @IDE :PyCharm @Author :张世航 @Date :2023/3/8 8:35 @Description :对参数进行正弦变换 """ import numpy as np def generate_sinusoids(dataset, sinusoid_degree): """ 变换方式 sin（x) :param dataset: 原始数据 :param sinusoid_degree:变换维度 :return: 变换后的参数 """ num_examples = dataset.shape[0] sinusoids = np.empty((num_examples, 0)) for degree in range(1, sinusoid_degree + 1): sinusoid_fatures = np.sin(degree * dataset) sinusoids = np.concatenate((sinusoids, sinusoid_fatures), axis=1) return sinusoids
预处理时候对数据进行归一化操作，让数据更适合训练
#!/usr/bin/env pytorch # -*- coding: UTF-8 -*- """ @Project :Follow-Tang-Yudi-machine-learning-algorithm-intensive-code-practice @File :normalize.py @IDE :PyCharm @Author :张世航 @Date :2023/3/7 16:02 @Description :归一化数据 """ import numpy as np def normalize(features): """ 特征归一化 :param features: 传入特征 :return: 归一化后的特征,特征均值,特征标准差 """ features_normalized = np.copy(features).astype(float) # 计算均值 features_mean = np.mean(features, 0) # 计算标准差 features_deviation = np.std(features, 0) # 标准化操作 if features.shape[0] > 1: features_normalized -= features_mean # 防止除以0 features_deviation[features_deviation == 0] = 1 features_normalized /= features_deviation return features_normalized, features_mean, features_deviation
搭建好网络模型后我们就可以进行模型的训练了
首先加载数据集数据
plotly.offline.init_notebook_mode() data_train = pd.read_csv("happiness_train_complete_prepared.csv", encoding="ISO-8859-1") data_test = pd.read_csv("happiness_test_complete_prepared.csv", encoding="ISO-8859-1") data_submit = pd.read_csv("happiness_submit.csv", encoding="ISO-8859-1") print(data_train.shape, data_test.shape, data_submit.shape) data_train1 = data_train.sample(frac=0.8) data_train2 = data_train.drop(data_train1.index) input_name = ["edu", "income"] output_name = "happiness" x_train1 = data_train1[input_name].values y_train1 = data_train1[[output_name]].values x_train2 = data_train2[input_name].values y_train2 = data_train2[[output_name]].values x_test = data_test[input_name].values
由于是个小练习，所以我们仅仅考虑了学历和收入两个维度来预测幸福度，由于测试集没有幸福度答案，所以我们将训练集分为两部分，80%作为训练集进行训练，20%作为测试集来判断模型训练的好坏
设置模型的超参数如下
num_epoch = 500
learning_rate = 0.02
polynomial_degree = 5
sinusoid_degree = 15
normalize_data = True
进行模型的训练并且绘制损失值
linear_model = LinearRegression(x_train1, y_train1, polynomial_degree, sinusoid_degree, normalize_data) theta, cost_history = linear_model.train(learning_rate, num_epoch) print("{}->{}".format(cost_history[0], cost_history[-1])) plt.plot(range(num_epoch), cost_history) plt.xlabel("epoch") plt.ylabel("cost") plt.title = "cost curve" plt.show()
使用模型对测试集进行预测并且获取得分
z_test = linear_model.predict(x_train2) sorce = 0.0 num_test = z_test.shape[0] for index, value in enumerate(z_test): sorce += (round(value[0]) - y_train2[index]) ** 2 rate = sorce / num_test print(rate)
预测下载好的测试集的幸福度并且保存到csv中
test_pridict = linear_model.predict(x_test)
result = list(test_pridict)
result = list(map(lambda x: round(x[0]), result))
data_submit["happiness"] = result
data_submit.to_csv("happiness.csv", index=False)
好了，我们可以运行一下试试了
得到如上图的损失曲线
打印出来的得分为0.65
那么我们将预测出的幸福度提交到阿里天池看看我们的小练习的得分吧
不出所料，仅仅考虑两个维度得分结果是0.7，比我们打印出来的得分还要差，如果感兴趣的同学，可以多考虑几个维度试试，不过需要注意下载的数据集有的列数据有空值，需要进行预处理，因为线性回归模型没办法处理空值的。
三、模型评估方法
1、sklearn工具包提供了机器学习的一些模型、损失函数以及数据集切分的方法
2、一般将数据集80%切分出来作为训练集，20%切分出来作为测试集
3、交叉验证：将训练集又分为几个部分，在训练中轮流取一部分作为内部测试集，其余作为训练集，进行验证，可以使用cross_val_score得到验证结果，也可以使用stratified kfold手动进行分割，自己计算结果
4、混淆矩阵（二分类问题）：测试结果可以列出下面矩阵
| 正类 | 负类 | |
| 被检测出 | TP | FP |
| 未被检测出 | FN | TN |
可以使用cross_val_predict获取交叉验证的估计值，也可以使用coufusion_matrixlai来获取混淆矩阵
精度的计算方式为TP/（TP+FP）表明模型的精准度，回归的计算方式为TP/（TP + FN）表明模型的泛化能力，可以使用内置函数precision_sorce来计算模型的精度，recall_score来获取模型的泛化能力
F1指标是综合考虑精度和回归能力的一个指标计算方法如下：
可以使用内置函数F1_score直接获取F1的得分
5、阈值
decision_function可以得到当前样本每一个分类的得分值。
一个模型设置的阈值越低，精度越低但是泛化能力越强，阈值越高，精度越高但是泛化能力会变弱
6、ROC曲线
precision_recall_curve可以得到不同阈值下的精度和回归值
可以计算正确预测值TPR = TP / （TP + FN），错误预测值FPR = FP / （FP+TN）
以TPR为x，FPR为y做的曲线叫做ROC曲线，ROC曲线所围成的面积越大说明模型就越好
roc_curve函数可以得到不同阈值下TPR和FPR的结果
四、线性回归实验
1、sklearn中，也内置了线性回归的模型LinearRegression
2、标准化的作用：取值范围不同的值会导致模型收敛速度变慢，进行标准化的操作有利于模型的收敛
3、学习率的作用：小的学习率收敛速度慢，大的学习率会产生抖动，所以一般先设置大的学习率，然后随着训练次数增加慢慢降低学习率
4、多项式回归：内置了PolonomialFeatures方法来将输入进行多项式展开增加参数的维度，但是维度太高会导致过拟合
5、metrics模块提供了方法进行评估模型，数据样本越大过拟合的概率也就越小
6、正则化：岭回归和lasso回归在损失函数中加入了正则项，可以消除过拟合
岭回归的损失函数为：
其中α为正则化惩罚力度
lasso回归损失函数为：
五、逻辑回归
1、逻辑回归是一个分类算法（二分类）
2、逻辑回归通过将线性回归结果映射成为概率值来进行分类问题的解决，这里提到一个sigmoid函数，函数公式如下：
自变量z是任意的实数，值域为[0,1]，将线性回归的结果通过sigmoid函数映射到[0，1]空间，就得到预测输入为这个类型的概率值，从而实现了线性回归到分类的转换
3、将线性回归带入sigmoid函数可以得到预测函数
那么第i个样本的概率值可以表示为如下方程式
可以得到似然函数
同样我们不关心似然函数的极值是多少，但是关注似然函数的极值点在哪里，对两边取对数，得到对数似然函数
引入J（θ）将对数似然函数转换为梯度下降任务
对J（θ）求偏导得到梯度
所以可以的写出参数θ的迭代方程式为
六、多分类逻辑回归实战
首先我们实现工具集sigmoid函数，很简单，代码如下
#!/usr/bin/env pytorch # -*- coding: UTF-8 -*- """ @Project :Follow-Tang-Yudi-machine-learning-algorithm-intensive-code-practice @File :sigmoid.py @IDE :PyCharm @Author :张世航 @Date :2023/3/15 15:57 @Description :sigmoid函数的定义 """ import math import matplotlib.pyplot as plt import numpy as np def sigmoid_function(z): """ sigmoid函数，把参数归一到0-1中 :param z:输入是x :return: 输出的f(x) """ fz = [] for num in z: fz.append(1 / (1 + math.exp(-num))) return fz def sigmoid(matrix): """ 矩阵的sigmoid方法 :param matrix: 输入的矩阵 :return: sigmoid后的矩阵 """ return 1 / (1 + np.exp(-matrix)) if __name__ == '__main__': z = np.arange(-10, 10, 0.01) fz = sigmoid_function(z) plt.title("sigmoid function") plt.xlabel("z") plt.ylabel("f(z)") plt.plot(z, fz) plt.show()
我们主要使用矩阵的sigmoid方法，写普通方法只是为了展示sigmoid函数的样子，运行一些就可以看到sigmoid函数的样子了
接下来我们搭建自己的逻辑回归模型
#!/usr/bin/env pytorch # -*- coding: UTF-8 -*- """ @Project :Follow-Tang-Yudi-machine-learning-algorithm-intensive-code-practice @File :logistic_regression_model.py @IDE :PyCharm @Author :张世航 @Date :2023/3/15 15:46 @Description :逻辑回归的模型 """ import numpy as np from scipy.optimize import minimize from util import prepare_train from util import sigmoid class LogisticRegression: def __init__(self, data, labels, polynomial_degree=0, sinusoid_degree=0, normalize_data=True): """ 1、对数据进行预处理操作 2、得到特征个数 3、初始化参数矩阵 :param data:原始数据 :param labels:数据的标签 :param polynomial_degree:多项式维度 :param sinusoid_degree:正弦维度 :param normalize_data:是否进行归一化 """ (data_processed, features_mean, features_deviation) = prepare_train.prepare_for_train(data, polynomial_degree, sinusoid_degree, normalize_data) self.data = data_processed self.labels = labels self.unique_labels = np.unique(labels) self.features_mean = features_mean self.features_deviation = features_deviation self.polynomial_degree = polynomial_degree self.sinusoid_degree = sinusoid_degree self.normalize_data = normalize_data num_features = self.data.shape[1] num_unique_labels = np.unique(labels).shape[0] self.theta = np.zeros((num_unique_labels, num_features)) def train(self, num_epoch=1000): """ 开始训练 :param num_epoch: 迭代次数 :return: 训练好的参数,损失值记录 """ cost_histories = [] num_features = self.data.shape[1] for label_index, label in enumerate(self.unique_labels): current_initial_theta = np.copy(self.theta[label_index].reshape(num_features, 1)) current_initial_label = (self.labels == label).astype(float) (current_theta, cost_history) = LogisticRegression.gradient_descent(self.data, current_initial_label, current_initial_theta, num_epoch) self.theta[label_index] = current_theta.T cost_histories.append(cost_history) return self.theta, cost_histories @staticmethod def gradient_descent(data, labels, current_initial_theta, num_epoch): """ 梯度下降算法 :param data:当前输入数据 :param labels: 当前的初始标签 :param current_initial_theta: 当前的初始参数 :param num_epoch: 迭代次数 :return: 损失值的记录, 训练后的参数 """ cost_history = [] num_features = data.shape[1] result = minimize( # 要优化的目标函数 lambda current_theta: LogisticRegression.cost_function(data, labels, current_theta.reshape(num_features, 1)), # 初始化的参数 current_initial_theta, # 优化策略 method="CG", # 梯度迭代的计算公式 jac=lambda current_theta: LogisticRegression.gradient_step(data, labels, current_theta.reshape(num_features, 1)), # 迭代次数 options={"maxiter" : num_epoch}, # 记录结果 callback=lambda current_theta: cost_history.append( LogisticRegression.cost_function(data, labels, current_theta.reshape(num_features, 1))) ) if not result.success: raise ArithmeticError("can not minimize cost" + result.message) optimized_theta = result.x.reshape(num_features, 1) return optimized_theta, cost_history @staticmethod def gradient_step(data, labels, theta): """ 梯度下降参数更新方法 :param data:传入的数据 :param labels: 传入的标签 :param theta: 初始化的参数 :return: 无 """ num_examples = data.shape[0] prediction = LogisticRegression.hypothesis(data, theta) label_diff = prediction - labels gradients = (1 / num_examples) * np.dot(data.T, label_diff) return gradients.T.flatten() @staticmethod def cost_function(data, labels, theta): """ 计算损失值函数 :param data:当前数据 :param labels:当前标签 :param theta:当前参数 :return:预测损失值 """ num_example = data.shape[0] predictions = LogisticRegression.hypothesis(data, theta) y_is_set_cost = np.dot(labels[labels == 1].T, np.log(predictions[labels == 1])) y_is_not_set_cost = np.dot(1 - labels[labels == 0].T, np.log(1 - predictions[labels == 0])) cost = (-1 / num_example) * (y_is_not_set_cost + y_is_set_cost) return cost @staticmethod def hypothesis(data, theta): """ 求预测值 :param data: 输入数据 :param theta: 当前参数 :return: 预测值 """ predictions = sigmoid.sigmoid(np.dot(data, theta)) return predictions def get_cost(self, data, labels): """ 获取损失值 :param data: 传入的数据 :param labels: 数据的标签 :return: 当前模型预测数据的损失值 """ (data_processed, _, _) = prepare_train.prepare_for_train(data, self.polynomial_degree, self.sinusoid_degree, self.normalize_data) return self.cost_function(data_processed, labels) def predict(self, data): """ 预测函数 :param data: 输入数据 :return: 预测的值 """ num_examples = data.shape[0] (data_processed, _, _) = prepare_train.prepare_for_train(data, self.polynomial_degree, self.sinusoid_degree, self.normalize_data) predictions = LogisticRegression.hypothesis(data_processed, self.theta.T) max_index = np.argmax(predictions, axis=1) class_prediction = np.empty(max_index.shape, dtype=object) for index, label in enumerate(self.unique_labels): class_prediction[max_index == index] = label return class_prediction.reshape((num_examples, 1))
我们选择著名的鸢尾花数据集作为我们模型的数据输入，为了方便展示选择花瓣的长和宽两个维度作为输入数据，打印数据分布看看
data = pd.read_csv("Iris数据集/iris.csv") iris_types = ["setosa", "versicolor", "virginica"] x_axis = "Petal.Length" y_axis = "Petal.Width" for iris_type in iris_types: plt.scatter(data[x_axis][data["Species"] == iris_type], data[y_axis][data["Species"] == iris_type], label=iris_type) plt.show()
我们构建自己的数据集和标签
num_examples = data.shape[0]
x_train = data[[x_axis, y_axis]].values.reshape((num_examples, 2))
y_train = data["Species"].values.reshape((num_examples, 1))
设置超参数
max_epoch = 1000
polynomial_degree = 0
sinusoid_degree = 0
训练数据打印损失值
model = LogisticRegression(x_train, y_train, polynomial_degree, sinusoid_degree, False) theta, cost_histories = model.train(max_epoch) labels = model.unique_labels plt.plot(range(len(cost_histories[0])), cost_histories[0], label=labels[0]) plt.plot(range(len(cost_histories[1])), cost_histories[1], label=labels[1]) plt.plot(range(len(cost_histories[2])), cost_histories[2], label=labels[2]) plt.show()
可以看出虽然对橙色那种类型的花损失值较多，但是总的来说都是收敛的
产生一些测试数据进行准确率验证
y_train_predictions = model.predict(x_train)
precision = np.sum(y_train_predictions == y_train) / y_train.shape[0] * 100
print("precision:", precision)
看起来准确率还不错
接下来打印决策边界看看
samples = 150 X = np.linspace(x_min, x_max, samples) Y = np.linspace(y_min, y_max, samples) Z_setosa = np.zeros((samples, samples)) Z_versicolor = np.zeros((samples, samples)) Z_virginica = np.zeros((samples, samples)) for x_index, x in enumerate(X): for y_index, y in enumerate(Y): data = np.array([[x, y]]) prediction = model.predict(data)[0][0] if prediction == "setosa": Z_setosa[x_index][y_index] = 1 elif prediction == "versicolor": Z_versicolor[x_index][y_index] = 1 elif prediction == "virginica": Z_virginica[x_index][y_index] = 1 for iris_type in iris_types: plt.scatter(x_train[(y_train == iris_type).flatten(), 0], x_train[(y_train == iris_type).flatten(), 1], label=iris_type) plt.contour(X, Y, Z_setosa) plt.contour(X, Y, Z_versicolor) plt.contour(X, Y, Z_virginica) plt.show()
看得出决策边界还是比较符合我们的样本数据的，有兴趣的同学可以试试四个维度的逻辑回归实验
七、逻辑回归实验分析
1、可以使用sklearn中的meshgrid函数产生对元组x1和x2生成全组合
2、softmax计算概率公式
本质为归一化，不过是加入了exp拉大了样本的差异
损失函数（交叉熵）
八、聚类算法
1、聚类算法是无监督学习，本质是把相似的东西分为一个一个簇
2、k-means算法：
2.1 k值：算法将样本分成的簇的数量，需要在训练前指定
2.2 质心：同簇各样本在各维度取均值
2.3 距离的度量：可以使用欧式距离和余弦相似度（需要先标准化）
2.4 优化目标：
2.5 工作流程
2.5.1 初始化k个点，计算各个点到中国点的距离，判断类别
2.5.2 将不同类别的点分别计算质心，作为新的k个点重新计算各个点与k个点之间的距离，再判断类别，进行循环迭代
2.6 优点：简单，快速，适合常规数据集
2.7 劣势：k值难确定，复杂度和样本呈线性关系，很难发现任意形状的簇，初值的影响很大
3、DBSCAN算法：
3.1 核心对象：若某个点的密度达到算法设定的阈值，那么它就是核心点
3.2 邻域距离的阈值：设定的半径
3.3 直接密度可达：若某点p在q的r领域内，且q是核心点，那么就称p和q直接密度可达
3.4 密度可达：若有一个点的序列q0、q1……qk对任意qi和qi-1是直接密度可达，那我们就称q0和qk密度可达，实际上就是直接密度可达的传播
3.5 密度相连：若从某核心点p出发点q和点k都是密度可达的，那么称q和k密度相连
3.6 边界点：属于某一类的非核心点，不能发展下线
3.7 噪声点：不属于任何一个类簇的点，从任何一个核心的都不可达
3.8 工作流程：
3.8.1 指定数据集，半径和密度阈值
3.8.2 随机选择一个为访问的对象p，如果所有对象都访问过则工作结束
3.8.3 如果p的半径内有大于阈值数量的点，创建一个簇C，将p放入C中；否则标记为噪声点，返回2
3.8.4 对p范围内的p'进行遍历，如果p'半径范围内也有核心点，那么把p'半径范围内的核心点也纳入遍历范围
3.8.5 如果p'不是任何簇的成员，将p'加入C，返回4
3.9 参数的选择：半径r，分析给定的数据集，寻找突变点；阈值，小一点，多次尝试
3.10 优势：不需要指定簇的个数，可以发现任意形状的簇，擅长寻找离群点，参数少
3.11 劣势：高维数据困难（可做降维），参数难选择（参数对结果影响很大），sklearn中效率很慢（数据削减策略）
九、k-means算法实战
先根据k-means的工作流程创建聚类模型
#!/usr/bin/env pytorch # -*- coding: UTF-8 -*- """ @Project :Follow-Tang-Yudi-machine-learning-algorithm-intensive-code-practice @File :k_means_model.py @IDE :PyCharm @Author :张世航 @Date :2023/3/21 15:17 @Description :k-means模型 """ import numpy as np class K_Means: def __init__(self, data, num_clusters): self.data = data self.num_clusters = num_clusters def train(self, num_epoch): """ 迭代模型 :param num_epoch: 迭代次数 :return:质心位置,样本最近的质心索引 """ # 1、随机选择k个中心点 centroids = K_Means.centroids_init(self.data, self.num_clusters) # 2、开始训练 num_examples = self.data.shape[0] closest_centroids_ids = np.empty((num_examples, 1)) for _ in range(num_epoch): # 3、得到当前样本点到k个中心点的距离，找到最近的 closest_centroids_ids = K_Means.centroids_find_closedt(self.data, centroids) # 4、更新质心 centroids = K_Means.centroids_compute(self.data, closest_centroids_ids, self.num_clusters) return centroids, closest_centroids_ids @staticmethod def centroids_init(data, num_clusters): """ 初始化中心点 :param data: 输入的数据 :param num_clusters: 要分成的簇的个数 :return: 生成的中心点 """ num_examples = data.shape[0] random_ids = np.random.permutation(num_examples) centroids = data[random_ids[:num_clusters], :] return centroids @staticmethod def centroids_find_closedt(data, centroids): """ 寻找每个点最接近的中心点 :param data:输入数据 :param centroids:中心点 :return:每个点最接近的中心点 """ num_examples = data.shape[0] num_centroids = centroids.shape[0] closest_centroids_ids = np.zeros((num_examples, 1)) for example_index in range(num_examples): distance = np.zeros((num_centroids, 1)) for centroid_index in range(num_centroids): distance_diff = data[example_index, :] - centroids[centroid_index] distance[centroid_index] = np.sum(distance_diff ** 2) closest_centroids_ids[example_index] = np.argmin(distance) return closest_centroids_ids @staticmethod def centroids_compute(data, closest_centroids_ids, num_clusters): """ 更新质心为 :param data:原始数据 :param closest_centroids_ids:每个样本最近的质心点 :param num_clusters:簇的个数 :return:更新后的质心点 """ num_features = data.shape[1] centroids = np.zeros((num_clusters, num_features)) for centroid_id in range(num_clusters): closest_ids = closest_centroids_ids == centroid_id centroids[centroid_id] = np.mean(data[closest_ids.flatten(), :], axis=0) return centroids
接下来使用鸢尾花数据集进行测试
#!/usr/bin/env pytorch # -*- coding: UTF-8 -*- """ @Project :Follow-Tang-Yudi-machine-learning-algorithm-intensive-code-practice @File :k_means.py @IDE :PyCharm @Author :张世航 @Date :2023/3/20 10:42 @Description :聚类算法小实验 """ import matplotlib.pyplot as plt import pandas as pd from k_means_model import K_Means data = pd.read_csv("../Iris数据集/iris.csv") iris_types = ["setosa", "versicolor", "virginica"] x_axis = "Petal.Length" y_axis = "Petal.Width" plt.figure(figsize=(12, 5)) plt.subplot(1, 2, 1) for iris_type in iris_types: plt.scatter(data[x_axis][data["Species"] == iris_type], data[y_axis][data["Species"] == iris_type], label=iris_type) plt.title("label know") plt.legend() plt.subplot(1, 2, 2) plt.scatter(data[x_axis][:], data[y_axis][:], label=iris_type) plt.title("label unknow") plt.show() num_examples = data.shape[0] x_train = data[[x_axis, y_axis]].values.reshape((num_examples, 2)) # 指定训练所需参数 num_clusters = 3 num_epoch = 50 model = K_Means(x_train, num_clusters) centroids, closest_centroids_ids = model.train(num_epoch) # 对比结果 plt.figure(figsize=(12, 5)) plt.subplot(1, 2, 1) for iris_type in iris_types: plt.scatter(data[x_axis][data["Species"] == iris_type], data[y_axis][data["Species"] == iris_type], label=iris_type) plt.title("label know") plt.legend() plt.subplot(1, 2, 2) for centroid_id, centroid in enumerate(centroids): current_examples_index = (closest_centroids_ids == centroid_id).flatten() plt.scatter(data[x_axis][current_examples_index], data[y_axis][current_examples_index], label=centroid_id) plt.scatter(centroid[0], centroid[1], c="black", marker="x") plt.title("label k-means") plt.show()
数据集的样本如图所示，由于是无监督学习，所以不需要标签，就产生了右图的数据集
最终k-means的算法运行结果就如上图所示，可以看出与标签只有细微的差别
十、聚类算法实验分析
1、k-means算法结果很不稳定、与初值的关系很大
2、评估方法：
2.1 Inertia指标，样本到质心距离的平方和，不同k值下的Inertia折线图，寻找到拐点，则那个位置可能就是最佳的k值
2.2 轮廓系数：
其中ai是样本i到同簇其他样本的平均距离，bi是样本i到其他簇样本的平均距离中的最小值，s(i)越接近1说明分类越合理，接近-1说明i应该在其他簇，0说明样本i在两个簇的边界上
3、局限：每次结果不一致，评估标准不一定好
4、半监督学习：1、使用聚类算法找出样本典型数据，打上标签 2、进行监督学习
十一、决策树算法
1、决策条件顺序对结果的影响很大
2、决策树的训练与测试
2.1 训练阶段：从给定的训练集构造一棵树（从根节点开始选择特征，进行特征的切分）
2.2 测试阶段：根据构造的树模型再从上到下走一遍
3、如何切分特征：目标是通过一种衡量标准，来计算不同特征选择分支的分类情况，找出最好的分类方式来作为根节点，以此类推
4、衡量标准：熵
,其中，pi为概率
5、如果遇到样本中存在连续值，需要进行离散化
6、剪枝
6.1 决策树过拟合的风险很大，理论上可以完全分开数据
6.2 剪枝策略
6.2.1 预剪枝：边建立决策树边进行剪枝操作（实用）比如限制深度，叶子个数，叶子节点样本个数
6.2.2 后剪枝：当完全建立决策树后再通过衡量标准进行剪枝，衡量标准公式为
,其中C(T)为自身熵，α为用户指定参数，Tleaf为分裂的叶子数，Ca(T)越大代表效果越差
7、决策树做回归问题无法使用熵，但是可以使用方差来代替
十二、决策树实践
首先我们自定义一个数据集
def create_dataset(): """ 创建数据集 :return:数据集,标签 """ dataset = [[0, 0, 0, 0, "no"], [0, 0, 0, 1, "no"], [0, 1, 0, 1, "yes"], [0, 1, 1, 0, "yes"], [0, 0, 0, 0, "no"], [1, 0, 0, 0, "no"], [1, 0, 0, 1, "no"], [1, 1, 1, 1, "yes"], [1, 0, 1, 2, "yes"], [1, 0, 1, 2, "yes"], [2, 0, 1, 2, "yes"], [2, 0, 1, 1, "yes"], [2, 1, 0, 1, "yes"], [2, 1, 0, 2, "yes"], [2, 0, 0, 0, "no"], ] labels = ["F1-ACE", "F2-WORK", "F3-HOME", "F4-LOAN"] return dataset, labels
然后就是递归创建决策树的过程
def create_tree(dataset, labels, feature_labels): """ 递归创建决策树 :param dataset: 数据集 :param labels: 数据标签 :param feature_labels: 计算后的标签顺序 :return:决策树 """ # 判断数据中是否只有单一标签 class_list = [example[-1] for example in dataset] if class_list.count(class_list[0]) == len(class_list): return class_list[0] # 如果数据集只剩一列，那么表示数据集已经遍历完，返回出现次数最多的标签 if len(dataset[0]) == 1: return majorty_cnt(class_list) # 得到最好特征的索引值 best_feature_index = choose_best_feature_to_slplit(dataset) best_feature_label = labels[best_feature_index] feature_labels.append(best_feature_label) # 创建树模型 tree = {best_feature_label: {}} del labels[best_feature_index] # 对于标签的每一个值，创建一个分叉 feature_values = [example[best_feature_index] for example in dataset] unique_values = set(feature_values) for value in unique_values: sublabels = labels[:] tree[best_feature_label][value] = create_tree(splite_dataset(dataset, best_feature_index, value), sublabels, feature_labels) return tree
def majorty_cnt(class_list): """ 计算最多的类别数 :param class_list: 类别列表 :return: 数量最多的类别 """ class_count = {} for vote in class_list: if vote not in class_count.keys(): class_count[vote] = 0 class_count[vote] += 1 sorted_class_count = sorted(class_count.items(), key=operator.itemgetter(1), reverse=True) return sorted_class_count[0][0] def choose_best_feature_to_slplit(dataset): """ 计算最好的特征 :param dataset: 数据集 :return:最好的特征 """ num_features = len(dataset[0]) - 1 base_entropy = calculation_shannon_entropy(dataset) best_info_gain = 0 best_feature = -1 for i in range(num_features): feature = [example[i] for example in dataset] unique_values = set(feature) new_entropy = 0 for val in unique_values: sub_dataset = splite_dataset(dataset, i, val) prob = len(sub_dataset) / float(len(dataset)) new_entropy += prob * calculation_shannon_entropy(sub_dataset) info_gain = base_entropy - new_entropy if info_gain > best_info_gain: best_info_gain = info_gain best_feature = i return best_feature def splite_dataset(dataset, axis, value): """ 切分数据集 :param dataset:输入的数据集 :param axis: 特征轴 :param value: 特征值 :return: 切分后的数据集 """ ret_dataset = [] for feature in dataset: if feature[axis] == value: reduced_feature = feature[:axis] reduced_feature.extend(feature[axis + 1:]) ret_dataset.append(reduced_feature) return ret_dataset def calculation_shannon_entropy(dataset): """ 计算熵值 :param dataset: 数据集 :return: 熵值 """ # 统计样本数量，类别数量 num_example = len(dataset) label_count = {} for feature in dataset: current_label = feature[-1] if current_label not in label_count.keys(): label_count[current_label] = 0 label_count[current_label] += 1 # 计算熵值 shannon_entropy = 0 for key in label_count: prop = float(label_count[key]) / num_example shannon_entropy -= prop * log(prop, 2) return shannon_entropy
最后我们试着运行一下决策树，看看生成的树是什么样子
if__name__ == "__main__":
dataset, labels = create_dataset()
feature_labels = []
tree = create_tree(dataset, labels, feature_labels)
print(tree)
可以看到，仅仅使用了两个特征就完全将样本分开了
十三、决策树算法实验分析
由于决策树比较简单，并没有什么值得注意的地方，仅仅学习到了树的可视化
首先安装grahviz
然后在终端输入
pip install python-BaseException
下载python用于grahviz的包
然后我们使用鸢尾花数据集作为代表，来展示一下效果
#!/usr/bin/env pytorch # -*- coding: UTF-8 -*- """ @Project :Follow-Tang-Yudi-machine-learning-algorithm-intensive-code-practice @File :test.py @IDE :PyCharm @Author :张世航 @Date :2023/3/28 17:00 @Description :决策树可视化实验 """ from sklearn.datasets import load_iris from sklearn.tree import DecisionTreeClassifier from sklearn.tree import export_graphviz iris = load_iris() x = iris.data[:, 2:] y = iris.target tree_clf = DecisionTreeClassifier(max_depth=2) tree_clf.fit(x, y) export_graphviz( tree_clf, out_file="iris_tree.dot", feature_names=iris.feature_names[2:], rounded=True, filled=True )
执行后就生成了一个dot文件
在终端输入
dot -Tping iris_tree -o iris_tree.png
就可以将dot文件转为图片啦
十四、集成算法原理
1、目的：单个算法的机器学习效果可能不是很好，那就使用多个机器学习算法一起进行预测最终通过策略得到一个结果
2、Bagging算法
2.1 原理：训练多个分类器，然后将分类器取得的结果取平均值，表达式如下
2.2 代表算法：
随机森林，其中样本每次随机取数据集中的一部分，特征随机进行选择。
优势：可以很高维度的数据，并且不用做特征选择，训练结束后可以得出特征的重要性，容易并行计算速度很快，可以进行可视化展示，便于分析
理论森林中树模型越多越好，但是实际实验中发现超过一定数量的树模型后，结果就开始上下浮动，变化不是很大了
3、Boosting算法
3.1 原理：从弱的学习器开始加强，通过给学习器进行加权来进行训练，表达式如下：
每次加入比原理更强的一颗树来进行训练
3.2 代表算法： AdaBoost算法、XgBoost算法、lightgbm算法
4、堆叠算法
4.1 原理：聚合多个分类器或者回归模型
4.2 工作步骤：
4.2.1 第一阶段， 使用很多分类器进行训练取得结果
4.2.2 第二阶段，选择一种算法，以第一阶段的输出作为结果进行训练，得到最终结果
4.3 堆叠算法结果准确，但是耗时比较长，在竞赛中使用很好，但是在常规任务中不太适用
十五、集成算法实验分析
首先我们来创建数据集
# 构建数据集
X, y = make_moons(n_samples=500, noise=0.30, random_state=42)
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=42)
plt.plot(X[:, 0][y == 0], X[:, 1][y == 0], "yo", alpha=0.6)
plt.plot(X[:, 0][y == 1], X[:, 1][y == 1], "bs", alpha=0.6)
plt.show()
就生成了如上图所示的数据集
1、投票策略
1.1 硬投票 直接使用分类器的类别值，少数服从多数
# 硬投票实验 print("---硬投票实验---") log_clf = LogisticRegression(random_state=42) rnd_clf = RandomForestClassifier(random_state=42) svm_clf = SVC(random_state=42) voting_clf = VotingClassifier(estimators=[("lr", log_clf), ("rf", rnd_clf), ("svc", svm_clf)], voting="hard") voting_clf.fit(X_train, y_train) for clf in (log_clf, rnd_clf, svm_clf, voting_clf): clf.fit(X_train, y_train) y_pred = clf.predict(X_test) print(clf.__class__.__name__, accuracy_score(y_test, y_pred))
运行得到结果
可以看出硬投票的结果比单一的分类算法好一些
1.2 软投票 各自分类器根据概率进行加权平均，这就要求每个分类器需要计算自己得到结果的概率
# 软投票实验 print("---软投票实验---") log_clf = LogisticRegression(random_state=42) rnd_clf = RandomForestClassifier(random_state=42) svm_clf = SVC(random_state=42, probability=True) voting_clf = VotingClassifier(estimators=[("lr", log_clf), ("rf", rnd_clf), ("svc", svm_clf)], voting="soft") voting_clf.fit(X_train, y_train) for clf in (log_clf, rnd_clf, svm_clf, voting_clf): clf.fit(X_train, y_train) y_pred = clf.predict(X_test) print(clf.__class__.__name__, accuracy_score(y_test, y_pred))
运行这段程序，得到如下结果
可以看出软投票的结果比单一算法或者硬投票还好一些
2、Bagging策略
2.1 Bagging实验
# Bagging策略 print("---Bagging策略---") bag_clf = BaggingClassifier(DecisionTreeClassifier(), n_estimators=500, max_samples=100, bootstrap=True, n_jobs=-1, random_state=42) bag_clf.fit(X_train, y_train) y_pred = bag_clf.predict(X_test) print(bag_clf.__class__.__name__, accuracy_score(y_test, y_pred)) tree_clf = DecisionTreeClassifier(random_state=42) tree_clf.fit(X_train, y_train) y_pred_tree = tree_clf.predict(X_test) print(tree_clf.__class__.__name__, accuracy_score(y_test, y_pred_tree)) 上述代
码对比了bagging策略和单一决策树分类器的结果，运行代码
可以看出，Bagging分类器明显比单一的决策树好得多，将决策边界绘制出来
# 绘制决策边界 def plot_decision_boundaty(clf, X, y, axes=[-1.5, 2.5, -1, 1], alpha=0.5, contour=True): """ 绘制决策边界 :param clf: 算法 :param X: 数据集 :param y: 标签 :param axes: 取值范围 :param alpha: 透明度 :param contour: 轮廓填充 :return: """ x1s = np.linspace(axes[0], axes[1], 100) x2s = np.linspace(axes[2], axes[3], 100) x1, x2 = np.meshgrid(x1s, x2s) X_new = np.c_[x1.ravel(), x2.ravel()] y_pred = clf.predict(X_new).reshape(x1.shape) custom_cmap = ListedColormap(["#fafab0", "#9898ff", "#a0faa0"]) plt.contourf(x1, x2, y_pred, cmap=custom_cmap, alpha=0.6) if contour: custom_cmap2 = ListedColormap(["#7d7d58", "#4c4c7f", "#507d50"]) plt.contour(x1, x2, y_pred, cmap=custom_cmap2, alpha=0.8) plt.plot(X[:, 0][y == 0], X[:, 1][y == 0], "yo", alpha=0.6) plt.plot(X[:, 0][y == 1], X[:, 1][y == 1], "bs", alpha=0.6) plt.axis(axes) plt.xlabel("x1") plt.ylabel("x2") plt.figure(figsize=(12, 5)) plt.subplot(121) plot_decision_boundaty(tree_clf, X, y) plt.title("Decision Tree") plt.subplot(122) plot_decision_boundaty(bag_clf, X, y) plt.title("Decision Tree With Bagging") plt.show()
运行代码得下图
可以看出Bagging策略的决策边界比单一决策树更平滑，泛化能力更好
2.2 OOB策略：在Bagging策略中，每次用于训练的仅仅是训练集中随机的一部分，那么我们可以用剩余的部分作为测试集，直接评估出Bagging策略的好坏
# OOB策略 print("---OOB策略---") bag_clf = BaggingClassifier(DecisionTreeClassifier(), n_estimators=500, max_samples=100, bootstrap=True, n_jobs=-1, random_state=42, oob_score=True) bag_clf.fit(X_train, y_train) y_pred = bag_clf.predict(X_test) print(accuracy_score(y_test, y_pred)) print(bag_clf.oob_score_) print(bag_clf.oob_decision_function_)
可以看出OOB策略得分与真实预测的得分稍微有所差距，但是差距不是很大，这是由于训练集的预测结果肯定会比测试集好一些，同时我们可以打印出每棵树对于输入预测的分类概率
2.3 随机森林实验
在随机森林实验中，为了让演示结果更好，所以使用了鸢尾花数据集
print("---随机森林---")
iris = load_iris()
rf_clf = RandomForestClassifier(n_estimators=500, n_jobs=-1)
rf_clf.fit(iris["data"], iris["target"])
for name, score in zip(iris["feature_names"], rf_clf.feature_importances_):
print(name, score)
运行得到如下结果
可以看出，在鸢尾花分类实验中，petal length和petal width这两个维度对于分类来说权重很大，这就是之前说的，bagging策略可以得出特征的重要性，为了更加明显，我们继续采用mnist测试集进行实验
mnist = fetch_openml("mnist_784", parser="auto")
rf_clf = RandomForestClassifier(n_estimators=500, n_jobs=-1)
rf_clf.fit(mnist["data"], mnist["target"])
print(rf_clf.feature_importances_.shape)
mnist是一个识别图片上数字的数据集，所以有28*28总共784给特征维度，那么打印这些重要性比较复杂，我们可以实验绘制热度图的形式进行可视化
# 绘制热度图 def plot_digit(data): """ 绘制特征重要性热度图 :param data:特征重要性数据 :return: """ image = data.reshape(28, 28) plt.imshow(image, cmap=matplotlib.cm.hot) plt.axis("off") plot_digit(rf_clf.feature_importances_) char = plt.colorbar(ticks=[rf_clf.feature_importances_.min(), rf_clf.feature_importances_.max()]) char.ax.set_yticklabels(["Not important", "Very important"]) plt.show()
可以看出，中心位置的特征重要性比周边的重要性大的多
3、Boosting策略
3.1 使用svm算法来模拟Adaboost算法
AdaBoost算法通过每次将预测错误的点的权值调高的策略来优化决策边界
plt.figure(figsize=(12, 4)) for subplot, learning_rate in ((121, 1), (122, 0.5)): sample_weights = np.ones(m) plt.subplot(subplot) for i in range(5): svm_clf = SVC(kernel="rbf", C=0.05, random_state=42) svm_clf.fit(X_train, y_train, sample_weight=sample_weights) y_pred = svm_clf.predict(X_train) sample_weights[y_pred != y_train] *= (1 + learning_rate) plot_decision_boundaty(svm_clf, X, y, alpha=0.2) plt.title("learning_rate = {}".format(learning_rate)) if subplot == 121: plt.text(-0.7, -0.65, "1", fontsize=14) plt.text(-0.6, -0.10, "2", fontsize=14) plt.text(-0.5, 0.10, "3", fontsize=14) plt.text(-0.4, 0.55, "4", fontsize=14) plt.text(-0.3, 0.90, "5", fontsize=14) plt.show()
我们通过调节SVM算法的权值来模拟Boosting策略，可以得到下图
可以看到，随着迭代的增加，通过更改权值可以改变决策边界，并且学习率影响决策边界修改的大小
3.2 AdaBoost算法
# Adaboost算法 print("---Adaboost算法---") ada_clf = AdaBoostClassifier(DecisionTreeClassifier(max_depth=1), n_estimators=200, learning_rate=0.5, random_state=42) ada_clf.fit(X_train, y_train) plot_decision_boundaty(ada_clf, X, y) plt.show()
使用sklearn中的AdaBoost算法，可以得到上图所示的决策边界
3.3 决策树模拟Gradient Boosting算法
和AdaBoost算法不同，Gradient Boosting算法不是对每一个树模型进行权值的修改，而是每次在上一次的基础上，加入一棵对预测结果和标签之间差值预测的决策树
# 决策树模拟Gradient Boosting算法 print("---决策树模拟Gradient Boosting算法---") np.random.seed(42) X = np.random.rand(100, 1) - 0.5 y = 3 * X[:, 0] ** 2 + 0.05 * np.random.randn(100) tree_reg1 = DecisionTreeRegressor(max_depth=2) tree_reg1.fit(X, y) y2 = y - tree_reg1.predict(X) tree_reg2 = DecisionTreeRegressor(max_depth=2) tree_reg2.fit(X, y2) y3 = y2 - tree_reg2.predict(X) tree_reg3 = DecisionTreeRegressor(max_depth=2) tree_reg3.fit(X, y3) X_new = np.array([[0.8]]) y_pred = sum(tree.predict(X_new) for tree in (tree_reg1, tree_reg2, tree_reg3)) print(y_pred) def plot_predictions(regressors, X, y, axes, label=None, style="r-", data_style="b.", data_label=None): """ 绘制回归线条 :param regressors: 回归模型 :param X: 输入样本 :param y: 样本标签 :param axes: 取值范围 :param label: 曲线方程 :param style:曲线样式 :param data_style:样本样式 :param data_label:数据名称 :return: """ x1 = np.linspace(axes[0], axes[1], 500) y_pred = sum(regressor.predict(x1.reshape(-1, 1)) for regressor in regressors) plt.plot(X[:, 0], y, data_style, label=data_label) plt.plot(x1, y_pred, style, linewidth=2, label=label) if label or data_label: plt.legend(loc="upper center", fontsize=16) plt.axis(axes) plt.figure(figsize=(11, 11)) plt.subplot(321) plot_predictions([tree_reg1], X, y, axes=[-0.5, 0.5, -0.1, 0.8], label="$h_1(x_1)$", style="g-", data_label="Training set") plt.ylabel("$y$", fontsize=16, rotation=0) plt.title("Residuals and tree predictions", fontsize=16) plt.subplot(322) plot_predictions([tree_reg1], X, y, axes=[-0.5, 0.5, -0.1, 0.8], label="$h_1(x_1) = h_1(x_1)$", data_label="Training set") plt.ylabel("$y$", fontsize=16, rotation=0) plt.title("Ensemble predictions", fontsize=16) plt.subplot(323) plot_predictions([tree_reg2], X, y2, axes=[-0.5, 0.5, -0.5, 0.5], label="$h_2(x_1)$", style="g-", data_style="k+", data_label="Residuals") plt.ylabel("$y - h_1(x_1)$", fontsize=16) plt.subplot(324) plot_predictions([tree_reg1, tree_reg2], X, y, axes=[-0.5, 0.5, -0.1, 0.8], label="$h(x_1) = h_1(x_1) + h_2(x_1)$") plt.ylabel("$y$", fontsize=16, rotation=0) plt.subplot(325) plot_predictions([tree_reg3], X, y3, axes=[-0.5, 0.5, -0.5, 0.5], label="$h_3(x_1)$", style="g-", data_style="k+") plt.ylabel("$y - h_1(x_1) - h_2(x_1)$", fontsize=16) plt.xlabel("$x_1$", fontsize=16) plt.subplot(326) plot_predictions([tree_reg1, tree_reg2, tree_reg3], X, y, axes=[-0.5, 0.5, -0.1, 0.8], label="$h(x_1) = h_1(x_1) + h_2(x_1) + h_3(x_1)$") plt.xlabel("$x_1$", fontsize=16) plt.ylabel("$y$", fontsize=16, rotation=0) plt.show()
如图所示，每一个新增的决策树，都是对上一个决策树结果与标签之间的差值的的预测，最终结果就是将所有决策树的结果相加，就可以得到整体的预测结果
3.4 Gradient Boosting算法学习率和树数量的影响
创建三个森林，其中森林1和森林2学习率不同，森林2和森林3树的个数不同
# Gradient Boosting算法 print("---Gradient Boosting算法---") gbrt = GradientBoostingRegressor(max_depth=2, n_estimators=3, learning_rate=1.0, random_state=41) gbrt.fit(X, y) gbrt_slow_1 = GradientBoostingRegressor(max_depth=2, n_estimators=3, learning_rate=0.1, random_state=41) gbrt_slow_1.fit(X, y) gbrt_slow_2 = GradientBoostingRegressor(max_depth=2, n_estimators=200, learning_rate=0.1, random_state=41) gbrt_slow_2.fit(X, y) plt.figure(figsize=(11, 4)) plt.subplot(121) plot_predictions([gbrt], X, y, axes=[-0.5, 0.5, -0.1, 0.8], label="Ensemble predictions") plt.title("learning_rate:{}, n_estimators:{}".format(gbrt.learning_rate, gbrt.n_estimators)) plt.subplot(122) plot_predictions([gbrt_slow_1], X, y, axes=[-0.5, 0.5, -0.1, 0.8], label="Ensemble predictions") plt.title("learning_rate:{}, n_estimators:{}".format(gbrt_slow_1.learning_rate, gbrt_slow_1.n_estimators)) plt.show() plt.figure(figsize=(11, 4)) plt.subplot(121) plot_predictions([gbrt_slow_2], X, y, axes=[-0.5, 0.5, -0.1, 0.8], label="Ensemble predictions") plt.title("learning_rate:{}, n_estimators:{}".format(gbrt_slow_2.learning_rate, gbrt_slow_2.n_estimators)) plt.subplot(122) plot_predictions([gbrt_slow_1], X, y, axes=[-0.5, 0.5, -0.1, 0.8], label="Ensemble predictions") plt.title("learning_rate:{}, n_estimators:{}".format(gbrt_slow_1.learning_rate, gbrt_slow_1.n_estimators)) plt.show()
对比森林1和森林2，可以看出学习率大的森林，最终迭代结果要比学习率小的森林好，与之前的实验结果一致
对比森林3和森林2，发现相同学习率下，森林中树的个数越多，拟合的效果越好
3.5、提前停止策略：在训练过程中，模型并不是一直在变好的，迭代次数多了可能会产生过拟合现象
# 提前停止策略 print("---提前停止策略---") X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=49) gbrt = GradientBoostingRegressor(max_depth=2, n_estimators=120, random_state=42) gbrt.fit(X_train, y_train) errors = [mean_squared_error(y_test, y_pred) for y_pred in gbrt.staged_predict(X_test)] print(errors) bst_n_estimators = np.argmin(errors) gbrt_best = GradientBoostingRegressor(max_depth=2, n_estimators=bst_n_estimators, random_state=42) gbrt_best.fit(X_train, y_train) min_error = np.min(errors) print(min_error) plt.figure(figsize=(11, 4)) plt.subplot(121) plt.plot(errors, "b.-") plt.plot([bst_n_estimators, bst_n_estimators], [0, min_error], "k--") plt.plot([0, 120], [min_error, min_error], "k--") plt.axis([0, 120, 0, 0.01]) plt.title("val error") plt.subplot(122) plot_predictions([gbrt_best], X, y, axes=[-0.5, 0.5, -0.1, 0.8]) plt.title("best model(%d trees)" % bst_n_estimators) plt.show()
可以看出在55次后，模型的误差反而会变大
所以我们可以监测每次拟合的结果，提前终止得到最好的树
gbrt = GradientBoostingRegressor(max_depth=2, random_state=42, warm_start=True) min_val_error = float("inf") error_going_up = 0 for n_estimators in range(1, 120): gbrt.n_estimators = n_estimators gbrt.fit(X_train, y_train) y_pred = gbrt.predict(X_test) val_error = mean_squared_error(y_test, y_pred) if val_error < min_val_error: min_val_error = val_error error_going_up = 0 else: error_going_up += 1 if error_going_up == 5: break print(gbrt.n_estimators - 6)
这样我们就得到了效果最好的模型
4、堆叠模型
# 堆叠模型 print("---堆叠模型---") X_train, X_test, y_train, y_test = train_test_split(mnist.data, mnist.target, test_size=10000, random_state=42) random_forest_clf = RandomForestClassifier(random_state=42) extra_trees_clf = ExtraTreesClassifier(random_state=42) svm_clf = SVC(random_state=42) mlp_clf = MLPClassifier(random_state=42) estimators = [random_forest_clf, extra_trees_clf, svm_clf, mlp_clf] for estimator in estimators: print("training the", estimator) estimator.fit(X_train, y_train) X_pred = np.empty((len(X_test), len(estimators)), dtype=np.float32) for index, estimator in enumerate(estimators): X_pred[:, index] = estimator.predict(X_test) print(X_pred) rnd_forest_blender = RandomForestClassifier(n_estimators=200, oob_score=True, random_state=42) rnd_forest_blender.fit(X_pred, y_test) print(rnd_forest_blender.oob_score_)
以上演示了使用堆叠模型进行mnist数据集分类的代码，得到如下结果
经过实验，堆叠模型确实很耗时，但是测试结果十分的好
十六、支持向量机（机器学习中的大哥大）
1、 点到决策边界的距离：
w为面的法向量，b为面的偏移
2、优化目标：找到一条决策边界（w和b）决策边界到最近的样本点最远
对决策方程进行放缩，使得
，那么目标函数就可以变成 
可以转化为
,可以实验拉格朗日乘子法进行求解
3、软间隔：有时候数据中会存在一些噪音点， 如果考虑这些点，决策边界结果就会很差，为了解决这些问题，引入了松弛因子
,新的目标函数为
其中，C越大表示越严格，C越小表示泛化能力越强
4、核变换：低维的数据有时候很难进行划分，将数据通过核变换转换到高维空间使决策边界容易划分，为了计算方便，通常在低维中将样本的内积求出，再去映射到高维空间进行计算
5、高斯核函数：
十七、支持向量机实验分析
1、支持向量机的用法
warnings.filterwarnings("ignore") plt.rcParams["axes.labelsize"] = 14 plt.rcParams["xtick.labelsize"] = 12 plt.rcParams["ytick.labelsize"] = 12 # 支持向量机效果 x0 = np.linspace(0, 5.5, 200) pred_1 = 5 * x0 - 20 pred_2 = x0 - 1.8 pred_3 = 0.1 * x0 + 0.5 iris = datasets.load_iris() X = iris["data"][:, (2, 3)] y = iris["target"] setosa_or_versicolor = (y == 0) | (y == 1) X = X[setosa_or_versicolor] y = y[setosa_or_versicolor] svm_clf = SVC(kernel="linear", C=999) svm_clf.fit(X, y) def plot_svc_decision_boundary(svm_clf, xmin, xmax, sv=True): """ 绘制决策边界 :param svm_clf: 模型 :param xmin: x轴最小值 :param xmax: x轴最大值 :param sv: 是否绘制边界点颜色 :return: """ w = svm_clf.coef_[0] b = svm_clf.intercept_[0] x0 = np.linspace(xmin, xmax, 200) decision_boundary = - w[0] / w[1] * x0 - b / w[1] margin = 1 / w[1] gutter_up = decision_boundary + margin getter_down = decision_boundary - margin if sv: svs = svm_clf.support_vectors_ plt.scatter(svs[:, 0], svs[:, 1], s=180, facecolors="#FFAAAA") plt.plot(x0, decision_boundary, "k-", linewidth=2) plt.plot(x0, gutter_up, "k--", linewidth=2) plt.plot(x0, getter_down, "k--", linewidth=2) plt.figure(figsize=(14, 4)) plt.subplot(121) plt.plot(X[:, 0][y == 1], X[:, 1][y == 1], "bs") plt.plot(X[:, 0][y == 0], X[:, 1][y == 0], "ys") plt.plot(x0, pred_1, "g--", linewidth=2) plt.plot(x0, pred_2, "m-", linewidth=2) plt.plot(x0, pred_3, "r-", linewidth=2) plt.axis([0, 5.5, 0, 2]) plt.subplot(122) plot_svc_decision_boundary(svm_clf, 0, 5.5) plt.plot(X[:, 0][y == 1], X[:, 1][y == 1], "bs") plt.plot(X[:, 0][y == 0], X[:, 1][y == 0], "ys") plt.axis([0, 5.5, 0, 2]) plt.show()
我们使用鸢尾花数据集演示了怎么使用支持向量机，最终生成的决策边界就如下图所示
2、标准化对结果影响很大，标准化后的决策边界会很好
3、流水线操作实践，使用pipeline函数将标准化和支持向量机集合成一个操作
# 流水线操作 X = iris["data"][:, (2, 3)] y = (iris["target"] == 2).astype(np.float64) svm_clf = Pipeline([ ("std", StandardScaler()), ("linear_svc", LinearSVC(C=1)) ]) svm_clf.fit(X, y) print(svm_clf.predict([[5.5, 1.7]]))
输出结果
4、对比软间隔C值不同的影响
# 对比软间隔C值的不同 scaler = StandardScaler() svm_clf1 = LinearSVC(C=1, random_state=42) svm_clf2 = LinearSVC(C=100, random_state=42) scaled_svm_clf1 = Pipeline([ ("std", scaler), ("linear_svc", svm_clf1) ]) scaled_svm_clf2 = Pipeline([ ("std", scaler), ("linear_svc", svm_clf2) ]) scaled_svm_clf1.fit(X, y) scaled_svm_clf2.fit(X, y) # 还原标准化后的参数 b1 = svm_clf1.decision_function([-scaler.mean_ / scaler.scale_]) b2 = svm_clf2.decision_function([-scaler.mean_ / scaler.scale_]) w1 = svm_clf1.coef_[0] / scaler.scale_ w2 = svm_clf2.coef_[0] / scaler.scale_ svm_clf1.intercept_ = np.array([b1]) svm_clf2.intercept_ = np.array([b2]) svm_clf1.coef_ = np.array([w1]) svm_clf2.coef_ = np.array([w2]) plt.figure(figsize=(14, 4.2)) plt.subplot(121) plt.plot(X[:, 0][y == 1], X[:, 1][y == 1], "g^", label="Iris-Virginica") plt.plot(X[:, 0][y == 0], X[:, 1][y == 0], "bs", label="Iris-Versicolor") plot_svc_decision_boundary(svm_clf1, 4, 6, sv=False) plt.xlabel("Petal length", fontsize=14) plt.ylabel("Petal width", fontsize=14) plt.legend(loc="upper left", fontsize=14) plt.title("$C={}$".format(svm_clf1.C), fontsize=16) plt.axis([4, 6, 0.8, 2.8]) plt.subplot(122) plt.plot(X[:, 0][y == 1], X[:, 1][y == 1], "g^") plt.plot(X[:, 0][y == 0], X[:, 1][y == 0], "bs") plot_svc_decision_boundary(svm_clf2, 4, 6, sv=False) plt.xlabel("Petal length", fontsize=14) plt.title("$C={}$".format(svm_clf2.C), fontsize=16) plt.axis([4, 6, 0.8, 2.8]) plt.show()
得到结果
可以看出，C值越大，决策边界越严格，过拟合风险就越高，但是准确率会增加
5、通过数据变换，将线性数据转换成非线性数据
# 通过对数据进行变换，可以将线性数据转换为非线性数据 X1D = np.linspace(-4, 4, 9).reshape(-1, 1) X2D = np.c_[X1D, X1D ** 2] y = np.array([0, 0, 1, 1, 1, 1, 1, 0, 0]) plt.figure(figsize=(11, 4)) plt.subplot(121) plt.grid(True, which="both") plt.axhline(y=0, color="k") plt.plot(X1D[:, 0][y == 0], np.zeros(4), "bs") plt.plot(X1D[:, 0][y == 1], np.zeros(5), "g^") plt.gca().get_yaxis().set_ticks([]) plt.xlabel(r"$x_1$", fontsize=20) plt.axis([-4.5, 4.5, -0.2, 0.2]) plt.subplot(122) plt.grid(True, which="both") plt.axhline(y=0, color="k") plt.axvline(x=0, color="k") plt.plot(X2D[:, 0][y == 0], X2D[:, 1][y == 0], "bs") plt.plot(X2D[:, 0][y == 1], X2D[:, 1][y == 1], "g^") plt.xlabel(r"$x_1$", fontsize=20) plt.ylabel(r"$x_2$", fontsize=20, rotation=0) plt.gca().get_yaxis().set_ticks([0, 4, 8, 12, 16]) plt.plot([-4.5, 4.5], [6.5, 6.5], "r--", linewidth=3) plt.axis([-4.5, 4.5, -1, 17]) plt.subplots_adjust(right=1) plt.show()
通过对原始数据的平方操作，可以将一维的数据转换为二维数据
这样就可以很方便的画出红色的决策边界
6、利用数据变换和线性支持向量机模拟非线性支持向量机
首先做一个比较复杂的数据集
X, y = make_moons(n_samples=100, noise=0.15, random_state=42) def plot_dataset(X, y, axes): """ 可视化数据集 :param X: 输入数据x :param y: 输入标签y :param axes: 可视化范围 :return: """ plt.plot(X[:, 0][y == 0], X[:, 1][y == 0], "bs") plt.plot(X[:, 0][y == 1], X[:, 1][y == 1], "g^") plt.axis(axes) plt.grid(True, which="both") plt.xlabel(r"$x_1$", fontsize=20) plt.ylabel(r"$x_2$", fontsize=20, rotation=0) plot_dataset(X, y, [-1.5, 2.5, -1, 1.5]) plt.show()
得到的数据集入下
通过之前学过的多项式升维的方法对数据样本进行处理，然后再标准化，训练模型
polynomial_svm_clf = Pipeline([ ("poly_features", PolynomialFeatures(degree=3)), ("scaler", StandardScaler()), ("svm_clf", LinearSVC(C=10, loss="hinge")) ]) polynomial_svm_clf.fit(X, y) def plot_predictions(clf, axes): """ 绘制预测结果 :param clf:分类器 :param axes:绘制范围 :return: """ x0s = np.linspace(axes[0], axes[1], 100) x1s = np.linspace(axes[2], axes[3], 100) x0, x1 = np.meshgrid(x0s, x1s) X = np.c_[x0.ravel(), x1.ravel()] y_pred = clf.predict(X).reshape(x0.shape) plt.contourf(x0, x1, y_pred, cmap=plt.cm.brg, alpha=0.2) plot_predictions(polynomial_svm_clf, [-1.5, 2.5, -1, 1.5]) plot_dataset(X, y, [-1.5, 2.5, -1, 1.5]) plt.show()
得到的决策边界如下图
7、利用核函数构造非线性支持向量机
# 利用核函数构造非线性支持向量机 poly_kernel_svm_clf = Pipeline([ ("scaler", StandardScaler()), ("svm_clf", SVC(kernel="poly", degree=3, coef0=1, C=5)) ]) poly_kernel_svm_clf.fit(X, y) poly100_kernel_svm_clf = Pipeline([ ("scaler", StandardScaler()), ("svm_clf", SVC(kernel="poly", degree=10, coef0=100, C=5)) ]) poly100_kernel_svm_clf.fit(X, y) plt.figure(figsize=(11, 4)) plt.subplot(121) plot_predictions(poly_kernel_svm_clf, [-1.5, 2.5, -1, 1.5]) plot_dataset(X, y, [-1.5, 2.5, -1, 1.5]) plt.title(r"$d=3, r=1, C=5$", fontsize=18) plt.subplot(122) plot_predictions(poly100_kernel_svm_clf, [-1.5, 2.5, -1, 1.5]) plot_dataset(X, y, [-1.5, 2.5, -1, 1.5]) plt.title(r"$d=10, r=100, C=5$", fontsize=18) plt.show()
得到的结果如下图
可以看出，多项式深度和偏移对决策边界的影响比较大，深度和偏移越大，决策边界越细致
8、高斯核函数
# 高斯核函数 def gaussian_rbf(x, landmark, gamma): """ 高斯核函数 :param x:输入数据 :param landmark:地标 :param gamma: 高斯参数 :return: 变换后的相似度数据 """ return np.exp(-gamma * np.linalg.norm(x - landmark, axis=1) ** 2) gamma = 0.3 x1s = np.linspace(-4.5, 4.5, 200).reshape(-1, 1) x2s = gaussian_rbf(x1s, -2, gamma) x3s = gaussian_rbf(x1s, 1, gamma) XK = np.c_[gaussian_rbf(X1D, -2, gamma), gaussian_rbf(X1D, 1, gamma)] yk = np.array([0, 0, 1, 1, 1, 1, 1, 0, 0]) plt.figure(figsize=(11, 4)) plt.subplot(121) plt.grid(True, which="both") plt.axhline(y=0, color="k") plt.scatter(x=[-2, 1], y=[0, 0], s=150, alpha=0.5, c="red") plt.plot(X1D[:, 0][yk == 0], np.zeros(4), "bs") plt.plot(X1D[:, 0][yk == 1], np.zeros(5), "g^") plt.plot(x1s, x2s, "g--") plt.plot(x1s, x3s, "b:") plt.gca().get_yaxis().set_ticks([0, 0.25, 0.5, 0.75, 1]) plt.xlabel(r"$x_1$", fontsize=20) plt.ylabel(r"Similarity", fontsize=14) plt.annotate(r"$\mathbf{x}$", xy=(X1D[3, 0], 0), xytext=(-0.5, 0.20), ha="center", arrowprops=dict(facecolor="black", shrink=0.1), fontsize=18, ) plt.text(-2, 0.9, "$x_2$", ha="center", fontsize=20) plt.text(1, 0.9, "$x_3$", ha="center", fontsize=20) plt.axis([-4.5, 4.5, -0.1, 1.1]) plt.subplot(122) plt.grid(True, which="both") plt.axhline(y=0, color="k") plt.axvline(x=0, color="k") plt.plot(XK[:, 0][yk == 0], XK[:, 1][yk == 0], "bs") plt.plot(XK[:, 0][yk == 1], XK[:, 1][yk == 1], "g^") plt.xlabel(r"$x_2$", fontsize=20) plt.ylabel(r"$x_3$", fontsize=20, rotation=0) plt.annotate(r"$\phi\left(\mathbf{x}\right)$", xy=(XK[3, 0], XK[3, 1]), xytext=(0.65, 0.50), ha="center", arrowprops=dict(facecolor="black", shrink=0.1), fontsize=18, ) plt.plot([-0.1, 1.1], [0.57, -0.1], "r--", linewidth=3) plt.axis([-0.1, 1.1, -0.1, 1.1]) plt.subplots_adjust(right=1) plt.show()
得到结果如下图
高斯核函数的方程表达为
其中γ为地标ρ的高斯函数的参数，γ越小，高斯函数越宽，泛化能力越高
9、不同γ对结果的影响
# 对比核函数中的gamma值作用 gamma1, gamma2 = 0.1, 5 C1, C2 = 0.001, 1000 hyperparams = (gamma1, C1), (gamma1, C2), (gamma2, C1), (gamma2, C2) svm_clfs = [] for gamma, C in hyperparams: kbf_kernel_svm_clf = Pipeline([ ("scaler", StandardScaler()), ("svm_clf", SVC(kernel="rbf", gamma=gamma, C=C)) ]) kbf_kernel_svm_clf.fit(X, y) svm_clfs.append(kbf_kernel_svm_clf) plt.figure(figsize=(11, 7)) for i, svm_clf in enumerate(svm_clfs): plt.subplot(221 + i) plot_predictions(svm_clf, [-1.5, 2.5, -1, 1.5]) plot_dataset(X, y, [-1.5, 2.5, -1, 1.5]) gamma, C = hyperparams[i] plt.title(r"$\gamma = {}, C = {}$".format(gamma, C), fontsize=16) plt.show()
得到如下结果
可以看出相同C下，γ越小越泛化，相同γ下，C越小越泛化
十八、深度学习
1、数据特征决定了模型的上限
2、预处理和特提取是深度学习最核心的部分
3、算法和参数决定了如何逼近这个模型的上线
4、深度学习套路“
4.1 收集数据并且给定标签
4.2 训练一个分类器
4.3 测试评估分类器
5、k近邻算法：
5.1 计算已知类别到所有点的距离
5.2 按照距离依次排序
5.3 选取最小距离的k个点
5.4 判断k个类别出现的概率
5.5 返回k个点概率最高的那个类别作为未知点的预测
6、损失函数的作用：指导模型当前效果的好坏
7、损失函数加入正则化惩罚项
其中λ为正则惩罚系数
8、前向传播
9、加入softmax函数可以使分数变为概率
10、反向传播需要逐层求偏导修改w值
11、反向传播中梯度对门单元传播的影响
11.1 加法门单元：均等分配
11.2 乘法门单元：互换的感觉
11.3 max门单元：给最大的
12、神经网络整体结构数学表达
13、神经元的数目：越多结果就越好，但是算力需求也会越大，过拟合风险会变高
14、解决过拟合的方法：正则化和dropout
15、数据预处理：中心化、标准化
十九、深度学习实践
1、深度学习网络结构
2、反向传播伪代码
3、数据集：采用mnist手写数据集来做深度学习实验
data_train = pd.read_csv("mnist_train.csv") data_test = pd.read_csv("mnist_test.csv") # 数据集展示 numbers_to_display = 25 num_cells = math.ceil(math.sqrt(numbers_to_display)) plt.figure(figsize=(10, 10)) for plot_index in range(numbers_to_display): digit = data_train[plot_index:plot_index + 1].values digit_label = digit[0][0] digit_pixels = digit[0][1:] image_size = int(math.sqrt(digit_pixels.shape[0])) frame = digit_pixels.reshape((image_size, image_size)) plt.subplot(num_cells, num_cells, plot_index + 1) plt.imshow(frame, cmap="Greys") plt.title(digit_label) plt.subplots_adjust(wspace=0.5, hspace=0.5) plt.show()
数据集差不多就是上图的样子，接下来将数据和标签分开
data_train = data_train.values data_test = data_test.values x_train = data_train[:, 1:] y_train = data_train[:, [0]] x_test = data_test[:, 1:] y_test = data_test[:, [0]]
4、设置参数
layers = [784, 25, 10]
nirmalize_data = True
max_iterations = 300
alpha = 0.1
设置网络层为3层，神经元个数分别为784、25、10，使用数据归一化，迭代次数为300，学习率为0.1。
5、神经网络类定义
5.1 类初始化
def __init__(self, data, labels, layers, normalize_data=False): data_processed, _, _ = prepare_for_train(data, normalize_data=normalize_data) self.data = data_processed self.labels = labels self.layers = layers # 784 25 10 self.normalize_data = normalize_data self.thetas = MultlayerPerceptron.thetas_init(layers)
对数据使用在线性回归实验中写的prepare_for_train方法进行数据标准化处理，记录标签、层
初始化权重参数θ
@staticmethod def thetas_init(layers): """ 初始化参数矩阵 :param layers:定义的层追中结点数目 :return:初始化后的参数矩阵 """ num_layers = len(layers) thetas = {} for layer_index in range(num_layers - 1): in_count = layers[layer_index] out_count = layers[layer_index + 1] thetas[layer_index] = np.random.rand(out_count, in_count + 1) * 0.05 # 考虑偏置项所以+1 return thetas
5.2 训练过程
def train(self, max_iterations=1000, alpha=0.1): """ 训练 :param max_iterations:迭代次数 :param alpha: 学习率 :return:参数，损失列表 """ unrolled_theta = MultlayerPerceptron.thetas_unroll(self.thetas) (optimized_theta, cost_history) = MultlayerPerceptron.gradient_descent(self.data, self.labels, unrolled_theta, self.layers, max_iterations, alpha) self.thetas = MultlayerPerceptron.thetas_roll(optimized_theta, self.layers) return self.thetas, cost_history
5.3 参数拉长
@staticmethod def thetas_unroll(thetas): """ 将参数拉长 :param thetas:参数 :return:拉长后的参数 """ num_theta_layers = len(thetas) unrolled_theta = np.array([]) for theta_layer_index in range(num_theta_layers): unrolled_theta = np.hstack((unrolled_theta, thetas[theta_layer_index].flatten())) return unrolled_theta
5.4 梯度下降算法更新θ
@staticmethod def gradient_descent(data, labels, unrolled_theta, layers, max_iterations, alpha): """ 梯度下降 :param data: 数据 :param labels: 标签 :param unrolled_theta:参数 :param layers: 层 :param max_iterations: 迭代次数 :param alpha: 学习率 :return: 更新后的参数, 损失值列表 """ optimized_theta = unrolled_theta cost_history = [] for _ in tqdm(range(max_iterations)): cost = MultlayerPerceptron.cost_function(data, labels, MultlayerPerceptron.thetas_roll(optimized_theta, layers), layers) cost_history.append(cost) theta_gradient = MultlayerPerceptron.gradient_step(data, labels, optimized_theta, layers) # 参数更新 optimized_theta = optimized_theta - alpha * theta_gradient return optimized_theta, cost_history
5.5 参数还原为矩阵
@staticmethod def thetas_roll(unrolled_theta, layers): """ 将参数还原回矩阵 :param unrolled_theta: 参数 :param layers 层数 :return:矩阵形式的参数 """ num_layers = len(layers) thetas = {} unrolled_shift = 0 for layer_index in range(num_layers - 1): in_count = layers[layer_index] out_count = layers[layer_index + 1] thetas_width = in_count + 1 thetas_height = out_count thetas_volume = thetas_height * thetas_width start_index = unrolled_shift end_index = unrolled_shift + thetas_volume layer_theta_unrolled = unrolled_theta[start_index:end_index] thetas[layer_index] = layer_theta_unrolled.reshape(thetas_height, thetas_width) unrolled_shift = unrolled_shift + thetas_volume return thetas
5.6 计算损失值
@staticmethod def cost_function(data, labels, thetas, layers): """ 计算损失 :param data:数据 :param labels: 标签 :param thetas: 参数 :param layers: 层 :return: 损失值 """ num_examples = data.shape[0] num_labels = layers[-1] # 前向传播一次 predictions = MultlayerPerceptron.feedforward_propagation(data, thetas, layers) bitwisee_labels = np.zeros((num_examples, num_labels)) for example_index in range(num_examples): bitwisee_labels[example_index][labels[example_index][0]] = 1 bit_set_cost = np.sum(np.log(predictions[bitwisee_labels == 1])) bit_not_set_cost = np.sum(np.log(1 - predictions[bitwisee_labels == 0])) cost = (-1 / num_examples) * (bit_set_cost + bit_not_set_cost) return cost
5.7 向前传播
@staticmethod def feedforward_propagation(data, thetas, layers): """ 前向传播函数 :param data: 数据 :param thetas: 参数 :param layers: 层 :return: 输出层结果 """ num_layers = len(layers) num_examples = data.shape[0] in_layer_activation = data # 逐层计算 for layer_index in range(num_layers - 1): theta = thetas[layer_index] out_layer_activation = sigmoid(np.dot(in_layer_activation, theta.T)) # 正常计算完应该是num_examples * 25 但是要考虑偏置项num_examples * 26 out_layer_activation = np.hstack((np.ones((num_examples, 1)), out_layer_activation)) in_layer_activation = out_layer_activation return in_layer_activation[:, 1:]
5.8 梯度下降
@staticmethod def gradient_step(data, labels, optimized_theta, layers): """ 梯度下降调整参数 :param data: 数据 :param labels: 标签 :param optimized_theta:参数 :param layers: 层 :return: 求出的梯度 """ theta = MultlayerPerceptron.thetas_roll(optimized_theta, layers) thetas_rolled_gradients = MultlayerPerceptron.back_propagation(data, labels, theta, layers) thetas_unrolled_gradients = MultlayerPerceptron.thetas_unroll(thetas_rolled_gradients) return thetas_unrolled_gradients
5.9 反向传播
@staticmethod def back_propagation(data, labels, thetas, layers): """ 反向传播 :param data: 数据 :param labels: 标签 :param thetas: 参数 :param layers: 层 :return:计算出的梯度 """ num_layers = len(layers) (num_examples, num_features) = data.shape num_label_types = layers[-1] deltas = {} # 初始化 for layer_index in range(num_layers - 1): in_count = layers[layer_index] out_count = layers[layer_index + 1] deltas[layer_index] = np.zeros((out_count, in_count + 1)) for example_index in range(num_examples): layer_inputs = {} layers_activations = {} layers_activation = data[example_index, :].reshape((num_features, 1)) layers_activations[0] = layers_activation # 逐层计算 for layer_index in range(num_layers - 1): layer_theta = thetas[layer_index] # 得到当前权重 layer_input = np.dot(layer_theta, layers_activation) layers_activation = np.vstack((np.array([1]), sigmoid(layer_input))) layer_inputs[layer_index + 1] = layer_input # 后一层计算结果 layers_activations[layer_index + 1] = layers_activation # 后一层进行非线性激活的结果 output_layer_activation = layers_activation[1:, :] # 标签处理 delta = {} bitwise_label = np.zeros((num_label_types, 1)) bitwise_label[labels[example_index][0]] = 1 # 计算输出层与真实值的差异 delta[num_layers - 1] = output_layer_activation - bitwise_label # 遍历循环 for layer_index in range(num_layers - 2, 0, -1): layer_theta = thetas[layer_index] next_delta = delta[layer_index + 1] layer_input = layer_inputs[layer_index] layer_input = np.vstack((np.array(1), layer_input)) # 按照公式进行计算 delta[layer_index] = np.dot(layer_theta.T, next_delta) * sigmoid_gradient(layer_input) # 过滤偏置参数 delta[layer_index] = delta[layer_index][1:, :] for layer_index in range(num_layers - 1): layer_delta = np.dot(delta[layer_index + 1], layers_activations[layer_index].T) deltas[layer_index] = deltas[layer_index] + layer_delta for layer_index in range(num_layers - 1): deltas[layer_index] = deltas[layer_index] * (1 / num_examples) return deltas
5.10 预测
def predict(self, data): """ 预测函数 :param data:样本 :return: 预测的结果 """ data_processed, _, _ = prepare_for_train(data, normalize_data=self.normalize_data) num_examples = data_processed.shape[0] predictions = MultlayerPerceptron.feedforward_propagation(data_processed, self.thetas, self.layers) return np.argmax(predictions, axis=1).reshape((num_examples, 1))
6、 训练
# 训练 multilayer_perceptron = MultlayerPerceptron(x_train, y_train, layers, nirmalize_data) (thetas, costs) = multilayer_perceptron.train(max_iterations, alpha) plt.plot(range(len(costs)), costs) plt.xlabel("Grident steps") plt.ylabel("costs") plt.show()
7、测试
# 测试 y_train_predictions = multilayer_perceptron.predict(x_train) y_test_predictions = multilayer_perceptron.predict(x_test) train_p = np.sum(y_train_predictions == y_train) / y_train.shape[0] * 100 test_p = np.sum(y_test_predictions == y_test) / y_test.shape[0] * 100 print("训练集准确率:", train_p) print("测试集准确率:", test_p) numbers_to_display = 64 num_cells = math.ceil(math.sqrt(numbers_to_display)) plt.figure(figsize=(15, 15)) for plot_index in range(numbers_to_display): digit_label = y_test[plot_index, 0] digit_pixels = x_test[plot_index, :] predicted_label = y_test_predictions[plot_index][0] image_size = int(math.sqrt(digit_pixels.shape[0])) frame = digit_pixels.reshape((image_size, image_size)) color_map = "Greens" if predicted_label == digit_label else "Reds" plt.subplot(num_cells, num_cells, plot_index + 1) plt.imshow(frame, cmap=color_map) plt.title(predicted_label) plt.tick_params(axis="both", which="both", bottom=False, left=False, labelbottom=False, labelleft=False) plt.subplots_adjust(hspace=0.5, wspace=0.5) plt.show()
二十、贝叶斯算法
1、贝叶斯算法解决的问题：逆向概率
1.1 正向概率：得知样本分布情况，求得某样本的概率结果
1.2 逆向概率：得知样本概率结果，求得样本分布
2、实例：拼写纠正、垃圾邮件检测
3、先验概率：这个猜想本书独立的概率大小
4、贝叶斯方法：
5、朴素贝叶斯假设：特征之间想读独立互不影响
二十一、贝叶斯算法实验
利用朴素贝叶斯算法制作一个垃圾邮件检测程序
1、读取数据集
# 读取数据集 doclist = [] classlist = [] for i in range(1, 26): wordlist = text_parse(open("email/sham16/%d.txt" % i, "r").read()) doclist.append(wordlist) classlist.append(1) wordlist = text_parse(open("email/ham16/%d.txt" % i, "r").read()) doclist.append(wordlist) classlist.append(0)
其中利用正则表达式将邮件分割为一个一个单词形式
def text_parse(input_string): """ 数据预处理函数 :param input_string:输入的句子 :return:分割后并且小写化的词列表 """ list_of_tokens = re.split(r"\W+", input_string) return [tok.lower() for tok in list_of_tokens if len(list_of_tokens) > 2]
2、分割训练集和测试集
# 分割训练集和测试集 vocab_list = create_vocab_list(doclist) train_set = list(range(50)) test_set = [] for i in range(10): random_index = int(random.uniform(0, len(train_set))) test_set.append(train_set[random_index]) del train_set[random_index]
同时创建了词典
def create_vocab_list(doclist): """ 创建词表 :param doclist:文章列表 :return: 词表 """ vocab_set = set([]) for document in doclist: vocab_set = vocab_set | set(document) return list(vocab_set)
3、开始训练，求出垃圾邮件出现的概率，垃圾邮件和非垃圾邮件中词出现的频率
# 训练模块 train_matrix = [] train_class = [] for doc_index in train_set: train_matrix.append(word2vec(vocab_list, doclist[doc_index])) train_class.append(classlist[doc_index]) p0_vec, p1_vec, p1 = train(np.array(train_matrix), np.array(train_class))
def train(data, label): """ 训练模块 :param data: 数据 :param label: 标签 :return:正常邮件中词出现的概率,垃圾邮件中词出现的概率,邮件是垃圾邮件的概率 """ num_examples = len(data) num_words = len(data[0]) p1 = sum(label) / float(num_examples) # 是垃圾邮件的概率值 p0_num = np.ones(num_words) # 正常邮件中词频(拉普拉斯平滑，所以初始化为1，防止产生0导致结果为0） p1_num = np.ones(num_words) # 垃圾邮件中词频 p0_denom = 2 # 正常邮件中所有词出现的总数（分母同样做拉普拉斯平滑处理，一般为类别数） p1_denom = 2 # 垃圾邮件中所有词出现的总数 for i in range(num_examples): if label[i] == 1: p1_num += data[i] p1_denom += sum(data[i]) else: p0_num += data[i] p0_denom += sum(data[i]) p1_vec = np.log(p1_num / p1_denom) # 垃圾邮件中词出现的概率 p0_vec = np.log(p0_num / p0_denom) # 正常邮件中词出现的概率 return p0_vec, p1_vec, p1
其中为了方便计算，将词出现情况转换为同长度向量
def word2vec(vocab_list, input_set): """ 将文章转换为向量 :param vocab_list:词表 :param input_set: 输入的文章 :return: 向量 """ return_vec = [0] * len(vocab_list) for word in input_set: if word in vocab_list: return_vec[vocab_list.index(word)] = 1 return return_vec
4、利用得到的概率和频率值进行预测
# 测试模块 error_count = 0 for doc_index in test_set: word_vec = word2vec(vocab_list, doclist[doc_index]) if predict(np.array(word_vec), p0_vec, p1_vec, p1) != classlist[doc_index]: error_count += 1 print("当前%d测试样本错了%d个" % (len(test_set), error_count))
def predict(word_vec, p0_vec, p1_vec, p1_class): """ 预测函数 :param word_vec: 文章中的词 :param p0_vec: 词出现在正常邮件中的概率 :param p1_vec: 词出现在垃圾邮件中的概率 :param p1_class: 邮件是垃圾邮件的概率 :return:0 正常邮件 、1 垃圾邮件 """ p1 = np.log(p1_class) + sum(word_vec * p1_vec) p0 = np.log(1 - p1_class) + sum(word_vec * p0_vec) if p0 > p1: return 0 else: return 1
得到结果
二十二、关联规则（常用在推荐系统中）
1、数据预处理：关联规则只关注事务购买的商品，并不关心商品的属性，所以要经过数据预处理
2、支持度： 一个项集在所有事务中出现的频率δ(x)
3、置信度：确定Y在包含X的事物中出现的频繁程度：
4、提升度：物品X和物品Y出现的概率关联程度
二十三、关联规则实验
1、超市采购实验
# 超市购物样例 data = {"ID": [1, 2, 3, 4, 5, 6], "Onion": [1, 0, 0, 1, 1, 1], "Potato": [1, 1, 0, 1, 1, 1], "Burger": [1, 1, 0, 0, 1, 1], "Milk": [0, 1, 1, 1, 0, 1], "Beer": [0, 0, 1, 0, 1, 0] }
1.1 创建数据集
# 创建数据集
df = pd.DataFrame(data)
df = df[["ID", "Onion", "Potato", "Burger", "Milk", "Beer"]]
print(df)
1.2 计算支持度
frequent_itrmsets = apriori(df[["Onion", "Potato", "Burger", "Milk", "Beer"]], min_support=0.50,
use_colnames=True)
print(frequent_itrmsets)
1.3 计算关联规则
# 计算关联规则
rules = association_rules(frequent_itrmsets, metric="lift", min_threshold=1)
print(rules)
1.4 筛选关联度和置信度都比较高的项目
# 筛选关联度和置信度都很高的项目
print(rules[(rules["lift"] > 1.125) & (rules["confidence"] > 0.8)])
得出结论洋葱和马铃薯、汉堡和马铃薯可以搭配着卖，如果客户买了洋葱和汉堡，那么他买马铃薯的概率也很高
2、超市采购实验2
2.1 pandas中的get_dummies函数可以将原始数据转换为one-hot编码，用于数据预处理
2.2 创建数据集
# 创建数据集 retail_shopping_basket = {"ID": [1, 2, 3, 4, 5, 6], "Basket": [["Beer", "Diaper", "Prretzels", "Chips", "Aspirin"], ["Diaper", "Beer", "Chips", "Lotion", "Juice", "BabyFood", "Milk"], ["Soda", "Chips", "Milk"], ["Soup", "Beer", "Diaper", "Milk", "IceCream"], ["Soda", "Coffee", "Milk", "Bread"], ["Beer", "Chips"]]} retail = pd.DataFrame(retail_shopping_basket) retail = retail[["ID", "Basket"]] pd.options.display.max_colwidth = 100 print(retail)
2.3 数据预处理
retail_id = retail.drop("Basket", axis=1) print(retail_id) retail_basket = retail.Basket.str.join(",") print(retail_basket) retail_basket = retail_basket.str.get_dummies(",") print(retail_basket) retail = retail_id.join(retail_basket) print(retail) frequent_itrmsets_2 = apriori(retail.drop("ID", axis=1), use_colnames=True)
2.4 计算支持度
# 计算支持度
frequent_itrmsets_2 = apriori(retail.drop("ID", axis=1), use_colnames=True)
print(frequent_itrmsets_2)
2.5 计算规则
print(association_rules(frequent_itrmsets_2, metric="lift"))
3、电影题材关联实验
3.1 读取数据集
# 电影题材关联实验 # 读取数据集 movies = pd.read_csv("ml-latest-small/movies.csv") print(movies.head(10)) movies_ohe = movies.drop("genres", axis=1).join(movies.genres.str.get_dummies()) print(movies_ohe.head()) print(movies_ohe.shape) movies_ohe.set_index(["movieId", "title"], inplace=True) print(movies_ohe.head())
3.2 计算支持度
# 计算支持度
frequent_itrmsets_movies = apriori(movies_ohe, use_colnames=True, min_support=0.025)
print(frequent_itrmsets_movies)
3.3 计算关联规则
rules_movies = association_rules(frequent_itrmsets_movies, metric="lift", min_threshold=1.25)
print(rules_movies)
print(rules_movies[(rules_movies.lift > 4)].sort_values(by=["lift"], ascending=False))
print(movies[(movies.genres.str.contains("Children")) & (~movies.genres.str.contains("Animation"))])
得到结论：儿童和动物关联性很大
4、手写关联规则计算实验
4.1 如果一个项集是非频繁的，那么它所有的超集也是非频繁的，利用这个可以进行剪枝加快计算
4.2 apriori函数计算过程
4.3 创建数据集
def load_dataset():
"""
加载数据集
:return:数据集
"""
return [[1, 3, 4], [2, 3, 5], [1, 2, 3, 5], [2, 5]]
4.4 创建1项集
def create_c1(data_set): """ 创建1项集 :param data_set:数据集 :return:1项集 """ c1 = [] for transaction in data_set: for item in transaction: if not [item] in c1: c1.append([item]) c1.sort() return list(map(frozenset, c1))
4.5 扫描项集并且筛选
def scan_data(data_set, CK, min_suuport): """ 扫描项集 :param data_set: 项集 :param CK: 标准 :param min_suuport: 最小支持度 :return:大于最小支持度的项集,项集支持度 """ ssCnt = {} for tid in data_set: for can in CK: if can.issubset(tid): if not can in ssCnt: ssCnt[can] = 1 else: ssCnt[can] += 1 num_items = float(len(list(data_set))) ret_list = [] support_data = {} for key in ssCnt: support = ssCnt[key] / num_items if support >= min_suuport: ret_list.insert(0, key) support_data[key] = support return ret_list, support_data
4.6 拼接项集
def apriori_gen(LK, k): """ 拼接项集 :param LK:上一步的项集 :param k:递归计数 :return:拼接好的项集 """ ret_list = [] len_LK = len(LK) for i in range(len_LK): for j in range(i + 1, len_LK): L1, L2 = list(LK[i])[:k - 2], list(LK[j])[:k - 2] if L1 == L2: ret_list.append(LK[i] | LK[j]) return ret_list
4.7 apriori算法主流程
def apriori(data_set, min_support=0.5): """ apriori算法主函数 :param data_set:数据集 :param min_support:最小支持度，默认0.5 :return:所有项集, 项集支持度 """ c1 = create_c1(data_set) L1, support_data = scan_data(data_set, c1, min_support) L = [L1] k = 2 while len(L[k - 2]) > 0: CK = apriori_gen(L[k - 2], k) LK, support = scan_data(data_set, CK, min_support) support_data.update(support) L.append(LK) k += 1 return L, support_data
4.8 生成关联规则函数
def gen_rules(L, support_data, min_conf=0.6): """ 生成规则函数 :param L: 所有项集 :param support_data: 项集支持度 :param min_conf: 阈值 :return: """ rule_list = [] for i in range(1, len(L)): for freq_set in L[i]: H1 = [frozenset([item]) for item in freq_set] rules_from_conseq(freq_set, H1, support_data, rule_list, min_conf)
4.9 计算规则结果
def rules_from_conseq(freq_set, H, support_data, rule_list, min_conf=0.6): """ 计算规则结果 :param freq_set:频繁项集 :param H: 频繁项集中一项 :param support_data: 支持度 :param rule_list: 规则列表 :param min_conf: 最小阈值 :return: """ m = len(H[0]) while len(freq_set) > m: H = cal_conf(freq_set, H, support_data, rule_list, min_conf) if len(H) > 1: H = apriori_gen(H, m + 1) m += 1 else: break
4.10 计算置信度函数
def cal_conf(freq_set, H, support_data, rule_list, min_conf=0.6): """ 计算置信度度的函数 :param freq_set:频繁项集 :param H: 频繁项集中一项 :param support_data: 支持度 :param rule_list: 规则列表 :param min_conf: 最小阈值 :return: 置信度大于阈值的项 """ prunedh = [] for conseq in H: conf = support_data[freq_set] / support_data[freq_set - conseq] if conf >= min_conf: print(freq_set - conseq, "-->", conseq, "conf:", conf) rule_list.append((freq_set - conseq, conseq, conf)) prunedh.append(conseq) return prunedh
4.11 总流程
data_set = load_dataset() L, support_data = apriori(data_set) i = 0 min_conf = 0.5 for freq in L: print("项数", i + 1, ":", freq) i += 1 gen_rules(L, support_data, min_conf)
二十四、词向量算法
1、将数据通过方法转换为向量，就可以把问题转换成计算向量相似度的问题
2、维度越高蕴含的信息也就越多，越准确
3、输入需要转换成向量表示，所以需要建立一张词表，词表初始是随机的，在算法学习中逐渐修改
4、训练数据的来源：随便一篇文章，通过滑动窗口来创建单词对
5、预测一个词知乎的词计算量过大，可以转换为输入两个词预测他们的关系
6、由于样例中都是正常的文章，所以关系都是正例，所以需要添加负样本
7、对比普通的神经网络，词向量算法需要在负反馈传播时候更新词表
二十五、词向量算法实践
1.创建日志目录
current_path = os.getcwd() parser = argparse.ArgumentParser() parser.add_argument( "--log_dir", type=str, default=os.path.join(current_path, "log"), help="The log directory for TensorBoard summaries." ) FLAGS, unparsed = parser.parse_known_args() if not os.path.exists(FLAGS.log_dir): os.makedirs(FLAGS.log_dir)
2.下载数据文本
url = "http://mattmahoney.net/dc/" def maybe_download(filename, expected_bytes): """ 下载数据集文件并且校验长度 :param filename: 文件名称 :param expected_bytes: 文件长度 :return: 文件句柄 """ local_filename = os.path.join(gettempdir(), filename) if not os.path.exists(local_filename): local_filename, _ = urllib.request.urlretrieve(url + filename, local_filename) statinfo = os.stat(local_filename) if statinfo.st_size == expected_bytes: print("下载完成") else: print(statinfo.st_size) raise Exception("无法下载" + local_filename + "请尝试在浏览器下载") return local_filename filename = maybe_download("text8.zip", 31344016) print(filename)
3.读取文本内容
def read_data(filename): """ 将压缩包解压为词表 :param filename:压缩包文件路径 :return: 词表数据 """ with zipfile.ZipFile(filename) as f: data = f.read(f.namelist()[0]).decode(encoding="utf-8").split() return data vocabulary = read_data(filename) print("Data size", len(vocabulary))
4.设置超参数
vocabulary_size = 50000 # 词表大小
feature_length = 100 # 特征长度
C = 3 # 窗口大小
K = 100 # 负样本大小
epoch = 1 # 迭代次数
batch_size = 128 # 取样大小
5.创建数据字典
def build_dataset(words, n_words): """ 创建数据集 :param words:词表 :param n_words: 关注词的数量 :return: """ count = [["UNK", -1]] count.extend(collections.Counter(words).most_common(n_words - 1)) dictionary = dict() for word, _ in count: dictionary[word] = len(dictionary) data = list() unk_count = 0 for word in words: index = dictionary.get(word, 0) if index == 0: unk_count += 1 data.append(index) count[0][1] = unk_count reversed_dictionary = dict(zip(dictionary.values(), dictionary.keys())) return data, count, dictionary, reversed_dictionary data, counts, dictionary, reversed_dictionary = build_dataset(vocabulary, vocabulary_size) del vocabulary print("Most common words(+UNK)", counts[:5]) print("Sample data", data[:10], [reversed_dictionary[i] for i in data[:10]])
6.创建数据集
class WordEmbeddingDataset(Dataset): def __init__(self, data, counts, dictionary, reversed_dictionary): super(WordEmbeddingDataset, self).__init__() self.data = torch.Tensor(data).long() word_counts = np.array([count[1] for count in counts], dtype=np.float32) word_freqs = word_counts / np.sum(word_counts) word_freqs = word_freqs ** (3. / 4.) self.word_freqs = torch.Tensor(word_freqs / np.sum(word_freqs)) self.dictionary = dictionary self.reversed_dictionary = reversed_dictionary def __len__(self): return len(self.data) def __getitem__(self, item): center_word = self.data[item] pos_indics = list(range(item - C, item)) + list(range(item + 1, item + C + 1)) pos_indics = [i % len(self.data) for i in pos_indics] pos_words = self.data[pos_indics] neg_words = torch.multinomial(self.word_freqs, K * pos_words.shape[0], True) return center_word, pos_words, neg_words dataset = WordEmbeddingDataset(data, counts, dictionary, reversed_dictionary) dataloader = torch.utils.data.DataLoader(dataset, batch_size=batch_size, shuffle=True)
7.编写模型
class EmbeddingModel(nn.Module): def __init__(self, vocabulary_size, feature_length): super(EmbeddingModel, self).__init__() self.vocabulary_size = vocabulary_size self.feature_length = feature_length initrange = 0.5 / self.feature_length self.out_embed = nn.Embedding(self.vocabulary_size, self.feature_length, sparse=False) self.out_embed.weight.data.uniform_(-initrange, initrange) self.in_embed = nn.Embedding(self.vocabulary_size, self.feature_length, sparse=False) self.in_embed.weight.data.uniform_(-initrange, initrange) def forward(self, input_labels, pos_labels, neg_labels): input_embedding = self.in_embed(input_labels) pos_embedding = self.out_embed(pos_labels) neg_embedding = self.out_embed(neg_labels) # 计算损失 input_embedding = input_embedding.unsqueeze(2) log_pos = torch.bmm(pos_embedding, input_embedding).squeeze() log_neg = torch.bmm(neg_embedding, input_embedding).squeeze() log_pos = F.logsigmoid(log_pos).sum(1) log_neg = F.logsigmoid(log_neg).sum(1) loss = log_pos + log_neg return -loss def input_embeddings(self): return self.in_embed.weight.data.cpu().numpy() model = EmbeddingModel(vocabulary_size, feature_length) if torch.cuda.is_available(): model = model.cuda()
8.开始训练
optimizer = torch.optim.SGD(model.parameters(), lr=0.2) writer = SummaryWriter("log") for e in range(epoch): for i, (input_labels, pos_labels, neg_labels) in enumerate(tqdm(dataloader, desc=str(e))): input_labels = input_labels.long() pos_labels = pos_labels.long() neg_labels = neg_labels.long() if torch.cuda.is_available(): input_labels = input_labels.cuda() pos_labels = pos_labels.cuda() neg_labels = neg_labels.cuda() optimizer.zero_grad() loss = torch.mean(model(input_labels, pos_labels, neg_labels)) writer.add_scalar("test_loss", loss, e * len(dataloader) + i) loss.backward() optimizer.step() embedding_weights = model.input_embeddings() np.save("feature_length-{}".format(feature_length), embedding_weights) torch.save(model.state_dict(), "embedding-{}.th".format(feature_length)) model.load_state_dict(torch.load("embedding-{}.th".format(feature_length))) embedding_weights = model.input_embeddings()
9.寻找词语相关词汇
def find_nearest(word): index = dictionary[word] embedding = embedding_weights[index] cos_dis = np.array([scipy.spatial.distance.cosine(e, embedding) for e in embedding_weights]) return [reversed_dictionary[i] for i in cos_dis.argsort()[:10]] for word in ["good", "fresh", "monster", "green", "like", "america", "chicago", "work", "computer", "language"]: print(word, find_nearest(word))
可以看出相关词汇都是一些高频词汇，可能是数据集表示很好导致的
或者还可以这样使用
man_idx = dictionary["man"] king_idx = dictionary["king"] woman_idx = dictionary["woman"] embedding = embedding_weights[woman_idx] - embedding_weights[man_idx] + embedding_weights[king_idx] cos_dis = np.array([scipy.spatial.distance.cosine(e, embedding) for e in embedding_weights]) for i in cos_dis.argsort()[:20]: print(reversed_dictionary[i])
找出是女性并且是王，理论上应该有queen女王，但是由于训练数据和参数的原因，这次实验并不是很成功
10.展示一下词向量的空间
meta_list = list(dictionary.keys())
writer.add_embedding(embedding_weights, metadata=meta_list)
writer.close()
