预测 NBA 胜率
我们在 Metis 的第二个项目叫做 Project Luther，它是关于网络抓取和线性回归的。
我是一个超级篮球迷，所以我立即想办法把这项运动和这个项目结合起来。经过几个小时的头脑风暴后，我决定尝试创建一个模型，在赛季开始前预测某个 NBA 球队的胜率。我的目标是尝试采用这种模式，看看是否有人可以利用它通过赌博赚钱。
数据采集和初始步骤
我想使用前几个赛季的球队统计数据和任何其他在赛季开始时可以得到的球队信息。我决定从网站 basketball-reference 获取我的统计数据和信息。这个网站有过去几十年每一季的大量统计数据。我最终获得了大约 40 年的各种统计数据，包括整体球队统计数据、教练信息和球员信息。我用 BeautifulSoup 从网站上搜集了这些信息。
一旦我有了一个熊猫数据框架中的数据，我立即分离出一部分数据用于最终测试。我选择搁置 2018 赛季，因为我想把我对那一年的预测与赌博数据结合起来。然后，我运行一个基线线性模型来获得一个起点。我选择的特征是上一个赛季的平均胜率，因为我觉得这可以很好地预测下一个赛季的成功。
微调我的模型
我的下一步是尝试运行我所有的功能在一起。这稍微增加了我的 R 值，但是我的模型现在有超过 30 个特征。我开始删除那些似乎无法改善我的模型的功能。事实证明，有相当多的功能，我可以删除，而不会降低我的 R 值非常多。
在缩小了我的特征列表后，我尝试添加多项式项，看看是否可以改进我的模型。大多数情况下没有，但有两个特定的多项式要素略微增加了 R。在我的要素工程中，我使用套索交叉验证来测试我的模型，并继续移除不需要的要素。
最终模型
最终，我的最终模型包含了套索正则化，并包含了球员的平均年龄、胜率、得分数、回归球员数和盖帽数等特征。
Metrics for final model
R 没有我希望的那么高，但是考虑到这个问题，这可能是可以预期的。均方根误差乍一看很棒，直到你记得我预测的是一个百分比。RMSE 相当于一个赛季 82 场比赛中的 9 场，这是相当可变的。该模型最重要的特征是胜算和玩家的平均年龄，这并不奇怪。
将模型应用于赌博
我的下一个问题是，这种模式可以用来赌博赢钱吗？我从一个赌博网站刮来了超额/欠额赌注。如果你不熟悉体育博彩，我将简要解释过度/不足投注。维加斯给每支球队在即将到来的赛季中的预期胜场数，然后你打赌这支球队会赢更多的比赛(多)或赢更少的比赛(少)。如果你最终是正确的，你就赢了。
我最初根据赌博信息测试了我对所有球队的预测，只有 50%是正确的，就像掷硬币一样。有趣的是，如果我只关注我的模型预测非常接近拉斯维加斯超额/不足数的八支球队，我的正确率为 87.5%。这看起来很棒，但在相信这个数字不是偶然之前，我想在其他年份进行测试。
模型错在哪里？
由于我的模型不擅长预测，我决定看看是否有具体的原因。下面的残差图表明，我的模型在实际胜率较低时预测过度，在实际胜率较高时预测不足。
我挑选了几个我的模型预测非常错误的例子，并发现有一致的原因。下表显示了球队，年份，我的预测和实际胜率的差异，以及我的预测如此错误的原因。
Examples of when model prediction was very different from actual winning percentage
总的来说，当我的模型严重高估时，球队因为某种原因失去了最好的球员。当我的模型严重低估时，球队在休赛期增加了最好的球员。虽然我的模型考虑了回归球员的数量，但它没有区分超级明星和替补球员。我相信如果我整合了这些信息，我的模型会是一个更好的预测器。
预测自行车共享用户的数量
用监督机器学习解决的经典回归问题
背景
这个项目始于我为机器学习课程( EE 660 )所做的最后一个班级项目。随着我对机器学习(ML)概念和实践的理解变得成熟，我彻底更新了它，并在这个过程中将代码从 Matlab 转移到 R 中。我使用的 R 脚本、数据集和图片都可以在我的 Github 上公开获取。事不宜迟，让我们开始吧。
Photo by Kelly Sikkema on Unsplash
定义问题和项目目标
自行车共享系统 是一项服务，用户可以付费或免费短期租赁/使用可供共享使用的自行车。目前，全世界有超过 500 个自行车共享项目。此类系统通常旨在通过在城市地区提供免费/负担得起的自行车短途出行，而不是机动车，来减少拥堵、噪音和空气污染。对于这样的系统，任何一天的用户数量都可能有很大的变化。预测每小时用户数量的能力可以允许监管这些系统的实体(企业/政府)以更高效和成本有效的方式管理它们。我们的目标是使用和优化机器学习模型，这些模型能够有效地**预测在任何给定的 1 小时时间段内将使用的共乘自行车的数量，**使用关于该时间/天的可用信息。
使用的数据集
我们使用的数据集来自加州大学欧文分校的机器学习库。数据集编制者使用了华盛顿首都自行车共享系统 2011 年和 2012 年两年历史日志中的部分信息。该信息可公开获得。编译器每小时和每天汇总数据(供任何有兴趣深入研究的人使用)。接下来，他们从这里提取并添加所有相应的天气和季节信息。
我们的数据集是一个 csv 文件(可在 my Github 上获得),包含 731 天中 17，379 个小时的信息，每个小时有 16 个特征(信息类别)。其特点是:
- 记录索引
- 日期
- 季节(1:春天，2:夏天，3:秋天，4:冬天)
- 年份(0: 2011，1:2012)
- 月份(1 至 12)
- 小时(0 到 23)
- 节假日:当天是否为节假日
- 工作日:一周中的某一天
- 工作日:如果一天既不是周末也不是节假日，则值为 1。否则为 0
- 天气情况:
— 1:晴，少云，局部多云，局部多云
— 2:雾+多云，雾+碎云，雾+少云，雾
— 3:小雪，小雨+雷雨+散云，小雨+散云
— 4:大雨+冰盘+雷雨+雾，雪+雾 - 以摄氏度为单位的归一化温度。数值分为 41(最大值)
- 归一化的感觉温度，以摄氏度为单位。数值分为 50(最大值)
- 归一化湿度。这些值分为 100(最大值)
- 归一化风速。这些值被分为 67(最大值)
- 临时用户计数
- 注册用户数
- 租赁自行车总数，包括休闲自行车和注册自行车
从一开始看，数据点远远超过特征的数量，这使得这是一个“瘦”数据集，被认为是 ML 的理想选择。
探索性数据分析
在开始用算法处理数据集之前，直观地探索它总是一个好主意。本项目我们将使用。使用 ggplot2 和gg extra包，我们可以快速绘制一些图表来研究可用功能如何影响自行车使用计数。现在我们来看一些图表。
Scatter plot- Adjusted Temperature vs Usage
Scatter plot- Temperature vs Usage
从上面的散点图可以看出，在大部分温度范围内，使用温度和调整温度之间存在正相关关系，线性拟合与最佳拟合曲线相差不远。这在直觉上应该是有意义的，因为人们不太可能在寒冷的天气骑车外出。对于最高温度，这似乎是数据的一个小的子集，在这条曲线上有一个下降。同样，这应该是有意义的，因为当外面太热时，用户也可能不愿意骑自行车。
Scatter plot- Humidity vs Usage
湿度和使用率之间似乎存在负相关，线性拟合非常接近所有数据的最佳曲线拟合(排除一些湿度非常低的异常值)。这可以用 DC 华盛顿州的气候来解释，那里的气候是出了名的潮湿。较高的湿度与较高的降雨机会相关。我们还可以假设，由于高湿度而增加的排汗会阻止人们在户外骑自行车。基于目前为止的想象，我们假设天气状况会影响自行车的使用，降雨会阻碍自行车的使用，这并不是不合理的。我们下面的图表在一定程度上支持了这个假设。从直方图(在 x 轴上)中注意到，晴天(天气等级 1)比阴天或雨天(天气等级 2 或 3)多得多。
Scatter plot- Weather Situation vs Usage
Scatter plot- Wind Speed vs Usage
然而，查看风速数据并不能给我们一个关于它如何影响使用的清晰解释。这两个因素之间的相关性充其量是微弱的。下面是到目前为止讨论的所有连续变量的相关矩阵，它为我们观察到的趋势增加了一些数字。
****注意:有趣的是，临时使用计数与连续变量更相关，并且与我们之前的假设更一致。如果我们想一想，这是有道理的，因为使用自行车上下班的注册用户不太可能受到不舒服的天气条件的影响。我们可以得出这样的结论:分别预测这两个计数，然后将计数相加得到总计数，这样更有意义。然而，当我尝试这样做时，我发现最终的预测不如我们简单地预测总计数得到的准确。因此，对于项目的其余部分，我们将忽略注册的和临时的计数，只将总计数作为输出。如果你愿意，你可以在我的 Github 上访问数据集和我的代码，并尝试包含这些变量，看看你是否能找到更好的结果。继续下一点。让我们来看看时间和日期是如何影响使用的。
从时间的影响来看，使用率最低的似乎是深夜(凌晨 4-5 点之间最低),高峰是上午 8-9 点和下午 5-7 点，这是 DC 的高峰时段。这种拟合远远不是线性的。但是，通过一些简单的数据操作(下一节将详细介绍)，我们可以根据到凌晨 4 点的时间距离来表示使用率，并找到某种程度上的线性拟合(见下文)。
注意:具有线性预测结果的特征是理想的，因为它减少了对复杂的非线性 ML 算法的需要。
在月份与使用情况图(如下)中也可以观察到类似的趋势，在夏季较温暖的月份使用率明显较高，而在一月份使用率最低。通过类似于上一个图的一些操作，该数据也可以用于表示基于到一月份的时间距离的使用情况。然而，这种相关性没有被操纵的时间图那么强。
********
最后，查看“年份”变量(如下)，可以看到使用率从第 1 年到第 2 年上升，这可能表明该系统越来越受欢迎。使用此变量时需要注意的一个重要事项是，尽管它在提供的年份范围(2011-12 年)内进行预测时可能有用，但该算法必须进行大量外推才能预测它对未来日期(2018 年及以后)的影响，这可能会使此变量成为不同时间段的不太可靠的预测器。
预处理:数据清理&特征工程
对于任何 ML 项目，预处理数据都是至关重要的一步。该过程通常是广泛遵循的预处理实践的积累，以及根据模型设计人员的判断对数据进行的特定案例调整。如果数据争论不正确或不充分，它可能导致训练不当的 ML 算法提供不准确或(充其量)次优的结果。古老的谚语“垃圾进来，垃圾出去”在这里很适用。这个过程的一个重要部分是特征工程，它包括将可用的数据特征转化为更有用的变量，帮助我们预测结果。让我们走一遍我们正在采取的步骤:
- 使用先验知识来删除不添加重要信息的特征，在这种情况下，这只是“索引”特征。
- 从日期中提取星期数。在所提供的格式中，日期本身不是我们的算法可以处理的。但是，从这个日期开始，我们可以提取周数(对于特定的年份)并使用该变量作为使用计数的预测值。
- 一个热编码，这是将非二进制分类特征(月、周数、小时、天气状况、季节、工作日)分割成多个二进制子特征的过程，其中每个子特征指示原始特征的某个类别是否为真(1 或 0)。如果可能的话，在不使数据集变得“胖”的情况下，最好将每个多分类特征分成多个二进制特征。由于我们有一个包含 17，000 多个数据点的大规模数据集，我们可以扩展这些功能，并且仍然有一个“瘦”数据集，这样我们就不会有过度拟合的风险。请注意，这确实显著增加了我们的数据集的维数，并增加了我们的计算时间，其系数大于 5X。由于我们在这个项目中只使用非复杂的回归算法，我们可以增加计算的复杂性，以提高预测的准确性。
- 修改循环变量值以表示距单个时间点的时间距离。正如我们在数据可视化部分看到的，对于循环变量(月/周-数字/小时)，试图找到变化值的精确线性模式可能很困难。尽管第 1 周中的数据点在时间上非常接近第 53 周中的数据点，但是与第 53 周相比，第 1 周的最佳线性模型拟合在第 1 周中可能具有非常不同的值。我们通过更改这些值来表示与固定时间点的时间距离来解决这个问题。根据我们在探索性数据分析部分的发现，将基本时间点设置在最少使用时间是合理的。因此，我们使用的时间距离是从凌晨 4 点开始计算小时，从一月中旬开始计算星期和月份。
说明时间序列的随机性
到目前为止，我们一直在处理时间序列数据的确定性方面。对于那些需要复习的人来说，确定性系统是一个在系统未来状态的发展中不涉及随机性的系统。换句话说，确定性模型总是从给定的起始条件或初始状态产生相同的输出。另一方面，如果一个系统是随机的，这个系统的一个或多个部分也具有随机性。我们观察到的大多数时间序列数据通常最好用确定性和随机成分的组合来建模。到目前为止，我们一直在处理我们的模型的确定性方面。现在让我们来看看建模的随机性。
将一个自回归模型整合到我们的特征中是一种解决系统中随机性的简单方法。这是基于这样一个假设，即任何给定小时的用户数取决于之前特定小时数的用户数。这种假设在时间序列数据中通常是有效的。
Ex:Y(t)= B0+B1 * Y(t-1)+B2 * Y(t-2)+…**
然而，为我们的模型选择正确数量的滞后值是很重要的。自相关是指一个时间序列中的观测值相互关联的方式，通过当前观测值(Yt)与当前观测值之前的 p 个观测值(Y(t-p))之间的简单相关性来衡量。R 中的自相关函数(ACF)告诉我们当前值和滞后值之间的自相关，并允许我们决定在模型中包含多少滞后值。
ACF results on the count per hour
不出所料，在一小时内的用户数和之前的 2 个滞后值之间确实存在强正相关，在一小时内的用户数和第 3 个滞后值之间存在中等正相关。因此，我们将 3 个滞后值添加到我们的数据集，作为新特征。
我们最终用 117 个特征从 17，377 个数据点中预测每小时自行车使用量。
应用机器学习算法
对于这个项目，我们将使用两个比较著名的回归算法:最大似然估计和最大后验概率。对于线性回归，我们感兴趣的是找到最佳的参数/权重， w ，使得给定我们的特征，X ，我们的预测结果，Y _ predict = Xw，尽可能接近真实结果， Y_test***
关于构造训练+测试集的注意事项: 在大多数回归问题中，当我们在查看横截面数据时，我们会在选择训练和测试集之前随机化数据点的顺序。然而，当我们处理时间序列数据时，我们必须进行时序检验，因为时间序列排序很重要。我们不能在我们的数据中间切出一块作为测试集，并在这部分之前和之后的数据上进行训练。我们需要在一组比测试数据更老的数据上进行训练。
R 和 Python 中的包允许我们用两三行代码方便地应用 ML 算法。虽然这种方法有其优势，但我认为理解 ML 算法背后的基本统计和概率概念也很重要，以便更好地理解我们的结果。与神经网络等更复杂的模型不同，我们将使用的算法更容易理解和解释。因此，我将非常简要地介绍一下他们每个人正在做什么，如果你愿意，你可以查看我在 Github 上的注释代码，一步一步地了解我们正在做什么。
1)最大似然估计
最大似然估计 (MLE)用于估计概率分布设置中的某个变量。假设我们有一个似然函数， P(D|w) 。这是拥有我们整个数据集的可能性， D ，给定某个符合数据集特征的 w ， X ，到数据集结果， y 。当我们对 w 进行 MLE 时，我们试图推断的参数是:
w _ MLE=argmax wP(D | w)****
换句话说，我们要找到最大化可能性的 w ， P(D|w)。 从这里开始，我们可以做一点线性代数，想出一个我们需要的 成本/损失函数 ，以便计算出最佳权重，然后使用导数最小化这个函数(回忆微积分)来寻找最佳权重。有了几个基本假设，我们可以使用下面公式中总结的一种叫做 普通最小二乘法(OLS) 的方法找到最佳参数/权重:
Ordinary Least Squares formula
使用 OLS，我们在称为训练集的数据子集上找到最佳参数。然后，我们在不同的独立数据子集上测试这些参数，称为测试集**，以查看我们的预测 yMLE 与实际输出 yTest 相比如何。******
我在这里有意加快了数学步骤，因为:a)有免费的在线资源比我更深入地解释了这些方法，b)在本文中，我们更关注这些算法的应用，而不是它们如何工作背后的统计数据。请注意，您需要对微积分、线性代数和概率分布有扎实的理解，才能彻底理解 MLE 或 MAP。
2)最大后验概率
除了 MLE，我们还将尝试另一种叫做最大后验概率(MAP) 的方法。地图由贝叶斯统计得出，贝叶斯统计是一个基于贝叶斯定理的统计领域:****
与基于频率统计的 MLE 不同，MAP 基于这样一种观点，即假设我们对分布有一些有用的先验知识。顾名思义，这适用于后验分布，而不仅仅是可能性。从上面的公式中，我们可以导出由以下公式定义的后验分布:
P(w | D)=(P(D | w) P(w))/P(D)*
假设 P(D) 或者我们的数据集的分布保持不变，我们可以得出结论:
P(w | D)∧P(D | w) P(w)*
第一部分， P(D|w) ，就是我们之前处理过的可能性项。并且我们假设 P(w) 遵循高斯分布，使得:
weights/parameters are assumed to follow a Gaussian distribution
由于我们实际上没有任何关于权重的先验信息，我们不知道权重分布应该遵循什么样的 m_w (均值)或τ平方(方差)。因此，我们使用嵌套的 for 循环来尝试这些术语的值的一千种不同组合，以在我们的验证集上测试我们的算法。验证集最好定义为训练集的子集，用于在最终测试集上测试参数之前微调参数。遵循我们为 MLE 概述的线性代数、微积分和概率步骤，我们发现每个训练集的由以下公式计算:****
在对验证集的参数进行微调后，我们使用测试集特征上的这些参数来预测地图结果，并将这些预测与 yTest 值进行比较。
评估结果和模型的性能
Scatter plot of MLE predictions for y VS real yTest values
Scatter plot of MAP predictions for y VS real yTest values
汇总统计评测:
MLE 运行时间:0.70 秒
中值预测误差:27.05
平均预测误差:56.53
R 平方值:0.65
图 运行时间:258.36 秒
中值预测误差:27.93
平均预测误差:57.85
R 平方值:0.63
一个有用的方法是将我们的模型的性能与天真的预测进行比较。简单预测是一种估计技术，其中上一期的值用作本期的预测，没有任何试图建立因果关系的调整。
朴素预测模型 中位预测误差:37
平均预测误差:69.9
R 平方值:0.65
Examining how MLE, MAP and Naive Forecast predictions vary from the real data using a sub-set of the Test data
讨论
从上面的图和值可以看出，MAP 和 MLE 给出了非常相似的结果，正如我们在没有任何关于我们试图导出的目标函数的先验信息的情况下使用 MAP 时所预期的那样。我们发现我们的预测误差中值约为 27，预测误差平均值约为 57。这明显优于朴素预测模型，朴素预测模型的中值预测误差为 37，平均预测误差为 69.9。请注意，用户计数的总体范围是从 0 到几乎 800。
两个模型和原始预测的 R 平方值约为 0.63-0.65。因此，当看这个指标时，我们的模型并不比天真的预测表现得更好。穿过散点图的最佳拟合线非常接近理想值，即穿过原点的梯度为 1 的线。请注意，我们的预测没有一个小于零，因为我们假设使用计数遵循泊松分布，并且我们使用了泊松回归模型。
进一步研究一下，让我们看看我们的测试数据的一个子集，它只有 100 个数据点。我们看到，当总使用计数较低时，MAP 和 MLE 都可以很好地预测用户计数，但当使用计数超过 450 时，就不那么好了。这些点导致具有大误差的异常值，这解释了为什么我们的平均误差比中值误差大得多。
考虑到我们的中值和均值误差比天真的预测低得多，我们可以得出结论，我们已经创建了预测每小时自行车共享用户数量的有效模型。当然还有进一步改进的空间，尤其是在用户数量较高的时段。因此，如果您对如何改进这些模型有任何建议，请随时告诉我。我们在这里的工作到此结束。感谢您的阅读。
如果你已经做到了这一步，我希望你能像我喜欢写这篇文章一样喜欢读这篇文章。我也希望你比以前更了解探索性的数据分析和机器学习。如果你认为这对你认识的人来说是有趣或有教育意义的，请与他们分享。如果你喜欢这篇文章，有什么想和我分享的，请随时评论或通过电子邮件联系我，地址是 nadir.nibras@gmail.com 或 https://www.linkedin.com/in/nadirnibras/。我致力于改进我的方法、分析或数据集，所以如果你有任何建议，请告诉我。如果你想关注我在数据科学方面的更多工作，请在 Medium 和 Linkedin 上关注我。
如果你是数据科学或其在理解世界中的应用的粉丝，请联系我。和其他统计爱好者聊天总是很有趣，我很乐意在项目上合作:)****
预测美国各县的用药过量死亡率
背景和动机
阿片类药物的流行已经成为这一代美国人的主要公共健康灾难之一。类似于烟草/吸烟或艾滋病毒/艾滋病对前几代人的影响，阿片类药物的流行似乎是这个时代的决定性公共健康危机。我想看看是否有可能建立一个模型来预测一个县一个县的阿片类药物相关死亡率，因为这种类型的模型可能会让我们深入了解在哪里以及如何进行干预。
我怀疑，美国受危机影响最大的地区的相当大的差异可能可以用人口或经济因素来解释，但我实际上更好奇的是，其他更容易修改的预测因素是否会变得重要。例如，常见的说法是，阿片类药物危机始于 20 世纪 90 年代医生开始更自由地开阿片类止痛药，部分原因是制药公司向医生保证这些止痛药成瘾性低，副作用少，但这两者都不是事实。但是阿片类药物处方率仍然是阿片类药物过量死亡的主要驱动因素吗？还是有其他变量是更强的预测因素，例如，海洛因或芬太尼等非法阿片类药物的量可能是更强的预测因素？
设计
该项目的目标是使用多元线性回归模型预测定量目标变量。该项目的设计相当简单:确定一个适当的目标变量，然后找到公开可用的数据集，这些数据集将提供我认为有用和适当的预测信息。然后，我计划使用这些独立变量/预测因子建立一个多元线性回归模型，然后进行特征选择，以查看哪些预测因子是最重要的。最后，我在维持测试数据集上测试了该模型，以查看该模型的性能如何。
数据
在探索了 CDC 对他们公共数据库的描述后，我决定使用 CDC WONDER 数据库来刮取我的目标变量，即 2016 年每个县因药物过量而导致的粗死亡率。关于使用这个作为目标变量的决定的一些评论:我最初想预测每个州的死亡率，但那只会给我 50 个记录。CDC MCD(多种死亡原因)数据库实际上允许您查询具体由于阿片类药物过量导致的死亡(您甚至可以更具体，例如查询由于海洛因、阿片类镇痛剂和其他类型的阿片类药物导致的死亡数量)，但对于个别类别来说，数字变得如此之小，以至于更容易获得由于药物过量导致的总体死亡率的实际数字，而无需指定药物类型。这是因为当每个县的总死亡人数在 10-20 之间时，CDC 会给你死亡人数，但不会计算粗略或经年龄调整的死亡率，因为他们声称这个小数字导致“不可靠”的统计数据。如果总死亡人数在 0-9 之间，那么他们根本不会给出一个数字，而是会声明结果是“被抑制的”，因为在一个小县的少量死亡人数，可能会识别出个别患者。
因此，最容易获得每个县因药物过量而导致的粗死亡率数据，可用的最新数据集是 2016 年。从统计学上来说，比较年龄调整死亡率会更准确，但 CDC 没有提供任何县的年龄调整死亡率，因为它给出了“不可靠”的结果(10-20 岁之间的总死亡人数)，而且没有每个年龄段死亡人数的具体细分，我无法直接计算年龄调整死亡率。看起来确实有间接计算年龄调整死亡率的方法，但我决定继续使用粗略死亡率，然后使用每个县的中位年龄作为预测值，作为一种粗略的方法来说明这一点。这意味着即使美国有大约 3000 个县，因为我最终有 1000 个记录(没错！)实际上可以计算出这些县的粗略死亡率。
Overdose-related crude mortality rate per 100,000 in 2016
就我的自变量而言，大致来说，它们分为 4 类:
人口统计学的
- 白种人百分比
- %高中毕业或以上(指至少高中毕业的人的百分比)
- 年龄中值
经济的
- 家庭收入中位数
- 贫困百分比
- 失业率%
医学的
- 该州的 PDMP 年龄(处方药监控计划)
- 每 100 人阿片类药物处方率
地理学的
- 美国人口普查地区部
人口和经济数据都是使用美国人口普查 API 查询 ACS(美国社区调查)5 年 2016 数据集获得的。确定要查询哪个美国人口普查数据集，然后确定如何构建查询，这无疑是一个学习过程。美国人口普查提供了如此多的信息，看起来似乎有可能获得你想要的精确、具体记录的极其细粒度的查询，但缺点是，一开始弄清楚如何提取甚至简单的数据是相当困难的！无论如何，在查看了我的预测值的配对图后，我意识到我得到的经济数据彼此之间具有相对较高的相关性，例如%贫困和%失业，尽管它们测量的不是完全相同的东西，但它们显然彼此相关。我认为这可能是导致共线性的一个大问题，但回想起来，我认为拥有所有这些变量实际上是没问题的——当我到达结果部分时，我将讨论我认为设计一个结合经济信息的变量可能对构建更好的模型更有用。
每 100 人的阿片类药物处方率来自 CDC，CDC 向表提供了过去十多年来每个县每 100 人的年度阿片类药物处方率。一些数字令人瞠目结舌；在弗吉尼亚州诺顿县，2014 年每 100 人有 563 份阿片类药物处方。没错，那个县的人均阿片类药物处方超过 5 张！许多县的数字没有提供，所以我最后用该州的平均比率来代表这些数字。我加载了 2014 年、2015 年和 2016 年每个县的处方率，当查看相关矩阵时，2014 年的处方率与死亡率具有最强的原始相关性，因此我最终使用 2014 年的阿片类药物处方率而不是其他年份作为我的预测器。由于每年的处方率非常相似，因此没有必要包括所有 3 年。
Opioid prescription rate per 100 people in 2014
处方药监控计划是州级计划(密苏里州除外！)收集和监测关于某些药物处方的信息，特别是阿片类药物和苯二氮卓类药物。不同的州在不同的年份实施了 PDMPs，在 2000 年代，许多州开始实施电子 PDMPs，以便医生、药剂师和其他机构可以访问这些数据。从处方提供者和药剂师的角度来看，PDMP 的主要目标之一是防止患者“逛医生”和去多个提供者那里收集大量阿片类药物或苯并咪唑类药物的处方。如果 PDMPs 真的对减少阿片类药物过量处方有很强的作用，那么希望这可以转化为较低的阿片类药物过量死亡率，因此人们可以假设，总的来说，在有较老的 PDMPs 的州，过量死亡率会较低。
最后，我认为地理区域也可能是粗死亡率的一个有用的预测指标，因为这个国家的一些地区比其他地区受到的打击更大。人们可能会争论如何细化，例如时区、区域、部门，甚至是州与州之间的划分(尽管对于州来说，问题在于州与州之间的大小差异很大)。我最终使用了美国人口普查的地区划分，因为这似乎是宏观/微观之间的一个很好的折衷，因为有 10 个地区。这最终被编码为伪变量，为了避免伪陷阱，我删除了“新英格兰”区域的伪变量。
算法
我们被分配使用的模型是多元线性回归模型。这可以使用 sklearn 或 statsmodels 直接实现，尽管 statsmodels 提供了更多关于模型的信息。考虑到所有的预测因素，我认为 R 平方在 0.3 左右。Damien 然后帮助我使用 sklearn 的多项式特性转换预测因子，以生成二次组合，然后将预测因子的数量从 16 个(8 个定量预测因子，8 个虚拟变量)增加到 53 个，包括定量预测因子的所有交互项。这将整个数据集的 R 平方增加到大约 0.43，这清楚地表明相互作用项确实提供了信息。
下一步是执行特征选择。我最后尝试了两种方法，带交叉验证的后向淘汰法，在每一轮中，每个交叉投票选择一个系数 p 值最高的变量，然后我去掉投票最多的变量。我选择的最终模型是平均验证误差最小的模型。我尝试的另一个策略是使用 sklearn 的 LassoCV 进行特征选择；我再次选择了验证误差最小的套索模型。
有趣的是，LassoCV 和逆向排除法都选择“年龄”作为最重要的预测因子，但之后的四个最重要的预测因子是不同的。此外，逆向排除法选择了几个相当重要的虚拟变量，而 lasso 仍然将它们包括在最终模型中，但它们不那么重要了。他们都给出了相似的 R2 值，大约 0.41 到 0.42，以及大约相同的验证 MSE，0.17 到 0.18。两个模型的最终维持测试数据集的测试 MSE 也大致相同；对于逆向淘汰模型，它是 0.181，对于套索模型，它是 0.184。
工具
为了获得预测值和目标变量数据，我对没有提供 API 的变量使用了 Selenium 和 BeautifulSoup，然后对他们的数据使用了美国人口普查 API。清理数据是在熊猫身上完成的，为了进行分析，我使用了 sklearn 和 statsmodels。
结果
我查看了两个模型中值最大的前 5 个正系数。如下图所示，它们非常相似。反向消除包括组合系数“家庭收入中位数 x 阿片类药物处方率”，而 lasso 包括“%贫困”，但除此之外，它们共有 4 个共同的系数。不出所料，阿片类药物的处方率是前五个积极的预测因素之一，尽管我没有包括任何关于非法药物使用的数据，如海洛因、芬太尼等。我不清楚非法药物使用或处方药是否是过量死亡的更强预测因素。我确实认为组合预测很有趣。家庭收入中值乘以失业率是一个强有力的正预测因子，比任何一个单独的特征都强，这有点奇怪，因为你会认为如果家庭收入中值上升，失业率就会下降，它们会相互抵消。但这可能不是一个真实的假设，因为如果有收入不平等程度高的县，富人变得更富，家庭收入中位数可能会上升(特别是在有收入的家庭)，但失业率也可能上升，因为社会经济光谱相反一端的人也可能同时失业。
我还在完整的数据集上重新调整了向后消除，然后将残差映射到 choropleth 图。我想看看这个模型是否有高估或低估粗死亡率的地理模式。上面的 choropleth 地图是使用 GeoPandas 创建的，geo pandas 是一个为处理空间和地理数据而构建的模块，并带有绘图方法以允许制图。然而，有另一个叫做 follow 的库，它允许你创建交互式地图，并把它们保存为 html 页面，所以我用 follow 创建了一个剩余地图。不幸的是，我不能将互动地图嵌入到一个中等帖子中，但如果你感兴趣，链接是这里。粉红色表示模型高估了该县的用药过量死亡率，而绿色表示模型低估了死亡率；颜色越深，残差越大。
在查看残差数据框架并按残差大小排序后，我可以看到具有最大正残差(意味着模型对真实粗死亡率的低估最严重)的县是马里兰州的巴尔的摩市。负残差最大的县(意味着模型的高估程度最差)是纽约州的萨拉托加县。
如你所见，该模型低估了巴尔的摩的真实过量死亡率 64 人/100K！在萨拉托加县，它被高估了一个较小的数量级，但我认为查看这两个异常值的每个要素的值是很有趣的，因为我认为它提供了一些关于为什么模型对这些县不太适用的见解。首先，我没有意识到美国人口普查局认为巴尔的摩市是自己的县，但它显然是这样的，与该国其他地区相比，巴尔的摩的白人人口比例非常低，请记住，白人比例在向后排除和 lasso 线性回归模型中是第二或第三重要的预测因素。2016 年贫困率为 18.3%，而全国贫困率为 12.7%。最后，巴尔的摩的中位年龄低于全国的中位年龄 38 岁，而且由于年龄是两个模型中最重要的预测因素，这也可能导致模型低估了用药过量的死亡率。
对于萨拉托加县，高中或更高学历的比例显然相当高，达到 94%。同样，失业率%和贫困率%都是全国平均水平；就贫困率而言，你可以看到它比 12.7%的全国贫困率低很多。
我认为，观察这些异常值告诉我，我的模型缺少某种整体经济福祉/繁荣或整体生活质量衡量指标，这可能没有在二次多项式回归的交互项中得到充分体现。巴尔的摩因其高犯罪率而臭名昭著，我认为非法阿片类药物肯定有可能是过量死亡的驱动因素，但该模型无法使用，因为我没有包括任何数据来源。在光谱的另一端，在查找了一些关于萨拉托加县的信息后，它显然是纽约州北部“技术谷”的中心，有大量的国际高科技和电子公司在那里设有校区，这也可以解释它的高教育率和低贫困率和失业率。创建包含更多经济变量(例如教育率 x 贫困率 x 失业率)的交互项可能会为模型带来更大的预测能力。
通过查看残差图，我注意到的另一件事是，该模型似乎倾向于低估城市县(或至少包括主要城市的县)的死亡率，而高估了非城市县的死亡率。我确实考虑过将美国人口普查的 2013 年城市/农村分类纳入我的模型，但人口普查提供的数据并不是哪些县是“城市”和哪些县是“农村”的一对一映射，当我想到这一点时，已经接近这个项目的截止日期，并且需要花费太长时间来弄清楚如何对信息进行编码。
结论
我从这一分析中得出的总体结论是，使用我选择的预测因子来解释用药过量死亡率的一些差异是可能的，但我制作的模型肯定不能解释所有的差异。模型的学习曲线表明，它们处于高偏差状态，我仍然对数据拟合不足，需要更多的特征。同样，我认为找到能提供不同国家非法阿片类药物流动信息的数据源将是有用的，我认为特征工程提出一个更全面的变量来捕捉整体生活质量或经济繁荣将是有用的。此外，该模型在城市和农村的表现之间似乎确实存在差异，因此添加这一功能也可能会增加该模型的预测能力。
我的大多数预测因素是经济或人口数据，但如果我要扩展模型，我认为真正有趣的是使用不同的公共卫生干预措施作为下一年过量死亡率的预测变量，看看这些干预措施的存在与否(或年龄)是否是死亡率的强负预测因素。许多州和县正在尝试不同的干预措施，例如发布纳洛酮命令，使用“活体控制法庭”，开设美沙酮诊所等。我试图通过在我的模型中使用 PDMP 的年龄来做到这一点，但它不是一个强有力的预测器。我认为这方面的一个问题是，我获得信息的 PDMP 网站给出了任何 PDMP 开始的年份，但提供商和药剂师可以实时查询的电子 PDM 直到 2000 年才真正开始，因此 PDM 可能至少在那之前无法影响向公众发布的阿片类药物的实际数量。因此，可能有必要进行更多的研究，为 PDMPs 找到一个更好的功能，而不仅仅是程序的年龄，例如，找出电子 PDMP 首次实施的年份，或者甚至获得一些可以表明每个 PDMP 使用了多少的方法。
最后，另一件有趣的事情是阿片类药物危机是如何随着时间的推移而演变的。我认为随着时间的推移，随着不同的驱动程序接管，不同的功能变得或多或少重要是完全可能的。为了进行这种分析，有必要为每一年建立不同的模型，并计算出一些预测因子的适当时滞；例如，长期纳洛酮命令的存在可能对一个县 1 或 2 年内的死亡率具有强烈的负相关性，但不一定是命令发布的那一年。在这个问题上肯定还有很多工作要做！
预测用户评论的极性
credits - freepick
*“自动化应用于高效的运作将会放大效率”——*比尔·盖茨
根据客户互动和在线评论管理软件 Podium 对 2005 名美国消费者的调查， 93%的消费者表示在线评论对他们的购买决定有影响。
这意味着，对于任何大批量的买卖业务来说，拥有一个强大的系统来对用户评论进行分类是至关重要的，因为评论是一种形式的社会证明。为了在规模上准确地实现这一点(同时考虑时间和成本节约)，添加自动化是有意义的。
在这篇博客文章中，我已经写了我们如何在内部使用哪些机器学习技术来创建这样一个自动化层的第一个版本。
问题
每天，电子商务网站都有数百万的产品列表上传到他们的平台上。这一数据变得更加复杂，因为成百上千的消费者最终会为几乎所有的商品写评论。
但是这个数据是有用的，因为评论被用来创建平均分数来评价产品。
现在，手动对用户评论进行分类既费钱又耗时。这也是为什么在小队，我们正在努力使给用户评论分配极性的过程自动化。目标是预测任何评论的极性，如正面、负面或中性/建议。
自动化情感分析的业务影响
通过包含一个简单的模型，该模型可以根据与评论相关的用户情绪自动对评论进行分类，它可以极大地帮助企业节省通常用于处理内部人工评论团队或外部 BPO 高管团队的成本。事实上，这样的模型实际上帮助团队变得更有生产力，因为现在他们只处理异常，而不是做端到端的猴子工作。
此外，评论的极性可以帮助准确地对产品评级，从而帮助新的在线消费者基于其他用户的喜欢/不喜欢推荐来过滤产品。
从数据集统计开始
我们通过收集各种产品的评论准备了一个用户评论数据集，并使用 Squad 的人在回路中的工作流来标记它们。我们在图 1 中可以看到，正面评价 42.04%，负面评价 45.02%，建议评价 12.94%。数据非常不准确，我们能提供建议的样本非常少。
Figure 1: The distribution of Positive, Negative and Suggestion reviews.
我们实验的不同机器学习管道
这就是机器学习的有趣之处。
通过删除所有重复条目和空白回答对数据进行预处理。然后，对于每个评论，计算 Tf-Idf(术语频率-逆文档频率)特征:
- sublinear_df 设置为 True 以使用对数形式的频率
- min_df 是一个单词必须出现在文档中才能保留的最小比例
- stop_words 设置为*【英语】删除所有常用代词(【a】、【the】*、…)以减少噪声特征的数量
Figure 2: Code snippet to build a Tf-Idf model.
然后，这些数据被分成 70-30 的比例，70%的数据用于训练，30%的数据用于测试我们的模型。有三种目标类别:积极、消极和建议。
我们尝试了两条管道:
**一、管道 1 — Tf-Idf + SVM(支持向量机)
*答:*使用 GridSearchCV 和各种参数进行了交叉验证，如图 3 所示。
B. 在 1.8Ghz 英特尔酷睿 I5 处理器(MacBook Air)和 8GB RAM 上训练这个模型大约需要一个小时。
Figure 3: Hyperparameters to select the best model.
**二世。管道 2 — Tf-Idf + NN(神经网络)
A. 我们做了一个定制的单层神经网络。其架构如下所述。
Figure 4: Code snippet to build a neural network.
B. 模型检查指针用于在每个历元之后保存最佳模型，早期停止用于在验证损失停止改善时停止训练模型。
Figure 5: Callbacks used while training.
C. 在我在 Pipeline 1 中使用的同一台 MacBook Air 机器上训练这个模型花了大约 20 分钟。
比较和结果
使用神经网络的第二种技术给出了相对更好的结果。我们评估了模型在某个置信阈值上的性能，即所有三个类别(负面、正面和建议)中发生的最大概率应该高于该值。置信度小于该阈值的样本被认为是不明确的，然后由 HIL 小队(回路中的人)系统验证。因此，自动化表示不需要人工验证的样本比例。
管道 1 (Tf-Idf + SVM)在 0.9 置信阈值下的结果:
管道 2 (Tf-Idf + NN)在置信度阈值为 0.9 时的结果:
结论
我们开发了一个用户评论分类模型。积极班和消极班的成绩略好于暗示班。从一般的观察中，可以得出这样的结论:暗示样本有时与阳性混淆，有时与阴性混淆。
建议类的性能可以通过进一步增加 Tf-Idf 特征向量的大小来提高；为每个类别制作 Tf-Idf 矢量器，并通过来自每个矢量器的每次检查，连接特征，然后训练分类器。
下一步是什么？
我们目前正在研究短文本增强技术，以增加用户评论数据集中倾斜类的比例，从而提高我们模型的准确性。此外，我们还在进行研究，通过使用基于方面的情感分析等技术来发现用户对产品特定功能的情感。
在我的下一篇文章中，我将写另一个用例，我们目前正在为我们的一个企业客户工作。在此之前，请查看—班参加乔尔考试 。
预测流行的纽约时报文章—质量评分算法
为什么有人会关心《纽约时报》的博客文章有多受欢迎？我认为你应该关心，因为你的领域有可转移的知识——如何监控你的库存质量。例如，我认为脸书的“库存”就是用户的帖子。你可能会喜欢看你最好的朋友的蜜月照片；你的前搭档就没这么多了。向用户提供低质量的库存会让脸书付出与实体店同样高昂的代价。
背景:我参加了在 edX 上举办的一个优秀的 MITx 课程 Analytics Edge，当时我正在研究这个数据集。我是 MOOCs(大规模开放在线课程)的超级粉丝，自 2014 年夏天以来，我一直在 edX 和 Coursera 上注册一门或多门课程。我非常尊重学习过程，并且不会发布代码和数据。
重建这个数据集是非常可能的，只需订阅纽约时报的费用。原始数据包括 2014 年初发表的 6500 篇文章，如果每篇文章收到 25 条以上的评论，这些文章就会被标记为受欢迎。这项任务是预测 2014 年秋季会流行的文章。提供了文章的发表日期、章节、标题和摘要。让我们从一些可视化的描述性分析开始，看看我们是否能做出任何推论。
首先，分离出受欢迎的文章并检查出版时间。
Popular NY Times Blog Articles by The Day of Week
一周内发表的文章似乎比周末发表的文章更受欢迎。这是术语影响的一个很好的例子。“受欢迎”是通过用户评论的数量来衡量的；这并不一定意味着纽约时报在周末的收视率会下降。
让我们更进一步，检查文章发布的时间。
Number of Popular Articles by the Hour of the Day of Posting
中午前后是热门文章发表的高峰，3 小时后的下午 3 点是另一个高峰。这可能是因为西海岸的作者发表文章的时间周期与他们东海岸的同事相似。或许这些文章的发表是为了配合各海岸工人的午休时间？
好的，那么文章发表的部分呢？
Popularity of Articles by Posted Section
一般来说，文章受欢迎的可能性是不受欢迎的 5 倍。然而，我们看到两个突出类别的文章更有可能受欢迎(即。其中绿色条比红色条长)。观点、健康和纵横字谜/游戏部分在用户中的受欢迎程度超出预期。
最后，我们来看一下文章内容。我生成了 6500 篇文章中使用最多的 100 个单词的单词云。
100 Most Popular Words in NY Times Blog Articles
你会注意到一些有政治倾向的词，比如奥巴马，共和党，民主党。直觉上，这是有道理的，因为 2014 年是美国中期选举年
从这四张图表中我们可以推断出什么？如果你写了一篇关于奥巴马医疗保健法的博客，并在美国东部时间周一中午 12 点发表，你的评论区会火起来。
这些都是描述性的，而不是预测性的分析。预测未来几个月哪些文章会流行呢？
好问题。我将在本周晚些时候的后续文章中讨论这个问题。
预测《纽约时报》评论的受欢迎程度(上)
大家好！我刚刚完成了多伦多瑞尔森大学 大数据认证的顶点项目 。这个项目很有挑战性，但是需要大量的学习。我花了近一个月的时间，从头到尾完成了它，并准备了报告。请按照我的 GitHub 上详细的分步 R 码。
第 1 部分将讨论问题陈述，并向您介绍数据集和不同的可视化方法，以便更好地理解数据。第 2 部分将更多地讨论预测建模，使用自然语言处理的文本到向量框架。
介绍
《纽约时报》(NYT)拥有庞大的读者群，在塑造公众舆论和对时事的看法以及设定公共话语基调方面发挥着重要作用，尤其是在美国。NYT 的文章评论区非常活跃，能够洞察读者对文章主题的看法。每条评论都可以收到其他读者的推荐。
NYT 版主面临的挑战
- 高达每篇文章 700 条评论NYT 版主一天内人工审核~ 12000 条评论。
- 版主需要根据预测的相关性和受欢迎程度，更快地做出筛选和分类评论的决定。
- 找到一种更简单的方法来分组相似的评论，并在读者中保持有用的对话。
Number of popular comments per article by News Desk
关键研究课题
- 广泛分析 NYT 的文章和评论
- 读者评论的流行度预测
方法
***** 这个项目的分步代码可以在我的** NYT NLP 顶点 GitHub 资源库中找到。*
步骤 1:数据收集
NYT 文章数据集:
- 数据集由九个部分组成。2017 年 1 月至 5 月和 2018 年 1 月至 4 月发表的文章的 csv 文件(可在网上获得，也可以使用 NYT API 抓取)。
- 共有 9335 篇不同的文章，包含 15 个变量。
NYT 评论数据集:
- 还有另外一组九个。包含对这些文章的评论集合的 csv 文件(可以在网上找到，也可以使用 NYT API 抓取)。
- 总计 2，176，364 条评论，包含 34 个变量。
** 项目中使用的所有数据文件都可以在这里找到。
步骤 2:数据清理和预处理
限制和减少
- 由于庞大的数据量和有限的计算资源，我决定将数据集限制在 14 个可用新闻服务台中的前 6 名。
- 我将某些“角色”特征转换为“因素”
- 我更改了一些特性的数据类型，尤其是将 UNIX 时间戳格式更改为 POSIXct。
- 我删除了分析不需要的特征。
Number of articles published per News Desk
Volume of comments per News Desk
步骤 3:情感倾向得分和计算
组织和清理文本:
- 使用 unnest_tokens 函数处理评论、摘录和文章标题中的文本。
- 在新数据帧的每一行中，文本主体被分割成记号(单个单词)。
- 我去掉了标点符号，把记号转换成小写。
语义取向分数确定:
- 使用基于词典的方法来区分文本取向和单词的语义取向。
- 适当的词典( BING )是从 r。
- 我用 BING 给每个标题、文章片段和文章评论分配了一个情感分数。
- 基于情感得分，每个标题、片段和评论被分为三个情感类别:负面、正面或中性。
步骤 4:提取和分析特征
除了数据中已经存在的变量之外，还衍生了一些其他特征:
- 针对每篇文章标题、摘录和评论正文的合计情感得分 (负面、正面、中性)
- 每条评论的总字数
- 每条评论的总字数
- 平均每句话字数 每个评论
- 时间特征: 添加评论的时间与文章发表日期/时间的时间差
- 星期几 添加注释时
- 添加注释的时间
将建议转换为二进制变量:
此外，为了使用分类方法构建预测模型(如我的下一篇文章中所讨论的)，目标变量“建议”(数字)被转换为可能值为 0 或 1 的二进制变量。
流行与非流行变量是从转换前的推荐变量的五个数字的汇总统计中得出的。其总体中值为 4；因此任何有< =3 张赞成票的评论都被标记为非热门**。**
0 = Non-Popular, 1 = Popular
The final data frame with 26 selected features
步骤 5:数据探索和可视化
在建立任何模型时，我们需要了解预测因素和响应变量之间的相关性。下面的可视化提供了对我的数据更好的理解，也研究了不同变量之间的关系。
Frequency of articles based on sentiment: positive, neutral or negative
Do certain words in article headlines elicit more comments?
Comment popularity across the top six News Desks
Correlation between number of comments and article word counts
Correlation between comment popularity and article sentiments
Correlation between comments and temporal features
Most commonly used words in comments
文本特征与评论流行度的相关性探究
进行 Anova 测试以确定响应变量和数值预测值之间相关性的统计显著性。对于所有三个预测因子，相关性在统计上是T5t，如下所示。
Number of words per sentence in comments vs. comment popularity
Number of sentences in comments vs. comment popularity
Article-comment time gap vs. comment popularity
数据探索摘要
- NYT 最受欢迎的新闻台是专栏**、国内、国外、华盛顿、商业和社论。版主在审核读者添加的评论时，可以专注于这些类别。**
- 在 2017 年和 2018 年，文章倾向于负面情绪多于正面情绪。这可以与美国和世界普遍的政治局势联系起来。
- 对于评论最多的文章以及这些评论中最常用的词，前 25 个术语相似，包括特朗普**、俄罗斯、难民、健康和秘密。**
- 根据评论热度分布，全国性文章是读者阅读和喜欢最多的文章**。这可以归因于美国在 2017-2018 年期间发生的政治变化。**
- 大多数评论是在上午和下午发表的。
- 周二、周四和周五是评论最活跃的日子，而最不活跃的日子是周末。
- 评论的受欢迎程度和一些衍生特征之间有很强的相关性:评论中每句话的平均字数，评论中的平均句子数，以及文章和评论发表之间的平均时间间隔。这意味着这些特征可以用来预测评论流行度。
结论
由于许多预测因素和反应变量之间没有很强的相关性，我选择了一种不同的方法来处理这个问题，我将在下一篇博客中讨论这个问题。
你对四张赞成票成为受欢迎程度的分界点有何看法？你认为提高或降低阈值会如何影响结果？请在评论区告诉我你的想法，并继续关注第 2 部分:-)
(首发@www.datacritics.com)
使用 EA 球员评分和 Tensorflow 预测足球比赛
Photo by Mario Klassen on Unsplash
快速警告
因为我用顶峰体育的收盘赔率作为我的目标，你可以说我只是在建模。然而，假设是到比赛开始的时候，Pinnacle 的线会移动足够多，因为人们打赌赔率应该非常接近真实的概率。有点像群体智慧效应。这是否是一个好的假设…？我要说的是，它不仅效果很好，而且达到了我希望的效果，我将在后面解释。
这个想法
去年我开发了一个预期目标(xG)模型——只是为了好玩和做一些 Python 编程。这在很大程度上受到了两个人的启发，他们是我在推特上关注的迈克尔·凯利和特德·克努特森。在我看来，xG 绝对是梦幻般的。如果你想知道更多的细节，你应该看看迈克尔的推特。然而，和所有的足球分析方法一样，它也有缺点。迈克尔又一次更深入地探讨了这些问题，但我将快速地列出它们，让大家了解为什么我决定做一些不同的事情。首先，它没有考虑参与游戏的玩家。例如，如果埃登·阿扎尔在切尔西缺阵(一个令人悲伤的事实)，你显然会期望切尔西的 xG 会更低，然而没有真正的方法来解释这一点。其次，在赛季开始时，当新球员到来时，xG 无法衡量这个新球员会对球队产生什么样的影响。最后，我还没有看到一个 xG 模型可以很好地适用于所有的五大联盟(这是我还没有尝试用我的模型做的事情，但它是我在某个时候要实现的一个扩展)。
我的 xG 模型也很晚才问世——它已经问世 4、5 年了，所以我的模型不是以前没有人做过的东西。我想尝试创造一些东西，a)可以预测足球比赛的结果，b)我以前从未见过，c)使用某种形式的球员数据。我可能大错特错，但我还没有看到很多人将神经网络应用于足球问题的例子。因此，考虑到以上所有因素，我决定在 Tensorflow 中写点东西，尝试根据一个球队的首发 11 人预测足球比赛的结果。
数据
最明显的玩家评级系统是 EA Sport 的 FIFA 游戏。它们每年都会更新，它们独立于联盟之外，而且它们足够准确(以至于我从来没有听到有人这样谈论它们)。我只使用了整体评分，但可能是其他评分的组合(防守，进攻等。)可能效果更好。我收集了过去 5 年来自 https://www.fifaindex.com/的的英超终极团队评分。我还需要知道哪些玩家开始了哪些游戏。对于这一点，【http://www.betstudy.com/】的是一个伟大的资源。它提供了过去 5 年英超联赛中所有比赛的 11 名首发球员的姓名、号码和国籍。来自 fifaindex 的名字与在 betstudy 上看到的名字并不完全匹配，所以我使用了名字、号码、球队和国籍的组合来匹配球员。一旦匹配完成，我有 22 个评分(11 个主场球员和 11 个客场球员)，在过去 5 年的英超联赛中的每场比赛都在 50 到 100 之间。这是一个好的开始。
What a vector would look like for a Chelsea line-up of Courtois; Moses, Rudiger, Alonso, Christensen, Azpilicueta; Fabregas, Kante; Hazard, Morata, Willian.
然而，仅仅使用这 22 个数字就意味着我们失去了潜在的重要信息。当糟糕的球队对阵优秀的球队时，他们往往会通过启动更多的防守球员来停止比赛——我希望我的模型能够考虑到这种变化。因此，我选择用一个大小为 18 的向量来模拟每个从 11 开始的值。在每个 18 维向量中，第一个分量是对守门员的评分，接下来的六个分量是对后卫的评分——如果只有四名后卫，其中两个分量为零。我在上面加了一张图，试图帮助想象这个 18 维向量在给定切尔西阵容的情况下会是什么样子。根据同样的方法，中场球员有七个组成部分，前锋有四个组成部分。输入的这种结构应该允许从地层中得出一些推论。例如，如果有六名防守球员在比赛，很可能这支球队正在寻找停巴士的机会，因此应该相应地计算赔率。基本上，团队的形成改变了向量的哪些分量被保留为 0，并允许神经网络从中进行推断。此外，当我们将这些数据传递给神经网络时，我们会将它们作为一个大型的 36 维向量进行传递——主队占据前 18 维，客场队占据最后 18 维。通过使用这种结构，网络也可以占据主场优势。
最后，模型的输出将是匹配的 1X2 赔率。我使用了顶尖体育赔率，这些赔率是由 http://www.football-data.co.uk/ T2 的约瑟夫·布赫达尔多年来收集的。
具有训练和验证损失的模型
The neural network architecture. Drawn using http://alexlenail.me/NN-SVG/index.html
收集并格式化所有数据后，必须确定网络的“内部”结构。我最初在 2013-2014、2014-2015、2015-2016 和 2016-2017 英超赛季训练该模型，以便我可以在 2017-2018 赛季回测它。在 1540 个游戏中，我保留了 50 个作为验证集。然后我训练这个模型，使用退出和提前停止使用多种不同的网络结构。实现最小验证误差的是具有两个隐藏层的网络，第一层具有 16 个节点，第二层具有 8 个节点。训练集和验证集的误差图如下所示，垂直线表示提前停止。我不会过多地讨论神经网络的细节，但如果有人对此感兴趣，可以在 GitHub(底部的链接)上找到所有代码。一旦我对模型进行了回溯测试(我在下面展示了回溯测试的一些结果)，我就对包括 2017-2018 赛季在内的所有数据重新训练了网络，为今年做好了准备。
A plot to show the training and validation losses (y-axis) against the number of iterations (x-axis). The vertical line represents the early stopping
回溯测试
为了进行回溯测试，我从一个 100 人的虚拟银行开始，并使用凯利标准来确定股份大小。我没有下赔率大于 3.2 的赌注，也没有下根据模型期望值小于 2%的赌注。在 2017-2018 赛季，该模型实现了 11%的投资回报率。这好得离谱好。我多次重新运行它，多次重新训练该模型，并不断得出相同的结果。所以我也深入研究了这些数字:
- 我赢了 50%的赌注
- 巅峰时期的平均赔率是 2.37
- 平均预测赔率为 2.01
- 赌注的平均价值为 7.23%，最高为 21.3%
这么大的投资回报率肯定是偶然的，但上面的观察似乎很有希望，总的来说，该模型将会盈利。我在我的资金下面绘制了一段时间的情况，以及赌注大小，如果有人特别感兴趣，我可以分享更多关于实际赌注的信息。
A plot to show the bank and stake sizes (y-axis) against the number of bets (x-axis).
2018/19 赛季
鉴于回溯测试的良好结果，我重新培训了包括 2017-2018 赛季在内的模型。使用这个新模型，我可以模拟整个赛季，但有几个注意事项。首先，我必须猜测每支球队最有可能的阵容(这很好，但这些显然会在整个赛季中发生变化)。其次，我不得不使用 2017-2018 年 EA 球员评级，这不是一个大问题，但有像萨拉赫这样的球员，他们的评级肯定会在 FIFA 19 中发生变化-我会用新的评级更新这个模拟，并在 FIFA 19 出来时重新运行模拟。此外，转会窗口仍然开放，这是有问题的，因为它看起来几乎肯定会有一些球员离开-像危险。无论如何，带着这个想法，我运行了这个赛季 100 万次(在我的笔记本电脑上大约需要 8 分钟)，并计算了平均积分、胜败和平局。然后，我计算了 1，000，000 个赛季中球队获得第一名、前 4 名和降级区的百分比。成绩在下面，看起来符合人们的普遍预期。
Predicted Premier League table averaged over 1,000,000 simulated seasons
进一步来说，我可以开始和团队混在一起。举个例子，假设哈扎德离开了，切尔西改打威利安。与哈扎德的 91 分相比，威利安的得分为 84 分。在模拟赛季中交换这一点，切尔西会下降大约 4 分！我们可以对曼联做同样的事情，把德基换成罗梅罗，他们平均下降 3 分。这些积分差异肯定与你的预期有关——没有哈扎德的切尔西是一支更差的球队，没有德基的曼联更差。(我还没有重新运行这个，包括利物浦所有的新签约——但这是我接下来要做的事情)。
总的来说，我对这种工作方式非常满意。这显然不是模拟足球比赛的完美方式。有很多因素需要考虑——教练、天气、疲劳、赛季末和赛季初等等。这个模型不可能捕捉到潜在的分布。然而，这很有趣，事实上像切尔西这样的球队在你拿走他们的明星球员后会丢分，这表明这实际上是按计划进行的。下一步肯定是看利物浦的预测如何变化，当我加入他们的新球员时。然后我想让它设置为自动使用 Smarkets API，并自动下注。最后，当我大学毕业时，我想在过去的 5 个赛季里，在所有 5 个欧洲顶级联赛中训练它，看看我是否能够产生一个“适合所有人的模式”。这也将允许我在世界杯和欧锦赛上使用这个模型——这将非常有趣。
编辑:我做了一个快捷的网络应用，如果你想自己尝试一下的话——https://pl-predictor-tensorflow.herokuapp.com
如果你在推特之外的地方发现了这个-https://twitter.com/BradleyGrantham。
如果有人对实际代码感兴趣，这些都在 GitHub 上-【https://github.com/BradleyGrantham/fifa-stats-scraper
使用神经网络预测概率分布
本帖原载于 Taboola 的工程博客
如果你最近一直在关注我们的科技博客，你可能已经注意到我们正在使用一种特殊类型的神经网络，称为混合密度网络(MDN)。MDN 不仅预测目标的期望值，还预测潜在的概率分布。
这篇博文将关注如何使用 Tensorflow 从头开始实现这样一个模型，包括解释、图表和一个包含全部源代码的 Jupyter 笔记本。
什么是 MDN，为什么它们有用？
现实生活中的数据是嘈杂的。虽然相当恼人，但这种噪音是有意义的，因为它提供了数据来源的更广泛的视角。根据输入的不同，目标值可能有不同程度的噪声，这可能会对我们对数据的理解产生重大影响。
这个用一个例子更好解释。假设以下二次函数:
给定 x 作为输入，我们有一个确定的输出 f(x) 。现在，让我们把这个函数变成一个更有趣(也更真实)的函数:我们将在 f(x) 中加入一些正态分布的噪声。这种噪音会随着 x 的增加而增加。我们将称新函数为 g(x) ，形式上等于g(x)= f(x)+𝜺(x)，其中𝜺 (x) 为正态随机变量。
让我们针对不同的 x 值对 g(x) 进行采样:
紫色线代表无噪声函数 f(x) ，我们可以很容易地看到添加噪声的增加。让我们检查一下 x = 1 和 x = 5 的情况。对于这两个值， f(x) = 4，so 4 是 g(x) 可以取的合理值。4.5 对于 g(x) 也是一个合理的预测吗？答案显然是否定的。虽然 4.5 似乎是 x = 5 的合理值，但当 x = 1 时，我们不能接受它作为 g(x)的有效值。如果我们的模型只是简单地学习预测y’(x)= f(x)，这些有价值的信息就会丢失。这里我们真正需要的是一个能够预测*y’(x)= g(x)*的模型。这正是 MDN 所做的。
MDN 的概念是由 Christopher Bishop 在 1994 年发明的。他的原始论文很好地解释了这个概念，但它可以追溯到史前时代，当时神经网络不像今天这样普遍。因此，它缺少实际实现的部分。这正是我们现在要做的。
我们开始吧
让我们从一个简单的神经网络开始，它只从有噪声的数据集中学习 f(x) 。我们将使用 3 个隐藏的密集层，每个层有 12 个节点，如下所示:
我们将使用均方差作为损失函数。让我们在 Tensorflow 中对此进行编码:
x = tf.placeholder(name='x',shape=(None,1),dtype=tf.float32)
layer = x
for _ in range(3):
layer = tf.layers.dense(inputs=layer, units=12, activation=tf.nn.tanh)
output = tf.layers.dense(inputs=layer, units=1)
训练后，输出如下所示:
我们看到网络成功地学习了 f(x) 。现在缺少的只是对噪声的估计。让我们修改网络来获得这条额外的信息。
走向 MDN
我们将继续使用我们刚刚设计的相同网络，但我们将改变两件事:
- 输出层将有两个节点，而不是一个，我们将这些节点命名为 mu 和 sigma
- 我们将使用不同的损失函数
现在我们的网络看起来像这样:
我们来编码一下:
x = tf.placeholder(name='x',shape=(None,1),dtype=tf.float32)
layer = x
for _ in range(3):
layer = tf.layers.dense(inputs=layer, units=12, activation=tf.nn.tanh)
mu = tf.layers.dense(inputs=layer, units=1)
sigma = tf.layers.dense(inputs=layer, units=1,activation=lambda x: tf.nn.elu(x) + 1)
让我们花一点时间来理解 sigma 的激活函数——记住，根据定义，任何分布的标准差都是非负数。指数线性单位 (ELU)，定义为:
产出 -1 作为其最低值，因此 ELU+1 将总是非负的。一定要去 ELU 吗？不，任何总是产生非负输出的函数都可以，例如 sigma 的绝对值。ELU 似乎做得更好。
接下来，我们需要调整我们的损失函数。让我们试着去理解我们现在到底在寻找什么。我们的新输出层为我们提供了正态分布的参数。这种分布应该能够描述采样产生的数据 g(x) 。我们如何衡量它？例如，我们可以从输出中创建一个正态分布，并最大化从中抽取目标值的概率。从数学上来说，我们希望最大化整个数据集的正态分布的概率密度函数(PDF)的值。同样，我们可以最小化 PDF 的负对数:
我们可以看到损失函数对于𝜇和𝝈.都是可微的您会惊讶于编码是多么容易:
dist = tf.distributions.Normal(loc=mu, scale=sigma)
loss = tf.reduce_mean(-dist.log_prob(y))
仅此而已。在训练模型之后，我们可以看到它确实在 g(x) 上拾取:
在上面的图中，蓝色的线和点代表用于生成数据的实际标准偏差和平均值，而红色的线和点代表由网络预测的不可见的 x 值的相同值。大获成功！
后续步骤
我们已经了解了如何预测简单的正态分布——但不用说，MDN 可以更通用——我们可以预测完全不同的分布，甚至几种不同分布的组合。您需要做的就是确保分布变量的每个参数都有一个输出节点，并验证分布的 PDF 是可微分的。
我只能鼓励你尝试和预测更复杂的分布——使用本文附带的笔记本，修改代码，或者在现实生活数据中尝试你的新技能！愿机会永远对你有利。
用深度学习预测职业棋手的棋步
Source: Pexels
所以我不擅长下棋。
我爸爸在我小的时候教我，但我想他是那种总是让孩子赢的爸爸。为了弥补这个世界上最受欢迎的游戏之一的技能缺乏，我做了任何数据科学爱好者都会做的事情:建立一个人工智能来击败我无法击败的人。不幸的是，它远不如 AlphaZero(甚至不如普通玩家)。但我想看看象棋引擎在没有强化学习的情况下如何运行，以及如何在网络上部署深度学习模型。
在这里对战吧！
获取数据
FICS 有一个数据库，里面有 3 亿个游戏，每个人下的棋，结果，以及参与者的评分。我下载了 2012 年所有至少有一个玩家 2000 ELO 以上的游戏。总共进行了大约 97000 场游戏，移动了 730 万次。获胜分布是:43000 个白色获胜，40000 个黑色获胜，14000 个平局。
极大极小
为了理解如何制作深度学习象棋人工智能，我必须首先理解传统象棋人工智能是如何编程的。进来的是极小极大。Minimax 是“最小化最大损失”的缩写，是博弈论中决定零和博弈应该如何进行的概念。
Minimax 通常与两个玩家一起使用，其中一个玩家是最大化者，另一个玩家是最小化者。代理人，或者任何使用这种算法获胜的人都认为他们是最大化者，而对手是最小化者。该算法还要求有一个董事会评估功能，这是一个谁是赢家的数字措施。这个数字介于-∞和∞之间。最大化者希望最大化这个值，而最小化者希望最小化这个值。这意味着当你，最大化者，面临两个移动之间的十字路口时，你会选择给你更高评价的那个，而最小化者会做相反的事情。这个游戏是在假设两个玩家都处于最佳状态并且没有人犯错的情况下进行的。
Source: GlobalSoftwareSupport
以上面的 GIF 为例。你，最大化者(圆圈)，有三个动作可以选择(从顶部开始)。你应该选择的移动直接取决于你的对手(方块)在移动后选择的移动。但是你的对手选择的移动直接取决于你在移动到移动后选择的移动，以此类推直到游戏结束。一直玩到游戏结束会消耗大量的计算资源和时间，所以你选择一个深度，在上面的例子中，是 2。如果最小化点(最左边的方块)选择向左移动，你有-1 和 1 可以选择。你选择 1 是因为它会给你最高分。如果极小值选择正确的移动，你选择 0，因为那更高。现在轮到最小化者了，他们选择 0，因为那更低。这出戏会一直演下去，直到所有的棋都走完了，或者你没有思考的时间了。对于我的象棋引擎，引擎假设白色是最大值，而黑色是最小值。如果引擎是白色的，算法决定哪个分支将给出最高的最低分，假设人类每次都选择最低分，反之亦然。为了获得更好的性能，该算法还可以与另一种算法结合使用: alpha-beta 剪枝。Alpha-beta 剪枝应用一个截断系统来决定是否应该向下搜索一个分支。
深度学习架构
我的研究始于 Erik Bernhardsson 关于象棋深度学习的伟大文章。他讲述了他如何采用传统方法让人工智能下棋，并将其转化为使用神经网络作为引擎。
第一步是将棋盘转换成输入层的数字形式。我借用了 Erik Bernhardsson 的编码策略，将棋盘与每个方格中的棋子进行热编码。这总共是一个 768 元素的数组(8 x 8 x 12，因为有 12 个元素)。
Bernhardsson 选择将输出层设为 1 表示白棋胜，设为-1 表示黑棋胜，设为 0 表示和棋胜。他认为游戏中的每个棋盘位置都与结果有关。每一个单独的位置都被训练成如果黑棋赢了就“偏向黑棋”，如果白棋赢了就“偏向白棋”。这允许网络返回一个介于-1 和 1 之间的值，这将告诉你该位置更有可能导致白棋获胜，还是黑棋获胜。
我想用一个稍微不同的评估函数来处理这个问题。网络能不能看不出是白棋还是黑棋赢了，但能看到哪一步棋会赢？首先，我尝试将 768 元件电路板表示成输出，其中一个位置是输入，下一个位置是输出。当然，这没有成功，因为这将它变成了一个多分类问题。这为引擎正确选择合法移动提供了太多的错误，因为输出层中的所有 768 个元素都可能是 1 或 0。所以我查阅了 Barak Oshri 和 Nishith Khandwala 的斯坦福论文使用卷积神经网络预测国际象棋中的棋步，看看他们是如何解决这个问题的。他们训练了 7 个神经网络，其中一个网络是棋子选择网络。这个网络决定了哪个方块最有可能被移走。其他六个网络是专门为每一件作品类型，并会决定一个特定的作品应该移动到哪里。如果棋子选择器选择了一个有棋子的方格，只有棋子神经网络会用最有可能移动到的方格作出响应。
我借鉴了他们两人的想法，做了两个卷积神经网络。第一个是从网络移动的*，它将接受 768 元素数组表示，并输出 pro 从哪个方块(在方块 0 和方块 63 之间)移动。第二个网络:移动到网络，将做同样的事情，除了输出层将是 pro 移动到的地方。我没有考虑谁赢了，因为我假设不管最终结果如何，训练数据中的所有移动都是相对最优的。*
我选择的架构是两个 128 卷积层和 2x2 滤波器，后面是两个 1024 神经元全连接层。我没有应用任何池，因为池提供了位置不变性。图片左上角的猫就像图片右下角的猫一样。然而，对于国际象棋来说，王卒的价值与车卒的价值是完全不同的。我对隐藏层的激活函数是 RELU，而我对最后一层应用了 softmax，所以我基本上得到了一个概率分布，所有方块的概率总和为 100%。
我的训练数据是训练集的 600 万个位置，剩余的 130 万个位置用于验证集。训练结束时，我从网络得到的移动的验证准确率为 34.8%，从移动到网络的验证准确率为 27.7%。这并不意味着 70%的时间它没有学会合法的动作，这只是意味着 AI 没有做出与验证数据中的职业球员相同的动作。相比之下，Oshri 和 Khandwala 的网络平均验证准确率为 37%。
将深度学习与极小极大相结合
因为这现在是一个分类问题，其中输出可以是 64 个类中的一个，这留下了很多出错的空间。训练数据(来自高水平玩家的游戏)的一个警告是，优秀的玩家很少能玩到将死。他们知道自己什么时候输了，通常不需要看完整场比赛。这种平衡数据的缺乏使得网络在接近尾声时严重混乱。它会选择要移动的车，并尝试对角移动它。该网络甚至会试图指挥反对派的作品，如果它是失败的(厚脸皮！).
为了解决这个问题，我按照概率对输出进行了排序。然后，我使用 python-chess 库获得了一个给定位置的所有合法移动的列表，并选择了具有最高合成概率的合法移动。最后，我应用了一个预测得分公式，对选择可能性较小的移动进行惩罚:400-(所选移动的指数之和)。合法移动在列表中的位置越靠后，它的预测得分就越低。例如，如果从网络移动的的第一个索引(索引 0)与从移动到网络的第一个索引相结合是合法的，那么预测得分是 400-(0+0)，这是最高可能得分:400。
我选择了 400 作为最大的预测分数，在和材料分数一起玩了这个数字之后。材料分数是一个数字，它会告诉你所做的移动是否会抓住一个棋子。根据捕捉到的棋子，这步棋的总得分会提高。我选择的材料值如下:
卒:10，骑士:500，主教:500，车:900，皇后:5000，国王:50000。
这尤其有助于结束游戏。在将死棋是第二个最有可能合法的棋步并且预测得分较低的情况下，国王的物质价值会超过它。棋子的分数如此之低，是因为网络在游戏初期考虑得很好，如果是战略行动，它会选择棋子。
然后，我将这些分数组合起来，返回给定任何潜在移动的棋盘评估。我通过深度为 3 的最小最大算法(带有 alpha-beta 修剪)输入这个，得到了一个有效的象棋引擎，它被将死了！
使用 Flask 和 Heroku 部署
我使用了 Youtube 上的 Bluefever 软件指南，展示了如何制作一个 javascript 象棋 UI，并通过向 flask 服务器发出 AJAX 请求来路由我的引擎。我使用 Heroku 将 python 脚本部署到 web 上，并将其连接到我的自定义域:【Sayonb.com。
结论
即使引擎的表现不如我希望的那样好，但我学到了很多关于 AI 的基础，将机器学习模型部署到 web，以及为什么 AlphaZero 不单独使用卷积神经网络来玩！
可以通过以下方式进行可能的改进:
- 利用二元模型 LSTM 将从移出的和从移至的网络组合成一个时间序列。这可能有助于将从移出的和从移至的决策放在一起，因为目前每个决策都是独立进行的。**
- 通过增加获取材料的位置来改进材料评估(在棋盘中心获取棋子比在边上获取更有利)
- 在使用神经网络预测分数和材料分数之间切换，而不是在每个节点都使用两者。这可以允许更高的最小最大搜索深度。
- 考虑边缘情况，例如:减少孤立自己棋子的可能性，增加将骑士放在棋盘中心附近的可能性
看看代码，或者在 GitHub repo 用自己的训练数据自己训练一个新的网络！
如果您有任何批评或担忧，您可以通过 LinkedIn 或我的个人网站联系我。
用机器学习预测 Reddit 评论投票数
在本文中，我们将使用 Python 和 scikit-learn 包来预测 Reddit 上一条评论的投票数。我们拟合了各种回归模型，并使用以下指标比较了它们的性能:
- r 来衡量拟合优度
- 在测试集上测量精确度的平均绝对误差(MAE)和均方根误差(RMSE)。
这篇文章基于来自这个 Github 库的工作。代码可以在这个笔记本里找到。
背景
Reddit 是一个流行的社交媒体网站。在这个网站上，用户在各种子主题中发布主题，如下图所示。
A thread in the “AskReddit” subreddit
用户可以对主题或其他评论进行评论。他们还可以给其他线程和评论投赞成票或反对票。
我们的目标是预测评论将获得的支持票数。
数据
该数据是一个 pickle 文件，包含发生在 2015 年 5 月的 1，205，039 行(评论)，存放在 google drive 上，可以使用此链接下载。
下面列出了将用于建模的目标变量和相关特征。它们可以分为几类。
目标变量
- 得分:评论的投票数
评论级别功能
- 镀金:评论上镀金标签(高级赞)的数量
- 区分:页面上的用户类型。“版主”、“管理员”或“用户”
- 争议性:一个布尔值，表示(1)或(0)评论是否有争议(热门评论的支持票数与反对票数接近)
- over_18 :线是否被标记为 NSFW
- time_lapse :评论和线程上第一条评论之间的时间，以秒为单位
- 小时评论:发布了一天中的小时评论
- 工作日:星期几发表评论
- is_flair :评论是否有 flair 文本(https://www . Reddit . com/r/help/comments/3 tbuml/whats _ a _ flair/)
- is_flair_css :是否有注释 flair 的 css 类
- 深度:线程中注释的深度(注释拥有的父注释数)
- no_of_linked_sr: 注释中提到的子条目数
- **链接网址数量:**评论中链接的网址数量
- 主观性:“我”的实例数
- is_edited : 评论是否被编辑
- 被引用:评论是否引用另一个评论
- no_quoted :评论中引用的次数
- senti_neg :负面情绪得分
- senti_neu :中性情绪得分
- senti_pos :正面情绪得分
- senti_comp :复合情绪得分
- 字数:评论的字数
父级功能
- time_since_parent :注释和父注释之间的时间，以秒为单位
- parent_score: 父评论的分数(如果评论没有父，则为 NaN)
- parent_cos_angle :注释与其父注释的嵌入之间的余弦相似度(https://nlp.stanford.edu/projects/glove/)
注释树根功能
- is_root :评论是否为根
- time _ since _ comment _ tree _ root:评论和评论树根之间的时间(秒)
- 评论树根得分:评论树根得分
线程级功能
- link_score :在线评论的 upvotes 开启
- upvote_ratio :帖子评论上的所有投票中，upvotes 的百分比为 on
- link_ups :线程上的投票数
- time_since_link :线程创建以来的时间(秒)
- no_past_comments :发表评论前帖子上的评论数
- score_till_now :该评论发布时的帖子得分
- title_cos_angle :注释与其线程标题嵌入之间的余弦相似度
- is_selftext :线程是否有 selftext
设置
让我们加载我们需要的所有库。
**import** **pandas** **as** **pd
import** **numpy** **as** **np
import** **matplotlib.pyplot** **as** **plt
import** **seaborn** **as** **sns****from** **sklearn.metrics** **import** mean_squared_error, r2_score, mean_absolute_error
**from** **sklearn.model_selection** **import** train_test_split
**from** **sklearn.preprocessing** **import** LabelBinarizer**from** **sklearn.dummy** **import** DummyRegressor
**from** **sklearn.linear_model** **import** LinearRegression
**from** **sklearn.linear_model** **import** LassoCV
**from** **sklearn.linear_model** **import** RidgeCV
**from** **sklearn.linear_model** **import** ElasticNetCV
**from** **sklearn.neighbors** **import** KNeighborsRegressor
**from** **sklearn.tree** **import** DecisionTreeRegressor
**from** **sklearn.ensemble** **import** RandomForestRegressor
**from** **sklearn.ensemble** **import** GradientBoostingRegressor**import** **warnings** warnings.filterwarnings('ignore')
我们还定义了一些与模型交互的函数。
**def** model_diagnostics(model, pr=True):
"""
Returns and prints the R-squared, RMSE and the MAE for a trained model
"""
y_predicted = model.predict(X_test)
r2 = r2_score(y_test, y_predicted)
mse = mean_squared_error(y_test, y_predicted)
mae = mean_absolute_error(y_test, y_predicted)
if pr:
print(f"R-Sq: **{r2:.4}**")
print(f"RMSE: {np.sqrt(mse)}")
print(f"MAE: **{mae}**")
**return** [r2,np.sqrt(mse),mae]**def** plot_residuals(y_test, y_predicted):
""""
Plots the distribution for actual and predicted values of the target variable. Also plots the distribution for the residuals
"""
fig, (ax0, ax1) = plt.subplots(nrows=1, ncols=2, sharey=**True**)
sns.distplot(y_test, ax=ax0, kde = **False**)
ax0.set(xlabel='Test scores')
sns.distplot(y_predicted, ax=ax1, kde = **False**)
ax1.set(xlabel="Predicted scores")
plt.show()
fig, ax2 = plt.subplots()
sns.distplot((y_test-y_predicted), ax = ax2,kde = **False**)
ax2.set(xlabel="Residuals")
plt.show()**def** y_test_vs_y_predicted(y_test,y_predicted):
"""
Produces a scatter plot for the actual and predicted values of the target variable
"""
fig, ax = plt.subplots()
ax.scatter(y_test, y_predicted)
ax.set_xlabel("Test Scores")
ax.set_ylim([-75, 1400])
ax.set_ylabel("Predicted Scores")
plt.show()**def** get_feature_importance(model):
"""
For fitted tree based models, get_feature_importance can be used to get the feature importance as a tidy output
"""
X_non_text = pd.get_dummies(df[cat_cols])
features = numeric_cols + bool_cols + list(X_non_text.columns)
feature_importance = dict(zip(features, model.feature_importances_))
**for** name, importance **in** sorted(feature_importance.items(), key=**lambda** x: x[1], reverse=**True**):
print(f"**{name:<30}**: **{importance:>6.2%}**")
print(f"**\n**Total importance: {sum(feature_importance.values()):.2%}")
**return** feature_importance
读入数据
df = pd.read_pickle('reddit_comments.pkl')
处理缺失值
数据有一些缺失值，通过插补或删除观察值来处理。以下各列中出现缺失值，原因如下:
- parent_score :一些评论没有父评论(估算)
- comment_tree_root_score 和time _ since _ comment _ tree _ root:一些评论是评论树的根(估算)
- parent_cosine，parent_euc，title_cosine，title_euc :部分评论缺少有手套单词嵌入的单词(dropped)。此外，一些评论没有父项(parent_cosine，parent_title imputed)
df = df[~df.title_cosine.isna()] # drop where parent/title_cosine is NaNparent_scrore_impute = df.parent_score.mode()[0] # impute with mode of parent_score column
comment_tree_root_score_impute = df.comment_tree_root_score.mode()[0] # impute with mode of comment_tree_root_score column
time_since_comment_tree_root_impute = df.time_since_comment_tree_root.mode()[0] # impute with mode of time_since_comment_tree_root column
parent_cosine_impute = 0
parent_euc_impute = 0
df.loc[df.parent_score.isna(), 'parent_score'] = parent_scrore_impute
df.loc[df.comment_tree_root_score.isna(), 'comment_tree_root_score'] = comment_tree_root_score_impute
df.loc[df.time_since_comment_tree_root.isna(), 'time_since_comment_tree_root'] = time_since_comment_tree_root_impute
df.loc[df.parent_cosine.isna(), 'parent_cosine'] = parent_cosine_impute
df.loc[df.parent_euc.isna(), 'parent_euc'] = parent_euc_impute
选择变量
在下一步中，我们定义在训练模型时使用哪些变量。我们为布尔变量、多类别变量和数字变量制作了一个列表。
bool_cols = ['over_18', 'is_edited', 'is_quoted', 'is_selftext']cat_cols = ['subreddit', 'distinguished', 'is_flair', 'is_flair_css','hour_of_comment', 'weekday']
numeric_cols = ['gilded', 'controversiality', 'upvote_ratio','time_since_link',
'depth', 'no_of_linked_sr', 'no_of_linked_urls', 'parent_score',
'comment_tree_root_score', 'time_since_comment_tree_root',
'subjectivity', 'senti_neg', 'senti_pos', 'senti_neu',
'senti_comp', 'no_quoted', 'time_since_parent', 'word_counts',
'no_of_past_comments', 'parent_cosine','parent_euc',
'title_cosine', 'title_euc', 'no_quoted','link_score']
使用我们的变量列表，我们可以为建模准备数据。下面的步骤使用 scikit-learn 的LabelBinarizer从分类列中生成虚拟变量，然后组合所有变量。
lb = LabelBinarizer()
cat = [lb.fit_transform(df[col]) **for** col **in** cat_cols]
bol = [df[col].astype('int') **for** col **in** bool_cols]
t = df.loc[:, numeric_cols].values
final = [t] + bol + cat
y = df.score.values
x = np.column_stack(tuple(final))
我们使用 80–20 的比例将数据分成训练集和测试集。
X_train, X_test, y_train, y_test = train_test_split(x, y, test_size=0.2, random_state=10)
建模
在本节中，我们使用 scikit-learn 来拟合 Reddit 数据上的模型。我们从基线模型开始，然后尝试用套索、岭和弹性网回归来改进结果。此外，我们尝试了 K 近邻、决策树、随机森林和梯度推进回归。
首先，让我们定义一个字典来存储模型诊断的结果。
model_performance_dict = dict()
线性回归模型
基线模型
我们用一个简单的模型来建立一个基线。这个模型总是预测平均票数。
baseline = DummyRegressor(strategy='mean')
baseline.fit(X_train,y_train)
model_performance_dict["Baseline"] = model_diagnostics(baseline)
线性回归
linear = LinearRegression()
linear.fit(X_train,y_train)
model_performance_dict["Linear Regression"] = model_diagnostics(linear)
套索回归
lasso = LassoCV(cv=30).fit(X_train, y_train)
model_performance_dict["Lasso Regression"] = model_diagnostics(lasso)
岭回归
ridge = RidgeCV(cv=10).fit(X_train, y_train)
model_performance_dict["Ridge Regression"] = model_diagnostics(ridge)
弹性网回归
elastic_net = ElasticNetCV(cv = 30).fit(X_train, y_train)
model_performance_dict["Elastic Net Regression"] = model_diagnostics(elastic_net)
非线性回归模型
K 近邻回归
knr = KNeighborsRegressor()
knr.fit(X_train, y_train)
model_performance_dict["KNN Regression"] = model_diagnostics(knr)
决策树回归
dt = DecisionTreeRegressor(min_samples_split=45, min_samples_leaf=45, random_state = 10)
dt.fit(X_train, y_train)
model_performance_dict["Decision Tree"] = model_diagnostics(dt)
随机森林回归
rf = RandomForestRegressor(n_jobs=-1, n_estimators=70, min_samples_leaf=10, random_state = 10)
rf.fit(X_train, y_train)
model_performance_dict["Random Forest"] = model_diagnostics(rf)
梯度推进回归
gbr = GradientBoostingRegressor(n_estimators=70, max_depth=5)
gbr.fit(X_train, y_train)
model_performance_dict["Gradient Boosting Regression"] = model_diagnostics(gbr)
模型比较
我们基于三个指标来比较这些模型:R、MAE 和 RMSE。为此，我们定义了下面的函数。
**def** model_comparison(model_performance_dict, sort_by = 'RMSE', metric = 'RMSE'):
Rsq_list = []
RMSE_list = []
MAE_list = []
**for** key **in** model_performance_dict.keys():
Rsq_list.append(model_performance_dict[key][0])
RMSE_list.append(model_performance_dict[key][1])
MAE_list.append(model_performance_dict[key][2])
props = pd.DataFrame([])
props["R-squared"] = Rsq_list
props["RMSE"] = RMSE_list
props["MAE"] = MAE_list
props.index = model_performance_dict.keys()
props = props.sort_values(by = sort_by)
fig, ax = plt.subplots(figsize = (12,6))
ax.bar(props.index, props[metric], color="blue")
plt.title(metric)
plt.xlabel('Model')
plt.xticks(rotation = 45)
plt.ylabel(metric)
让我们使用这个函数来比较基于每个指标的模型。
model_comparison(model_performance_dict, sort_by = 'R-squared', metric = 'R-squared')
model_comparison(model_performance_dict, sort_by = 'R-squared', metric = 'MAE')
model_comparison(model_performance_dict, sort_by = 'R-squared', metric = 'RMSE')
解释结果
考虑到性能和训练时间，随机森林模型是一个合理的选择。平均绝对误差约为 9.7，这意味着平均而言，模型估计误差约为 9.7。让我们看一些图，了解更多关于模型性能的信息。
y_predicted = rf.predict(X_test)
plot_residuals(y_test,y_predicted)
比较测试分数和预测分数的直方图，我们注意到，当目标变量较小时，模型往往会高估目标变量。此外，模型从不预测目标变量会远大于 2000。在目标变量较大的少数情况下，结果似乎有偏差。大多数意见只有少量的支持票，但 model 预计这些意见会得到更多的支持票。然而，当评论有极端数量的支持票时，模型会低估它。
残差的这种分布表明，合乎逻辑的下一步是探索堆叠模型的结果。堆叠是一种集合技术(如随机森林、梯度推进等。)这通常可以提高性能。我们将首先拟合一个分类器来预测向上投票的数量(具有像“少”、“一些”、“许多”这样的类),并且结果将被用作回归模型中的附加预测器。这种方法有可能减少错误并提高拟合度，因为除了我们的原始信息之外，回归模型还会有关于评论数量的提示，以帮助它进行预测。
基于树的模型也允许我们量化他们使用的特性的重要性。
rf_importances = get_feature_importance(rf)
最不重要的特征是不同子数据的指示变量。由于该数据仅包括五个最受欢迎且相当普通的子主题(食品、世界新闻、电影、科学和游戏)的评论，我们不认为这些功能非常重要。此外，还有许多不太重要或不重要的评论。这些特征可以被移除。这有助于避免过度拟合，并减少训练模型所需的时间。
五个最重要的特性是描述评论所在线程或评论父线程的特性。我们可能会预料到这一点，因为受欢迎和有趋势的内容会显示给更多的用户，所以接近有很多支持票的内容的评论也更有可能获得很多支持票。
还需要注意的是，许多非常重要的特性都缺少值。因此，对缺失值的处理方式进行更深入的分析可以提高模型的性能(例如，当我们删除评论树根时，parent score 是最重要的特性，大约占 25%)。使用平均值、中值的插值或使用简单线性回归的预测也值得测试。
结论
在本文中，我们概述了一个使用 scikit-learn python 库预测 Reddit 评论投票的机器学习工作流。我们比较了线性和非线性回归模型的性能，发现随机森林回归是最佳选择。
在快速检查了这个模型的残差之后，我们看到了许多改进的空间。该项目的下一步可能包括:
- 使用较少的特征拟合模型，并将它们的性能与原始模型进行比较
- 分析缺失值及其对模型性能的影响
- 堆叠模型以提高性能
本文基于一个项目，该项目最初由 Adam Reevesman、Gokul Krishna Guruswamy、乐海、Maximillian Alfaro 和 Prakhar Agrawal 在旧金山大学数据科学硕士的机器学习入门课程中完成。相关的工作可以在这个 Github 仓库中找到，本文的代码可以在这个笔记本中找到。
我将很高兴收到以上任何反馈。在 areevesman@gmail.com，你可以通过 LinkedIn 或者电子邮件找到我。
用人口普查收入数据预测学校表现
Source: Pixabay
我小时候，我们家经常搬家。出于这个原因，我的父母一直在考虑找一个合适的社区居住。有几个因素影响了这个决定，但对于忧心忡忡的父母来说，有一个主导因素:一所好学校。
弗雷泽研究所是一家加拿大研究和教育机构，每年发布加拿大学校排名。这些评级是基于安大略省标准化评估 EQAO 的平均表现。
在这个数据科学项目中，我们将根据收入、位置、学校董事会以及学校是小学还是中学来分析多伦多的学校表现数据。在那之后，我们将看看我们是否可以预测 2017 年这些学校的评级。
获取学校评级数据
为了获得安大略省学校的评级，我们可以下载弗雷泽研究所发布的 pdf，然后将相关页面复制并粘贴到一个文本文件中，这样我们就可以使用 python 轻松地逐行扫描。小学和中学以单独的 pdf 格式给出，因此我们最终得到两个文本文件，其中共有 3810 所学校，格式如下:
Ajax [Public] Ajax OSSLT count: 363
ESL (%): 4.7 Special needs (%): 21.8
Actual rating vs predicted based 2016–17 Last 5 Years
on parents’ avg. inc. of $ n/a: n/a Rank: 596/747 n/a
Academic Performance 2013 2014 2015 2016 2017 Trend
Avg. level Gr 9 Math (Acad) 2.7 2.6 n/a 2.8 2.6 n/a
Avg. level Gr 9 Math (Apld) 1.8 1.8 n/a 1.9 1.9 n/a
OSSLT passed (%)-FTE 81.1 80.5 n/a 76.7 75.4 n/a
OSSLT passed (%)-PE 51.3 49.3 n/a 49.0 37.3 n/a
Tests below standard (%) 31.9 34.4 n/a 31.2 34.4 n/a
Gender gap (level)-Math M 0.1 M 0.1 n/a F 0.2 F 0.1 n/a
Gender gap OSSLT F 8.4 F 5.2 n/a F 14.1 F 20.2 n/a
Gr 9 tests not written (%) 0.7 1.5 n/a 0.0 0.7 n/a
Overall rating out of 10 5.5 5.0 n/a 5.4 4.7 n/a
第一行显示学校名称、学校所在的城市以及参加测试的学生人数。与我们最相关的最后一行显示了学校的评分(满分 10 分)。
我们可以将这些信息与来自多伦多数据目录的学校位置数据库联系起来。数据库为我们提供了以下信息:
- 学校名称
- 学校董事会(天主教、公立或私立)
- 位置(地址、纬度、经度和自治市)
- 政府特定标识符
使用 Python，我们可以编写一个脚本来解析这个文本文件，提取评级，并将这些评级作为列添加到多伦多学校数据库中。这样做只是将正则表达式应用到正确的行上。
import pandas as pdfile = 'Elementary Schools 2017.txt'
file2 = 'Secondary Schools 2017.txt'
years = [2013, 2014, 2015, 2016, 2017]
df = pd.read_excel('SCHOOL.xlsx')with open(file) as f:
for line in f:
if '[' in line:
city = re.search('] (.*?) ', line).group(1)
if city == 'Toronto':
school_is_in_toronto = True
school_name = re.search('(.*) \[', line)
school_name = school_name.group(1).upper()
if 'Overall' in line:
if school_is_in_toronto:
ratings = re.search('10 (.*)', line).group(1)
ratings = ratings.split()[0:-1] row = df['NAME'].str.startswith(school_name) if file == 'Elementary Schools 2017.txt':
df.loc[row, 'SCHOOL_TYPE'] = 'ELEMENTARY'
else:
df.loc[row, 'SCHOOL_TYPE'] = 'SECONDARY' for i in range(len(ratings)):
column = str(years[i]) + ' RATING'
df.loc[row, column] = ratings[i]school_is_in_toronto = False
我们首先浏览文件，找到有方括号的行(唯一有方括号的行是有学校名称和城市信息的行)。接下来，我们检查城市是否在多伦多，如果是，我们存储学校的名称。在代码的第二部分，我们检查单词“总体”(只出现在最后一行)并提取评级。使用我们从文本文件中记录的学校名称，我们使用“name”列将其与多伦多学校数据库进行匹配。因为小学和中学存储在两个文本文件中，所以我们也可以添加这个特性。最后，我们添加了 5 列，代表每年的评级。
获取收入数据
告诉我们一个地区的家庭收入中值的分类普查数据可以很容易地从人口普查地图的 API 中下载。但是，普查地点数据是按传播区域 (DA)排序的。要将 DA 转换成地理位置，我们可以下载 StatsCan 提供的文件，它告诉我们加拿大每个传播区域的中心纬度和经度。
在合并了熊猫的两个文件和收集了多伦多的数据后，我们现在有了一个包含 3676 个传播区域的表格，包括它们的经度和纬度，以及它们的家庭收入中值。
income_at_latlongs_df = pd.merge(cnsmpr_df, statscan_df,
how='inner',on=['DA_ID'])
print(income_at_latlongs_df.head())
为了将收入信息与评级信息联系起来，我们可以找到与每所学校距离最小的 DA，并假设它是该邻域的中值收入。
为了实现这一点，我们可以使用库 GeoPy( pip install geopy)来寻找 vincenty 距离(球面上两点之间的距离)。
from geopy.distance import vincentyincomes = []
for school_coordinate in school_coordinates:
minimum_distance = 9999
for i, DA_coordinate in enumerate(DA_coordinates):
distance = vincenty(school_coordinate, DA_coordinate)
if distance < minimum_distance:
minimum_distance = distance
row = income_table.index[i] income_for_school_area = float(income_at_latlongs_df.loc
[row, 'Median_Household_Income'])
incomes.append(income_for_school_area)df['INCOME'] = incomes
将这些收入添加到学校位置数据框架后，我们现在有了完整的数据。
数据可视化
首先，我们可以使用图书馆 gmplot ( pip install gmplot)制作一个散点图，显示每个学校的位置。
该图显示，多伦多中部和东北部的学校评级较高(较亮的点表示评级较高)。通过查阅由人口普查局提供的收入热图，我们可以看到有可能存在相关性。
我们来看看多伦多给出的评分分布。
from statistics import stdev, mean
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(color_codes=True)ratings = df['2017 RATING']sns.distplot(ratings)
plt.title('Distribution of 2017 Ratings in Toronto')
plt.xlabel('School Rating Out of 10')
plt.ylabel('Proportion of Schools')
plt.show()print(mean(ratings), stdev(ratings))
数据大致呈正态分布，均值为 6.35，标准差为 1.73。
既然小学和中学的数据我们都有，那就来看看两种类型的区别吧。
fig = sns.boxplot(x=df['SCHOOL_TYPE'], y=ratings)
plt.title('Distribution of 2017 Ratings by School Type')
plt.xlabel('School Type')
plt.ylabel('School Rating Out of 10')
labels = ['Secondary', 'Elementary']
fig.set_xticklabels(labels)
plt.show()
在这里我们看到，小学的评级通常高于中学。
同样，也可以绘制每个学校董事会评级的箱线图。
不出所料，私立学校排名最高，而公立和天主教学校排名相近。
我们也可以做一个多伦多每个自治市的收视率箱线图。
除了约克市的四分位间距低于其他城市外，其他城市的评级似乎相似。稍后，我们可以测试这种差异在统计上是否显著。
为了恰当地将学校的评级形象化为居民年收入中值的函数，我们应该将收入以 2 万美元为间隔放入垃圾箱。
highest_bin = 200000
num_bins = 10
intervals = highest_bin / num_bins
values = []
column = 'INCOME'
for _, row in df.iterrows():
bin = row[column]/intervals
bin = min(bin, num_bins)
bin = (math.floor(bin))
values.append(bin)df[column+'_BINS'] = values
绘制数据的箱线图得出:
这张图暗示了在一个街区的中值收入和离那个街区最近的学校的评级之间有某种关联。
统计显著性
即使我们在上面的方框图中直观地看到了差异，这也不足以证明这种差异不仅仅是偶然的。我们必须证明这种差异具有统计学意义。为了找到统计意义，我们可以进行单向方差分析。通过方差分析，我们陈述了零假设，即类别中的样本是从具有相同均值的总体中抽取的。例如，对于学校董事会列，我们假设 _catholic = _public = _private。另一个假设是，样本是从拥有不同均值的人群中抽取的，因此会受到他们所属群体的影响。
from scipy import statscolumns_to_analyze = ['MUN', 'SCHOOL_BOARD', 'SCHOOL_TYPE']
for column in columns_to_analyze:
grouped_dfs = []
for group in df.groupby(column).groups:
grouped_df = df.groupby(column).get_group(group)
grouped_df = grouped_df.reset_index()['2017 RATING']
grouped_dfs.append(list(grouped_df.dropna()))
F, p = stats.f_oneway(*grouped_dfs)
print(f'{column}: {p: 0.2e}')MUN: 3.70e-05
SCHOOL_BOARD: 9.25e-07
SCHOOL_TYPE: 1.33e-05
每一列的方差分析给出了远低于 0.05 的 p 值，因此我们可以说差异具有统计学意义。我们拒绝每列的零假设，即各组的总体均值相同。
特征工程
我们可以使用回归技术来看看学校的评级是否可以使用我们创建的特征来预测。
但是，在我们将特征放入模型之前，必须填充缺失值，并且所有特征必须被热编码。为了处理缺失值，我们用已收集的评分平均值替换任何缺失值。我们应该确保平均值中不包括 2017 年的评级，因为这是我们试图预测的值。我们可以在遍历前面的文本文件时填充缺失的值。
def fill_missing_values(ratings_with_na):ratings_with_na = [np.nan if rating == 'n/a'
else float(rating)
for rating in ratings_with_na]if all([rating == 'n/a' for rating in ratings_with_na[0:-1]]):
return ratings_with_namean_rating = round(np.nanmean(ratings_with_na[0:-1]), 1)filled_missing_values = []
for rating in ratings_with_na:
if rating in [np.nan]:
filled_missing_values.append(mean_rating)
else:
filled_missing_values.append(rating)return filled_missing_values
该函数采用一个包含 5 个评级的列表，并将所有“不适用”的评级转换为 np.nan，因此我们可以应用 np.nanmean()来计算平均值。然后，我们检查前 4 个评分是否缺失，如果缺失，我们将返回包含缺失值的列表(稍后我们将删除这些行，因为它们非常少)。最后，我们找到 2013-2016 年间的平均评级，并返回一个没有缺失值的新评级列表。
为了对每个分类值进行热编码，我们可以使用 pd.get_dummies()，它将每个分类值转换成多个列。根据数据点是否属于该类别，每列的值可能为 0 或 1。
ordinal_columns = df[['INCOME_BINS', '2013 RATING', '2014 RATING',
'2015 RATING', '2016 RATING', '2017 RATING']]
categorical_columns = ['MUN', 'SCHOOL_TYPE', 'SCHOOL_BOARD']
categorical_columns = pd.get_dummies(df[categorical_columns])
ohe_df = pd.concat([categorical_columns, ordinal_columns], axis=1)
包含重命名列的新数据框架如下所示。
既然所有缺失值都已填充，分类列也已热编码，现在让我们绘制皮尔逊相关矩阵，找出我们应该使用的特征。
from matplotlib.colors import LinearSegmentedColormapcorr_sorted = abs(ohe_df.corr()['2017 Rating']).sort_values()
sorted_df = ohe_df[list(corr_sorted.index)]
corr = round(sorted_df.corr(), 2)min_color = 'white'
max_color = (0.03137254, 0.18823529411, 0.41960784313, 1)
cmap = LinearSegmentedColormap.from_list("", [max_color,
min_color,
max_color])
fig = sns.heatmap(corr, annot=True, cmap=cmap,
xticklabels=corr.columns.values,
yticklabels=corr.columns.values,
cbar=False)
plt.xticks(rotation=0)
fig.xaxis.set_tick_params(labelsize=8)
fig.yaxis.set_tick_params(labelsize=8)plt.show()
对相关矩阵进行排序，颜色相对于 0 对称绘制，因此我们可以看到 0.19 的相关性与-0.19 的相关性一样显著。正如假设的那样，每年的学校评级与其他年份的评级密切相关，而收入和私立学校的地位与评级有一定的相关性。
预言；预测；预告
首先，我们必须建立指标来评估我们的模型是否有效。在测量回归性能时，我们将使用平均绝对误差和均方根误差。为了衡量我们的机器学习模型是否值得努力，我们将创建一个基线。例如，如果我们每次都预测平均值，并获得比任何机器学习模型更低的均方误差/平均绝对误差，那么我们的模型实际上是无用的。如果我们每次只预测平均评分(6.35)，以下是衡量标准。
y_predicted = mean(df_y)
rmse = (mean((y_predicted - df_y)**2))**0.5
mae = mean(abs(y_predicted - df_y))
print(f'RMSE: {rmse: 0.4f}, MAE: {mae: 0.4f}')RMSE: 1.7311, MAE: 1.3864
如果我们对所有情况的预测值都是 6.35，那么平均来说，预测值将与真实值相差约 1.39 个点。
现在，我们可以设置模型来预测 2017 年给出的评分！为了了解在没有往年信息的情况下能否预测 2017 年的评级，我们将使用学校类型、学校董事会、自治市和收入特征。下面是随机森林、梯度推进、线性回归和留一交叉验证的支持向量机的设置代码。
from sklearn.model_selection import cross_validate, LeaveOneOut
from sklearn.ensemble import RandomForestRegressor,
GradientBoostingRegressor
from sklearn.linear_model import LinearRegression
from sklearn.svm import SVRdef LOOCV(ohe_df, use_previous_years=False):
models = {
'RFR': RandomForestRegressor(n_estimators=50, random_state=0),
'GBR': GradientBoostingRegressor(max_depth=1, random_state=0),
'LIR': LinearRegression(),
'SVR': SVR(kernel='linear')
}
if use_previous_years is False:
ordinal_columns = df[['INCOME_BINS', '2017 RATING']]
ohe_df = pd.concat([categorical_columns,
ordinal_columns], axis=1) df_x = ohe_df.iloc[:, :-1]
df_y = ohe_df.iloc[:, -1] scoring = ['neg_mean_squared_error', 'neg_mean_absolute_error']
loo = LeaveOneOut.get_n_splits(df_x, df_y) for name, model in models.items():
scores = cross_validate(model, df_x, df_y, cv=loo,
scoring=scoring)
rmse = (-1*mean(scores['test_neg_mean_squared_error']))**0.5
mae = -1*mean(scores['test_neg_mean_absolute_error'])
print(f'{name} RMSE: {rmse: 0.4f}, MAE: {mae: 0.4f}')LOOCV(ohe_df)RFR RMSE: 1.6488, MAE: 1.2841
GBR RMSE: 1.5320, MAE: 1.2210
LIR RMSE: 1.5255, MAE: 1.2210
SVR RMSE: 1.5447, MAE: 1.2354
我们可以看到，线性回归的平均绝对误差为 1.22，RMSE 为 1.5255。这是不坏的表现，因为这意味着我们的预测平均在真实值的 1.22 点以内，这比我们每次预测平均值要好 0.16 点或 12%。
为了获得更多的信息，我们来看看与往年的信息有哪些误差。
LOOCV(models, df_x, df_y, use_previous_years=True)RFR RMSE: 1.0381, MAE: 0.7906
GBR RMSE: 0.9796, MAE: 0.7443
LIR RMSE: 0.9513, MAE: 0.7215
SVR RMSE: 0.9520, MAE: 0.7200
SVM 和线性回归这次效果最好，平均绝对误差为 0.72，比我们每次预测平均值的结果好 0.76 点或 48%！
结论
在本文中，我们汇总了来自多伦多数据目录、StatsCan 和弗雷泽研究所的数据。在数据被转换成可用的形式后，它被用来可视化一所学校的评级(满分为 10 分)与各种特征(如学校董事会和中等邻里收入)之间的相关性。发现有相关性，有了这些信息就有可能预测一个学校的评级平均 1.22 分。
数据从质量上说的是众所周知的:贫困可以预测许多特征，而教育是一个非常重要的特征。
接下来的步骤可能包括:
- 将其他城市纳入数据
- 利用学校雇用的教师的信息
- 获取除收入以外的家庭人口统计数据，看看它如何影响学校的表现
文章中使用的所有代码都可以在我的 GitHub 上找到。如果你对这篇文章有任何反馈或批评，你可以通过 LinkedIn 联系我。
预测亚马逊产品评论的评分
分类器算法
介绍
在机器学习中，分类用于根据包含类别事先已知的观察值的训练数据集，将新观察值分类到特定的集合/类别中。最常见的例子是电子邮件的“垃圾邮件”或“非垃圾邮件”类别。在电子商务中，可以使用分类器算法来基于单词对评论的情感进行分类。语言中的特定单词根据它们的积极或消极情绪被预先分类。
分类是监督学习的一个实例。训练集已正确识别观察值。分类器算法用于基于与训练数据集的相似性和/或距离从未分类的无监督数据中创建聚类/集合。
分析背后的统计数据
在简单分类器模型中，正和负数据点的简单计数将定义总体正或负集合。这就有问题了。例如，在句子中的单词的情况下，“Great”和“Good”都是肯定词。但是“很棒”比“好”的影响更大。我们需要训练我们的模型来衡量确定的数据点。
分类器模型边界可以是将“阳性”和“阴性”结果分开的简单线，也可以是将多个组分开的更复杂的超平面。
可以使用模型的误差或准确度以及“假阳性”和“假阴性”来观察模型的有效性。
使用的工具
- 计算机编程语言
- GraphLab
- s 帧(类似于熊猫数据帧)
数据加载
我们准备对亚马逊产品评论进行分类，了解正面或负面评论。亚马逊有 5 星评级。我们将用它来比较我们的预测。
*> product_reviews=graphlab.SFrame(‘amazon_baby.gl/’)**> product_reviews.head()*
*> product_reviews.show()*
创建一个单词字典，将使用我的模型进行分类。
*> product_reviews[‘wordcount’] = graphlab.text_analytics.count_words(product_reviews[‘review’])*
选择一个特定的产品来预测评论的情绪。Vulli Shopie 是一款适合婴儿长牙的长颈鹿玩具。过滤产品的所有评论。
*> vs_reviews = product_reviews[product_reviews[‘name’]==’Vulli Sophie the Giraffe Teether’]**> vs_reviews[‘rating’].show(view=’Categorical’)*
使用分类器模型的数据分析
过滤掉中间等级为 3 的。高于 4 的评级意味着它是正面评级。
*> product_reviews = product_reviews[product_reviews[‘rating’] !=3]**> product_reviews[‘binrating’] = product_reviews[‘rating’] >= 4*
创建培训和测试数据集。
*> train_data, test_data = product_reviews.random_split(0.8, seed=0)*
构建一个情感分类器来识别评论是正面的还是负面的情感。逻辑分类器模型将使用来自训练数据的单词(单词计数列)和评级(binrating)来开发预测目标(binrating)的模型。
*> sentiment_model = graphlab.logistic_classifier.create (train_data,* *target=’binrating’,* *features=[‘wordcount’],* *validation_set=test_data)*
评估测试数据的模型，并使用 ROC(受试者工作特性)曲线显示。
*> sentiment_model.evaluate(test_data, metric=’roc_curve’)*
结果
*> sentiment_model.show(view=’Evaluation’)*
91%的准确率就不错了。
使用我们建立的模型，预测 Vullie Sophie 长颈鹿玩具的每个评论情绪。
*> vs_reviews[‘predicted_sentiment_by_model’]=
sentiment_model.predict(vs_reviews, output_type=’probability’)*
我们的模式运行得非常好。对于更高数量的情感(接近 1)，我们可以观察到亚马逊产品星级为 5。
我们可以根据模型预测的情绪来查看最积极和最消极的评论。
*> vs_reviews=vs_reviews.sort(‘predicted_sentiment_by_model’, ascending=False)**> vs_reviews[0][‘review’]*
苏菲，哦，苏菲，你的时候到了。我的孙女，维奥莱特 5 个月大，开始长牙了。小索菲给瓦奥莱特带来了多大的快乐。索菲是由非常柔韧的橡胶制成的，结实但不坚韧。对瓦奥莱特来说，把索菲扭成前所未闻的姿势，让她把索菲送进嘴里是相当容易的。小鼻子和小蹄子完全适合小嘴巴，流口水是有目的的。索菲身上的油漆是食品级的。苏菲 1961 年出生于法国。制造商想知道为什么没有婴儿用品，于是在圣索菲节用最好的不含邻苯二甲酸盐的橡胶制作了索菲，这个名字就这样诞生了。从那时起，数以百万计的苏菲出现在这个世界上。她很柔软，婴儿的小手很容易抓住。紫罗兰特别喜欢苏菲那凹凸不平的头和角。苏菲的脖子很长，很容易抓住和扭动。她有可爱的，相当大的斑点，吸引了紫罗兰的注意。索菲有快乐的小吱吱声，带来了紫高兴的尖叫。她能让苏菲吱吱叫，这带来了很多快乐。索菲光滑的皮肤抚慰着维奥莱特的小牙龈。索菲有 7 英寸高，是婴儿抱和爱的正确尺寸。众所周知，婴儿抓到的第一件东西是放进嘴里的——拥有一个能刺激所有感官并有助于解决出牙问题的玩具是多么美妙。苏菲很小，可以放进任何尺寸的口袋或包里。索菲是从几个月到一岁婴儿的最佳选择。听到发现索菲不可抗拒的婴儿发出的咯咯笑声是多么美妙啊。苏菲万岁！强烈推荐。pris rob 12–11–09 "
最负面的评论。
*> vs_reviews[-1][‘review’] # most negative*
“我的儿子(现在 2.5 岁)喜欢他的索菲，我每次去参加婴儿送礼会都会买一个。现在，我的女儿(6 个月)就在今天差点被它呛到，我再也不会给她吃了。如果我不在听得见的范围内，那可能是致命的。她发出的奇怪声音引起了我的注意，当我走到她面前，发现她弯曲的前腿被塞进了她的喉咙，她的脸变成了紫色/蓝色时，我慌了。我把它拉出来，她在尖叫前吐了整个地毯。真不敢相信我对这个玩具的看法怎么就从必须要变成不能用了。请不要忽视任何关于窒息危险的评论，它们并没有过分夸张！”
用 Python 中的随机森林分类器预测恒星、星系和类星体
来自斯隆数字巡天数据集的见解
The SDSS website banner!
最近，在我寻找有趣(物理学相关)数据集的过程中，我偶然发现了 Kaggle 上的斯隆数字巡天(SDSS)数据集，以及 Faraz Rahman 的辉煌 Kaggle 内核，它使用 R 中的支持向量机预测了不同类型的天文物体(恒星、星系和类星体)。然而，由于 R 带回了一些可怕的记忆和训练 SVM 需要大量的计算工作，我决定使用 sci 来尝试一下这个分类问题
数据概述
标签
那么恒星、星系和类星体到底是什么？如果你在开始这个项目之前问我，我可能无法回答(真为我感到羞耻)。幸运的是，法拉兹的笔记本简明扼要地总结了它们是什么:
- 一个 星系 是一个由恒星、恒星残骸、星际气体、尘埃和暗物质组成的引力束缚系统。星系根据它们的视觉形态分为椭圆形、螺旋形和不规则形。许多星系被认为在其活动中心有超大质量黑洞。
- 恒星 是一种天文物体，由一个发光的等离子球体组成，通过自身重力聚集在一起。离地球最近的恒星是太阳。
- 一颗 类星体 又称准恒星天体，是一颗极其明亮的活动星系核(AGN)。类星体辐射的能量是巨大的。最强大的类星体的亮度超过 1041 瓦，是普通大星系如银河系的数千倍。
安德鲁·邦克教授关于恒星和星系(S26) 的短选项课程的页面也是一个很好的资源。
特征
细节的细节可以在 Kaggle 数据集概述上找到。更重要的功能总结如下:
- ra，dec —分别为赤经和赤纬
- u，g，r，I，z-过滤波段(又称光度系统或天文星等)
- 运行、重新运行、camcol、字段—图像中字段的描述符(即 2048 x 1489 像素)
- 红移——由于天文物体的运动，波长增加
- 车牌——车牌号码
- mjd —修正儒略历观测日期
- 光纤 id —光纤 id
探索性分析
这里的分析遵循了法拉兹的分析。我会让观想自己说话。
加载库和数据
*import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inlinedf = pd.read_csv("skyserver.csv")*
数据描述
自然，我从df.head()、df.describe()和df.info()开始。df.info()的输出如下所示:
*<class 'pandas.core.frame.DataFrame'>
RangeIndex: 10000 entries, 0 to 9999
Data columns (total 18 columns):
objid 10000 non-null float64
ra 10000 non-null float64
dec 10000 non-null float64
u 10000 non-null float64
g 10000 non-null float64
r 10000 non-null float64
i 10000 non-null float64
z 10000 non-null float64
run 10000 non-null int64
rerun 10000 non-null int64
camcol 10000 non-null int64
field 10000 non-null int64
specobjid 10000 non-null float64
class 10000 non-null object
redshift 10000 non-null float64
plate 10000 non-null int64
mjd 10000 non-null int64
fiberid 10000 non-null int64
dtypes: float64(10), int64(7), object(1)
memory usage: 1.4+ MB*
没有一个条目是 NaN，这是维护良好的数据集所期望的。清洁不是必须的。
唯一条目
nunique()方法返回包含每列唯一条目数量的 Series 对象。
*df.nunique().to_frame().transpose()*
每个天文实体的出现
然后我在班级栏上运行了value_counts()。
*occurrences = df['class'].value_counts().to_frame().rename(index=str, columns={'class': 'Occurrences'})occurrences*
我们看到大多数条目不是星系就是恒星。只有 8.5%的条目被归类为类星体。
密度分布图
使用核密度估计(kde)，我绘制了各种特征的(平滑)密度分布。
*featuredf = df.drop(['class','objid'], axis=1)
featurecols = list(featuredf)
astrObjs = df['class'].unique()colours = ['indigo', '#FF69B4', 'cyan']plt.figure(figsize=(15,10))
for i in range(len(featurecols)):
plt.subplot(4, 4, i+1)
for j in range(len(astrObjs)):
sns.distplot(df[df['class']==astrObjs[j]][featurecols[i]], hist = False, kde = True, color = colours[j], kde_kws = {'shade': True, 'linewidth': 3}, label = astrObjs[j])
plt.legend()
plt.title('Density Plot')
plt.xlabel(featurecols[i])
plt.ylabel('Density')
plt.tight_layout()*
还为每个类别绘制了滤波器频带密度。
*filterbands = pd.concat([df.iloc[:,3:8], df['class']],axis=1)plt.figure(figsize=(15,5))
plt.suptitle('Density Plots')
sns.set_style("ticks")
for i in range(len(astrObjs)):
plt.subplot(1, 3, i+1)
for j in range(len(featurecols2)):
sns.distplot(df[df['class']==astrObjs[i]][featurecols2[j]], hist = False, kde = True, kde_kws = {'shade': True, 'linewidth': 3}, label = featurecols2[j])
plt.legend()
plt.xlabel(astrObjs[i])
plt.ylabel('Density')
plt.tight_layout()*
附加可视化
为了完整起见，我包括一个 3D 图，与原始笔记本相同。最初的意图似乎是确定 SVM 的线性核是否有效(如果我错了，请纠正我)。底部有很多群集，我取了红移的对数(忽略误差)以使可视化更清晰。
*from mpl_toolkits.mplot3d import Axes3Dfig = plt.figure(figsize=(5,5))
ax = Axes3D(fig)for obj in astrObjs:
luminous = df[df['class'] == obj]
ax.scatter(luminous['ra'], luminous['dec'], np.log10(luminous['redshift']))ax.set_xlabel('ra')
ax.set_ylabel('dec')
ax.set_zlabel('log redshift')ax.view_init(elev = 0, azim=45)plt.show()*
构建随机森林分类器
训练集和测试集分离
传统的列车测试分割可以通过以下方式完成:
*from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifierx_train, x_test, y_train, y_test = train_test_split(features, labels, test_size=0.3, random_state=123, stratify=labels)clf = RandomForestClassifier()*
训练复杂性
当我最初试图用sklearn.svm.linearSVC训练我的数据时，我的笔记本电脑开始严重过热。训练时间复杂度一般在 O(mn)到 O(mn)之间，其中 m 是特征数，n 是观测数，如 Jessica Mick 此处所解释。另一方面，growing CART(分类和回归树)的训练复杂度为 O(mn logn)和 O(mn)，这里的解释为(随机森林是 CART 的系综)。由于时间、耐心和手头的咖啡有限，我决定换成随机森林模型。
事后看来，我可以做的一件事是将我的数据缩放到[-1，1]来加速 SVM(甚至随机森林)，正如 Shelby Matlock 在同一篇文章中提到的。这样我也能得到更稳定的预测结果。
*from sklearn.preprocessing import MinMaxScaler
scaling = MinMaxScaler(feature_range=(-1,1)).fit(x_train)
x_train_scaled = scaling.transform(x_train)
x_test_scaled = scaling.transform(x_test)*
超参数优化
对于超参数调谐，我发现这个和这个相当方便。我们首先实例化一个随机森林，并查看可用超参数的默认值。美化打印get_params()方法:
*from pprint import pprint
pprint(clf.get_params())*
这给出了:
*{'bootstrap': True,
'class_weight': None,
'criterion': 'gini',
'max_depth': None,
'max_features': 'auto',
'max_leaf_nodes': None,
'min_impurity_decrease': 0.0,
'min_impurity_split': None,
'min_samples_leaf': 1,
'min_samples_split': 2,
'min_weight_fraction_leaf': 0.0,
'n_estimators': 10,
'n_jobs': None,
'oob_score': False,
'random_state': None,
'verbose': 0,
'warm_start': False}*
我决定关注的超参数是:
n_estimators(森林中的树木数量)max_features(max。节点分割中使用的特征数量，通常为。<数据集中的特征数量)max_depth(最大。每个决策树中的层数)min_samples_split(敏。分割节点之前节点中的数据点数)min_samples_leaf(最小。节点中允许的数据点数)criterion(用于确定决策树停止标准的度量)
使用随机搜索进行调整
为了缩小搜索范围，我首先进行了随机搜索交叉验证。在这里，我使用 k = 10 倍交叉验证(cv=10)、跨越 100 个不同组合(n_iter=100)并同时使用所有可用内核(n_jobs=-1)来执行参数的随机搜索。随机搜索随机选择特征的组合，而不是遍历每个可能的组合。更高的n_iter和cv分别导致更多的组合和更小的过拟合可能性。
*from sklearn.model_selection import RandomizedSearchCVhyperparameters = {'max_features':[None, 'auto', 'sqrt', 'log2'],
'max_depth':[None, 1, 5, 10, 15, 20],
'min_samples_leaf': [1, 2, 4],
'min_samples_split': [2, 5, 10],
'n_estimators': [int(x) for x in np.linspace(start = 10, stop = 100, num = 10)],
'criterion': ['gini', 'entropy']}rf_random = RandomizedSearchCV(clf, hyperparameters, n_iter = 100, cv = 10, verbose=2, random_state=123, n_jobs = -1)rf_random.fit(x_train, y_train)*
一大堆东西出现了。为了获得最佳参数，我调用了:
*rf_random.best_params_*
这给出了:
*{'n_estimators': 100,
'min_samples_split': 5,
'min_samples_leaf': 2,
'max_features': None,
'max_depth': 15,
'criterion': 'entropy'}*
使用网格搜索进行调整
我现在可以指定更小范围的超参数来关注。GridSearchCV 非常适合超参数的微调。
*from sklearn.model_selection import GridSearchCVhyperparameters = {'max_features':[None],
'max_depth':[14, 15, 16],
'min_samples_leaf': [1, 2, 3],
'min_samples_split': [4, 5, 6],
'n_estimators': [90, 100, 110],
'criterion': ['entropy']}rf_grid = GridSearchCV(clf, hyperparameters, cv = 10, n_jobs = -1, verbose = 2)
rf_grid.fit(x_train, y_train)*
这花了我大约 50 分钟。我打了电话:
*rf_grid.best_params_*
这返回了:
*{'criterion': 'entropy',
'max_depth': 14,
'max_features': None,
'min_samples_leaf': 2,
'min_samples_split': 5,
'n_estimators': 100}*
训练分类器
我最后用最佳超参数更新了分类器。
*clf.set_params(criterion = 'entropy', max_features = None, max_depth = 14, min_samples_leaf = 2, min_samples_split = 5, n_estimators = 100)*
测试和评估
然后，我在测试集上测试了更新后的分类器，并根据几个指标对其进行了评估。
准确性得分和 F1 得分
注意，accuracy_score是指校正预测的分数，f1_score是精度(分类器不将阴性样本标记为阳性的能力)和召回率(分类器找到所有阳性样本的能力)的加权平均值。scikit-learn 文档最好地解释了这些概念:
Formula for accuracy_score
**
Formula for f1_score, where beta=1
sklearn.metrics有这些现成的。列表f1_score中分数的顺序对应于类的编码方式，这可以通过使用分类器的.classes_属性来访问。
*from sklearn.metrics import accuracy_score, f1_scoresortedlabels = clf.classes_
accscore = accuracy_score(y_test, y_pred)
f1score = f1_score(y_test, y_pred, average = None)print(accscore)
for i in range:
print((sortedlabels[i],f1score[i]), end=" ")*
这将返回非常令人满意的分数。
*0.99
('GALAXY', 0.9900265957446809) ('QSO', 0.9596774193548387) ('STAR', 0.9959935897435898)*
混淆矩阵
A 混淆矩阵 C 有矩阵元素 C_(i，j)对应于已知在 I 组但预测在 j 组的观测值的数量，换句话说，对角线元素代表正确的预测，而非对角线元素代表错误标记。我们的目标是得到混淆矩阵的大对角线值 C(I，I)。
我们看一下分类报告和混淆矩阵。
*from sklearn.metrics import classification_report, confusion_matrixcm = confusion_matrix(y_test, y_pred, sortedlabels)print(classification_report(y_test, y_pred))
print(cm)*
这将返回:
*precision recall f1-score support
GALAXY 0.98 0.99 0.98 1499
QSO 0.95 0.89 0.92 255
STAR 0.99 1.00 1.00 1246
micro avg 0.98 0.98 0.98 3000
macro avg 0.97 0.96 0.97 3000
weighted avg 0.98 0.98 0.98 3000
[[1481 11 7]
[ 29 226 0]
[ 1 1 1244]]*
可视化混淆矩阵
我一般觉得有用(而且好看！)来绘制混淆矩阵。在我选择的 seaborn 调色板中，较深的颜色意味着更多的条目。
*cm = pd.DataFrame(cm, index=sortedlabels, columns=sortedlabels)sns.set(font_scale=1.2)
sns.heatmap(cm, linewidths=0.5, cmap=sns.light_palette((1, 0.2, 0.6),n_colors=10000), annot=True)
plt.xlabel('Predicted')
plt.ylabel('True')*
如果你有反馈/建设性的批评，请在下面随意评论！
用 Python 预测股票价格
用 100 行代码
投资股票市场过去需要大量资金和经纪人，经纪人会从你的收入中提成。然后 Robinhood 颠覆了这个行业，让你只需投资 1 美元，完全避开经纪人。Robinhood 和类似的应用向任何拥有联网设备的人开放了投资，并让非投资者有机会从最新的科技初创企业中获利。
“space gray iPhone X turned on” by rawpixel on Unsplash
然而，让我们这些不是经济学家或会计师的人自由投资“最热门”或“趋势”股票并不总是最好的财务决策。
成千上万的公司使用软件来预测股票市场的走势，以帮助他们的投资决策。普通的罗宾汉用户没有这个功能。原始预测算法，如时间序列线性回归，可以通过利用 python 包如 scikit-learn 和iexfinance来进行时间序列预测。
这个程序将从网上抓取给定数量的股票，预测它们在一定天数内的价格，并向用户发送短信，通知他们可能适合查看和投资的股票。
设置
为了创建一个程序来预测股票在一定天数内的价值，我们需要使用一些非常有用的 python 包。您需要安装以下软件包:
如果您还没有这些包，您可以通过pip install PACKAGE或者通过克隆 git 库来安装它们。
下面是一个用 pip 安装 numpy 的例子
pip install numpy
用 git
git clone [https://github.com/numpy/numpy](https://github.com/numpy/numpy)
cd numpy
python setup.py install
现在打开您最喜欢的文本编辑器，创建一个新的 python 文件。从导入以下包开始
import numpy as np
from datetime import datetime
import smtplib
import time
from selenium import webdriver#For Prediction
from sklearn.linear_model import LinearRegression
from sklearn import preprocessing, cross_validation, svm#For Stock Data
from iexfinance import Stock
from iexfinance import get_historical_data
注意:python 附带了 datetime、time 和 smtplib 包
为了清理雅虎股票筛选程序，你还需要安装 Chromedriver 来正确使用 Selenium。在这里可以找到和
获得股票
使用 Selenium 包，我们可以从雅虎股票筛选程序中筛选股票的代码缩写。
首先，创建一个带有参数n的函数getStocks，其中 n 是我们希望检索的股票数量。
def getStocks(n):
在函数中创建你的 chrome 驱动程序，然后使用driver.get(url)来检索你想要的网页。我们将导航至https://finance . Yahoo . com/screener/predefined/aggressive _ small _ caps？offset=0 & count=202 这将显示 200 只列在“激进小型股”类别中的股票。如果你去 https://finance.yahoo.com/screener，你会看到雅虎提供的所有筛选类别的列表。然后，您可以根据自己的喜好更改 URL。
#Navigating to the Yahoo stock screenerdriver = webdriver.Chrome(‘**PATH TO CHROME DRIVER**’)url = “https://finance.yahoo.com/screener/predefined/aggressive_small_caps?offset=0&count=202"driver.get(url)
确保将下载 chromedriver 的路径添加到粗体代码所在的位置。
您现在需要创建一个列表来保存 ticker 值stock_list = []。
接下来，我们需要找到 ticker 元素的 XPath，以便能够抓取它们。转到 screener URL 并在您的 web 浏览器中打开开发者工具(Command+Option+I/Control+Shift+I或 F12 for Windows)。
单击“选择元素”按钮
单击股票并检查其属性
最后，复制第一个 ticker 的 XPath，HTML 元素应该如下所示
<a href=”/quote/RAD?p=RAD” title=”Rite Aid Corporation” class=”Fw(b)” data-reactid=”79">RAD</a>
XPath 应该看起来像这样
//*[[@id](http://twitter.com/id)=”scr-res-table”]/div[2]/table/tbody/tr[**1**]/td[1]/a
如果您检查第一个下面的 ticker 属性，您会注意到 XPath 完全相同，除了上面代码中加粗的 1 对于每个 ticker 增加 1。所以第 57 个 ticker XPath 值是
//*[[@id](http://twitter.com/id)=”scr-res-table”]/div[2]/table/tbody/tr[**57**]/td[1]/a
这对我们帮助很大。我们可以简单地创建一个for循环，每次运行时递增该值，并将跑马灯的值存储到我们的stock_list中。
stock_list = []
n += 1for i in range(1, n):
ticker = driver.find_element_by_xpath(‘//*[@id = “scr-res-table”]/div[2]/table/tbody/tr[‘ + str(i) + ‘]/td[1]/a’)stock_list.append(ticker.text)
n是我们的函数getStocks(n)将检索的股票数量。我们必须增加 1，因为 Python 是 0 索引的。然后我们使用值i 来修改每个 ticker 属性的 XPath。
使用driver.quit()退出网络浏览器。我们现在有了所有的股票价格，并准备好预测股票。
我们将创建一个函数来预测下一节中的股票，但现在我们可以创建另一个for循环来遍历列表中的所有股票值，并预测每个股票的价格。
#Using the stock list to predict the future price of the stock a specificed amount of daysfor i in stock_list:
try:
predictData(i, 5)
except:
print("Stock: " + i + " was not predicted")
用 try 和 except 块处理代码(以防我们的股票包不能识别股票代码值)。
预测股票
创建一个新函数predictData，它采用参数stock和days(其中 days 是我们希望预测未来股票的天数)。从 2017 年 1 月 1 日到现在，我们将使用大约 2 年的数据进行预测(尽管您可以使用任何您想要的数据)。设置start = datetime(2017, 1, 1)和end = datetime.now()。然后使用 iexfinance 函数获取给定股票df = get_historical_data(stock, start=start, end=end, output_format=’pandas’)的历史数据。
然后将历史数据导出到一个. csv 文件，为预测创建一个新的虚拟列，并设置forecast_time = int(days)
start = datetime(2017, 1, 1)
end = datetime.now()#Outputting the Historical data into a .csv for later use
df = get_historical_data(stock, start=start, end=end, output_format='pandas')csv_name = ('Exports/' + stock + '_Export.csv')
df.to_csv(csv_name)
df['prediction'] = df['close'].shift(-1)
df.dropna(inplace=True)
forecast_time = int(days)
使用 numpy 操作数组，然后预处理这些值，并创建 X 和 Y 训练和测试值。对于这个预测，我们将使用 test_size 的0.5，这个值给了我最准确的结果。
X = np.array(df.drop(['prediction'], 1))
Y = np.array(df['prediction'])
X = preprocessing.scale(X)
X_prediction = X[-forecast_time:]X_train, X_test, Y_train, Y_test = cross_validation.train_test_split(X, Y, test_size=0.5)
最后，对数据进行线性回归。创建变量clf = LinearRegression()，拟合 X 和 Y 训练数据，并将 X 值预测存储在变量prediction中。
#Performing the Regression on the training dataclf = LinearRegression()
clf.fit(X_train, Y_train)
prediction = (clf.predict(X_prediction))
在下一节中，我们将定义通过 SMS 发送股票预测的函数sendMessage。在predictData函数中添加一个if语句，该语句存储一个字符串作为输出，并调用sendMessage函数，向其传递参数output。
变量output可以包含您认为有用的任何信息。我让它告诉我股票名称，1 天预测和 5 天预测。
#Sending the SMS if the predicted price of the stock is at least 1 greater than the previous closing pricelast_row = df.tail(1)
if (float(prediction[4]) > (float(last_row['close']))):output = ("\n\nStock:" + str(stock) + "\nPrior Close:\n" + str(last_row['close']) + "\n\nPrediction in 1 Day: " + str(prediction[0]) + "\nPrediction in 5 Days: " + str(prediction[4]))sendMessage(output)
发送消息
创建一个将output作为参数的函数sendMessage。要发送短信，我们将使用smtplib包，这样我们就可以通过电子邮件发送短信。
将您的电子邮件用户名、密码和接收号码存储为变量。我的手机运营商是威瑞森，所以我用的是@vtext 域名。多亏了这个网站，这里有一些受欢迎的电话公司的分机号码。
- number@txt.att.net(手机短信)，number@mms.att.net(彩信)
- T-Mobile: number@tmomail.net(短信和彩信)
- number@vtext.com(手机短信)，number@vzwpix.com(彩信)
- sprint:number @ messaging . sprint PCs . com(短信)，number@pm.sprint.com(彩信)
- 维珍手机:number@vmobl.com(短信)，number@vmpix.com(彩信)
def sendMessage(output):username = "EMAIL"
password = "PASSWORD"
vtext = "PHONENUMBER@vtext.com"
使用以下行发送带有正确消息的 SMS
message = output
msg = """From: %s To: %s %s""" % (username, vtext, message)server = smtplib.SMTP('smtp.gmail.com', 587)
server.starttls()
server.login(username, password)
server.sendmail(username, vtext, msg)
server.quit()
运行程序
最后，创建一个 main 方法来运行程序。我们将把预测的股票数量设置为 200。
if __name__ == '__main__':
getStocks(200)
结论
仅在 10 只股票上运行预测，实际 1 日价格和 1 日预测价格之间的平均百分比误差为 9.02% ，其中 5 日百分比误差是令人惊讶的 5.90% 。这意味着，平均而言，5 天的预测仅比实际价格低0.14 美元。
这些结果可以归因于样本量小，但无论如何，它们都是有希望的，可以在你投资股票时提供很大的帮助。
[## 使用 Tweepy 在 Python 中创建一个 Twitter 机器人
Twitter 大约有 15%是由机器人组成的，我想尝试一下。我谷歌了如何创建一个推特机器人…
medium.freecodecamp.org](https://medium.freecodecamp.org/creating-a-twitter-bot-in-python-with-tweepy-ac524157a607) [## 推特数据分析
在使用 Tweepy 和 Python 以及这段代码创建了免费的 Wtr bot 之后，我想用一种方式来看看 Twitter 用户是如何…
medium.co](https://medium.com/@lucaskohorst/twitter-data-analysis-d5fb7d42ebc7) [## Python 中的新闻倾向性排名
我最近在《华盛顿邮报》上读到一篇文章，题为“根据……从自由派到保守派对媒体进行排名。”
towardsdatascience.com](/ranking-news-bias-in-python-e9bb5d1ba93f)
Kaggle 内核:预测学生的成绩
教育数据挖掘是指将数据挖掘应用于教育数据集。在大多数情况下，EDM 类似于普通的数据挖掘。然而，我们应该考虑教育数据集的具体特征。数据通常有多个层次。此外，许多研究是在一个机构进行的，这使得大多数研究不太容易概括。如果你想了解更多关于 EDM 的内容，你可以看看 Romero 和 Ventura 的论文。两位作者都撰写了关于该主题的多个文献综述，这些综述提供了该领域随时间推移的发展的清晰视图。EDM 中与预测相关的重要主题有:预测注册人数、预测学生表现和预测流失。
在 Kaggle 上，我发现了这个学生成绩的数据集。在分析中，我观察了各种可视化效果，并比较了基于树的机器学习算法对学生成绩的预测。尽管数据集很小，但我们能够达到几乎 82%的准确率。根据想象，女性在这项测试中的得分似乎比男性高。在线环境中被访问资源的数量是最终成绩最重要的预测因素。
您可以通过点击此链接查看 Kaggle 笔记本，也可以点击下面的静态笔记本。
预测坦桑尼亚水泵的运行状况
水是生命和环境的基础；它在经济和社会发展活动中发挥着核心作用。水涉及人类生活的所有领域，包括:家庭、牲畜、渔业、野生动物、工业和能源、娱乐和其他社会经济活动。它通过加强粮食安全、家庭卫生和环境，在减贫方面发挥着关键作用。安全和清洁的水的供应提高了生活水平，而水的不足造成了严重的健康风险，并导致生活水平和预期寿命的下降。
坦桑尼亚的主要淡水来源包括:湖泊、河流、小溪、水坝和地下水。然而，这些并没有在全国范围内得到很好的分配。一些地区既缺乏地表水资源，也缺乏地下水资源。
人口增长和城市化对可用水的数量和质量造成严重压力。当前和未来人类生活和环境的可持续性主要取决于适当的水资源管理。
Image Courtesy: flickr user christophercjensen
首先分析坦桑尼亚现有供水点的功能状态。在此分析中，我们将研究坦桑尼亚水泵的数据集，以预测供水点的运行状况。这给了我们分析的灵感。
通过预测功能正常但需要维修的水泵，降低了坦桑尼亚水利部的总成本。这可以改善水泵的维护操作，并确保坦桑尼亚各地的社区都能获得清洁的饮用水。
首先，让我们看看现有的数据。
训练集值 —训练集的独立变量。'https://S3 . amazonaws . com/driven data/data/7/public/4910797 b-ee55-40a 7-8668-10 EFD 5 C1 b 960 . CSV’
训练集标签 —训练集值中每一行的因变量(status_group)。'https://S3 . amazonaws . com/driven data/data/7/public/0 BF 8 BC 6 e-30d 0-4c 50-956 a-603 fc 693d 966 . CSV’
导入数据
training_values = pd.read_csv('[https://s3.amazonaws.com/drivendata/data/7/public/4910797b-ee55-40a7-8668-10efd5c1b960.csv](https://s3.amazonaws.com/drivendata/data/7/public/4910797b-ee55-40a7-8668-10efd5c1b960.csv)', encoding = 'iso-8859-1')training_lables = pd.read_csv('[https://s3.amazonaws.com/drivendata/data/7/public/0bf8bc6e-30d0-4c50-956a-603fc693d966.csv](https://s3.amazonaws.com/drivendata/data/7/public/0bf8bc6e-30d0-4c50-956a-603fc693d966.csv)', encoding = 'iso-8859-1')
数据集中的要素
list(training_values.columns.values)
我们有 40 个预测变量。在数据集中的 40 个特征中，我们有 31 个分类变量、7 个数值变量和 2 个日期变量。
['id', 'amount_tsh', 'date_recorded', 'funder', 'gps_height', 'installer', 'longitude', 'latitude', 'wpt_name', 'num_private',
'basin', 'subvillage', 'region', 'region_code', 'district_code',
'lga', 'ward', 'population', 'public_meeting', 'recorded_by',
'scheme_management', 'scheme_name', 'permit', 'construction_year',
'extraction_type', 'extraction_type_group', 'extraction_type_class',
'management', 'management_group', 'payment', 'payment_type',
'water_quality', 'quality_group', 'quantity', 'quantity_group',
'source', 'source_type', 'source_class', 'waterpoint_type',
'waterpoint_type_group']
水泵是这三种类型中的任何一种
- 功能性的，
- 不起作用，或
- 功能正常但需要修理。
training_lables.head(5)
现在，让我们检查 status_group 中泵的每种功能类型的计数，以便了解泵状态的功能场景。
training_lables.status_group.value_counts()functional 32259
non functional 22824
functional needs repair 4317
Name: status_group, dtype: int64
根据上述结果，我们可以说，有 54.31%的功能泵、38.42%的非功能泵和 7.27%的功能泵需要修理。
根据上面的数字，我们可以粗略估计，如果我们从数据中随机抽取一个泵作为功能泵，有 54.31%的可能性。
删除具有相似数据表示的特征
(提取类型，提取类型组，提取类型类)，(付款，付款类型)，
(水质，质量组)，(水源，水源类)，(子村，区域，区域代码，地区代码，lga，ward)，以及(水点类型，水点类型组)的特征组都包含不同粒度的数据的相似表示。因此，在训练过程中，我们冒着过度拟合数据的风险，在我们的分析中包括所有的特征，这些特征可以被丢弃。
- 可以删除 id，因为它对于每个实例都是唯一的。
- num_private 大约 99%为零
数据清理和分析
从数据来看，一些基于人类直觉的特征似乎是有区别的。amount_tsh(供水点可用的水量)、gps_height、basin、installer、population、scheme_management、施工年份、抽水类型、管理组、水质、付款类型、水源和供水点类型看起来对识别水泵状态非常重要。
我们的特征中有空值，需要更新以更好地训练我们的模型。
在清除空值的过程中，我们发现一些特征具有很高的 arity。因此，对于具有高 arity 的特征，让我们基于频率保留前 10 个值，并将所有剩余的值分配给第 11 个合成值作为“其他”。
将全部数据分成 10 组后，我们创建了一个数据透视表，将数据分组到子组中，该子组包含特定特性的每个合成值在每个条件下的泵数量，我们还创建了一个热图，以图形方式查看数据。例如，查看名为 basin 的特征，其实例广泛分布在所有类别中(这可能对我们的模型有影响)。
piv_df= training_df[['basin','status_group','status_values']]
piv_table = piv_df.pivot_table(index='basin',
columns='status_group', aggfunc='count')
piv_table
Heat Map constructed using the Pivot table created in the previous step
在丢弃不需要的特征并清理数据集后，我们剩下 20 个特征可用于训练我们的模型。我们将用于培训的功能有:
['amount_tsh', 'days_since_recorded', 'funder', 'gps_height',
'installer', 'basin', 'subvillage', 'population', 'public_meeting',
'scheme_management', 'scheme_name', 'permit', 'construction_year',
'extraction_type', 'payment_type', 'water_quality',
'quantity_group', 'source_type', 'source_class', 'waterpoint_type',
'waterpoint_type_group']
为分类器获取数据的一种常见方法是将可用数据分成两个集合，一个训练集和一个测试集。分类模型建立在训练集上，并应用于测试集。测试集从未被模型看到过，因此当模型应用于看不见的数据时，最终的性能将是一个很好的指导。
通常，选择的比例是 70%用于训练集，30%用于测试。这个想法是，更多的训练数据是一件好事，因为它使分类模型更好，而更多的测试数据使误差估计更准确。
X_train, X_test, y_train, y_test = sklearn.model_selection.train_test_split(training, test, test_size = 0.3, random_state = 0)
对于分类，我们将使用集成方法，集成学习使用几种用不同学习算法构建的机器学习模型来提高预测的准确性。它基本上是对几个模型的预测进行平均，从而得到一个更好的预测。我将在这个分析中使用的一些集成方法是随机森林，额外树和梯度推进。
我们将尝试不同的方法，选择最适合我们数据的方法。
随机森林分类器
rfc = sklearn.ensemble.RandomForestClassifier(n_estimators=1000,
min_samples_split=6,
criterion='gini',
max_features='auto',
oob_score=True,
random_state=1,
n_jobs=-1)
决策树:
dtc = DecisionTreeClassifier(criterion='gini',
max_depth = 10,
max_features = 'auto',
random_state = 1,
splitter = 'best')
额外的树分类器:
ETC = ExtraTreesClassifier(n_estimators=1000, min_samples_split=10)
梯度增强:
param_grid = {'learning_rate': [0.075,0.07],
'max_depth': [6,7],
'min_samples_leaf': [7,8],
'max_features': [1.0],
'n_estimators':[100, 200]}gbc = GridSearchCV(estimator=GradientBoostingClassifier(),
param_grid=param_grid,
n_jobs=-1)
每个分类器的验证分数如下。
对于每个分类器，我们构建了混淆矩阵来检查每个值中的真阳性和假阳性。让我们显示在随机森林分类器和额外树分类器的情况下表现良好的分类器的混淆矩阵。
随机森林分类器的混淆矩阵
额外树分类器的混淆矩阵
综合考虑训练数据集和测试数据集的准确率、F-1、查准率和查全率，随机森林分类器表现良好，根据以上证据，我们可以选择随机森林作为我们的分类模型。
当我们选择模型时，让我们检查每个特性的重要性，看看哪些特性影响了我们的模型。
['amount_tsh', 'days_since_recorded', 'funder', 'gps_height',
'installer', 'basin', 'subvillage', 'population', 'public_meeting',
'scheme_management', 'scheme_name', 'permit', 'construction_year',
'extraction_type', 'payment_type', 'water_quality',
'quantity_group', 'source_type', 'source_class', 'waterpoint_type',
'waterpoint_type_group']
在所有特性中，数量组、子村庄和 gps 高度似乎是最重要的。通过观察特征的重要性，我们还可以推断出水泵的面积和特定区域的可用水量将对水泵的功能状态产生更大的影响。
最终结果:
分析让我们相信数据集很难分类。对分类进行优先排序，以确定需要维修的泵。一个标准泵的价格从 100 美元到 2000 美元不等。安装这种泵需要钻井，费用可能在 1000 美元到 3000 美元之间。另一方面，维护水泵只需要花费几十美元，这可以帮助坦桑尼亚的水行业获得数百万美元。
未来工作:
- 我们可以从外部来源收集准确的人口数据，并将其添加到我们的数据集中，并检查特定区域的人口是否影响泵的功能状态，这可以帮助他们预测该区域并在该区域采取更多措施。
- 我们可以通过使用历史数据来确定水泵的寿命，这可以帮助坦桑尼亚水利部在水泵无法使用之前对其进行维修。
- 为了提高预测的性能和执行速度，我们可以使用极限梯度提升方法，在 python 中我们可以使用 XGBoost 库来处理它。
项目地点:
T3【https://github.com/ksdkalluri/identifying_faulty_pumps】T5
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
