与二重身战斗
原文：https://towardsdatascience.com/fighting-doppelgängers-2fc28762e169
如何去除数据中的邪恶双胞胎，减少特征空间
摘要
给定一个包含许多变量的大型数据集，其中一些变量可能代表相同的现象，或者彼此带来重叠的信息。面对这个问题的一个常见的预处理策略是分析相关性。在本文中，我描述了一种策略，用于扫描线性相关性网络，搜索要保留的最佳变量，并丢弃多余的变量。此外，我提供了一个 R 包来轻松地执行这项工作。
当处理由记录机械事件的传感器产生的数据时，包含成百上千个变量的大型数据集是常见的。在这些情况下，许多变量可以作为预测某些目标度量的候选变量。然而，尤其是在工业环境中，数据可以包括完全线性相关或非常相关的变量。
例如，传感器可以从相同的过程中提取几个特征，作为相同基的线性变换(如一组记录的总和、它们的平均值等)。).在其他情况下，有真正不同的测量，但本质上相关，或代表同一现象的两个相反方面(想象一种化学混合物的两个互补元素)。
非常相关的变量是多余的，通常它们不会带来额外的信息。而且，完全因变量会让一些机器学习算法崩溃。一般来说，这些“分身”使模型构建过程变得复杂，增加了计算量，在某些情况下会产生过度拟合。
所以，问题是:如何修剪数据集以去除邪恶的二重身？
设置相关性的阈值
我们用来识别重复或近似重复变量的有力武器是线性相关指数。
当两个变量的相关性为~|1|时，保留其中一个变量从纯统计学的角度来看是无关紧要的。
当相关性低于|1|时，情况变得复杂，但是仍然相关，并且基于理论的选择是不切实际的。在这些情况下，我们需要一个标准来确定相关性是否相关。一种策略是对相关值设置阈值。
给定一个阈值，我们可以关注强关系来实施修剪策略。
修剪策略
修剪的关键思想是丢弃尽可能多的变量，保留尽可能多的信息。
想象一种情况，其中三个变量 A、B 和 C 有很强的关系。b 与 A 和 C 都具有超过阈值的相关性，而 A 和 C 彼此具有低于阈值的相关性。下图恢复了这些关系，将每个变量表示为网络的一个节点，如果两个变量的相关性高于设定的阈值，则连接这两个变量。
图片来自作者——在 CC0 许可下发布
面对这种情况，修剪的方法可以是去掉 B ，保留相关性较低的 A 和 C 。然而，还有一种更为节俭的方法:只保留 B ，去掉 A 和 C 。 B 可以被认为是网络的代表，它可以作为整个社区的代言人。
因此，正如您所看到的，选择并不是唯一的。
为了更好地理解这两种方法之间的差异，我引入了一个具有更复杂关系的数据集。
一个更复杂的例子
下图是七个模拟变量( alpha 、 beta 、 gamma 、 delta 、 epsilon 、 zeta 、 eta )得到的网络关联。变量在图中的空间分布是根据它们之间关系的强度:相邻的变量之间的相关性更强。仅表示高于阈值|0.7| 的关系。
源数据在我的 R 包分身中包含的数据帧测量中可用。
图片来自作者——在 CC0 许可下发布
我手工构建了上面的图，但是如果你是一个 R 用户，可以考虑使用包[qgraph](http://sachaepskamp.com/qgraph/examples)将一个关联矩阵可视化为一个网络。如果你是 Python 用户，我建议看一下这篇关于 TDS 的文章。从现在开始，我使用 R 语言来构建例子。
我在下面报告完整的相关矩阵。
r <- cor(doppelganger::measures)
diag(r) <- NA
round(r, 3)
alpha beta gamma delta epsilon zeta eta
alpha NA 0.129 0.833 0.815 0.715 0.207 0.645
beta 0.129 NA 0.152 0.108 0.056 0.770 0.097
gamma 0.833 0.152 NA 0.984 0.620 0.193 0.702
delta 0.815 0.108 0.984 NA 0.588 0.167 0.711
epsilon 0.715 0.056 0.620 0.588 NA 0.072 0.519
zeta 0.207 0.770 0.193 0.167 0.072 NA 0.137
eta 0.645 0.097 0.702 0.711 0.519 0.137 NA
该图显示了两个社区的存在:一个是对 zeta - beta 的社区，另一个包括剩余的变量。与前一个天真的例子相比，情况更加混乱。如何进行？
区分变量的优先级
我们可以根据变量的中心度对变量进行排序，中心度意为一个节点与其他节点的连接强度。最中心的节点对应于网络中最具代表性的变量。
我们可以将中心度计算为与变量相关的相关性绝对值的平均值。下面的 R 代码计算每个节点的中心性，并按降序排列。
sort(colMeans(abs(r), na.rm=TRUE), decreasing=TRUE)
gamma delta alpha eta epsilon zeta beta
0.5805754 0.5620824 0.5572697 0.4685961 0.4282394 0.2575334 0.2188217
根据中心度的排序允许我们在选择保留什么和丢弃什么时对变量进行优先级排序。我们面前有两条路:
- 保留最中心的变量，以它们的代表性的名义(中心性准则)；
- 保留最外围的变量，以它们独立的名义(外围性判据)。
无论如何，我们可以应用相同的修剪程序，但是:
在“中心性”情况下，我们按照中心性程度向量以降序扫描变量，而在“边缘性”情况下，我们按照升序扫描变量。
我们开始保留排名中的第一个变量，去掉它的二重身。然后，让我们保留第二个变量，去掉它的二重身。等等，直到结束。很简单，不是吗？
算法应用
现在，将这个过程应用到上面显示的例子中。我描述了遵循中心性标准的操作。
**第一步。**最中心的变量是伽马。因此，我们“开启”伽马和“关闭”德尔塔、阿尔法和埃塔，它们与伽马具有过阈值相关性。在下图中，“关闭”的变量用浅灰色标记，它们的修正用虚线表示。
**第二步。**在伽玛之后，按照优先级顺序，第一个未被丢弃的变量是ε。然而，ε与其余变量的相关性没有超过阈值。它的复制品已经被丢弃，所以不需要任何操作。
**步骤三。**以下变量为ζ，与 beta 有很强的相关性。所以，我们去掉了贝塔，保留了泽塔。
图片来自作者——在 CC0 许可下发布
在迭代结束时，我们将继续使用伽玛、ε和泽塔。三个变量而不是七个。
中心性对边缘性
中心性标准利用了变量总体中孤立社区的存在。从理论的角度来看，当我们在社区内部有强连接而在社区之间有弱连接时，这种方法发挥了最大的作用。相反，外围性标准可以在包含一些子群的更同质的网络中发挥最大作用。
中心性标准倾向于保留更少但更相关的变量。边缘性标准倾向于保留更多的变量，但相关性较低，然而，它们不一定是每个亚组的最具代表性的。
R 的二重身套装
使用 R 包二重身只用一行代码就可以完成前面的步骤。
加载包和数据集:
library(doppelganger)
measures # comes with the package
# A tibble: 100 × 7
alpha beta gamma delta epsilon zeta eta
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 0.0185 -0.564 -0.526 -0.368 -0.934 -0.878 0.0913
2 0.120 0.00233 -0.414 -0.515 -0.193 -0.464 -0.324
3 -0.767 0.0641 -1.51 -1.26 -1.25 -0.0825 -1.70
4 -2.12 0.105 -2.38 -2.37 -1.75 -0.323 -0.426
5 -1.48 -0.126 -0.681 -0.340 -1.87 0.0678 -0.791
6 1.69 0.268 2.39 2.59 1.06 0.115 1.44
7 1.26 -1.17 0.290 0.231 1.93 -0.683 -0.253
8 1.23 1.40 1.10 0.997 0.717 1.07 0.265
9 0.295 -0.0263 -0.0221 -0.115 -0.125 0.179 -0.978
10 0.157 0.268 0.215 0.454 0.531 1.25 -0.0989
# … with 90 more rows
再现上述变量的修剪过程:
dg <- doppelganger(measures, priority="c", threshold=0.7)
同名函数doppelganger创建一个包含一些槽的列表对象，包括keep和drop，分别收集变量的名字来保存和删除。
dg$keep
[1] "gamma" "epsilon" "zeta"
如果您需要检查特定变量的相关性，该包提供了函数inspect_corrs:
inspect_corrs(dg, "gamma")
# A tibble: 6 × 3
variable is_alias cor
<chr> <lgl> <dbl>
1 delta TRUE 0.984
2 alpha TRUE 0.833
3 eta TRUE 0.702
4 epsilon FALSE 0.620
5 zeta FALSE 0.193
6 beta FALSE 0.152
非缺失数据的加权
软件包doppelganger通过用可用于指数计算的观察值数量对每个相关性进行加权来计算中心度。这种策略在存在缺失数据的情况下是有效的，因为它允许“奖励”具有许多与其他变量的填充记录相匹配的填充记录的变量。
我们希望保留最好的变量，但也保留尽可能多的记录！
对候选预测器池的修剪可以有效地帮助建立一个简单有效的模型，尤其是用于推理目的。最佳实践是将选择程序建立在潜在预测者的理论假设或预测能力的基础上。然而，通常，由于缺乏信息或数据集的维度，这是不切实际的。在这些情况下，一个不太精确的方法——就像本文中描述的方法——可以派上用场。
最后，记住什么是精确和回忆，不要在面试中害怕这些问题
一个直观的指南，一步一步地绘制混淆矩阵，并推导出精确度、召回率和 F1 分数的定义
预览。作者图片
理解混淆矩阵、精确度和召回率的概念对于评估任何分类算法的质量是最重要的。而且虽然这在现实任务中并不那么重要(首先，你会在 90%的情况下使用 F1-score，其次，你会花 10 秒钟谷歌一下来刷新你的记忆)，但是理解这些东西还是很有用的。另外，这个问题在数据科学和机器学习职位面试中极其常见。
我(以及我认识的一些更有经验的人)仍然需要花几秒钟来记住什么是什么。因此，我向你展示我的一步一步的指导，如何记住如何绘制混淆矩阵，并从中推导出精确度、召回率和 F1 分数的定义。
如果你已经知道什么是精确或者回忆，但是仍然经常在定义上感到困惑，你应该读读这篇文章。
这篇文章并不解释什么是精确或召回。这篇文章解释了如何记住什么是精确或回忆。
一步一步画困惑矩阵
当我们开始画混淆矩阵时，让我们画一个 2x2 的矩阵。
我们记得在某个地方一定有一个实际的类，和一个预测的类。取决于你如何排列它们，你的矩阵可能会相对于主对角线镜像。请记住，你在其他来源看到的矩阵可能与我画的不同。
为了避免混淆，让我们将预测的类放在的上面(使用直觉，单词Predirect 中的字母 P 是由一个竖条组成的，上面有一个圆圈**)。**
之后，让第一个元素(实际类的情况下为上，预测类的情况下为左)为负例，第二个元素(实际类的情况下为下，预测类的情况下为右)为正例。这里我用的直觉是，我们是从左到右，从上到下，从一个较小的数向一个较大的数移动(而负类是 0，正类是 1)。
画出混淆矩阵。第一步。作者图片
第二步也很简单——让我们填充主对角线上的元素。这些是我们的分类器正确识别的元素——真阳性和真阴性(真在这里意味着我们的分类器是对)。
画出混淆矩阵。第二步。作者图片
- 真阳性(TP) 是被分类为阳性且实际为阳性的元素——预测阳性和实际阳性。**
- 真阴性(TN) 是已经被归类为阴性的元素，实际为阴性——预测阴性和实际阴性。**
之后，我们还要在矩阵中放置假阳性和假阴性(假这里指我们的分类器是错的)。这就是问题开始出现的地方。
为了记住假阳性和假阴性的位置，记住在这些预测中分类器是错误的**。这意味着，尽管它们不正确，但它们是预测的值，因此它们对应于预测的类。**
画出混淆矩阵。第三步。作者图片
让我们给出正式的定义:
- 假阳性(FP) (也称为 I 型错误)是一个已经被归类为阳性但实际为阴性的元素——预测阳性，实际阴性。**
- 假阴性(FN) (又称 II 型错误 ) 是已经被归类为阴性但实际为阳性的元素——预测阴性，实际阳性。**
为了快速理解它们之间的区别，不至于混淆，只要记住分类器的预测是包含在题目中的，并且这个预测是错误的就够了。比如假正。它是正的，因为我们的分类器是这样预测的，但它是错误的。所以预测类是正，而实际类不是正的，所以负。
为了不混淆 TN、TP、FN 和 FP，问两个问题就足够了:
- 分类器预测了什么？负或正。
- 分类器正确吗？真或假。
两个答案都包含在问题中。
如何永远不要混淆 TN，TP，FN，FP？作者图片
现在我们有一个现成的混淆矩阵，这已经是一半了！
混乱矩阵。作者图片
如何停止恐惧，开始理解精密和回忆
精确度和召回率公式都很简单。它是真阳性与真阳性加上某个值的比率。分母(真阳性加上某个值)是实际阳性或预测阳性。那有什么意义？让我们弄清楚。
精度公式
记住精度公式中的所有项都包含字母 P ( P 精度— P ，很容易记忆)。
所以精度是 TP 与 (TP + FP) 的比值，其中 (TP + FP) 是预测阳性。
精确公式。作者图片
回忆公式
另一方面，召回是 TP 与 (TP + FN) 的比值，其中 (TP + FN) 是实际阳性。
回忆公式。作者图片
这些到底有什么意义？
我不认为有必要给出奇怪的定义来解释一切，但只会使其复杂化(比如回忆是一种检测概率)。我不否认它们是有用的，但它们往往弊大于利。最后，你可以很容易地发现他们在阅读其他文章。
我会尽可能直观地解释一切。我建议你在阅读下面两条语句时，保持精确度，并在眼前回忆公式。
- 如果 FP = 0，则精度等于 1。
所以精度高关系到****假阳性率低或者说I 型错误量小。**** - 如果 FN = 0，召回等于 1。
所以高召回率与低假阴性率或少量 II 型错误有关。
如何记住这些？高精度—很少的假阳性—很少的 I 型错误。记住精度总是精度-召回对中的第一个(I 类型错误)；我们已经记住的是——精度公式中的所有项都包含字母 P (很少有误正)。**
这是让您理解何时使用特定指标的基础。
- 如果 FP 价格高，用精度(预测借款人会还贷，虽然他破产了；“好”的借款人是正类)；
- 而如果 FN 价格高，用召回(预测病人是健康的，虽然他生病了；病是正类)。
这些定义取决于对类的选择(哪个类是负的，哪个是正的)——所以要时刻注意这个。例如，交换上面例子中的类，看看会发生什么。
f1-分数
F1-score 是一种将精确度和召回率结合成一个指标的方法。我们需要这个，因为只用一个数字评估质量很容易，而不是两个。我们不能使用简单平均值，因为它不考虑极小的值。
例如，如果我们让 F1-score 是精度和召回率的平均值，则相同的值得到:
精度= 0.5 和召回率= 0.5 的分类器 1 和
精度= 1 和召回率= 0 的分类器 2。
两个分类器都将得到 F1 分数= 0.5，尽管很明显第一个分类器要好得多。当然，这是个玩具例子，但离真相不远。
事实上，我们想要这样一个表达式，如果两个元素中至少有一个比一个小得多，那么这个表达式就是小的。这样的表达式称为加权平均值，或调和平均值。****
f1-得分公式。作者图片
这个表达式的分子是乘积，这是合乎逻辑的。如果其中一项比一项小得多，那么它们的乘积也将分别比一项小得多。需要乘以 2，以便 F1 分数像所有常用指标一样从 0 变为 1。
请记住，还有其他分类指标。比如非常受欢迎的AUC-ROC—ROC(受试者工作特性)曲线下的区域(你可以在这里了解一下)或者不太受欢迎的 AUC-PRC — 精确召回曲线下的区域。
结论
让我总结一下我们下面列出的所有概念:
- 真阳性(TP) 是被分类为阳性且实际为阳性的元素——预测阳性，实际阳性；****
- 真阴性(TN) 是已经被归类为阴性的元素，实际上是阴性的——预测阴性和实际阴性；****
- 假阳性(FP) (也称为 I 型错误)是已经被分类为阳性但实际为阴性的元素——预测阳性，实际阴性；****
- 假阴性(FN) (也称 II 型错误 ) 是已经被归类为阴性但实际为阳性的元素——预测阴性，实际阳性；****
- 高精度意味着低 FP 率(少数 I 型错误)；
- 高召回意味着低 FN 率(少数 II 型错误)；
- F1-score 是一个将精度和召回率结合成一个数字的好方法。
总结。作者图片
可能对你有用的资源
感谢您的阅读！
- 我希望这些材料对你有用。在 Medium 上关注我以获得更多类似的文章。
- 如果您有任何问题或意见，我将很高兴得到任何反馈。在评论中问我，或者通过 LinkedIn 或 Twitter 联系。
- 为了支持我作为一名作家，并获得数以千计的其他媒体文章，使用我的推荐链接获得媒体会员资格(不收取额外费用)。
使用 Upgini 查找机器学习模型的信息特征
用于丰富特性的 Python 库
Kushagra Kevat 在 Unsplash 上拍摄的照片
特征工程是构建机器学习模型的重要组成部分。我们通常会花费大量时间来寻找新的功能，为机器学习模型增加信息能力。
找到一个有助于预测目标变量的特征并不是一件容易的事情。我们最终可能会创建出对模型没有或增加很少价值的特性。另一方面，如果我们最终创建了一个增强模型的准确性和健壮性的特征。找到这样的特征几乎总是比超参数优化更有帮助。
没有免费的午餐，所以我们可能需要搜索大量的数据源或执行广泛的探索性数据分析，以便能够找到将产生重大影响的功能。
Upgini 正是为此目的而创建的开源 Python 库。它所做的是搜索成千上万的公共和社区数据源，然后创建将提高模型性能的要素。
特征工程可以被认为是一个两步任务。第一步是找到数据(即特征)，第二步是检查它与问题的相关性。Upgini 完成了这两个步骤。它不仅发现特征，而且评估它们。
Upgini 的 enricher 对象返回现有功能和新功能的 SHAP 值，以便我们了解新功能有多有用。SHAP 值经常用于机器学习模型的可解释性和可解释性。
让我们看一个例子来演示如何使用 Upgini 来丰富特性。第一步是安装 Upgini，可通过 pip 完成，如下所示:
pip install upgini
我们将使用 OpenML 网站上有 CCO 许可的 intuit 数据集。该任务是电子邮件响应预测，因此数据集包含一些要素和一个二进制目标变量。我们将使用 Pandas 来读取和清理数据集。
import pandas as pddata = pd.read_csv("csv_result-intuit.csv")data.head()
(图片由作者提供)
res1 是目标变量，取值为 0 和 1。数据集需要清理一下。
- id 功能是多余的，因此我们将删除它们。
- 邮政编码 99999 无效，因此我们也将删除这些行。
data = data.drop(["id","id.1"], axis=1)data = data[data.zip != 99999]
搜索关键字
Upgini 通过搜索关键字找到特征。它目前支持以下功能:
- 日期/日期时间
- 邮政编码
- 国家
- 电话号码
- 散列电子邮件
- IP 地址
当创建一个 enricher 对象时，我们需要指出包含搜索关键字的列，在我们的例子中是邮政编码。需要提及的重要一点是，在使用邮政编码时，我们还需要定义国家。因此，我们将创建一个列来显示这些邮政编码所在的国家，即美国。
data["country"] = "US"
数据清理和预处理
有些行具有重复的特征值，但目标值不同。我们还将使用 Pandas 的删除重复项功能来删除它们。
cols = list(data.columns)
cols.remove("res1")data = data.drop_duplicates(subset=cols)
原始数据集包含 20000 个观察值(即行)，就目标值而言，这是非常不平衡的。大约 95%的目标值为 0，这意味着没有响应。处理不平衡数据集有不同的方法，但这不是本文的主题。因此，我们将从数据集中提取一个不那么不平衡的样本。
data_1 = data[data.res1==1]
data_0 = data[data.res1==0].sample(n=4000)data = pd.concat([data_0, data_1], axis=0, ignore_index=True)data["res1"].value_counts(normalize=True)**Output**
0 0.807591
1 0.192409
现在% 80.76%的目标值是 0。
列车测试分离
下一步是将数据分成训练集和测试集。我们将通过使用 Scikit-learn 的列车测试分割功能来执行此操作。
from sklearn.model_selection import train_test_splitX = data.drop(columns=["res1"])
y = data["res1"]X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.25, random_state=42
)
使用 Upgini 查找功能
我们现在将创建一个 Upgini 的 enricher 对象。
from upgini import FeaturesEnricher, SearchKeyenricher = FeaturesEnricher(
search_keys = {
"zip": SearchKey.POSTAL_CODE,
"country": SearchKey.COUNTRY
},
api_key = "WRITE YOUR API KEY HERE"
)
搜索键是通过将列名和搜索键类型写成 Python 字典中的键值来定义的。另一个参数是 api 密匙，你可以通过注册 Upgini 免费获得。
下一步是使训练集适合刚刚创建的 enricher 对象。
enricher.fit(X_train,
y_train,
eval_set=[(X_test, y_test)]
)
(图片由作者提供)
发现了 3 个新特性。我们还可以看到现有要素的 SHAP 值和 SHAP 值。其他特性的 SHAP 值为 0，所以我没有把它们包含在上面的截图中。
衡量绩效
最后一步是衡量新特性对模型性能的影响。我将使用一个随机的森林分类器来完成这个任务，但是你可以随意使用任何你喜欢的算法。让我们首先创建模型。
from sklearn.ensemble import RandomForestClassifierclf = RandomForestClassifier(
n_estimators=100,
max_depth=3,
random_state=0
)
我们将使用 enricher 对象的计算度量方法来评估新特性。
enricher.calculate_metrics(
X_train, y_train,
eval_set = [(X_test, y_test)],
estimator = clf,
scoring = "accuracy"
)
(图片由作者提供)
性能改进显示在提升列中。训练集的精度提高了%6.8，而测试集的精度提高了%1.2。考虑到获得这些特性是多么简单，这实际上是一个很大的改进。
结论
在某些任务中，1%的精度改进可能至关重要。考虑一个有数百个位置的零售连锁店。需求预测增加 1%可以在库存规划、供应链运营和其他运营工作中节省数千美元。
正如我们在上面的例子中所看到的，用 Upgini 搜索新特性只需要几个简单的步骤。因此，考虑到它的潜力和易用性，它绝对值得一试。
你可以成为 媒介会员 解锁我的全部写作权限，外加其余媒介。如果你已经是了，别忘了订阅https://sonery.medium.com/subscribe如果你想在我发表新文章时收到电子邮件。
感谢您的阅读。如果您有任何反馈，请告诉我。
使用 Python 找到两个角度之间的差异
原文：https://towardsdatascience.com/find-the-difference-between-two-angles-using-python-ac553a1ef5c0
因为在角度从 355 度变为 5 度之前，这都是一种乐趣和游戏
在我的业务领域(海上风力)，我花了很多时间分析风力数据，但直到最近我才尝试使用 python 来分析风向随时间的变化。所以直到最近我才意识到我不知道该怎么做。
看起来很简单，如果方向从 350 度变到 320 度，从另一个中减去一个，得到绝对值，这就是 30 度的变化。但是等一下，等一下…
从 350 度到 5 度的变化是 15 度或 345 度的变化，这取决于变化发生的方向。同样，从 250 度到 200 度的变化可以是 50 度，也可以是 450 度。如果您正在处理大型数据集(或者只是通常重视您的时间)，打开巨大的文件，手动检查每个文件并编辑正确的结果是不可行的。你和你的电脑都可能会认输，所以需要一个自动化的解决方案。
在整篇文章中，我将使用风向作为例子，但这适用于任何使用角度或弧度的东西，并寻求随时间的变化。
目录
1。必要的假设
2。代码
3。解释
必要的假设
如前所述，当方向通过北方，或零度时，我们的简单方法就失败了。将度数转换成弧度也不能解决问题(相信我，我试过了)。相反，我们需要对方向行为进行假设。
如果方向从 350 度变为 5 度，那就是 15 度的变化，只有如果，**是顺时针方向。**即如果风向从 350 度变化到 351 度、352 度……359，0，1 … 5 度。也有可能风向已经从 350 度变成了 349 度，348 度…6、5 度或逆时针方向。这两个选项如下所示。
角度方向可以逆时针(红色)或顺时针(蓝色)改变。作者提供的图片-未按比例绘制。
这两种情况实际上是否可能取决于数据类型和数据粒度。例如，一艘船的航向有可能在三秒钟内改变 190 度吗？极端天气会导致极端风向改变吗？
为了使建议的解决方案起作用，我们必须假设不可能在给定的周期内方向改变超过 180 度。即，如果数据是一分钟粒度的时间序列，从分钟 1 到分钟 2，超过 180 度的真正方向改变是不可能的。
如果这对您的数据集不起作用，请格外小心！考虑你拥有的可以证明或反驳这一理论的任何其他数据点，如果可能的话，获取更高频率的数据，以便这一假设可以成立。
代码
这不是一个完美的 pythonic 解决方案，但是对于我的用例来说，它完成了工作。下面的片段使用了熊猫图书馆，关于它的文档可以在这里找到。
**# find diff between min and max wind direction**df['delta'] = df.max_winddirection - df.min_winddirection
**# check for if difference minimun < max - 180, if so, use different equation to get degree change
**
df.loc[df.min_winddirection < df.max_winddirection - 180, 'delta'] = 360 - df.max_winddirection + df.min_winddirection
解释
以 350 度到 5 度为例。
- 上面的代码首先为所有行分配基本的 delta 方程“最大角度-最小角度”。这将应用 350 作为该行的增量。
- 代码检查是否满足这些“穿过北方”情况的标准。这将检查是否:最小角度<最大角度—180”。在这种情况下，5 小于(350–180 ),因此为真。
- 该代码应用了一个新的公式，当δ为真时为*“360-最大角度+最小角度”。*导致 15 度。
如果不满足“穿过北方”标准，将毫无问题地应用基本最大值减去最小值公式。
希望这对你有所帮助。如果你有一个更干净的替代解决方案，我很乐意从中学习，请随时在评论中提出来。
用 Python 求矩阵的逆矩阵
原文：https://towardsdatascience.com/find-the-inverse-of-a-matrix-using-python-3aeb05b48308
如何在 Python 中使用行运算求矩阵的逆矩阵
介绍
本文跟随Python 中的高斯消去算法。介绍了一种利用行归约求逆矩阵的方法。
查看下面的文章，了解高斯消去法的必要介绍。
https://levelup.gitconnected.com/gaussian-elimination-algorithm-in-python-4e90cb3a0fd9 
照片由 Linus Mimietz 在 Unsplash 上拍摄
背景
矩阵代数基础|第二部分呈现逆矩阵。回想一下并非所有矩阵都是可逆的。 A 必须是平方 ( n×n )且有一个非零行列式。
本质上，矩阵的乘以它的逆得到单位矩阵**， I ，如等式 1 所示。**
等式 1 —计算矩阵的逆矩阵(图片由作者提供)
以等式 2 中的 3×3 矩阵 A 为例。
等式 2 —矩阵 A(图片由作者提供)
等式 3 等同于等式 1，变量替换为。
等式 3 —参数替换(图片由作者提供)
下面计算的 结果 是未知数 A⁻ 。
算法
从等式 3 的分量创建一个扩充的矩阵。这个新矩阵包含 A 串联的 列式与 I ，如等式 4 所示。
等式 4 —扩充矩阵(图片由作者提供)
通过对增广矩阵应用行运算获得逆矩阵。
对增广矩阵进行 高斯消去 型过程，得到 中的 A 降排梯队形式*(rref)*同时跃迁 I 为 A ⁻ 。
总而言之:
- 将 A 转换成rref。于是，A 变成了*单位矩阵。*****
a .转换到单位矩阵(图片由作者提供)
- 对I 运行上面的等价行运算，就变成了逆矩阵 A ⁻。
单位矩阵转换为逆矩阵(图片由作者提供)
例子
图 1 描述了逐步 操作改变扩充矩阵的前三列以实现 rref 的必要步骤。**
图 1 —实现简化的行梯队形式的行操作(图片由作者提供)
结果在意料之中。 A 成为单位矩阵，而 I 变换成以前未知的逆矩阵。**
Python 实现
有了 Python 中的编程的高斯消去算法，代码只需要微小的修改就可以得到的逆*。***
使用 Gist 1 将等式 2 中的 A 定义为 NumPy 数组。
要点 1 —使用 Numpy 定义 A
同样，实例化一个新的变量 I ，与 A 相同的方形*。***
要点 2 —定义单位矩阵
使用 NumPy 的列连接操作创建增广矩阵*，如要点 3 所示。***
要点 3 —创建身份矩阵
在不考虑某些边缘情况的情况下，下面在要点 4 中提供的代码是获得 逆 所必需的行操作的幼稚实现。**
要点 4 —在 Python 中寻找逆矩阵
与高斯消去算法相比，对代码的主要修改是代替的终止于行梯队形式，操作继续到达缩减行梯队形式。****
因此，不是只迭代枢轴下面的，而是枢轴上面的行也被遍历和操作*。***
执行脚本返回与图 1 中相同的答案。**
图 2-逆矩阵 Python 控制台输出
结果
给定任意数量的的可逆矩阵和任意大小的*，上述算法均适用。增加矩阵的大小也是可能的。***
Gist 5 提供了在 NumPy 中创建一个随机方阵的代码。**
要点 5 —大型随机矩阵
比较自定义算法与的运行时与 NumPy 当量突出了速度差。Gist 6 中的代码是一个简单的方法来记录时间*。***
要点 6 —比较方法运行时间
NumPy 比 T21 快一秒来反转矩阵。这个巨大的时差只会随着矩阵维度的扩大而增加。
结论
这篇文章概述了矩阵代数中使用的一个基本方法来计算一个矩阵的逆矩阵。**
使用概述的理论矩阵代数方法和等效的 Python 代码来理解操作如何工作。然而，Python 中的 NumPy 等库经过优化，可以高效地解密逆矩阵。在实践中，使用健壮的、维护良好的数学库。**
如果你对 Python、工程和数据科学感兴趣，可以看看我的其他文章。
***https://medium.com/@andrewdaviesul/membership
在 Gist 7 中找到所有 Python 代码。
要点 7——完成 Python 代码以找到矩阵的逆***
参考
[1] 工程师矩阵代数 — 杰弗里·查斯诺夫
求正定矩阵的最小拉伸方向
Python 中图像形成和相机校准综合教程系列的第 4 部分
在进入摄像机校准之前，让我们讨论一下正定矩阵和一般矩阵的一些性质。我们将在校准摄像机时使用这些属性。
特征值和特征向量
设𝐴是形状为𝑚×𝑛矩阵，𝑥⃗是大小为𝑛.的向量那么𝑥⃗就是𝐴的特征向量如果𝐴𝑥⃗ = 𝜆𝑥⃗，其中𝜆是一个标量叫做𝑥⃗.的特征值特征值和特征向量的性质取决于矩阵𝐴.
假设矩阵𝐴有𝑛独立特征向量，那么𝐴可以因式分解为:
其中，𝑆是以特征向量为列的特征向量矩阵，λ是特征值矩阵，它是沿对角线具有特征值的对角矩阵。注意，𝑆在这里是可逆的，这只有在列是独立的情况下才有可能，这意味着特征向量应该是独立的。
这里𝑥⃗1，𝑥⃗2，…，𝑥⃗𝑛是𝐴的特征向量；而𝜆1，𝜆2，…，𝜆𝑛是它们对应的特征值。
对称矩阵
𝐴是一个对称矩阵如果𝐴=𝐴⊺，其中𝐴⊺是𝐴.的转置我们将使用符号⊺来表示转置。对称矩阵将总是具有实特征值，并且特征向量可以被选择为正交的。
我们看到:
由于对称矩阵的特征向量可以选择正交，我们可以将特征向量矩阵 s 表示为𝑄，其中𝑄的列彼此正交。同样，𝑄的逆也是它的转置，因为它是正交的。因此，我们可以将𝐴分解为:
正定矩阵
如果在一个对称矩阵中所有的特征值都是正的，那么这个矩阵叫做正定矩阵。如果𝐴是正定矩阵，𝜆1，𝜆2，𝜆3…是𝐴的特征值，那么𝜆𝑖 > 0 和𝜆𝑖 ∈ 𝐑对于 i = 1，2，3，…
椭圆体
正定矩阵有一个有趣的性质:如果𝐴是正定矩阵，那么𝑥⃗⊺𝐴𝑥⃗ = 1 代表𝐑𝐧中以零为中心的椭球。𝐑𝐧是 n 维实空间，𝑥⃗ ∈ 𝐑𝐧.
这个椭球的主轴由矩阵𝐴.的特征值和特征向量决定
让我们看一个例子。
例子
设矩阵𝐴 = [[5，4]，[4，5]]。𝐴的特征值是 1 和 9，它们都是正的。因此 A 是正定矩阵。接下来，让我们计算𝑥⃗⊺𝐴𝑥⃗.这里是𝑥2].的𝑥⃗=[𝑥1
所以我们得到了𝑥⃗⊺𝐴𝑥⃗ = 5𝑥1 +8𝑥1𝑥2+5𝑥2。等等，那看起来不像椭圆。椭圆的方程具有以下形式:
那么是哪里出了问题呢？让我们用𝐴.的特征向量基来表示向量𝑥⃗本质上，我们需要在𝑥⃗上执行基变换的改变，并且我们在系列的第 2 部分中讨论了基变换的改变。基矩阵的变化是特征向量矩阵𝑄.的逆由于𝑄是一个标准正交矩阵，𝑄的逆是𝑄⊺，这是我们的基变换矩阵。
让𝑦⃗ = 𝑄⊺𝑥⃗.𝑦⃗是向量𝑥⃗ wrt 的特征向量基。现在让我们看看能否得到一个椭球体:
consider 𝑥⊺𝐴𝑥
= 𝑥⊺(𝑄Λ𝑄⊺)𝑥 (𝐴 = 𝑄Λ𝑄⊺)
= (𝑥⊺𝑄)Λ(𝑄⊺𝑥) (Matrix Multiplication is associative)
but, 𝑦 = 𝑄⊺𝑥 and 𝑦⊺ = 𝑥⊺𝑄
substituting y in the equation, we get:
𝑥⊺𝐴𝑥 = 𝑦⊺Λ𝑦
让我们看看𝑦⊺λ𝑦是什么样子的:
这里𝑦1，𝑦2，…，𝑦𝑛是矢量𝑦⃗.的分量或系数我们可以看到，𝑦⊺λ𝑦= 1 看起来像一个椭球的方程。我们可以将这个等式改写为:
这个方程明确地表示了一个椭球体，其中 1/√𝜆1，1/√𝜆2，… 1/√𝜆𝑛是主轴长度的一半。
回到我们的例子，𝐴 = [[5，4]，[4，5]]，𝐴的特征值是 1 和 9。将它们代入上述方程，我们得到 y1 + 9y2 = 1，这就是我们的椭圆方程。让我们画出这个:
+= 1
我们可以看到主轴的长度是 0.66 和 2，也就是 2 * 1/√𝜆1 和 2 * 1/√𝜆2.
接下来，让我们从特征向量基转换回标准的欧几里德基。转换到我们乘以𝑄⊺的特征向量基。为了从特征向量基转换回来，我们乘以𝑄⊺的逆，也就是𝑄.
𝑦⃗ = 𝑄⊺𝑥⃗ ⟹ 𝑥⃗ = 𝑄𝑦⃗
𝐴的特征向量的归一化形式为[-1/√2，1/√2]和[1/√2，1/√2]。所以𝑄的特征向量矩阵是:
请注意，这些列是正交的。我们将上面椭圆上的每一点乘以𝑄，然后像这样绘制变换后的点:
注意特征向量是如何沿着椭圆的主轴的。我们可以说拉伸的方向是沿着矩阵的特征向量，这些特征向量是椭圆的主轴。如果一个点沿着具有最大特征值的特征向量，它将被拉伸到最大，如果它沿着具有最小特征值的特征向量，它将被拉伸到最小。
正定矩阵的这个性质有一个重要的应用，我们将在相机校准中使用。
在此之前，有一个重要的定理你需要知道:如果𝐴是任意一个矩形矩阵，那么𝐴⊺𝐴就是一个正定矩阵。
设𝐴是一个形状为𝑚 x 𝑛的矩形矩阵，我们想找出|𝐴𝑥⃗|的最小值，其中𝑥⃗是任何大小为𝑛的向量，|𝐴𝑥⃗|代表𝐴𝑥⃗.的大小让我们看看如何计算:
We want to find the minimum value of |𝐴𝑥⃗|,
which is the same as finding the minimum value of |𝐴𝑥⃗|².
Now, |𝐴𝑥⃗|² can be written as (𝐴𝑥⃗)⊺(𝐴𝑥⃗),
which in turn can be written as 𝑥⃗⊺𝐴⊺𝐴𝑥⃗.
Now, 𝑥⃗⊺𝐴⊺𝐴𝑥⃗ = 𝑥⃗⊺(𝐴⊺𝐴)𝑥⃗ (associative property of Matrices),
But remember, 𝐴⊺𝐴 is a positive definite matrix.
⟹ 𝑥⃗⊺(𝐴⊺𝐴)𝑥⃗ = 1 represents an ellipsoid with the
eigenvectors of 𝐴⊺𝐴 along the principal axes.
So the direction of minimum stretching will be along
the eigenvector of 𝐴⊺𝐴 with the smallest eigenvalue.
因此，如果𝑥⃗沿着𝐴⊺𝐴的具有最小特征值的特征向量，那么|𝐴𝑥⃗|将是最小的。
结论
我希望你喜欢这篇文章。这篇文章的代码可以在这里找到。在的下一篇文章中，我们将执行摄像机校准。如果你有任何疑问或问题，请在下面的评论中告诉我。
参考
图像制作者名单
本文中的所有图片和数字，除非在标题中明确提及其来源，否则均由作者提供。
查找 ARIMA 模型的顺序
原文：https://towardsdatascience.com/find-the-order-of-arima-models-b4d99d474e7a
了解并找到时间序列基本建模的最佳参数
图片由 @m_____me
ARIMA 是开始单变量时间序列实验的最佳模型之一。它提供了最先进的性能，特别是在小数据集的情况下，深度学习模型尚未处于最佳状态。
这是一个简单的，可解释的模型，但 ARIMA 是一个参数模型。这意味着在拟合模型之前需要设置特定的参数。事实上，模型的自回归、移动平均和平稳部分分别需要参数 p、q 和 d
在本文中，您将了解这些参数是什么，以及如何使用 ACF 、 PACF 图和 BIC/AIC 信息标准找到最佳参数
了解 ARIMA 车型
首先，我们需要对这些模型有个基本的了解。
时间序列中“线性回归”或“经典模型”的等价物是 ar 和 MA。这些基本模型通常表现良好。它们是可解释的、轻量级的，是开始任何数据科学项目的良好基础。
ARMA 模型有一个严重的限制。它假设输入数据是稳定的(时间序列的均值和方差不依赖于时间序列被观察的时间)。如果您的数据是非平稳的，您可以使用 AR I MA 模型，其中 I 代表“集成的”。积分部分说明了“微分”，这是一种去除趋势并使时间序列平稳的方法。
1)自回归模型
订单 p 的 AR 是一个回归自身 p 过去值的模型。换句话说，系列 Vt 的当前值可以解释为 p 个过去值与随机误差的线性组合。我们假设下一个值与前 p 个值相关。用于预测下一个值的先前输入的数量被称为“顺序”。我们通常称之为“p”。这些用来计算下一个值的先验值被称为“滞后变量”。
作为滞后变量线性组合的 Vt 方程为:
AR§:Vt = ao+a1V { t-1 }+…+apV { t-p }+Nt。
当我们拟合 AR 模型时，我们在[0，p]中估计 I 的{ai}。 Nt 是一个白噪声过程(我们称之为“白噪声过程”是一系列均值为 0、方差恒定且误差不相关的随机值)和 ao 均值。
比如，当 p=3 时，为了预测下一个时间戳(和时间序列均值 a0)，我们代入最后三个值，分别乘以系数 a1，a2，a3。
您可以看到如何使用下面的代码片段模拟订单 p 的 AR 模型:
该模型将以下内容作为输入:
- 随机变量的实现列表 V
rand_v_list = np.random.normal(size=N)
- p 系数列表{ao，…，ap}。
模拟的 AR 模型是将 V 的 p 个连续值乘以其 p 个系数来创建一个时间序列。
我们的目标是:给定一个模拟的 AR 时间序列，如何找到 p，用来生成它的参数？
2)移动平均模型
移动平均模型利用观测值和应用于滞后观测值的移动平均的残差之间的相关性。换句话说，马将预测值建模为到 q 步的过去误差项的线性组合。误差或残差定义为 res =预测值-真值。
它利用回归方法中以前的预测误差来预测未来的观测值。在这种情况下，每个未来观测值都是以前预测误差的加权平均值。
MA(q):Et = B0+B1 * Et-1+…+BP * Et-q+Nt
我们估计{bi}在[0，q]中，Nt 是白噪声(均值为 0 且方差有限的不相关随机变量的集合)。
如果残差和预测值之间仍然存在相关性，这意味着我们仍然可以从误差中获得信息。实际上，我们可以认为误差是随机的。MA 过程只是来自白噪声过程的值 E 的加权组合。
您可以通过下面的代码片段看到该模型是如何工作的:
该模型采用与 ar 流程相似的输入。
3) ARMA
自回归移动平均 ARMA (p，q)结合了 AR §和 MA (q)过程，考虑了时间序列的滞后或过去值以及误差项。它是平稳时间序列最有效的线性模型[1]。
4) ARIMA
当输入数据不稳定时，用一个自回归综合移动平均代替 ARMA 模型。
在 ARIMA，积分部分将时间序列“固定化”[2]。当积分的阶数为 d=0 时，ARIMA 的行为类似于 ARMA 模型。当 d=1 时，模型将从 t 处观察到的值中减去 t-1 值。当数据中存在长期趋势时，您可以添加微分。术语综合指的是建模的时间序列必须不同多少次才能产生平稳性。
当你使用 ARIMA 模型时，不要手动进行微分，你可以直接在非平稳时间序列上拟合 ARIMA 模型。
图片由 @tezos
分析
现在我们了解了这些过程是如何工作的，以及它们的顺序意味着什么，让我们看看如何在实践中找到参数。为了分别找到模型的 AR、integrated 和 MA 部分的参数 p、d 和 q，我们可以使用:
- ACF 和 PACF 剧情。
- 领域知识
- 各种适合的度量标准 (BIC，AIC)
我们提醒 p 是模型中滞后观测值的数量(也称为滞后阶数)，d 是原始观测值的差分次数，q 是移动平均窗口的大小。
我们将在下一节介绍如何使用 ACF、PACF 图、BIC 和 AIC 标准。
自相关函数
时间序列过程的一个重要方面是自相关。自相关是一种统计属性，当一个时间序列与其自身的先前版本或滞后版本线性相关时，就会出现自相关。它还用于检测时间序列中可能的季节性。
我们使用自相关函数来评估时间序列中的依赖程度，并选择适当的模型(MA、AR 或 ARIMA)。
自相关函数( ACF)和偏自相关函数 (PACF)可以针对任何时间序列(不仅仅是平稳的)进行计算。
在实践中，我们使用这两个图的组合来确定 ARMA 过程的顺序。
我如何知道我的时间序列遵循哪个过程？
通用流程:
- ARIMA(0，0，0)是一个白噪声模型。(Vt=Nt)
- ARIMA(0，1，0)是随机游走(Vt-V{t-1}=c+Nt，c 是均值)
- ARIMA(0，1，1)指数平滑法
但是除了常见的模式之外，我们使用 PACF 和 ACF 图来可视化要识别的模式:
过程识别表。作者图片
在下图中，我们展示了阶为 1 的 AR、MA 或 ARMA 过程中 ACF 和 PACF 的行为。如果你的时间序列是 AR 或 MA 或 ARMA，你可以很容易地使用这个表来识别任何订单。在 AR 或 MA 过程的情况下，您可以计算图表中突然下降的峰值数量。在 ACF 和 PACF 中，0 阶滞后总是等于 1(该值与其自身完全相关)，因此我们在第一个峰值后开始计数。
1) ACF
我们使用自相关函数 ACF 来了解时间序列与其滞后值的相关程度。
我们使用 ACF 来标识 ARMA 模型的 MA(q) 部分。
事实上，我们确定哪个滞后变量与产出有很强的相关性。这些滞后变量对时间序列的下一个值有很强的预测能力，应该用于预测。
在实践中，我们计算与输出值有显著自相关的滞后数，以找到 MA 模型的 q 值。
让我们针对 500 个模拟随机观察值，模拟顺序为 q=1、q=3、q=5 和 q=7 的四个 MA 过程。
绘制这四个不同的模拟数据:
顺序递增的不同 MA 过程。作者图片
现在，让我们用 ACF 和 PACF 图绘制 MA(q=5)模拟数据。
在右图中，我们可以看到 PACF 图上的阻尼下降(正的和负的)。即使在滞后 10 年之后，滞后仍然很明显。而在左侧的 ACF 图上，在 1+5 个显著滞后之后，所有随后的滞后都不显著。我们可以看到一个并购过程。
为了找到 q，我们可以计算 ACF 图上显著滞后的数量(蓝带之外)。除了第一个滞后对应于值与其自身的相关性之外，我们可以看到与序列的接下来 5 个滞后的相关性。
ACF 和 PACF 用于 MA(q=5)。我们可以在 ACF 中读取 5 个显著或“高”峰，左图。作者图片
2) PACF
直觉
偏自相关函数(PACF) 表示在我们考虑一些其他变量集合的值的假设下，两个变量之间的相关性。在回归中，这种部分相关性可以通过关联两个不同回归的残差来找到。
我们使用 PACF 来寻找残差与下一个滞后的相关性。如果我们发现残差和下一个滞后中隐藏的信息具有合理的相关性，我们将使用下一个滞后作为建模的特征。
换句话说，目标时间戳和先前时间戳之间的自相关值由两个值之间的直接相关性加上间接相关性组成。我们认为这些间接相关性是观察值相关性的线性函数。
我们可以使用 PACF 来计算 AR 流程中显著滞后的数量。
实验
首先让我们模拟 p=1，p=3，p=5，p=7 的移动平均 AR。
订单可变的不同 AR 流程。作者图片
我们可以在上面的图中看到，随着 p 的增加，我们也增加了时间序列的振荡性质。当时间序列增长过快时，自回归分量会产生一种自然趋势，使其回到平均值。现在它在平均值 0 附近振荡，但也可能是任何值。因为时间序列既取负值也取正值，所以 ACF 和 PACF 图也取正值和负值。
现在我们可以画出这些 AR 过程的 ACF 和 PACF 函数，并读出显著滞后的数目作为 p 的值
ACF 和 PACF 用于 AR(p=7)。我们可以在右边的 PACF 图上看到七个重要的峰值。图片由作者提供。
3) AIC/BIC 标准
绘制 ACF/PACF 对于识别 AR 和 MA 过程是有效的。但是对于 ARIMA 进程，更常见的是使用 auto_arima 函数。Auto arima 是一种强力方法，它尝试不同的 p 和 q 值，同时最小化两个标准:AIC 和 BIC。
评估正则化拟合优度的最常见度量是:
- 贝叶斯信息准则(BIC)
- 赤池信息标准(AIC)。
这些度量提供了考虑模型复杂性的模型性能的度量。AIC 和 BIC 将一个反映模型与数据拟合程度的术语与一个根据参数数量按比例惩罚模型的术语结合起来[3]。
作为一种正则化技术，我们希望根据模型中参数的数量进行惩罚。事实上，p 和 q 越大，你用来预测下一个值的延迟就越多，你就越有可能过度拟合你的数据。
自动 ARIMA 过程旨在确定一个
**ARIMA**模型的最佳参数，确定一个单一的拟合 ARIMA 模型。[…]为了找到最佳模型，自动 ARIMA 为给定的
**information_criterion**(分别为‘AIC’，‘aicc’，‘BIC’，‘hqic’，‘OOB’)(Akaike 信息标准、修正的 aka ike 信息标准、贝叶斯信息标准、Hannan-Quinn 信息标准或“出袋”–用于验证评分)之一进行优化，并返回使值最小化的 ARIMA。
在实践中，我们使用现成的工具 auto_arima 从包 pmdarima 中自动找到时间序列的顺序。
让我们尝试 auto_arima 来找出我们模拟的几个 ma 过程的顺序:
q=1
Performing stepwise search to minimize aic
ARIMA(0,0,0)(0,0,0)[0] intercept : AIC=1605.796, Time=0.01 sec
ARIMA(1,0,0)(0,0,0)[0] intercept : AIC=1493.552, Time=0.03 sec
ARIMA(0,0,1)(0,0,0)[0] intercept : AIC=1461.981, Time=0.03 sec
ARIMA(0,0,0)(0,0,0)[0] : AIC=1604.553, Time=0.01 sec
ARIMA(1,0,1)(0,0,0)[0] intercept : AIC=1463.723, Time=0.05 sec
ARIMA(0,0,2)(0,0,0)[0] intercept : AIC=1463.755, Time=0.05 sec
ARIMA(1,0,2)(0,0,0)[0] intercept : AIC=1465.600, Time=0.13 sec
ARIMA(0,0,1)(0,0,0)[0] : AIC=1460.398, Time=0.02 sec
ARIMA(1,0,1)(0,0,0)[0] : AIC=1462.121, Time=0.03 sec
ARIMA(0,0,2)(0,0,0)[0] : AIC=1462.155, Time=0.02 sec
ARIMA(1,0,0)(0,0,0)[0] : AIC=1491.861, Time=0.01 sec
ARIMA(1,0,2)(0,0,0)[0] : AIC=1463.988, Time=0.06 sec
Best model: ARIMA(0,0,1)(0,0,0)[0]
Total fit time: 0.468 seconds
Optimal order for is: (0, 0, 1)
q=3
Performing stepwise search to minimize aic
ARIMA(0,0,0)(0,0,0)[0] intercept : AIC=1702.731, Time=0.01 sec
ARIMA(1,0,0)(0,0,0)[0] intercept : AIC=1570.816, Time=0.03 sec
ARIMA(0,0,1)(0,0,0)[0] intercept : AIC=1628.147, Time=0.04 sec
ARIMA(0,0,0)(0,0,0)[0] : AIC=1701.862, Time=0.01 sec
ARIMA(2,0,0)(0,0,0)[0] intercept : AIC=1528.848, Time=0.04 sec
ARIMA(3,0,0)(0,0,0)[0] intercept : AIC=1519.618, Time=0.06 sec
ARIMA(4,0,0)(0,0,0)[0] intercept : AIC=1485.096, Time=0.06 sec
ARIMA(4,0,1)(0,0,0)[0] intercept : AIC=1484.876, Time=0.11 sec
ARIMA(3,0,1)(0,0,0)[0] intercept : AIC=1509.277, Time=0.13 sec
ARIMA(4,0,2)(0,0,0)[0] intercept : AIC=1464.510, Time=0.17 sec
ARIMA(3,0,2)(0,0,0)[0] intercept : AIC=1465.074, Time=0.15 sec
ARIMA(4,0,3)(0,0,0)[0] intercept : AIC=1465.187, Time=0.28 sec
ARIMA(3,0,3)(0,0,0)[0] intercept : AIC=1464.135, Time=0.20 sec
ARIMA(2,0,3)(0,0,0)[0] intercept : AIC=1462.726, Time=0.23 sec
ARIMA(1,0,3)(0,0,0)[0] intercept : AIC=1462.045, Time=0.17 sec
ARIMA(0,0,3)(0,0,0)[0] intercept : AIC=1460.299, Time=0.09 sec
ARIMA(0,0,2)(0,0,0)[0] intercept : AIC=1507.915, Time=0.06 sec
ARIMA(0,0,4)(0,0,0)[0] intercept : AIC=1462.121, Time=0.09 sec
ARIMA(1,0,2)(0,0,0)[0] intercept : AIC=1467.963, Time=0.08 sec
ARIMA(1,0,4)(0,0,0)[0] intercept : AIC=1463.941, Time=0.23 sec
ARIMA(0,0,3)(0,0,0)[0] : AIC=1458.689, Time=0.12 sec
ARIMA(0,0,2)(0,0,0)[0] : AIC=1506.487, Time=0.03 sec
ARIMA(1,0,3)(0,0,0)[0] : AIC=1460.415, Time=0.11 sec
ARIMA(0,0,4)(0,0,0)[0] : AIC=1460.498, Time=0.04 sec
ARIMA(1,0,2)(0,0,0)[0] : AIC=1466.278, Time=0.07 sec
ARIMA(1,0,4)(0,0,0)[0] : AIC=1462.305, Time=0.11 sec
Best model: ARIMA(0,0,3)(0,0,0)[0]
Total fit time: 2.717 seconds
Optimal order for is: (0, 0, 3)
q=5
Performing stepwise search to minimize aic
ARIMA(0,0,0)(0,0,0)[0] intercept : AIC=1659.497, Time=0.01 sec
ARIMA(1,0,0)(0,0,0)[0] intercept : AIC=1570.804, Time=0.03 sec
ARIMA(0,0,1)(0,0,0)[0] intercept : AIC=1613.884, Time=0.03 sec
ARIMA(0,0,0)(0,0,0)[0] : AIC=1658.949, Time=0.01 sec
ARIMA(2,0,0)(0,0,0)[0] intercept : AIC=1495.855, Time=0.04 sec
ARIMA(3,0,0)(0,0,0)[0] intercept : AIC=1482.804, Time=0.05 sec
ARIMA(4,0,0)(0,0,0)[0] intercept : AIC=1484.509, Time=0.07 sec
ARIMA(3,0,1)(0,0,0)[0] intercept : AIC=1484.564, Time=0.11 sec
ARIMA(2,0,1)(0,0,0)[0] intercept : AIC=1484.926, Time=0.07 sec
ARIMA(4,0,1)(0,0,0)[0] intercept : AIC=1486.509, Time=0.15 sec
ARIMA(3,0,0)(0,0,0)[0] : AIC=1481.204, Time=0.03 sec
ARIMA(2,0,0)(0,0,0)[0] : AIC=1494.160, Time=0.02 sec
ARIMA(4,0,0)(0,0,0)[0] : AIC=1482.892, Time=0.03 sec
ARIMA(3,0,1)(0,0,0)[0] : AIC=1482.953, Time=0.05 sec
ARIMA(2,0,1)(0,0,0)[0] : AIC=1483.270, Time=0.03 sec
ARIMA(4,0,1)(0,0,0)[0] : AIC=1484.892, Time=0.08 sec
Best model: ARIMA(3,0,0)(0,0,0)[0]
Total fit time: 0.824 seconds
Optimal order for is: (3, 0, 0)
q=7
Performing stepwise search to minimize aic
ARIMA(0,0,0)(0,0,0)[0] intercept : AIC=2171.867, Time=0.01 sec
ARIMA(1,0,0)(0,0,0)[0] intercept : AIC=1789.289, Time=0.03 sec
ARIMA(0,0,1)(0,0,0)[0] intercept : AIC=1931.174, Time=0.04 sec
ARIMA(0,0,0)(0,0,0)[0] : AIC=2172.420, Time=0.01 sec
ARIMA(2,0,0)(0,0,0)[0] intercept : AIC=1788.083, Time=0.04 sec
ARIMA(3,0,0)(0,0,0)[0] intercept : AIC=1779.499, Time=0.07 sec
ARIMA(4,0,0)(0,0,0)[0] intercept : AIC=1778.438, Time=0.07 sec
ARIMA(4,0,1)(0,0,0)[0] intercept : AIC=1773.792, Time=0.26 sec
ARIMA(3,0,1)(0,0,0)[0] intercept : AIC=1780.497, Time=0.10 sec
ARIMA(4,0,2)(0,0,0)[0] intercept : AIC=1695.057, Time=0.32 sec
ARIMA(3,0,2)(0,0,0)[0] intercept : AIC=1738.073, Time=0.35 sec
ARIMA(4,0,3)(0,0,0)[0] intercept : AIC=1691.378, Time=0.51 sec
ARIMA(3,0,3)(0,0,0)[0] intercept : AIC=1711.992, Time=0.25 sec
ARIMA(4,0,4)(0,0,0)[0] intercept : AIC=1694.119, Time=0.72 sec
ARIMA(3,0,4)(0,0,0)[0] intercept : AIC=1701.593, Time=0.27 sec
ARIMA(4,0,3)(0,0,0)[0] : AIC=1689.749, Time=0.27 sec
ARIMA(3,0,3)(0,0,0)[0] : AIC=1710.347, Time=0.19 sec
ARIMA(4,0,2)(0,0,0)[0] : AIC=1693.396, Time=0.15 sec
ARIMA(4,0,4)(0,0,0)[0] : AIC=1692.698, Time=0.49 sec
ARIMA(3,0,2)(0,0,0)[0] : AIC=1736.557, Time=0.16 sec
ARIMA(3,0,4)(0,0,0)[0] : AIC=1699.989, Time=0.16 sec
Best model: ARIMA(4,0,3)(0,0,0)[0]
Total fit time: 4.481 seconds
Optimal order for is: (4, 0, 3)
我们对订单 1、3、5 和 7 的 ma 过程使用 auto arima。对于小订单 q=1 和 q=3，Auto_arima 准确地识别 ma 过程及其订单，但是对于订单 q=5 和 q=7，它混合了 AR 和 MA。
结论
当您开始时间序列分析时，从可能满足用例需求的简单模型开始是一个好的实践。ARIMA 模型简单而透明，你可以推导出严格的统计特性。它们在小数据集上表现良好，构建和重新训练都很便宜。
如果您需要使用它们，您需要了解它们是如何工作的，并明确地设置它们的参数。本文为您提供了充分自信地调整模型顺序的技术。
本文笔记本可用 此处 。
参考资料:
[1]使用 Python 进行时间序列分析:从基础到前沿技术https://www . Amazon . com/Hands-Time-Analysis-Python-Techniques/DP/1484259912(非附属)
[2]https://people . cs . Pitt . edu/~ Milos/courses/cs 3750/lectures/class 16 . pdf
[3][https://www . science direct . com/topics/psychology/Bayesian-information-criterion #:~:text = The % 20 a kaike % 20 information % 20 criterion % 20(AIC，to % 20 its % 20 number % 20 of % 20 parameters。](https://www.sciencedirect.com/topics/psychology/bayesian-information-criterion#:~:text=The%20Akaike%20information%20criterion%20(AIC,to%20its%20number%20of%20parameters.)
查找数据库中排名靠前且最慢的查询
原文：https://towardsdatascience.com/find-the-top-n-most-expensive-queries-48e46d8e9752
找到降低数据库处理速度的瓶颈查询
在这篇文章中，我们将寻找慢速查询(图片由 Nicolai Dürbaum 在 unsplash 上提供)
当您的数据库增长时，很容易失去对它正在执行的所有进程的跟踪。跟踪阻碍操作的缓慢查询不是很好吗？
在本文中，我们将创建一个查询，为您提供分析和优化数据库所需的所有工具:它选择有问题的查询，提供这些查询的相关信息，并为您提供查找和改进它们的方法。阅读本文后，您将能够揭示前 n 个查询以及每个查询的信息:
- 在哪个数据库和对象中查找查询(例如在
Coffee数据库中的FreshPot存储过程中) - 实际的查询本身
- 查询计划
- 创建和最后执行的时间
- 执行次数
- 返回的行数的总数、最小值、最大值和平均值、运行时间和 CPU 时间
注意，本文分析了 SQL Server 数据库，但是类似的功能也存在于 Postgres 数据库中。查看本文中关于 Postgres 数据库中所有查询的统计数据。
搜寻慢速查询
我们将编写一个查询，提供哪些查询速度慢的信息，并为我们提供分析查询的工具。那我们可以
- 检测:哪些查询有问题？
- 丰富:增加更多信息，让我们了解缓慢的原因
- 分析:在哪里可以找到查询，问题出在哪里？
- 优化:我们如何进一步改进我们的查询和数据库？
完整查询可在本文底部或这里 。注意，我们将要使用的一些视图需要VIEW DATABASE STATE的许可。
1.查找有问题的查询
在第一部分中，我们使用 CTE 来选择我们感兴趣的查询。我们想从sys.dm_exec_query_stats获取一些数据；这将跟踪 SQL Server 中缓存查询计划的性能统计信息。
这里的目标是选择我们感兴趣的记录。在本例中，我们选择的记录是新的(自上次执行后不到 30 天)和频繁使用的**(在过去 30 天内超过 100 次)。我们不介意很少执行的慢速查询。接下来，我们按照平均 CPU 时间对查询进行排序，然后返回前 10 个最慢的查询。**
当然，我们可以根据许多标准进行筛选:
- 执行时间变化很大的记录(例如，最小和最大 CPU 时间相差很大)
- 运行时间和 CPU 时间(CPU 实际执行的时间)相差很大的记录；可能查询浪费了很多等待的时间)
- 按返回的 avg_rows 或总 cpu 时间排序。
有些虫子伪装得很好(图片由 id23 在 Unsplash 上拍摄)
2.提供有关查询的更多信息
在上一部分中，我们选择了频繁执行且执行时间很长的查询；这些是我们的问题查询。这一部分的目标是尽可能多地收集关于这些查询的相关信息。让我们看一下代码
上面的代码很长，但它只做了几件事:
- 使用第 34 行的视图
sys.dm_exec_sql_text添加实际的 SQL 语句(并清除第 21 到 31 行的语句) - 使用视图
sys.dm_exec_query_plan添加查询计划(第 35 行) - 使用视图
sys.dm_exec_plan_attributes获取我们稍后需要的数据库 id 和对象 id。在第 37 到 41 行中，我们将 dbid 和 objectid 记录从行转换为列，这样我们可以更容易地在第 36 行交叉应用它们。
这是输出(带有匿名声明文本):
我们的数据库分析语句的输出(不是所有的列都适合，但最重要的列在这里，图片由作者提供)
3.分析我们的问题查询
我们现在可以开始分析查询；以下是一些建议:
- 检查 cpu 时间和运行时间之间的差异；可能查询等了很多？
- 检查最小和最大 cpu 和执行时间之间的差异；也许将某些作业安排在晚上服务器有更多可用内存的时候？
- 分析语句文本(格式化 vscode 中的 sql
cntrl-shift-p→format document with→sql - 单击并签出 query_plan
一旦您发现了问题或改进查询的机会，我们需要找到查询的位置。我们已经有了数据库名称；让我们找到查询所在的对象的名称(例如，在存储过程中)。我们可以通过在下面的查询中插入db_name和objectid 来做到这一点:
SELECT * FROM CoffeeDB.sys.objects WHERE object_id = 808389949
这将告诉您我们正在寻找的对象的类型和名称。例如:
现在我们可以在 CoffeeDB 中搜索一个名为‘My _ Stored _ Procedure’的存储过程，并加快查询速度！
4.进一步优化
随着所有拥塞的消失，我们的数据库进程再次自由流动(图片由 Alexander Schimmeck 在 unsplash 上提供)
结论
通过这篇文章，您将完全有能力优化您的数据库！我们已经了解了代码的工作原理以及如何使用它，现在是时候将它应用到您的情况中了。查看 完整查询此处 。
如果你有建议/澄清，请评论，以便我可以改进这篇文章。同时，看看我的其他关于各种编程相关主题的文章，比如:
- 了解索引如何加快查询速度
- 跟踪 Postgres 中所有查询的统计数据以防止瓶颈
- Docker 中的 Postgres 入门
- 删除进入另一个表格
- 更新到另一个标签页 le
- 在一条语句中插入、删除和更新
- 更新选择一批记录
- 插入唯一表格
编码快乐！
—迈克
页（page 的缩写）学生:比如我正在做的事情？跟我来！
**https://mikehuls.medium.com/membership **
发现数学背后的美
原文：https://towardsdatascience.com/finding-beauty-behind-all-the-math-844c8e39472d
作者聚焦
Hannah Roos 如何利用她的数据科学、商业和心理学背景帮助人们解决数据问题
在 Author Spotlight 系列中，TDS 编辑与我们社区的成员谈论他们在数据科学领域的职业道路、他们的写作以及他们的灵感来源。今天，我们很高兴与 汉娜·鲁斯 分享我们的对话。
Hannah Roos是一位热衷于数据科学的心理学家。受她在管理咨询和行为科学领域工作的启发，她撰写了动手编程教程，并深入研究了 R 和 Python。她被积极弥合循证研究和商业世界之间的鸿沟的使命所驱使。
你进入数据科学的切入点是什么？
在我的心理学学士学位期间，我们被告知对研究问题采取合理的方法的重要性，特别是在“软”科学中，这让我对大多数学生最害怕的部分感兴趣:统计和假设检验。
我意识到，我们的大脑可能无法直观地理解这些事情，但这实际上激发了我的雄心:我想了解如何通过数据解决问题。
你如何处理这个领域中更难的话题？
我喜欢数据科学具有挑战性的方面，因为一旦你经过长时间的深思熟虑和努力工作后掌握了一项任务，你会有一种巨大的成就感，你最终会看到所有数学背后的美。
在我的研究结束时，我担任视觉认知实验的研究助理。我所有的同事似乎都更擅长编程和复杂的数据分析:这可能是我旅程中特别具有挑战性的部分。似乎他们可以进入一个充满逻辑思想的平行世界，而且他们有合适的词汇来辩论它们。但在那个阶段，这个世界不知何故对我关闭了。也许我根本就不是为这个而生的？
但是我不再关注我是否“足够聪明”，而是继续学习变得越来越熟练。到目前为止，我已经完成了 60 多门关于 Datacamp 的课程，扩展了我关于 R、Python 和 SQL 的知识，以及统计学和机器学习中的常见概念。
读你的文章，你会感觉到你的心理学背景经常影响你的工作。
是啊，那是真的！作为一名心理学家，我总是试图通过数据来找出是什么让人们滴答作响，这不是一项简单的任务。
我对发生在心理学和商业世界之间的实际问题有一种感觉，这些真的激发了我的写作。例如，我最近写了一篇关于使用 AI 预测公司员工流动率的文章。当然，你可以建立一个机器学习算法，将员工分为 1)面临离职风险的员工和 0)可能留下的员工。但老实说，影响一个人在特定公司工作的动机有很多方面，这使得很难超越样本内预测进行概括——例如，与你的主管的关系、市场动态、学习机会……所有这些都可能在相对较短的时间内发生变化。那么，你如何在实践中准确预测如此复杂的决策呢？
虽然我相信在现实世界中使用数据来推动决策，但我认为我们也需要使用我们的常识和领域知识。这就是我们如何自信地说数据能告诉我们关于未来的什么，以及它不能告诉我们什么。
通过数据工作，你最想解决哪类问题？
我喜欢专注于数据的项目，这些项目完美地满足了我作为一名科学家以及广大公众的好奇心。总的来说，我想更深入地研究一下心理测试，也就是说，使用统计学来改进测量我们无法直接测量的特征(例如，性格)的方法。这有可能告诉你一些关于人们如何在一个非常深刻的层面上不同的事情。
目前，我正在分析 2020 年 COVID 爆发早期的心理影响数据——这就是为什么我的下一篇文章将包括时间序列分析和 Python 编码示例的方面。从更广泛的角度来看，我对跨时间的分析感兴趣，以调查变化和发展。
是什么促使你开始公开写你的工作和其他数据相关的话题？
我开始写关于 p 值和进行假设检验的真正含义。这是因为我意识到我已经错了很长一段时间，并且感觉学生和专业人士对这个概念仍然有很多误解。我想分享“啊哈！”我终于和其他同样在挣扎的人有了这样的时刻:如果我能写一篇文章，让人们更容易理解，会怎么样？
我喜欢写与数据相关的话题，因为从事我热爱的事情令人惊叹，但我也想生活在一个可以免费分享知识的世界里。希望我能以这种方式为这样一个世界做出贡献。
作为一名作家，看起来你已经有意识地选择了追求彻底性和深度，而不是产量。你的文章背后的流程是什么样子的？
我想对一个问题提供一个完整和整体的观点，而不是匆忙地通过编程部分。通常，我会从一个我想更深入理解的概念开始——比如说贝叶斯统计或结构方程建模。然后，我开始阅读一些科学论文，并寻找一个公共数据集，我可以用它来进行并行实验。
一旦我觉得我已经从我的学习项目中获得了一些关键的见解，一个好的用例，以及一个其他人也可以从中受益的有效的编码方法，我实际上就开始写作了。为了使它更容易，我从一个包含我想提出的关键点的普通 word 文件上的项目符号开始，然后一步一步地用一个连贯的文本替换它。
你有什么实用的建议给那些想把自己的工作写下来的数据科学家吗？
对于初学者，我会说你不需要用你的文章证明什么。也许只是探索一个有趣的话题，并写下你的见解。想象一下，你会把这些告诉你的朋友，他实际上并不太关心数据科学，但仍然礼貌地听你说。问问你自己:我能激发一点兴趣吗？如果是这样，你可能有一个好故事要讲。如果没有，你至少有你自己的发现之旅。
数据科学是一个不断变化的领域，你希望在未来几个月或几年看到什么特别的变化吗？
我希望看到数据科学变得更加包容，让每个人都可以使用，而不是询问来自行业和学术界的专家，我认为每个人都应该有机会成为这样的专家，并从数据中获得自己的见解。
希望像“人工智能”这样的术语不会引发一种不祥的黑匣子的想法，这种黑匣子利用我们的个人信息来赚取更多的利润——在理想的未来，来自更广泛公众的人们最终会理解算法如何工作，以评估它们在多大程度上可能是有益的。
对于更广阔的技术世界，我希望看到单调乏味的工作流程实现更多自动化，这样人们就可以更多地关注对他们真正重要的事情。我相信，一些科技初创公司已经有了与数据相关的伟大创新，能够彻底改变我们的生活方式，让我们成为更加自主的个体——如果以深思熟虑的方式设计和使用的话。
您可以通过媒体跟踪 Hannah，了解她的最新作品；为了体验 Hannah 过去的 TDS 深潜和教程，这里有一些亮点:
- 用《R》(2022 年 1 月 29 分钟)中的拉万揭开人格的潜在结构
- 智能处理 R 中的缺失数据(2021 年 2 月 15 分钟)
- 数据评书为 R 中的调查数据(2020 年 12 月 8 分钟)
- 比较好做的贝叶斯？(2021 年 4 月 30 分)
想和广大观众分享一些你自己的作品吗？我们很乐意收到你的来信。
这个 Q & A 是为了长度和清晰而稍加编辑的。
在图像中查找聚类
原文：https://towardsdatascience.com/finding-clusters-in-an-image-96fe5eefc784
不仅仅是在表格数据中查找聚类
Alex Shuper 在 Unsplash 上拍摄的照片
我们也很容易在数据的行和列中找到聚类。但是如何在图像中找到集群呢？让我用最近在卡塔尔举办的 2022 年世界杯的一张图片来说明这个问题。
这里展示的是我在巴西对塞尔维亚的比赛中拍摄的照片。这张照片是在理查森的惊人进球之前拍摄的。
自己拍的 2022 世界杯足球赛巴西 vs 塞尔维亚的照片(图片由作者提供)
你能在照片中观察到任何集群吗？视觉上，您可以找到一群人，如下所示。
2022 年世界杯足球赛巴西对塞尔维亚的照片，显示了集群，由本人拍摄(图片由作者提供)
现在让我告诉你如何自动识别集群的技术步骤。
1.检测图像中的对象
第一步是检测图像中的所有对象。检测什么对象将基于聚类的目的。这里我们的目标是检测人群。检测照片中的人物可以在 YOLOv5 上完成。使用 YOLOv5 进行物体检测的结果如下所示。
用 YOLOv5 进行物体检测(图片由作者提供)
YOLOv5 还提供了额外的统计数据，如对象检测的数量。在这张照片中，有 17 个物体(人)。下面是一段代码
import yolov5
from sklearn.cluster import DBSCAN
from sklearn.preprocessing import StandardScaler
import numpy as np
import pandas as pd
# load pretrained model
model = yolov5.load('yolov5s.pt')
# set model parameters
model.conf = 0.25 # NMS confidence threshold
model.iou = 0.45 # NMS IoU threshold
model.agnostic = False # NMS class-agnostic
model.multi_label = False # NMS multiple labels per box
model.max_det = 1000 # maximum number of detections per image
# set image
img = 'brazil_vs_serbian.png'
# inference with test time augmentation
results = model(img, augment=True)
# parse results
predictions = results.pred[0]
boxes = predictions[:, :4] # x1, y1, x2, y2
scores = predictions[:, 4]
categories = predictions[:, 5]
2.识别物体的 X 和 Y 位置
下一步是识别物体的 X 和 Y 位置。YOLOv5 结果存储图像中每个对象的位置。位置存储为 xmin、ymin、xmax、ymax。我们可以取 xmin 和 xmax，以及 ymin 和 ymax 的平均值，求出物体的 x，y 位置。
对象的坐标(图片由作者提供)
可以将坐标转换为数据框。有一个调整很重要，它与反转 Y 位置有关。通常，图像以 Y 轴反转的方式存储。这意味着 Y=0 在图像的顶部。在计算过程中，我们可以通过减去图像的高度(这里是 2253)来还原位置。高度也可以自动确定，但是为了简化代码片段，我直接输入了高度。
df = results.pandas().xyxy[0]
df['x'] = (df['xmin'] + df['xmax'])/2.0
df['y'] = (2253-(df['ymin'] + df['ymax'])/2.0) #2253 is image height
结果在这里显示为数据帧。
位置为 x，y 的数据框(图片由作者提供)
3.使用基于密度的聚类来查找聚类。
现在我们已经准备好寻找星团了。找到图像中的聚类意味着找到与图像的其他部分相比有更多对象的区域。本质上，这意味着根据图像中对象的密度来寻找聚类。一种有效的算法是基于密度的聚类(DBSCAN)。下面显示的是散点图中绘制的结果。
DSBSCAN 结果(图片由作者提供)
DBSCAN 算法将绿点识别为密集区域的一部分。下面是 DBSCAN 算法的代码片段。
#Density clustering
X = StandardScaler().fit_transform(df[['x','y']])
clustering = DBSCAN(eps=0.5, min_samples=5).fit(X)
df['c'] = clustering.labels_
4.分析结果
现在让我们对比照片分析散点图上的结果。
分析结果(图片由作者提供)
瞧啊。散点图看起来像是足球比赛的可视化数字表示！
结论
在图像中寻找聚类是一项有趣的技术，也有许多用例。在博客中，您看到了体育分析的一个用例。此外，还有各种其他用例，如识别公共空间中的过度拥挤或识别购物者在购物中心花费更多时间的区域。
希望你喜欢这个博客！
请用我的推荐链接加入 Medium 。
https://pranay-dave9.medium.com/membership
请订阅每当我发布一个新的故事时，请及时通知我。
https://pranay-dave9.medium.com/subscribe
额外资源
网站(全球资讯网的主机站)
你可以访问我的网站，这是一个学习数据科学的无代码平台。【https://experiencedatascience.com】T5T6
Youtube 频道
这是我的 YouTube 频道
https://www.youtube.com/c/DataScienceDemonstrated的链接
在 Python 中查找两个纬度和经度之间的距离
地理数据
在 Python 中查找两个纬度和经度之间的距离
地理坐标的特征工程距离
为模型准备数据时，有时查找两个位置之间的距离可能会很有用。这篇文章展示了如何在 Python 中从两个地点的纬度和经度找到两个地点之间的最短球面距离和旅行距离。
照片由 Unsplash 上的 Dariusz Sankowski 拍摄
🌐地理坐标
我们可以根据地理坐标定位地球上的任何地方。一个位置的地理坐标由它的纬度和经度位置组成。
作者图片
📍纬度
纬度是北极和南极之间垂直位置的量度。假想的水平纬线称为纬线。赤道是一条特殊的纬线，位于 0°纬度，介于南北极之间。
作者图片
📍经度
经度是水平位置的度量。假想的垂直经线称为经线。本初子午线是一条特殊的子午线，位于 0°经度。谈到时区，经度也很重要。
在本初子午线的对面，反子午线位于经度 180 度。
纬线像一个环，而经线像半个环。
📦设置
我们将导入这些库，并在澳大利亚墨尔本设置两个示例位置坐标:
import numpy as np
import pandas as pd
from math import radians, cos, sin, asin, acos, sqrt, pifrom geopy import distance
from geopy.geocoders import Nominatim
import osmnx as ox
import networkx as nxlat1, lon1 = -37.82120, 144.96441 # location 1
lat2, lon2 = -37.88465, 145.08727 # location 2
安装 [osmnx](https://osmnx.readthedocs.io/en/stable/#installation)可能比较繁琐。按照本教程设置环境的一个简单方法是使用 [Google 合作实验室](https://research.google.com › colaboratory):首先，创建一个新的笔记本；二、用!pip install osmnx安装库；第三，重启:从顶部菜单进入运行时>重启运行时>环境准备就绪！
🔵最短球面距离
地球的赤道半径是 6378 公里，极地半径是 6356 公里，所以地球不是一个完美的球体。然而，假设地球是球形的，使我们能够容易地找到近似的距离，这在某些应用中是令人满意的。在本节中，我们将使用哈弗辛公式从地理坐标中找出两个位置之间的球面距离。让我们先熟悉一下哈弗辛函数。
哈弗辛函数如下:
使用哈弗辛公式可以计算出圆心角的哈弗辛，等于球面距离除以球面半径:
我们可以使用哈弗辛函数的第一个定义来转换这个公式，并重新排列它，使d在左边:
现在，是时候把这个翻译成 Python 代码了。有两点需要强调:首先，纬度和经度是以度为单位的，所以在我们将它们代入公式之前，我们必须将它们转换成弧度。其次，我们将使用全球平均值 6371 公里作为球形地球的半径。
def calculate_spherical_distance(lat1, lon1, lat2, lon2, r=6371):
# Convert degrees to radians
coordinates = lat1, lon1, lat2, lon2
# radians(c) is same as c*pi/180
phi1, lambda1, phi2, lambda2 = [
radians(c) for c in coordinates
]
# Apply the haversine formula
a = (np.square(sin((phi2-phi1)/2)) + cos(phi1) * cos(phi2) *
np.square(sin((lambda2-lambda1)/2)))
d = 2*r*asin(np.sqrt(a))
return dprint(f"{calculate_spherical_distance(lat1, lon1, lat2, lon2):.4f} km")
或者，我们可以使用带余弦的哈弗辛函数的第二个定义，并重新排列等式来表示d:
这可以用 Python 表达如下:
def calculate_spherical_distance(lat1, lon1, lat2, lon2, r=6371):
# Convert degrees to radians
coordinates = lat1, lon1, lat2, lon2
phi1, lambda1, phi2, lambda2 = [
radians(c) for c in coordinates
]
# Apply the haversine formula
d = r*acos(cos(phi2-phi1) - cos(phi1) * cos(phi2) *
(1-cos(lambda2-lambda1)))
return dprint(f"{calculate_spherical_distance(lat1, lon1, lat2, lon2):.4f} km")
更实际的是，我们可以使用geopy包在一行代码中获得球面距离:
print(f"{distance.great_circle((lat1, lon1), (lat2, lon2)).km:.4f} km")
另外，用geopy包很容易找到其他距离。例如，我们可以这样基于椭球地球假设得到距离:distance.distance((lat1, lon1), (lat2, lon2)).km。有不同的椭球模型可用，前面的函数使用了WGS-84模型，这里有一个替代语法:distance.geodesic((lat1, lon1), (lat2, lon2), ellipsoid=’WGS-84').km。如果你想了解更多关于这个图书馆的信息，请查看它在计算距离上的资源。
🚗最短旅行距离
在这一节中，我们将看看如何使用 OpenStreetMap 和[OSMnx](https://github.com/gboeing/osmnx)包找到最短的旅行距离。我们将从绘制城市网络图开始:
mel_graph = ox.graph_from_place(
'Melbourne, Australia', network_type='drive', simplify=True
)
ox.plot_graph(mel_graph)
这段代码可能需要一段时间来运行。我们用network_type='drive'来获取行驶距离。其他网络类型也可用。例如，如果我们在步行距离之后，那么我们将代码调整为network_type='walk'。
现在，我们可以使用图表找到行驶距离:
orig_node = ox.distance.nearest_nodes(mel_graph, lon1, lat1)
target_node = ox.distance.nearest_nodes(mel_graph, lon2, lat2)
nx.shortest_path_length(G=mel_graph, source=orig_node, target=target_node, weight='length')
从位置 1 到位置 2 的最短行驶距离为 15，086.094 米。值得注意的是，从位置 2 到位置 1 的距离不一定与从位置 1 到位置 2 的距离相同。
我们可以创建一个计算距离的函数:
def calculate_driving_distance(lat1, lon1, lat2, lon2):
# Get city and country name
geolocator = Nominatim(user_agent="geoapiExercises")
location = geolocator.reverse(f"{lat1}, {lon1}")
address = location.raw['address']
area = f"{address['city']}, {address['country']}" # Get graph for the city
graph = ox.graph_from_place(area, network_type='drive',
simplify=True) # Find shortest driving distance
orig_node = ox.distance.nearest_nodes(graph, lon1, lat1)
target_node = ox.distance.nearest_nodes(graph, lon2, lat2)
length = nx.shortest_path_length(G=graph, source=orig_node,
target=target_node, weight='length')
return length / 1000 # convert from m to kmsprint(f"{calculate_driving_distance(lat1, lon1, lat2, lon2):.2f} km")
就是这样！如果你想了解更多关于这个库的信息，请查看 OSMnx 用户参考和 OSMnx 示例。
您想访问更多这样的内容吗？媒体会员可以无限制地访问媒体上的任何文章。如果您使用 我的推荐链接成为会员，您的一部分会费将直接用于支持我。
谢谢你看我的帖子。如果你感兴趣， 以下是我的一些帖子的链接:
◼️ 用这些技巧充实你的 Jupyter 笔记本
◼️ 用这些技巧整理你的 Jupyter 笔记本
◼️ 有用的 IPython 魔法命令
◼️Python 虚拟数据科学环境简介
◼️git 数据科学简介
◼️ 你会发现有用的 python 中的简单数据可视化
seaborn(python)
◼️️给熊猫用户的 5 个建议
◼️️ 在熊猫中编写 5 个常见的 SQL 查询
再见🏃💨
使用无监督学习寻找平面布置图上的“m”细节
使用 opencv、pytesseract、scikit-learn 和 pandas
斯文·米克在 Unsplash 上的照片
由于可以在搜索中应用个性化过滤器，在网上浏览一处房产会带来无数好处。然而，有一个非常重要的过滤器被遗漏了(至少对于像[ 1 ，[ 2 ]这样的英国房产搜索网站来说是这样的):房产的平方英尺分割。
有很多情况下，我会发现一个伟大的财产(大小，位置，价格等。)直到我看到它的平面图细节。如果我们可以根据房间的平方米来过滤属性结果，那将会节省我们很多时间。因此，我决定做一个实验，看看如何提取给定房地产广告平面图上的所有“m”细节。
介绍
在这篇文章中，我阅读了一个在线平面图图像并对其进行了处理，然后解释了为什么我使用无监督学习算法将相关的文本块分组在一起，演示了进一步的清理并测试了几个不同平面图图像的结果。 在结论部分 ，我提到了我遇到的局限性和任何改进的空间(隐喻性的)。
步骤:
**第一步:**读取并处理输入的平面图图像(opencv)
**第二步:**检测字符及其位置(pytesserac)
**第三步:**聚类字符(scikit-learn)
**第四步:**处理聚类(熊猫)
测试:
导入库
我们需要以下内容:
from matplotlib import pyplot as plt
import cv2
import numpy as npimport pandas as pd
pd.set_option('display.max_rows', 500)
pd.set_option('display.max_columns', 500)
pd.set_option('display.max_colwidth', 500)%matplotlib inline
plt.rcParams["figure.figsize"] = (40,20)import pytesseract
pytesseract.pytesseract.tesseract_cmd = r'/usr/local/bin/tesseract/'from sklearn.preprocessing import StandardScaler
from sklearn.cluster import DBSCAN
步骤 1:读取和处理输入图像
我们可以很容易地阅读如下平面图图像:
img = cv2.imread('Dummy Floorplanner Input Image.png')plt.imshow(img)
按作者输入图像
我们正在处理输入图像，这样，当输入到 pytesseract 时，它可以揭示更多关于文本块的信息。请阅读下文了解更多关于阈值的信息。
https://pyimagesearch.com/2021/04/28/opencv-thresholding-cv2-threshold/
经过处理后，我们的输入图像如下所示:
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
blur = cv2.GaussianBlur(gray, (3,3), 0)
thresh = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY_INV + cv2.THRESH_OTSU)[1]
plt.imshow(gray)
plt.imshow(blur)
plt.imshow(thresh)
按作者分类的图像(灰色、模糊、阈值)
**第二步:**检测字符及其位置
Pytesseract 从左上到右下开始读取图像像素(就像看书一样)。因此，如果一个房间的信息与另一个房间的信息在同一水平层上，那么两个房间的文本将被连接到一个房间，并且没有逻辑方法来正确解析完整的字符串，如下所示(餐厅和卧室的名称和尺寸被连接)。
print(pytesseract.image_to_string(thresh))
作者图片
幸运的是，还有另一个 pytesseract 函数(image_to_boxes)可以检测每个字符及其在图像上的位置。
df_img = pd.DataFrame([x.split(' ') for x in pytesseract.image_to_boxes(thresh).split('\n')], columns=['char', 'left', 'top', 'width', 'height', 'other'])df_img = df_img[ ~ df_img['left'].isnull()]# dropping whitespace characters like
# [',' '.' '/' '~' '"' "'" ':' '°' '-' '|' '=' '%' '”']df_img = df_img[ ~ df_img['char'].str.contains(r'[^\w\s]')].reset_index(drop=True)df_img[['left', 'top', 'width', 'height']] = df_img[['left', 'top', 'width', 'height']].astype(int)df_img.head(20)
作者图片
让我们根据这些字符的 x-y 坐标来绘制它们。
fig, ax = plt.subplots()ax.scatter(df_img['left'].tolist(), df_img['top'].tolist())for i, txt in enumerate(df_img['char'].tolist()):
ax.annotate(txt, (df_img['left'].tolist()[i],
df_img['top'].tolist()[i]),
textcoords='data',
fontsize=28)
作者图片
基本上，我们需要对近距离内的数据点(字符)进行分组。
**第三步:**聚类字符
如果需要，请阅读下面的帖子来了解 DBSCAN 是如何工作的。
使用 DBSCAN，我们能够将相关的字符块组合在一起。当该聚类算法预测值为-1 时，该数据点被视为异常值。
X = StandardScaler().fit_transform(df_img[['left', 'top']].values)db = DBSCAN(eps=0.19, min_samples=10)db.fit(X)y_pred = db.fit_predict(X)
plt.figure(figsize=(10,6))
plt.scatter(X[:,0], X[:,1],c=y_pred, cmap='Paired')
plt.title("Clusters determined by DBSCAN")df_img['cluster'] = pd.Series(y_pred)
df_img.groupby(['cluster'])['char'].apply(lambda x: ' '.join(x)).reset_index()
作者提供的图片
从上面的图像中我们可以看到，聚类算法在对平面图图像上的字符进行分组时工作得非常好。然而，它仍然需要一些进一步的清洁和处理。
步骤 4:处理集群
通常，在平面图上，我们有位置的类型(卧室、厨房、花园等。)，其尺寸以米为单位，最后以英尺为单位。有了这些信息，我们可以假设我们第一次看到一个聚类的数字是我们可以连接文本的最低水平。换句话说，在一个群中，如果在一个数字第一次出现的级别之下有一个数字，这可能是一个平方英尺信息。让我们从数据框中删除所有平方英尺(不是米)的行。
df_cc = df_img.copy().reset_index(drop=True)for cluster_no in df_cc['cluster'].unique():
index_char_top_list = []
# if the data point is not an outlier
if cluster_no!=-1:
index_char_top_list = [(index, char, top) for index, char, top in
zip(df_cc[(df_cc['cluster']==cluster_no)].index,
df_cc[(df_cc['cluster']==cluster_no)]['char'].values,
df_cc[(df_cc['cluster']==cluster_no)]['top'].values)ifchar.isdigit()] if index_char_top_list:
df_cc = df_cc[~ ((df_cc['cluster']==cluster_no) & (df_cc['top'] <= ( index_char_top_list[0][2] - 5 ))) ]df_img.shape[0], df_cc.shape[0]
# (149 rows went down to 104 rows)
让我们创建一个函数，它可以将一串数字解析成宽度、长度维度。
def dimension_splitter(input_text):
input_text_len = len(input_text) if input_text_len%2==0:
split_text_by = int(input_text_len/2)
else:
split_text_by = int(input_text_len/2+0.5)
dim1 = input_text[:split_text_by] dim2 = input_text[split_text_by:]
dim1 = float('{}.{}'.format(dim1[:-2], dim1[-2:])) dim2 = float('{}.{}'.format(dim2[:-2], dim2[-2:])) return dim1, dim2
让我们对清理后的聚类数据框进行分组，并引入维度列。
df_cc_grouped = df_cc.groupby(['cluster'])['char'].apply(lambda x: ' '.join(x)).reset_index(name='text')df_cc_grouped['text_digits'] = df_cc_grouped.apply(lambda x: ''.join([y for y in x['text'] if y.isdigit()]), axis=1)df_cc_grouped['text_digits_len'] = df_cc_grouped.apply(lambda x: len([y for y in x['text'] if y.isdigit()]), axis=1)df_cc_grouped = df_cc_grouped[(df_cc_grouped['cluster']!=-1) &
(df_cc_grouped['text_digits_len']>=5)].reset_index(drop=True)df_cc_grouped['room'] = df_cc_grouped.apply(
lambda x: x['text'][:[x.isdigit() for x in x['text']].index(True)].strip()
, axis=1)df_cc_grouped['length'] = df_cc_grouped.apply(lambda x: dimension_splitter(x['text_digits'])[0]
, axis=1)df_cc_grouped['width'] = df_cc_grouped.apply(lambda x: dimension_splitter(x['text_digits'])[1]
, axis=1)df_cc_grouped['surface'] = np.round(df_cc_grouped['length'] * df_cc_grouped['width'], 2)df_cc_grouped
作者提供的图片
我们可以正确地提取所有的平方米信息。让我们看看我们的逻辑是否适用于其他图像。
试验
在第一项测试中，在检测步骤中，pytesseract 误读了“接待室”中的一些字母，并在检测到 Utility 之前多加了一个“a”。然而，在此之后，无监督学习如预期那样工作，并正确地揭示了所有信息。
img_test1 = cv2.imread('Dummy Floorplanner Test 1.png')plt.imshow(img_test1)df_cc_grouped_test1
作者提供的图片
执行第 2 号测试时没有出现任何问题。
img_test2 = cv2.imread('Dummy Floorplanner Test 2.png')plt.imshow(img_test2)df_cc_grouped_test2
作者提供的图片
尽管平面布局不同(纵向或横向)，我们能够在大多数地方提取平方米信息。
结论
在这篇文章中，我演示了如何使用 OCR 库(pytesseract)来检测房地产平面图上的字母和数字，然后使用无监督学习算法(DBSCAN)来收集相关信息，以找到每个房间的平方米信息。
如果 pytesseract 可以正确地检测字符，我们可以毫无问题地找到房产的平方米分割。我们利用了这样一个事实，即英国的房产平面图具有几乎相同的结构:(1)房间的类型，下面的尺寸(2)米，下面的尺寸(3)英尺。但是，该脚本需要进一步开发，以处理具有不同结构的平面图图像，或者对缺失/未检测到的数据进行一些科学猜测。
感谢阅读！
通过分析寻找适合商业地产的产品市场
为房地产营销和分析的直觉增加数据
本文代表我个人的想法、观点和分析，与任何其他组织无关。
为租户匹配商业空间对双方来说都是一项具有挑战性的工作。房东通常会对企业做出不确定的长期承诺，反过来，企业也会承担在租赁期内支付租金的风险。
一个空间的价值通常取决于它的位置。在许多地区，你可以开车或步行去感受一个地区。你还看到其他什么企业在运营？谁倾向于走在街上？虽然作为一个开始很有帮助，但这是一个有限的时间点样本，并且错过了更大的画面——该领域与其他选项相比如何？
历史上，数据被汇总到共同认可的地理位置进行分析，作为展示位置的一种手段。我们直观地了解一个城市或邮政编码。不过，这些地区很大，对于决定长期租赁的租户来说，对周围地区更有信心有助于做出更好的决定。
今天，我们可以超越传统的邮政编码或城市界限，但仍然可以开发直观的数据工具来帮助企业做出这些关键决策。我们可以问类似*这样的问题:“有多少人住在离一处房产 15 分钟步行路程的范围内？他们的人口构成是怎样的？”*或“附近还有哪些其他类型的企业？”。
具有多种方法和时间选项的仪表板示例
我很幸运生活在一个快速发展和成长的地区。每天，我都会路过处于不同发展阶段的建筑。我注意到一栋建筑的开发商在 LinkedIn 上发布了一个帖子，寻找商业租户，以占据一个混合用途开发项目的底层。
北卡罗来纳州罗利市中心附近新开发项目的物业营销
我想知道——谁最适合这份工作？我大致了解这个领域，但是我可以引入哪些数据来揭示新的机会或见解呢？
免责声明:这不考虑许可限制、空间成本、建筑限制或其他相关信息。这只是一个分析框架，用于确定一处房产对潜在租户的价值。
邻里
以下是一些关于格伦伍德南区的基本信息:
- 它位于罗利市中心附近，大型豪华公寓建筑发展迅速
- 它最近的名声是作为一个聚会社区，俱乐部和酒吧散布在这个小区域
- 考虑到最近的发展、许可申请和价格上涨，未来 3-5 年可能会转向更高档的区域
除了高端办公空间和众多的酒吧和餐厅之外，它非常适合步行，步行即可到达多所大学。
房产
新开发项目位于 Glenwood South 的北端，交通繁忙。它位于住宅和商业混合区的正中央，附近步行即可到达大学、高中和小学，此外还有高速公路入口(在下图中用黄色标出；不完全是高速公路，但也不是普通的街道)。
来自谷歌地图的带有注释的酒店位置图片
鉴于该物业的位置和周边地区，我的观点是，它将从步行交通中受益最多，其次(尽管仍然很有价值)是位于高能见度街道上。我怀疑——像许多城市地区一样——将会有有限的停车位，限制了在特定时间可以停靠的汽车数量。
所以，知道了以上内容，这让我们不禁要问…有多少人住在离该地点步行的距离内，他们的人口构成是什么？该位置周围还有哪些企业可能会让某类租户受益或受到抑制？
分析
这就是有趣的地方——至少如果你是像我一样的分析呆子。
让我们假设在一个城市地区，某人将步行 20 分钟到达一个位置。这是地图上 20 分钟步行半径的样子。给定这个半径，我们对该区域发生的活动了解多少？
距离该位置约 20 分钟的步行半径；作者图片
数据显示，在这段步行时间内，大约有 10，500 人和 25，000 个工作岗位。考虑到数据报告和今天之间的时间差，我估计这两个数字今天会更高。这些数字明显高于县基准，表明这是一个密集的混合使用区(“生活-工作-娱乐”，俗话说)。
在那些生活在 20 分钟步行范围内的人中，大约 50%在 18 岁至 34 岁之间(相比之下，该县为 24%)，40%的工作是办公室工作(相比之下，该县为 30%)，这意味着“年轻专业人士”可能是人口的一个适当特征。
在人口和办公室工作之间，每天有足够多的人在 20 分钟的步行距离内填满当地的 NHL 竞技场
该地区非常适合步行，满分为 100 分，得分约为 85 分，典型家庭收入约为 73，000 美元。家庭收入低于县平均水平，然而，典型的家庭规模也更小，这表明人均可支配资金更高。
估计的覆盖区域数据，距离酒店大约 20 分钟的步行时间
在步行范围内大约有 100 家有执照的餐馆，其中大部分在南边，因为北边是住宅区。
步行半径范围内获得许可的餐饮机构的大致位置；作者图片
步行区内约有 16 个被视为“豪华”的公寓区——更不用说额外的公寓大楼了——还有同一个街区的多块土地，作为奖励，它们可能在租约有效期内开发。
最后，来自北卡罗来纳州的汽车交通数据表明，在给定的一天，大约有 11，500 辆汽车经过这个地方。这是 2020 年的数据，受到 COVID 的影响，因此，我的观点是，今天这一数据已经大幅增加。
这意味着什么
因此，综上所述，我们可以说这个地区非常适合步行和城市化，周边人口年轻、富裕且不断增长。在汽车交通、当地人口和附近的工作岗位之间，它可能有一个每天至少 40，000 人的可寻址市场(或者，用当地人可以欣赏的术语来说，足以填满附近 UNC 教堂山的篮球体育场两倍)。
最高价值的解决方案是一个企业将受益于(1)周围地区的人口和工人；(2)在(非常)邻近的小学接送孩子的家长，以及(3)在附近步行距离内上高中的 2000 多名学生。
哪些类型的企业可以从这一目标受众和市场规模中受益？
银行在寻找市场份额的过程中一直在这个三角地带进行着激烈的竞争，因此资本相对较低的租户(相比之下，比如说，必须建造一个商业厨房)可能是一家寻找高知名度分支机构的银行。他们可以通过步行，加上路过的道路交通带来的隐性品牌效应，接触到年轻、高收入的客户群。
高端健身精品店可能会有机会，但另一个需要考虑的问题是，所有的豪华公寓楼都将有健身房，一些写字楼也将有健身房，因此潜在市场可能会比步行距离内的人口/员工的预期要小。
酒吧或餐厅/快速服务餐厅(QSR)是一种可能性，尽管如地图所示，许多已经存在，更多的正在计划中，作为新开发的一部分，在附近步行即可到达，因此它需要是一个独特的概念，尚未在更广泛的区域内捕获。
像 CVS 或 Walgreens 这样的药店可能会根据空间的大小而发挥作用，因为它们都有大约一英里远的位置，这在城市环境中可能会错过希望立即方便的人的市场机会。
我要提到的最后一个想法是类似于 Amazon Go 商店的东西，这无疑是一个很长的机会/非常低的可能性，但我认为具有很高的价值。它的价值跨越年龄组，它将受益于混合使用驱动交通通过。该位置位于步行交通的中心，可用于取包裹或用作微型配送中心。除了所讨论的步行可达性之外，一个企业还可以通过 10 分钟的车程到达约 120，000 人，释放出更多的价值。它的快速进出模式将有利于(可能的)低可用停车位。
这只是现代数据工具和资源可用能力的一个示例。我们选择步行时间，因为它与位置相关，但任何驾驶时间或径向距离(即 2 英里的圆圈)也是选项。餐厅被引用是因为它们是空间的常见用途，但如果医疗保健客户是目标受众，它也可以很容易地被牙科诊所或医生取代。
如今，数据可以为商业房地产代理商带来的本地市场知识和专业知识增加一层，从而为决策带来更多的说服力。找到正确的位置从来都不容易，但是增加相关的、有针对性的信息的数据访问可以显著增加决策的价值。
对这类材料和分析感兴趣？ 在 LinkedIn 上联系我，或者在 jordan@jordanbean.com 给我发邮件。
使用 TF-IDF 和 Python 查找相关文章
原文：https://towardsdatascience.com/finding-related-articles-with-tf-idf-and-python-d6e1cd10f735
如何找到与 TF-IDF 相关的文章？用 Python 实现 TF-IDF 算法。
马丁·范·登·霍维尔在 Unsplash 上的照片
正如我在关于 总结文字 的文章中提到的，我收集了过去十年的新闻文章。总结完这些文章后，我面临的下一个挑战是找到相关的文章。如果我有一篇文章，我如何找到与这篇文章相关的文章？立即导致问题什么是相关文章？
经过初步研究，选择 TF-IDF 作为算法。这是一个复杂性有限的旧算法，因此可以创建一个实现，并从中学习如何找到类似的文章。通常，有更好的实现可用，例如在 scikit-learn 中，但是我更喜欢为了学习而构建自己的实现。请注意，完全实现需要一些时间。
TF-IDF 算法
那么，什么是 TF-IDF 呢？我们分解缩写的时候，说的是词频——逆文档频，一口。术语频率是指特定术语(词)在文本中出现的相对次数；
词频(图片由作者提供)
文档 d 中的 t 项的 TF 等于文档 d 中出现的 t 项的频率除以文档 d 中所有*t’*项的频率之和。示例:
'The green man is walking around the green house'
TF(man) = 1 / 9 = 0.11
TF(green) = 2 / 9 = 0.22
这是对一个单词在文档中的重要性的度量，根据文档的长度进行校正。
IDF 是对单词在整个语料库(分析的所有文档的集合)中的重要性的度量。这个想法是，如果一个单词出现在很多文档中，它不会增加很多信息。标准功能是
反向文档频率(作者图片)
IDF 的计算方法是将文档数( N )除以一个术语出现的文档总数，并取其对数值。如果一个术语在 10.000 个文档中出现 10 次，则 IDF 等于 3。在同一组文档中出现 100 次的术语的 IDF 值为 2。当一个术语出现在所有文档中时，IDF 值等于 0.0。对数值用于减少 IDF 可能具有的大范围值。
最后，一项的 TF-IDF 值等于 TF 乘以 IDF:
TF-IDF 公式(图片由作者提供)
上面的公式是 TF 和 IDF 的标准公式。更多变种可以在 TF-IDF 的维基百科页面找到。
为所有文档的语料库中的每个单词计算 TF-IDF 值，包括术语出现次数为 0 的文档。作为例子，下面两句话
0: 'the man walked around the green house',
1: 'the children sat around the fire'
产生以下 TF-IDF 值:
计算值(图片由作者提供)
对于所有唯一的单词，为每个句子计算 TF-IDF 值。在以下情况下，该值为 0.0
- 由于 TF 为 0.0，特定单词不会出现在句子中(例如，第二句中的“green”)
- 一个特定的词出现在所有的句子中，例如“大约”,因为 1 的对数值等于 0.0。在这种情况下，所有句子的 TF-IDF 值都是 0.0
计算距离
既然计算了 TF-IDF 值，那么可以通过计算这些值之间的余弦相似度来确定两个文档之间的距离。这是通过使用每个单词的 TF-IDF 值(上面数据帧中的一行)作为向量来完成的。
余弦距离是 N 维空间中从原点到两点的直线之间的角度差。在二维空间中形象化是相对简单的。上面的例子是 9 维的，超出了我的想象。在二维空间中:
余弦差的定义(图片由作者提供)
余弦差可以使用欧几里得点积公式计算:
计算余弦差(图片由作者提供)
它是两个向量的点积除以向量长度的乘积。当两个向量平行时，余弦差为 1.0，当它们正交时，余弦差为 0.0。
该功能需要应用于语料库中的所有文档组合。从 A 到 B 的差等于从 B 到 A 的差，从 A 到 A 的差是 1.0。下表显示了距离矩阵在对角线上的镜像:
差 A-B 等于 B-A(图片由作者提供)
完成所有这些计算后，是时候开始实现了！
实施 TF-IDF
干得好，你挺过了理论背景！但是现在的问题是如何实现这个功能。目标是计算几个文档之间的距离。反过来，我们需要 TF-IDF 来计算数据，我们需要 TF 和 IDF 来计算 TF-IDF。
计算 TF 和 IDF 需要对文档进行一些处理。我们需要 TF 的每个文档中每个单词的出现次数和文档中单词的总数，以及 IDF 的每个单词出现在文档中的数量和文档的数量。这些都需要相同的通用文档解析，因此增加了一个预处理步骤:
class DocumentDistance:
distancespre = {}
def add_documents(self, docs: typing.List[str],
) -> None:
"""
Calculate the distance between the documents
:param documents: list of documents
:return: None
"""
word_occ_per_doc, doc_lens, doc_per_word = self.pre_proces_data(docs)
# Calculate TF values
tfs = []
for i in range(len(word_occ_per_doc)):
tfs.append(self.compute_tf(word_occ_per_doc[i], doc_lens[i]))
# Calculate IDF values
idfs = self.compute_idf(doc_per_word, len(docs))
# Calculate TF-IDF values
tfidfs = []
for i in range(len(tfs)):
tfidfs.append(self.compute_tfidf(tfs[i], idfs))
# Calculate distances
self.calculate_distances(tfidfs)
定义了一个类，它将所有文档之间的距离存储为一个二维数组。这个数组将由calculate_distances方法填充。但是首先，我们将创建 pre_process_data 方法，该方法将解析列表中的所有文档，并返回每个文档的单词出现次数、文档长度以及每个单词在文档中出现的次数:
from nltk import word_tokenize
class DocumentDistance:
...
def pre_proces_data(self,
documents: typing.List[str]
) -> (typing.Dict[int, typing.Dict[str, int]],
typing.Dict[int, int],
typing.Dict[str, int]):
"""
Pre proces the documents
Translate a dictionary of ID's and sentences to two dictionaries:
- bag_of_words: dictionary with IDs and a list with the words in the text
- word_occurences: dictionary with IDs and per document word counts
:param documents:
:return: dictionary with word count per document, dictionary with sentence lengths
"""
# 1\. Tokenize sentences and determine the complete set of unique words
bag_of_words = []
unique_words = set()
for doc in documents:
doc = self.clean_text(doc)
words = word_tokenize(doc)
words = self.clean_word_list(words)
bag_of_words.append(words)
unique_words = unique_words.union(set(words))
# 2\. Determine word occurences in each sentence for all words
word_occurences = []
sentence_lengths = []
for words in bag_of_words:
now = dict.fromkeys(unique_words, 0)
for word in words:
now[word] += 1
word_occurences.append(now)
sentence_lengths.append(len(words))
# 3\. Count documents per word
doc_count_per_word = dict.fromkeys(word_occurences[0], 0)
# Travese all documents and words
# If a word is present in a document, the doc_count_per_word value of
# the word is increased
for document in word_occurences():
for word, val in document.items():
if val > 0:
doc_count_per_word[word] += 1
return word_occurences, sentence_lengths, doc_count_per_word
第一部分通过在 NLTK 包的word_tokenize、部分的帮助下对文档进行标记，将文档分割成单词。在标记之前，会进行一些文档清理，比如将文档转换成小写(稍后讨论)。在clean_word_list method中清理词表，目前为空。每个文档的令牌列表存储在bag_of_words变量中，这是一个每个文档有一个条目的列表，包含一个令牌列表。在这个循环中，创建了一组所有唯一的单词。这个集合unique_words包含了语料库中出现的所有唯一单词。
步骤 2 为所有文档确定所有唯一单词的出现次数。对于每个文档(在bag of words上的循环)，为所有唯一的单词(dict.fromkeys(…))创建一个值为 0(零)的字典。然后，代码遍历文档中的所有单词，并将字典值增加 1(一)。该词典被添加到word_occurences中，为所有文档创建一个词典列表，其字数和文档长度存储在sentence_lengths中。
最后一步，步骤 3，统计每个单词在文档中出现的次数。首先，创建一个列表doc_count_per word，其中每个唯一单词出现 0 次。当文档中特定单词的单词计数大于零时，该单词存在于文档中。
预处理的结果是:
Documents:
'the man walked around the green house'
'the children sat around the fire'
word_occurences:
{
0: {'green': 1, 'fire': 0, 'house': 1, 'sat': 0, 'around': 1, 'man': 1,
'the': 2, 'children': 0, 'walked': 1},
1: {'green': 0, 'fire': 1, 'house': 0, 'sat': 1, 'around': 1, 'man': 0,
'the': 2, 'children': 1, 'walked': 0}
}
sentence_lengths:
[7, 6]
doc_count_per_word:
{'around': 2, 'green': 1, 'house': 1, 'walked': 1, 'sat': 1,
'children': 1, 'man': 1, 'fire': 1, 'the': 2}
有了这些数据集，可以相对直接地计算 TF 和 IDF:
def compute_tf(self,
wordcount: typing.Dict[str, int],
words: typing.List[str]
) -> typing.Dict[str, float]:
"""
Calculates the Term Frequency (TF)
:param wordcount: dictionary with mapping from word to count
:param words: list of words in the sentence
:return: dictionary mapping word to its frequency
"""
tf_dict = {}
sentencelength = len(wordcount)
for word, count in wordcount.items():
tf_dict[word] = float(count) / sentencelength
return tf_dict
def compute_idf(self,
doc_count_per_word: typing.List[typing.Dict[str, int]],
no_documents: int
) -> typing.Dict[str, int]:
"""
Calculates the inverse data frequency (IDF)
:param doc_count_per_word: dictionary with all documents. A document is a dictionary of TF
:param no_documents: number of documents
:return: IDF value for all words
"""
idf_dict = {}
for word, val in doc_count_per_word.items():
idf_dict[word] = math.log(float(no_documents) / val)
return idf_dict
def compute_tfidf(self,
tfs: typing.Dict[str, float],
idfs: typing.Dict[str, float]
) -> typing.Dict[str, float]:
"""
Calculte the TF-IDF score for all words for a document
:param tfs: TFS value per word
:param idfs: Dictionary with the IDF value for all words
:return: TF-IDF values for all words
"""
tfidf = {}
for word, val in tfs.items():
tfidf[word] = val * idfs[word]
return tfidf
对每个文档调用compute_tf方法。每个单词的 TF 值是通过将出现次数除以句子长度来计算的(计算结果被强制转换为浮点类型)。
使用每个单词的文档计数和文档总数来调用compute_idf。讨论的公式适用于这些值。
最后，通过将 TF 值乘以相应的 IDF 值来计算每个句子每个单词的 TF-IDF 值。
tfs:
[{'sat': 0.00, 'green': 0.14, 'walked': 0.14, 'the': 0.28,
'around': 0.14, 'children': 0.00, 'fire': 0.00, 'man': 0.14,
'house': 0.14},
{'sat': 0.16, 'green': 0.00, 'walked': 0.00, 'the': 0.33,
'around': 0.16, 'children': 0.16, 'fire': 0.16, 'man': 0.00,
'house': 0.00}]
idfs:
{'sat': 0.69, 'green': 0.69, 'walked': 0.69, 'the': 0.00,
'around': 0.00, 'children': 0.69, 'fire': 0.69, 'man': 0.69,
'house': 0.69}
tfidfs:
[{'sat': 0.00, 'green': 0.09, 'walked': 0.09, 'the': 0.00,
'around': 0.00, 'children': 0.00, 'fire': 0.00, 'man': 0.09,
'house': 0.09},
{'sat': 0.11, 'green': 0.00, 'walked': 0.00, 'the': 0.00,
'around': 0.00, 'children': 0.11, 'fire': 0.11, 'man': 0.00,
'house': 0.0}
]
现在我们有了 TF-IDFS 值，就可以计算不同文档之间的距离了:
def normalize(self,
vector: typing.Dict[str, float]
) -> float:
"""
Normalize the dictionary entries (first level)
:param tfidfs: dictiory of dictionarys
:return: dictionary with normalized values
"""
sumsq = 0
for i in range(len(vector)):
sumsq += pow(vector[i], 2)
return math.sqrt(sumsq)
def calculate_distances(self,
tfidfs: typing.List[typing.Dict[str, float]]
) -> None:
"""
Calculate the distances between all elements in tfidfs
:param tfidfs: The dictionary of dictionaries
:return: None
"""
vectors = []
# Extract arrays of numbers
for tfidf in tfidfs:
vectors.append(list(tfidf.values()))
self.distances = [[0.0] * len(vectors) for _ in range(len(vectors))]
for key_1 in range(len(vectors)):
for key_2 in range(len(vectors)):
distance = np.dot(vectors[key_1], vectors[key_2]) / \
(self.normalize(vectors[key_1])* self.normalize(vectors[key_2]))
self.distances[key_1][key_2] = distance
如本文第一部分所述，两个向量之间的余弦距离是通过将向量的点积除以归一化向量的积来计算的。
为了计算点积，使用了来自 numpy 库的dot方法。点积通过对成对相乘的值求和来计算:
计算矩阵的点积(图片由作者提供)
向量的标准化值等于原点和向量所标识的点之间的直线长度。这等于每个维度的平方和的平方根:
标准化矢量(图片由作者提供)
矢量的归一化是通过normalize()方法实现的。所有矢量组合之间的距离用calculate_distances()方法计算。所有距离都存储在该类的distances属性中。在计算距离之前，该变量被初始化为 0.0 值的 N x N 矩阵。
一个例子:
Sentences:
'the man walked around the green house'
'the children sat around the fire'
'a man set a green house on fire'
Distances:
[1.00, 0.28, 0.11]
[0.28, 1.00, 0.03]
[0.11, 0.03, 1.00]
请注意，值为 1.0 时距离最小，值为 0.0 时距离最大。由于两个句子没有共同的单词，所以它们的距离是 0.0。一句话和它本身的距离是 1.0。
性能改进
有了由三个文档组成的语料库，每个文档由一个小句子组成，距离可以快速计算出来。当文档变成文章并且数量增加到数百时，这个时间增加得很快。下图显示了独特词的数量与新闻文章数量的函数关系。
独特单词的数量(图片由作者提供)
当处理 1000 篇文章时，唯一单词的数量增加到几乎 15000 个。这意味着描述文档的每个向量有 15000 个条目。一个点乘需要 2.25 亿次乘法和 2.25 亿次加法。每矢量。归一化一个向量也是 2.25 亿次乘法(计算平方)，2.25 亿次加法和一个平方根。因此，两个向量之间的距离计算是 7.75 亿次乘法，7.75 亿次加法，1 次平方根，1 次除法和 1 次加法。所有这 100 万次来填充整个距离数组。可以想象这需要一些时间…
那么如何才能减少工作量呢？黄金法则是尽可能少花时间。所以让我们看看我们能做些什么来优化。
距离计算
在第一部分中已经提到，距离矩阵包含重复值。A 到 B 的距离等于 B 到 A 的距离，我们只需要计算一次，分两个地方加到矩阵里。
其次，矩阵的对角线是 1.0，因为 A 和 A 之间的距离总是 1.0。这里不需要计算。这两步将使计算的距离减半。
第三，对于每次计算，所有向量都被归一化。这是多余的。我们可以预先计算所有的归一化，并将其存储在数组中以备再次使用。这将把距离计算减少 2/3，对于每个矢量组合，只需要计算点。
这些改进将所需时间减少了 85%!
def normalize(self,
tfidfs: typing.List[typing.Dict[str, float]]
) -> typing.Dict[int, float]:
"""
Normalize the dictionary entries (first level)
:param tfidfs: dictiory of dictionarys
:return: dictionary with normalized values
"""
norms = []
for i in range(len(tfidfs)):
vals = list(tfidfs[i].values())
sumsq = 0
for i in range(len(vals)):
sumsq += pow(vals[i], 2)
norms.append(math.sqrt(sumsq))
return norms
def calculate_distances(self,
tfidfs: typing.List[typing.Dict[str, float]]
) -> None:
"""
Calculate the distances between all elements in tfidfs
:param tfidfs: The dictionary of dictionaries
:return: None
"""
norms = self.normalize(tfidfs)
vectors = []
# Extract arrays of numbers
for tfidf in tfidfs:
vectors.append(list(tfidf.values()))
self.distances = [[1.0] * len(vectors) for _ in range(len(vectors))]
for key_1 in range(len(vectors)):
for key_2 in range(key_1 + 1, len(vectors)):
distance = np.dot(vectors[key_1], vectors[key_2]) / (norms[key_1] * norms[key_2])
self.distances[key_1][key_2] = distance
self.distances[key_2][key_1] = distance
对代码的一些小的改变导致计算时间减少了 85%。但是有更多的可能性
减少向量大小
距离计算是代码中最昂贵的部分。上述更改大大减少了计算时间。但是还有一个方面会极大地影响性能，那就是文档向量的大小。对于 1.000 篇文章，每个文档向量是 15.000 个元素。如果我们能够移除这个向量大小，计算时间将会受益。
那么我们能做些什么来减少向量的大小呢？我们如何找到可以去除的没有影响的元素？查看 IDF 函数，它显示停止字对距离计算没有影响。它们出现在每个文档中，因此所有向量中的 IDF 值和 TF-IDF 将为 0.0。将所有这些零相乘仍然需要时间，所以第一步是从唯一单词列表中删除停用词。
NLTK 工具包中提供了每种语言的停用词。工具包提供了停用词列表，可以从唯一的词中过滤掉这些停用词。
第二，只在一个文档中出现的单词对于所有其他文档总是乘以零。我们可以安全地从向量中移除这个单词。
最后，最后一步是使用词干。单词列表将包含“house”和“houses”等单词。通过将这些合并到 house，同样的意思，我们可以进一步减少单词列表并提高整个算法的质量。它将不再把“house”和“houses”视为不同的词，而是一个词。NLTK 提供了几个词干分析器，支持多种语言，其中使用了 雪球词干分析器 。
词干化也减少了所有动词的基数。像‘like’，‘liking’，‘liked’等词都将被简化为‘like’这个词干。
方法reduce_word_list已经被引入(但为空),所以现在我们可以用它来应用这些规则:
from nltk.stem.snowball import SnowballStemmer
from nltk.corpus import stopwords
stop_words = set(stopwords.words("dutch"))
stemmer = SnowballStemmer("dutch")
def clean_word_list(self, words: typing.List[str]) -> typing.List[str]:
"""
Clean a worlist
:param words: Original wordlist
:return: Cleaned wordlist
"""
words = [x for x in words if x not in self.stop_words and
len(x) > 2 and not x.isnumeric()]
words = [self.stemmer.stem(plural) for plural in words]
return list(set(words))
第二个规则是去除只存在于一个文档中的单词，通过在第三个步骤之后增加第四个步骤，该规则被应用于预处理功能的修改版本中:
# 4\. Drop words appearing in one or all document
words_to_drop = []
no_docs = int(len(documents) * .90)
for word, cnt in doc_count_per_word.items():
if cnt == 1 or cnt > no_docs:
words_to_drop.append(word)
for word in words_to_drop:
doc_count_per_word.pop(word)
for sent in word_occurences:
sent.pop(word)
第一个循环收集出现一次或总是被收集的单词。在这个循环中不可能修改单词列表，因为它会改变用于迭代的集合。第二步需要从文档计数器和每个文档中删除( pop )多余的文字。这是一个相对昂贵的操作，但向量大小的减少是值得的。
对所有这些向量缩减步骤的效果感到好奇吗？让我们看看:
独特单词的数量(图片由作者提供)
蓝线是唯一单词的原始数量。灰线是减少的字数。因此，在 1.000 篇文章中，字数从 15.000 减少到 5.000。每个向量乘法减少 66%。
代码用 cProfile 进行分析，结果用 snake_viz 可视化。原始代码的执行时间:
剖析原始代码(作者截图)
添加文档需要 31.2 秒，其中 29.6 秒用于计算距离，0.4 秒用于预处理数据。
29.6 calculate_distance(...)
0.6 compute_tf(...)
0.4 pre_proces_data(...)
0.4 compute_tfidf(...)
0.1 compute_idf(...)
优化版本的执行时间:
剖析优化版本(作者截图)
添加文档所需的时间现在为 5.39 秒(was 31.2)，其中 4.32 秒(was 29.6)用于计算距离，0.9 秒用于处理数据(was 0.4)。
所以预处理需要双倍的时间，增加了 0.5 秒的执行时间，但是距离计算只需要原来时间的 14%。额外的预处理时间很容易收回。
4.2 calculate_distance(...)
0.9 pre_proces_data(...)
0.1 compute_tf(...)
0.0 compute_tfidf(...)
0.0 compute_idf(...)
开心吗？
对结果满意吗？老实说，没有。我真的很喜欢实现所有计算文档距离的逻辑。看到可能的性能改进是令人满意的。所以从算法开发的角度来看，我很高兴。
但是应用该算法并没有产生很好的效果。例如，关于寒冷天气的新闻文章似乎与一篇关于伊朗妇女权利的文章更相关，而不是与一篇关于有足够的天然冰可以滑冰的文章更相关。我以为会是另一种情况…
这并不像是浪费了几个小时的快乐编程时间，但是一定会进行一些额外的研究来找到更好的实现。尽管如此，我还是希望您学到了一些关于用 Python 实现算法和寻找提高性能的方法的知识。
最后的话
一如既往，完整的代码可以在我的 Github 上找到。在这里你可以找到本文中讨论的实现，以及一个扩展版本，支持多种语言和TF 和 IDF 的多种实现，如 TF-IDF 的维基百科页面上所讨论的。
我希望你喜欢这篇文章。要获得更多灵感，请查看我的其他文章:
如果你喜欢这个故事，请点击关注按钮！
免责声明:本文包含的观点和看法仅归作者所有。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
