验证人群健康慢性指标之间的关系
糖尿病、慢性肾病和心血管疾病之间的联系
在的最后一个故事中，我们开始研究美国疾病控制和预防中心(CDC)15 年的慢性病数据集。探索性数据分析始于理解数据的列和行，以及与进一步分析相关的内容。
在这篇文章中，我们将深入理解这 40 万行和 17 个主题类别，这需要将数据帧转换成数据透视表汇总和可视化的格式。在看了前面的 df_new.head()之后，我认为我们不会使用下面的列，这导致了 15 列的较小集合:
df_new = df_new.drop(columns=['YearEnd','LocationDesc','DataSource','DataValue','DataValueFootnoteSymbol','DatavalueFootnote','LowConfidenceLimit','HighConfidenceLimit','GeoLocation'])<class 'pandas.core.frame.DataFrame'>
RangeIndex: 403984 entries, 0 to 403983
Data columns (total 15 columns):
YearStart 403984 non-null int64
LocationAbbr 403984 non-null object
Topic 403984 non-null object
Question 403984 non-null object
DataValueUnit 374119 non-null object
DataValueType 403984 non-null object
DataValueAlt 273666 non-null float64
StratificationCategory1 403984 non-null object
Stratification1 403984 non-null object
LocationID 403984 non-null int64
TopicID 403984 non-null object
QuestionID 403984 non-null object
DataValueTypeID 403984 non-null object
StratificationCategoryID1 403984 non-null object
StratificationID1 403984 non-null object
dtypes: float64(1), int64(2), object(12)
表 1:按主题、问题和值汇总
使用数据透视表(很像 Excel)，让我们总结每个问题及其相关值。为了通知每个问题和值，我们将从每个问题的类别(主题)、问题 ID、问题、数据单元(数据值单元)以及与稍后的分层相关的数据类型(数据值类型)创建多个索引。Pandas Pivot_table 需要数值，默认为平均值(numpy.mean)。此外，我已经删除了 LocationID 和 YearStart 列，因为我以后不会使用它们，并且为了可读性，将其四舍五入为两位小数:
df_QD = df_new.pivot_table(index=['Topic','QuestionID','Question','DataValueUnit','DataValueType'],columns=None,dropna=True)df_QD.drop(columns=['LocationID','YearStart']).round(2).head(25)
Table 1. Pivot Table showing the Question and Values
表 2:使用带有分层、问题和值的 Groupby】
除了 pivot_table，我还可以使用 groupby()创建一个类似的汇总表。我想在生成的数据帧中展示更多关于分层的信息。
df_unittype_table1 = df_new.groupby(['Topic','QuestionID','Question','StratificationID1','Stratification1','DataValueUnit','DataValueType']).mean().round(2)df_unittype_table1.drop(columns=['YearStart','LocationID'])
Table 2. Using groupby() to summarize stratification
表 3:按问题和地点汇总
一种有趣的方法是用数据值按每个状态(LocationAbbr)显示数据。下面的数据透视表包括额外的索引，这些索引为诸如单位和类型之类的数字提供了上下文。
df_new_qloc = df_new.pivot_table(values='DataValueAlt',index=['Topic','QuestionID','Question','DataValueUnit','DataValueType'], columns='LocationAbbr',aggfunc='mean',dropna=True).round(2)
Table 3. Pivot Table of Question and Location
可视化所有指标之间的相关性
现在，我们已经看到这些表格按问题、分层和位置呈现数据值。那还是很多信息啊！我们先来了解一下问题指标本身说明了什么。每对问题指示器之间有什么关系？这就是数据可视化将描绘出一幅理解整体关系的图画的地方。
最初，我使用表 3 中的 pivot_table，但是轴标签太长了，很难阅读。稍微倒退一下，我们将基于新的列问题创建一个不同的 pivot_table (df_new_qloc2 ),取问题的前 37 个字符:
df_new['QuestionAbbr'] = df_new['Question'].str[:37]df_new_qloc2 = df_new.pivot_table(values='DataValueAlt',index=['Topic','QuestionID','QuestionAbbr'], columns='LocationAbbr',aggfunc='mean',dropna=True).round(2)
通过转置 df_new_qloc2 数据透视表，我们准备应用关联方法。corr()并将数据可视化。使用来自 seaborn 文档的代码，我们可以绘制一个相关矩阵热图。这使每一对指标可视化，以了解正相关对位于何处。由于相关矩阵会在对角线的上方和下方产生重复的视图，为了简单起见，我们将通过使用 np.zeros_like()创建一个零数组并使用 triu _ indices _ from()返回对角线上方区域的索引来屏蔽热图的上半部分。使用 matplotlib，我们设置了一个大的图形尺寸，以允许在 jupyter 笔记本上放大。cmap 为图形设置色彩映射表(更多选项此处)。最后，sns.heatmap()使用遮罩呈现热图。
# Table 4 - Using Table 3 with question indicators as columns
new_qloc2_corr = df_new_qloc2.transpose().corr()# Generate a mask for the upper triangle
mask = np.zeros_like(new_qloc2_corr, dtype=np.bool)
mask[np.triu_indices_from(mask)] = True# Set up the matplotlib figure
f, ax = plt.subplots(figsize=(150, 150))
f.suptitle('All Chronic Indicators Correlation Heatmap', x=0.4,y=0.85,fontsize=150)
ax.tick_params(labelsize=50)# Generate a custom diverging colormap
cmap = sns.diverging_palette(190, 10, as_cmap=True)# Draw the heatmap with the mask and correct aspect ratio
sns.heatmap(new_qloc2_corr, mask=mask, cmap=cmap, vmax=.3, center=0,square=True, linewidths=2, cbar_kws={"shrink": .3})
Figure 1. All Indicators Correlation Heatmap (Higher Correlation in Pink)
即使做了一些调整，信息也是密集的，可以更方便用户阅读(我们将在下一篇文章中解决这个问题)。尽管如此，热图轴显示了主题、问题 ID 和缩写问题 QuestionAbbr。我们可以看到一些粉红色的区域正相关程度更高，绿色的区域负相关程度更高。从视觉上看，心血管疾病(CVD)和癌症(CAN)的某些领域似乎具有更高的相关性。看到这个演示很有帮助，但最终，我们希望通过创建一个顶级相关性列表来获得更细粒度的外观。
表五。指标对顶部关联表
仅从 df_new 中取出 QuestionID、LocationAbbr 和 DataValueAlt，我创建了 new_qloc2_1_corr 表。使用线性代数，我们可以通过将 new_qloc2_1_corr 相关数据帧和与先前数据帧形状相同但值为 1 的数组的下半部分相乘来对数组进行矢量化。的。stack()将数组放入一个列表中。sort_values()将列表按降序排序。Top_corr 由大约 36k 项组成。
top_corr = (new_qloc2_1_corr * np.tril(np.ones(new_qloc2_1_corr.shape),-1)).stack().sort_values(by=['DataValueAlt'],ascending=False)
通过几个额外的步骤来重置索引，更改列名，并根据 QuestionID 添加几个新列，我们就有了一个包含顶部相关对的数据框架。由于主题有许多重复对，我将删除重复项，不包括主题类别的相同指标对，这样会产生大约 33k 个项目:
top_corr_pair_dd = []for row in range(len(top_corr_pair)): if top_corr_pair.loc[row]['Topic1'] != top_corr_pair.loc[row].
['Topic2']:
top_corr_pair_dd.append(top_corr_pair.loc[row])
Table 5. Top Correlation by Topic Descending
在这一点上，我们看到了指标中出现的主题。表 5 向我们展示了基于主题类别的特定模式。例如，慢性肾病(CKD)和心血管疾病(CVD)中的那些问卷(QID1 和 QID2 列)；慢性肾病和糖尿病患者(DIA)；CVD 和 DIA 内的那些；总体条件内的指标(OVC)和先前指标；最后，患有 DIA 和 CVD 慢性阻塞性肺病(COP)的患者。
总的来说，我们看到 CVD、CKD 和 DIA 等指标之间的关系，这有意义吗？
心血管疾病包括许多与心脏相关的疾病，包括心律失常(不正常的不规则心脏跳动)、高血压(对动脉壁的高作用力)、心脏骤停(功能突然丧失)、冠状动脉疾病(主要血管受损)和充血性心力衰竭(心脏泵血不畅的慢性疾病)。这是一种非常常见的疾病，源于血管堵塞，影响着四分之一的美国人。据 CDC 称，增加这种疾病的风险因素包括饮食、肥胖、糖尿病、过度饮酒和缺乏锻炼。
慢性肾脏疾病是肾脏及其功能丧失导致衰竭的长期疾病，也称为终末期肾脏疾病(ESRD)。这也是一种常见的疾病，影响着七分之一的美国人。根据梅奥诊所的说法，增加这种疾病的风险因素包括糖尿病、高血压、吸烟、心血管疾病、老年和某些种族背景。在肾病患者中，心血管疾病倾向于“诊断不足和治疗不足”(S，Hernandez GT。慢性肾脏疾病和心血管疾病之间的联系。 J 肾病学。2014;3(3):99–104.).根据美国肾脏基金会的说法，当肾脏的功能没有达到最佳状态时，这需要心脏和心血管系统更加努力地工作。
糖尿病与体内葡萄糖失调有关，包括一些疾病，如 1 型、2 型和妊娠糖尿病。虽然 1 型糖尿病的原因不明，如家庭背景和环境，但妊娠期糖尿病通常伴随着妊娠而来。二型糖尿病与上述疾病有相似的风险因素，包括缺乏锻炼、年龄、家庭背景、高血压、体重、多囊卵巢综合征和胆固醇/甘油三酯水平异常(来源:梅奥诊所)。患有糖尿病会导致各种并发症，包括心血管和慢性肾脏疾病。
答案是肯定的！该数据集与文献一致，表明心血管疾病、慢性肾病和糖尿病之间存在高度相关性。
在本系列的最后一篇博文中，我们将继续深入了解以下内容:
- 首要条件是什么？
- 具体的顶级指标还有哪些亮点？
- 按年份和分层我们能看到什么趋势？
希望你喜欢读到这里，并随时给我一些想法和建议！
相对与绝对:如何进行成分数据分析？第二部分
这是我早先关于成分数据分析的文章的延续，在那篇文章中，我展示了将成分数据视为绝对数据而不是相对数据的缺陷。在这篇文章中，我将总结我们可以用来正确分析成分数据的技术，并用具体的例子展示使用 RNA-Seq 数据。
有两种主要策略用于处理成分数据，特别是 NGS 数据:
1。归一化以取回绝对计数
2。成分数据分析(CoDA)方法，在样品参考范围内转换数据(例如:ALR、CLR)
绝对计数标准化
当需要跨样本比较时，这是 NGS 数据预处理中最广泛使用的技术。相对读取计数被“标准化”为总读取深度，以“恢复”绝对计数。然而，当 RNA 或细胞的总绝对量或相关生物材料的量在样品间显著变化时，这不能恢复绝对计数。这通常会给分析人员带来错误的安全感，并导致将这些“标准化”样品视为绝对计数。在样本间进行比较时，这可能会导致错误的结论。让我们用模拟数据来证明这一点。
模拟细节
这里，我模拟了 100 个基因的数据，其中
a. 5 个基因在对照和实验条件之间具有 1 的真实对数倍数变化(在选择下接近耐受或不生长)；
b . 2 个基因在实验条件下具有与> 1 相同的真实对数倍数变化(在选择下耐受并表现出生长)；
c . 2 个基因在实验条件下具有与< 1 相同的真实对数倍数变化(不耐受或耐受)，
我模拟了其余 91 个基因的不同比例发生变化的 5 种不同情况。在发生变化的基因中，约 90%被耗尽，约 10%被富集。
其他基因的缺失/富集影响相对计数值和读取深度标准化计数，即使总读取深度固定在 200K 读取
读取深度标准化(RDN)计数
- **差异表达或丰度:**即使所有读数具有相同的总深度(计数之和)，使用读数深度标准化计数(RDN 计数)计算的基因对数倍数变化(LFC)与真实的对数倍数变化相比有所偏移(见下图 1)。有趣的是，转变的方向并不总是可以根据基因改变的比例来预测的。例如，当大约 70%的基因改变时，使用 RDN 计数计算的 LFC 与真实的 LFC 相比下移。另一方面，当 90%的基因改变时，使用 RDN 计数计算的 LFC 与真实 LFC 相比向上移动。这是因为前一种情况下的绝对真实计数高于后一种情况。一般来说，我们无法预测或估计样本的真实总绝对计数。
Fig 1: Comparing True Log-Fold Changes to Log-Fold Changes Calculated using RDN Counts
**2。基因之间的相关性:**看事物如何在相对计数和绝对计数之间进行比较。，我计算了所有 5 个样本中非恒定基因的相关性(每个样本都有 0.1、0.2、0.4、0.6、0.9 个比例的变化基因)。在 200K 读取深度下，我使用了真实计数和使用聚酯模拟计数数据的相对计数。
Fig 2: Comparing True Correlations Between Genes to Correlations Calculated using RDN Counts
从上图可以看出，使用 RDN 计数计算的一些相关系数与真实相关系数有显著差异，存在负偏差。
上面的两个例子显示了使用 RDN 计数来估计基因之间的差异表达或相关性的缺陷。当试图从相对组成数据中恢复绝对计数时，应始终使用加标控制，而不是使用 RDN 计数。接下来我们将展示这一点
峰值标准化计数
为了真正校正绝对计数的变化，我们需要在测序步骤之前以相同的丰度(数量)添加到所有样品中的掺入对照或基因。这样做可以将所有样本标准化到相同的总丰度级别，并使比较正确。这仅在数据由于测序而关闭时有效(因为我们是在测序之前添加掺入物)，如果约束是生物学的或者发生在测序步骤的上游，这将没有帮助。在这种情况下，我们需要在这个约束步骤之前添加尖峰，但是由于添加尖峰的物理和生物限制，这并不总是可能的。
让我们用我们的数据来看看这是如何工作的。在我们的数据中，我们有 92 个不同的对照或掺入基因，它们具有真正的绝对丰度。让我们使用这些来“标准化”数据，从而使所有样本达到相同的绝对计数范围。
- **差异表达或丰度:**下面的图 3 类似于图 1，但是用掺入标准化计数代替 RDN 计数。该图添加了人工抖动(噪声)来显示所有数据，但真实数据都位于对角线上。这表明了插队的力量。如果刚好在导致数据闭合或约束的步骤之前添加尖峰信号，则适当设计的尖峰信号可以恢复绝对计数(达到常数乘法因子)，但这并不总是可能的。
Fig 3: Comparing True Log-Fold Changes to Log-Fold Changes Calculated using Spike-in Normalized Counts
2.**基因之间的相关性:**查看基因之间的相关性，我们看到使用尖峰标准化计数计算的系数可以恢复真实的系数。下图 4:
Fig 4: Comparing True Correlations Between Genes to Correlations Calculated using Spike-in Normalized Counts
看来我们找到了解决问题的方法。我们所要做的就是添加一些控制，我们很好！不幸的是，没那么快。在这个模拟的案例中，组成数据的闭合来源是测序，我们能够在模拟测序数据之前添加一些控制。在现实生活的数据生成过程中，闭包的来源可以出现在提取 DNA/RNA 的通常复杂的工作流程中的任何地方。此外，生物系统本身可能是固有组成的(例如，细胞产生 RNA 的能力是有限的)，在这种情况下，细胞外引入的尖峰信号不能恢复真实的绝对计数。
成分数据分析(CoDA)方法
一种替代掺入归一化的方法是使用 CoDA 方法，这种方法通常根据样本内参考转换计数数据。加法对数变换(ALR)和中心对数变换(CLR)是一些常用的尾波变换的例子。这些方法最初是由 John Aitchison 在 1986 年提出的。核心思想是组件相对于另一个参考组件的对数比变换可以被视为任何其他无约束数据。这将数据从原始的单纯形空间(如我们在第一部分中的三元图)转换到欧几里得空间。这使得我们可以对这些数据使用所有经典的分析技术。
注意:这些技术并不像前一节中的“规范化”方法那样要求开放数据。这些技术也适用于所有数据，无论是相对数据还是绝对数据。另一点需要注意的是，使用尖峰信号进行归一化与使用加性对数比(ALR)变换是一样的。使用一般 ALR 变换的好处是，即使我们没有样品间丰度恒定的掺入，它也是适用的。一般 ALR 变换的缺点是，我们需要正确选择基准电压源，以理解数据并回答相关问题。现在让我们使用和以前一样的数据集来更详细地看看 CoDA 方法。
**1。差异表达或丰度:**有许多方法可以发现处理前后成分数据的变化。令人惊讶的是，这些方法中有许多来自微生物组文献，而基因表达文献主要依赖于传统方法，如 DESeq2 和 EdgeR ，这些方法没有明确考虑数据的组成性质。DESeq2 和 EdgeR 隐含地假设绝对丰度不会因处理而改变。这相当于使用 CoDA 方法的中心对数比(CLR)转换。这种转化使用基因或组分的几何平均值作为参照，因此所有结果都必须根据几何平均值来解释。在这个阶段，很容易将这种假设解释为基因的几何平均值在对照组和治疗组之间没有变化。也许几何平均数改变了，也许没有，除了来自测序的相对计数之外，没有正交信息，没有办法确切知道。DESeq2 和其他差异表达工具的大多数用户都落入了这个陷阱，并认为算法调用的任何重大变化都意味着绝对计数的重大变化。相反，这些只是相对于所有组件的几何平均值的显著变化。
有一些新兴的方法将统计学的严谨性应用于成分数据中的 DA。最流行的方法是 ALDEx2 和 ANCOM 。这些方法的主要原理是依靠相对于参考成分的转换相对数据的对数比测试，并仔细解释这些结果。这些方法的主要问题是只能根据所选的参比来解释结果，并且没有提供如何选择参比的指导。朱利亚诺·克鲁兹给我介绍了一种更近的方法，这种方法使用了差异排名(DR)并提出了一种更合理的选择参考的方法。这是我将在这里简要使用的，并希望在将来的某个帖子中深入到运行这些算法的血淋淋的细节。
DR 的主要思想是选择一些随机参考成分来计算治疗和对照中所有成分的对数比。在下一步中，这些成分按照它们在处理和控制条件之间的差δ(对数比)的顺序排列。使用已知相对计数计算的δ(对数比)值的等级顺序应该与使用未知真实绝对计数计算的δ(对数比)值的等级相同。例如，对于 90%的基因差异表达的情况，我在下面显示了使用相对计数计算的δ(对数比)值与使用绝对计数计算的δ(对数比)值:
Fig 5: Δ(log-ratio) values Calculated using Absolute vs. Relative Counts
如您所见，δ(对数比)值的大小因我们使用相对计数还是绝对计数而异，但δ(对数比)值的等级保持不变。这并不意味着排名靠前的基因在治疗中的数量高于对照组，而排名靠后的基因数量较低。可能发生的情况是，与对照条件相比，在治疗条件下，排名靠前的基因已经耗尽了绝对计数，但是在治疗条件下，排名靠后的基因的耗尽甚至更糟。简而言之，我们不能说治疗和病情之间的绝对读数有什么变化。
我现在将选择排名第一的基因作为我的参考，并使用这个新的参考再次计算δ(对数比)值。
Fig 6: Δ(log-ratio) Values Calculated using the Top-Ranking Gene as Reference
从该图中，我们可以使用 0.5 的任意截止值，并选择超过该值的任何基因作为我们的潜在 DA 基因进行进一步测试。当然，如果我们想测试更多的基因，我们可以放宽截止日期。
另一个避免选择参照的建议是在人群中建立某种阳性或阴性对照。假设，我们知道一个在治疗条件下绝对丰度会增加的基因，那么我们可以使用这个基因作为计算对数比的天然参考，并对δ(对数比)值进行排序。任何大于 1 的对数比意味着该基因优于阳性对照，而小于 1 的对数比意味着比阳性对照差。更好的是，有 2 个对照来限制效应大小，并参照这两个基因来解释对数比。
在我的模拟中，我每次复制只有一个样本，因此无法进行任何统计分析。在以后的文章中，我将为每个条件生成多个重复，并使用 ALDEx2、ANCOM 和 DR 算法来测试它们的灵敏度和特异性。
2.**基因之间的相关性:**如本系列第 1 部分所示，相关性不是亚组成连贯的，因此不遵循 CoDA 的原则之一。简而言之，任何两个基因之间的相关系数取决于数据中存在的其他成分或基因。Aitchison 首先提出使用对数比转换值的方差(VLR)来估计两个变量的相关性。例如，为了计算跨 n 个样本的特征或基因 g 和 h 之间的相关性，我们将使用:
VLR 是子成分相干的，因此不会导致伪相关。使用 VLR 的主要问题是，即使当基因 g 和 h 完全相关时，它等于 0，但当基因完全独立时，它没有上限。由于这个比例问题，很难将 VLR 的一个基因对与 VLR 的另一个基因对进行比较。基于 VLR 提出了几种方法/度量来估计组合之间的依赖关系，最著名的是 SparCC、SPIEC-EASI 和比例。在这篇博客中，我只是稍微详细地回顾一下比例。所有这些方法都试图使用 VLR 来导出类似于相关系数的度量，因此可以在不同的成分对之间进行比较。
基于洛弗尔等人的工作，在 R 包 propr 中提出了三个基于比例的度量标准。艾尔。和奎恩等人。艾尔。这些指标是根据对数转换数据计算的。定义见 propr 包。下面的 Ai 是指数据中基因或成分“I”的对数转换值。Ai 可以是绝对计数或相对计数，这些定义仍然适用。
- phi(φ)= var(Ai-Aj)/var(Ai)
- ρ(⍴)= var(ai-aj)/(var(ai)+var(aj))
- phis(φs)= var(Ai-Aj)/var(Ai+Aj)
与传统相关系数最接近的度量是 rho，范围从-1 到 1。phi 是无界的，可以从 0 变化到 Inf，phis 是 phi 的对称变体，是 rho 的单调函数。在博客的其余部分，我将重点关注 rho。
a .使用绝对计数 :如果我们有一个尖峰控制，我们可以从相对计数恢复绝对计数。由于我们已经有可用的峰值数据，我将使用峰值转换后的数据计算 rho 值，即使用峰值 TPM 计数作为归一化计数的基因 1 的 A = log(TPM _ counts/Spike _ TPM _ counts)。这将恢复原始绝对计数。现在，我们可以使用上面的等式计算ρ值。我将绝对计数和 rho 值之间的相关性绘制如下:
Fig 7: Correlations of True Absolute Counts Between Genes vs. Rho Values Calculated using Spike-Normalized Data
从该图中可以看出，使用比例关系，我们可以捕获真实绝对计数之间的大部分原始相关性。当然，这是一个虚构的例子，我们有一个很好的 spike-in 来检索绝对计数。即使有了这个人为的例子，我们仍然可以看到真实的相关性和使用峰值归一化计数计算的比例值之间的一些差异。这是由于基于比例的度量的计算方式使得它们对对数变换值的方差估计极其敏感。这里我们只有 5 个样本来计算方差，在大多数情况下，前 3 个样本具有相同的值。我怀疑这导致了计算分解指标的公式。我得多琢磨琢磨这个。这一假设的证据是，如果我们只寻找在至少 4 个不同样本中具有不同值的分量，则相关值和比例度量匹配得相当好，如下所示:
Fig 8: Correlations of True Absolute Counts Between Genes vs. Rho Values Calculated using Spike-Normalized Data: Only for Genes with At Least 4 Distinct Values out of 5 Samples
一般而言，rho 和其他基于比例的度量具有良好的精度和较差的召回率，并且具有更多样本给出了对方差的更好估计，因此给出了 rho 值的更好估计。此外，引导通常用于建立一个重要关系的界限。例如，在上面的图中，重要“rho”值的临界值可能是 0.75。
b .使用相对计数 :这个世界复杂得不公平，不幸的是，我们身边通常没有一个好的替代物可供我们使用😢。因此，我们必须使用相对数据，或者更具体地说，加性对数转换(ALR)相对数据。或者，如果我们确信样本间计数的几何平均值不会改变，我们可以使用中心对数变换(CLR)(我们知道这不适用于此处的模拟数据)。本质上，在这种情况下，我们能做的最好的事情就是计算相关数据之间的关系。因此，让我们比较相对数据(相对于选择的参考基因)的 rho 值与真实绝对计数之间的相关性。下图显示了使用 2 个随机选择的参考基因计算的相对计数之间的相关性:
Fig 9: Correlations of True Absolute Counts Between Genes vs. Rho Values Calculated using 2 randomly chosen reference genes
从上面的图中可以明显看出，参考基因的选择对我们称之为成比例的基因对产生了巨大的影响。我们无法事先知道如何选择比例值的截止值，以确保精确度和召回率之间的良好平衡。Aitchison 提出的 VLR，他的公式在这篇文章的前面已经展示过了，不会有这样的问题，因为当我们取一个基因对中的两个基因的对数比时，参考被抵消了，但是我们不能用 VLR 来比较基因对，因为它没有一个可解释的尺度。不幸的是，我们又回到了起点。这在数学上用 Erb 和 Notredame 来解释。另一种选择是不去考虑比例本身，而是计算不同条件下的微分比例。要做到这一点，我们只需要比较不同条件下的基因对，VLR 就足够了，我们不必担心缩放。这篇文章已经太长了，这可能是另一篇文章的讨论。此外，现在我已经看到了基于比例的度量的问题，其他方法，如 SparCC 似乎是一个不错的选择，至少尝试一下。再说一次，也许是为了另一篇文章。总的来说，对于成分数据，基于比例的度量提供了一种有趣的替代基于相关性的度量，但在缺乏良好的加标控制的情况下解释这些数据是很棘手的。
外卖食品
总之，在没有加标控制的情况下，没有从相对数据中恢复绝对计数的灵丹妙药。然而，学习 CoDA 方法帮助我更好地理解解释成分数据的细微差别。我认为这类似于贝叶斯统计和频率统计，前者迫使从业者以先验的形式列出所有的假设。在 CoDA 方法的情况下，选择参考迫使用户警惕对数据的错误解释。
你可以在这里找到这篇文章和之前的文章的所有代码和数据:https://github.com/kimoyerr/CODA
相对与绝对:用模拟理解成分数据
我在工作中经常被要求做的一件事是比较不同样本的 ngs 数据。这些可能是复制，也可能来自完全不同的环境或处理步骤。我们都直觉地理解这些数据总是受到测序深度的限制，因此需要仔细地得出结论。直到我看到 @ WarrenMcG 的这篇论文的预印本，我才缺乏分析这些数据的词汇或技术。这开始了几天对成分数据的深入挖掘，我生成了一些合成数据，以便更好地理解其中的细微差别。在这篇文章中，我试图展示我所学到的东西，并希望有助于传播将组成数据视为常规无约束数据的危险。
什么是成分数据
成分数据是指数字构成整体比例的任何数据。例如，你将一枚硬币抛 10 次，得到的正面和反面的数量就是成分数据。这些数据的主要特征是它们是“封闭的”或受限的，如果我告诉你我们通过投掷双面硬币得到了 6 个正面，你就知道有 4 个反面，而我没有给你提供任何进一步的信息。从概率的角度来说，这意味着数据中的所有成分并不是相互独立的。因此，组成数据不存在于欧几里得空间中，而是存在于一个称为单形的特殊约束空间中。对于那里的工程师来说，每次绘制多相平衡图时，你可能会遇到单纯形图，就像这里使用 R 库三元绘制的水的多相平衡图一样:
本质上，一个特定的固定量的水可以在不同的条件下存在于不同的相中，但是蒸汽、液态水和冰的总量应该总是等于我们开始时的固定量。成分数据存在于比我们想象的更多的地方。一些例子包括投票数据，食品标签上的营养信息，你呼吸的空气中不同气体的比例，甚至你屏幕上像素的 RGB 值。
将成分数据视为无约束数据的问题
根据 Aitchison ，一般来说，任何成分数据分析方法都应满足 4 个不同的标准。传统的分析通常不能满足所有这些要求，因此会导致错误的结果。我将用模拟数据来说服我自己，希望你能理解传统分析方法的这些缺陷。我使用absimseq生成一些虚拟的 RNA-seq 之类的数据。简而言之，我有 4 个来自控制条件的样本和 4 个来自处理条件的样本。每个样本由 100 个天然基因和 92 个刺突(对照)基因组成，样本中所有基因的计数总计约 20K 个读数。与对照条件相比，模拟处理条件具有 20%的差异表达的天然基因(用于描述细胞/组织内基因拷贝数变化的奇特生物学术语)。差异表达的确切大小与我们在此试图证明的无关。既然我们已经建立了数据，让我们继续讨论传统分析方法在处理成分数据时通常无法考虑的 4 个标准。
你可以在这里找到我使用的所有数据和代码:https://github.com/kimoyerr/CODA.git
比例不变性
组成数据的一个共同属性是它们可以用不同大小的向量来表示。空气由 78%的氮气、21%的氧气和 1%的剩余体积组成。在向量符号中，这可以写成[78，21，1]。就体积的百万分率(ppmv)而言，这也可以指定为[780840、209460 和 9780]。尽管这些符号非常不同，但它们代表了相同的物理样本。然而，传统的分析方法，如欧几里德距离和层次聚类，会被这种表示上的差异所迷惑。让我们用模拟的 RNA-Seq 数据来证明这一点。在 R 中，执行以下操作:
load('sim_counts_matrix_100.rda')
orig.dist <- dist(t(counts_matrix))orig.dendro <- as.dendrogram(hclust(d = dist(t(counts_matrix))))
# Create dendro
dendro.plot <- ggdendrogram(data = orig.dendro, rotate = TRUE)# Preview the plot
print(dendro.plot)
Clustering of original data
正如所料，我们看到四个对照样品(样品 1 至 4)聚集在一起，四个处理样品(样品 5 至 8)聚集在一起。
现在让我们缩放一些不同的样本，如下所示。例如，如果您对一些样品进行了比平时更高深度的测序，但仍想一起分析所有样品，就会发生这种情况。
# Scale the samples differently
scaled.counts_matrix = counts_matrix %*% diag(c(1,1,5,5,1,1,5,5))
scaled.dist <- dist(t(scaled.counts_matrix))scaled.dendro <- as.dendrogram(hclust(d = dist(t(scaled.counts_matrix))))
# Create dendro
scaled.dendro.plot <- ggdendrogram(data = scaled.dendro, rotate = TRUE)# Preview the plot
print(scaled.dendro.plot)
Clustering of scaled data
新的树状图现在显示了聚集在一起的缩放样本(样本 3、4、7 和 8 ),这是不正确的，因为唯一改变的是这些样本的读取总数，但它们的单个组件相对于彼此是相同的。
扰动不变性
这与我们的分析结论不应依赖于用来表示成分数据的单位的观察结果有关。回到空气成分的例子，我们可以混合使用百分比和 ppmv 值来表示数据，比如:[78%、21%和 9780ppmv]。这仍然代表着相同的物理样本，我们应该能够将它与另一个类似空气的样本进行比较，在这个新的坐标系和旧的坐标系中，并得出相同的结论。让我们再次使用模拟数据来计算新的扰动坐标空间中样本之间的距离:
#Multiply each row (gene) with a scalar value
perturbed.counts_matrix = counts_matrix * c(seq(1,192,1))
colnames(perturbed.counts_matrix) = colnames(counts_matrix)
perturbed.dist <- dist(t(perturbed.counts_matrix))perturbed.dendro <- as.dendrogram(hclust(d = dist(t(perturbed.counts_matrix))))
# Create dendro
perturbed.dendro.plot <- ggdendrogram(data = perturbed.dendro, rotate = TRUE)# Preview the plot
print(perturbed.dendro.plot)
Clustering of perturbed data
从上面的图可以清楚地看出，与比例不变性问题相比，这是一个小得多的问题。来自对照和处理条件的样本聚集在一起，但是样本之间的距离不同。这可能会导致错误的结论，具体取决于所提的问题。在大多数 RNA-seq 或 NGS 数据中，在不同条件下比较样品时，扰动方差不是大问题。
亚位置相干性
这一特性确保了无论包含哪些数据成分，分析都会得出相同的结论。例如，在我们的空气示例中，我们可以删除代表 1%体积的组件。现在样本在 2 维空间([78.5%，21.5%])中表示，而不是在原始的 3 维空间([78%，21%，1%])。在这个低维空间中的任何分析都不应不同于在原始三维空间中的比较。从下一代测序(NGS)数据的角度来看，这是一个重要的属性，在这种情况下，我们并不总是能够保证始终找到数据中的所有成分。这可能有很多原因，有时是有意的，有时是无意的，不可避免的。
为了在传统分析方法中模拟这种特性(或缺乏这种特性)，我从由 100 个基因+ 92 个对照组成的原始数据集生成了由前 50 个基因+ 92 个对照组成的第二个数据集。较小数据集中的 50 个基因具有与原始数据集中相同的绝对(无约束)值。然后，通过使所有读取的总和为大约 20K(如在原始数据集中一样),这些不受约束的数据被封闭(受约束)。
计算两个数据集中 50 个常见基因之间的相关系数，然后进行比较:
load('sim_counts_matrix_100.rda')
counts.all <- counts_matrix
# Load the sub-compositional data made up of only the first 50 genes (features) + 92 controls from the original data of 100 genes (features) + 92 controls
load('sim_counts_matrix_50.rda')
counts.sub.comp <- counts_matrix# Get the correlation between the 50 common genes
cor.all <- as.vector(cor(t(counts.all[1:50,])))
cor.sub.comp <- as.vector(cor(t(counts.sub.comp[1:50,])))
tmp <- as.data.frame(cbind(cor.all,cor.sub.comp))
names(tmp) <- c('correlation_all', 'correlation_sub_comp')
tmp$abs.diff <- as.factor(ifelse(abs(tmp$correlation_all - tmp$correlation_sub_comp)>0.5,1,0))# Plot
ggplot(tmp,aes(correlation_all,correlation_sub_comp, color=abs.diff)) + geom_point(size=2) + th + scale_colour_manual(values = c("1" = "Red", "0" = "Blue")) + theme(legend.position = "none")
Differences in the correlation between genes (features) depending on what other genes are present in the data
在上面的图中，我们可以看到两个数据集中的大多数相关系数是相似的，有些系数有显著差异(以红色显示)。根据所提问题的不同，这些差异可能会导致结论的显著差异。这可能是基因网络分析中的一个主要问题。
置换不变性
这一特性意味着不管原始数据的顺序如何，分析方法都应该得出相同的结论。我能想到的大多数分析方法一般都遵循这个特性。请让我知道，如果你知道任何违反这一点的分析方法，我会在帖子中包括它们。
如何正确分析成分数据
现在，我们确信传统的分析方法并不总是满足分析成分数据的所有 4 个重要属性，我们能做什么？三十多年来，Aitchison 做了一些令人难以置信的工作，想出了更好的技术。根据我处理 NGS 数据的经验，我感觉他的建议在很大程度上被当成了耳旁风。希望这篇博文能够说服其他人在分析成分数据时更加谨慎。在下一篇文章中，我将展示一些分析成分数据的方法，再次强调 NGS 数据和应用。
https://quotecites.com/quote/Hermann_Hesse_378
不懈追求 1%:分析部门的角色。
当我与希望购买我的团队服务的营销人员交谈时，或者他们希望建立自己的分析团队时，我经常告诉他们，一个合适的分析团队就是要在他们的营销组织中安装智能。这并不意味着适用于每个营销部门——有些部门很小，没有专家也能很好地分析他们的数据——但是一旦一家公司达到他们需要专业化的规模，你就有几个方向可以选择。在本帖中，我们将探讨营销人员可以投资的三个主要领域，以提高他们的营销智能能力。
智能:进化，而不是革命。
智力是一个有趣的词。很美，也很有争议。智能是一种选择的工具，它引导人类主宰食物链，保护我们的家庭，发展文明，并创造我们茁壮成长的地方来做刺激我们大脑的事情。(作者注:在这篇博文中，我主要关注智力的积极方面。)
那么什么是智能，为什么它对可持续、可实现的业务增长如此重要？智能是利用数据获取有关营销环境的知识并利用这些知识进行调整的能力。是动作驱动；它导致适应；这是任何企业为了保持竞争优势都必须履行的基本职能。这也不是我们(作为决策者)能够在一夜之间实现的。智能依赖于准确的信息才能智能。一个根据错误信息采取行动的营销部门，最终会因为缺乏真正的环境知识而受到伤害。
智能是利用数据获取和应用知识的能力。这是每个企业保持竞争优势所需要的。
强大的情报能力始于干净的数据。我们可以信赖的数据和捕捉营销环境重要方面的数据。对于我领导的团队来说，干净有用的数据来自强大的数据运营团队的规划、协调和工程。
强大的情报能力始于干净的数据。
数据操作:获取干净的数据。
从强大的分析能力开始的稳步增长之旅始于良好的数据运营实践。在我工作的团队中，数据运营可以是数据工程师、产品经理和信息技术团队代表的组合。
如果营销组织很小，也可以简单到拥有合适的工具和基础设施。
可靠数据运营的目标是识别一个组织拥有的数据、我们需要的数据，并创建一个计划，以建立一个关于公司、客户和该组织所处的营销环境的可信的真实来源。
要了解您是否需要投资数据运营，请问自己这样一个问题:如果我正在查看的数据与我的想法和决策相矛盾，我会改变主意吗？我有信心数据是好的吗？如果答案是否定的，那么让我们通过看一看我们的数据操作来开始智能进化。
分析和统计:在不确定性下改变我的想法的科学
如果数据操作正常，下一个主要投资是分析师和统计学家，他们可以使用方法来创建关于我们收集的数据的有用见解。这些人有各种各样的口味，所以让我把他们放在两个桶里，冒犯了他们中的大多数。
第一桶:数据挖掘专家/故事讲述者/灵感来源。
第二桶:统计学家/智囊团/风险——秀逗魔导士。
这种过于简化的做法非常有用，因为它有助于我们理解分析团队中的一些基本角色。这两个组拥有独特的技能，您将在非常具体的用例中求助于它们。
数据挖掘专家正在探索数据，以帮助激励营销决策者。他们在我们已经知道的东西中找到洞察力。他们的工具是性能报告、数据可视化和分析。这种风格的分析师充满好奇，是一个伟大的故事讲述者。他们有着卓越的商业头脑，只是在他们小心翼翼的洞察力发展过程中黯然失色。这是您在监控营销绩效时最先接触到的角色之一。
统计学家与数据挖掘者相关，但拥有不同的技能。他们专注于具有不确定性的决策，并将这种不确定性的风险降至最低。他们会用数学来帮助你理解你可能采取的不同决策路径中的风险。利用大脑互动体的最佳时机是当你有关于你要做的决定的极好的数据时——但是你有几个选择。
数据科学家和技术专家:构建认知解决方案
一旦你掌握了数据挖掘，并通过统计将风险最小化，现在是时候创建商业解决方案，利用你正在收集的数据的力量。这将把您的分析能力扩展到您提供的产品和服务中——注入它以解决特定的问题，或者用更智能/个人/相关的产品取悦您的客户。
数据科学家和软件工程师—这两个角色可以帮助您构建数据驱动的解决方案。
认知解决方案是能够创造和应用知识的商业解决方案。
数据科学家/软件工程师组合非常强大。如果你能把这两个角色放在同一个房间里——由可靠的数据操作支持——那么你就有了一个正在形成的业务解决方案发电站。
仍然需要一些关于如何发展你的营销智慧的指导。在下面的评论区给我一个欢呼吧；或者查看一下我的其他博客，寻找一些灵感。
用空间和自然语言处理重温复仇者联盟 3：无限战争
发现地球上最强大的英雄的台词中的顶级名词、动词、实体和文本相似性
经过漫长的一年等待，复仇者联盟 4：终局之战终于来了。我和你一样，世界上大多数人都会在第一天赶到电影院看电影，体验复仇者联盟如何拯救世界，结束一个十年的故事。为了平静我的神经和缓解等待，我想重温以前的电影，无限战争，但不同的和互动的。而且，因为我是一个数据专家，当然，它必须涉及数据和一些流行词汇。
答案？自然语言处理，简称 NLP。
使用 spaCy ，一个 NLP Python 开源库，旨在帮助我们处理和理解大量文本，我分析了电影的脚本，以研究以下概念:
- 电影中的十大动词、名词、副词和形容词。
- 由特定角色说出的顶级动词和名词。
- 电影中的前 30 个命名实体。
- 每个字符对的台词之间的相似性，例如，托尔和灭霸的台词之间的相似性。
在本文中，我将讨论并展示我的发现，同时用代码解释我是如何用 spaCy 做到这一点的。
你对代码和技术词汇不感兴趣吗？今天是你的幸运日！我想说的是，我在这里使用的词汇和术语大多是非技术和用户友好的，所以即使你没有自然语言处理、人工智能、机器学习或在这里插入流行语的经验，你也应该能够掌握我想要传达的主要思想和概念。所以，请随意忽略这些代码:)
The Mad Titan. Credits: Marvel
处理数据
用于实验的数据或文本语料库——在 NLP 中通常被称为——是电影的剧本，可在这个链接获得。然而，在使用这些数据之前，我必须对其进行清理。因此，我删除了一些不必要的东西，如描述一个动作或场景的评论，如*“[灭霸粉碎宇宙魔方，露出蓝色的太空石……]”，以及说这句话的角色的名字(实际上，这个名字是用来知道谁说了什么，但不是用于分析的实际语料库的一部分)。此外，作为空间数据处理步骤的一部分，我忽略了标记为停止词的术语，换句话说，即常用词，例如【我】【你】【安】。另外，我只使用了引理，这是每个单词的标准形式。例如，动词 "talks】、、、、、【talking】、是同一个词位的形式，其引理是、【talk】、*。
要在 spaCy 中处理一段文本，首先，我们需要加载我们的语言模型，然后在文本语料库上调用该模型。结果是一个Doc对象，一个保存已处理文本的对象。
Creating a Doc object in spaCy
现在我们有了一个干净的经过处理的语料库，是时候开始了！
十大动词、名词、副词和形容词
仅仅通过看动词就能知道电影的整体动作或情节是什么吗？本文的第一张图解决了这个问题。
“I know”, “you think” are some of the most common phrases
【知道】【走】【来】【得到】【想】【告诉】【杀死】【需要】【停止】**【想要】。从中我们可以推断出什么？因为我看过这部电影几次——这也暗示我有偏见——我愿意根据这些动词得出结论，复仇者联盟 3：无限战争是关于知道、思考和研究如何去阻止某事或某人。
这就是我们如何获得带有空间的动词:
描述动词的单词，即副词呢？
“I seriously don’t know how you fit your head into that helmet” — Doctor Strange
对于一部关于阻止一个人摧毁半个宇宙的电影来说，口语副词中有很多实证主义——“对”、“完全正确”和“更好”等词就是例子。
所以，我们知道了动作，以及它们是如何被描述的，现在是时候看看名词了。
“You will pay for his life with yours. Thanos will have that stone.” — Proxima Midnight
看到第一个结果是*【stones】并不奇怪，毕竟电影讲的就是他们。在第二个位置，我们有术语“生命”，，这是灭霸想要摧毁的东西，接下来是“时间”，恰恰是复仇者联盟没有的(注意:“时间”*也可能是因为提到了时间石)。
最后，我将用形容词或者描述名词的单词来结束这一节。与副词类似，我们有一些表达肯定的术语，如*、、、【对】、，还有一些表达肯定的术语，如、、、【肯定】、*。
“I’m sorry, little one.” — Thanos
特定人物提到的顶级动词和名词
之前，我们看到了电影中提到的最常见的动词和名词有哪些。虽然这些知识让我们对电影的整体感觉和情节有所了解，但它并没有说太多关于我们角色的个人冒险经历。因此，我应用了与查找前十名动词和名词相同的过程，但是是在字符级别上。
由于电影中有许多角色，我只选择了一些实际上说了合理数量台词的角色，加上一些我最喜欢的:)。这些角色是托尼·史塔克、奇异博士、加莫拉、托尔、火箭、彼得·奎尔(星际领主)、乌木·莫和灭霸。抱歉，队长，你没入选。然而，在我的 GitHub(链接在文章末尾)中，你可以找到一个文件夹，里面有每个角色的图形。
接下来的图片显示了这些人物使用的顶级名词。
What’s with Quill calling Drax so much?
我发现很奇怪，甚至令人耳目一新的是，在大多数情况下，我们亲爱的英雄们使用的顶级名词是提及船员。比如托尼说了*【小子】*九次(指蜘蛛侠)，火箭叫了奎尔(星主)三次，而奎尔本身就叫了(更像是冲着)德拉克斯七次。
通过进一步的考察，我们可以推断出对每个角色来说什么似乎是最重要的。就钢铁侠而言，数据表明地球对他来说很有价值。与他相似的是加莫拉，他总是在思考更高的目标——、、、、、【行星】、——并最终为此付出代价。奇异博士还有另一个目的——保护他的宝石——他反复提到了这个目的。然后，还有托尔，他对灭霸有个人恩怨，说了八次他的名字，还有一个新的兔子最好的朋友。最后，还有疯狂的泰坦本身，灭霸，他无法停止思考收集无限宝石，或者他的女儿。
虽然名词富有表现力和意义，但动词却不是这样。正如你将在下一张图中看到的，动词不像名词那样丰富多彩。像*、、、【想要】、、、这样的词占据了大部分的榜首。然而，有一个角色可能拥有整个语料库中最独特的动词:乌木·莫。如果你不记得，乌木的血盆大口，是灭霸的头号心腹。就像他是一个好仆人一样，他的目标是——除了得到时间之石——宣扬他主人的使命，这个任务他用了诸如、和“欢呼”这样的词来完成*蠕变。
“Hear me, and rejoice. You have had the privilege of being saved by the Great Titan…” — Ebony Maw
此外，以下是格鲁特的热门名词。
“I am Groot”
命名实体
到目前为止，我们已经探索了我们的英雄和反派在这部史诗电影中使用的最常见的动词、名词、副词和形容词。然而，为了完全解释我们一直在仔细研究的所有这些词，我们需要一些上下文，即命名实体。
我引用 spaCy 网站上的话，一个命名实体是一个*“被赋予名称的真实世界的物体——例如，一个人，一个国家，一个产品或者一本书的名字。”*所以，知道实体，就意味着，意识到人物所谈论的事物。在 spaCy 包中，实体有一个预测标签，它将实体分类为许多类型中的一种，如人、产品、艺术作品等等(https://spacy.io/api/annotation#named-entities)，为我们提供了额外的粒度级别，这对于进一步分类它们很有用。不幸的是，出于简单的原因，我不会使用实体类，而只是实体本身。
这些是排名前 30 的实体。
“MAYEFA YA HU” is the chant of Wakanda’s Jabari warriors
首先，我们有灭霸，考虑到这部电影全是关于他的，这一点也不奇怪。跟随他的是他的女儿加莫拉，影片的中心人物之一。然后在第三个位置，我们有格鲁特(我肯定我没有解释为什么)，其次是托尼，其他复仇者，以及一些地方，如纽约，阿斯加德和瓦坎达[永远]。除了英雄和地点，两块*【六】*(见 14 号实体)无限之石——时间和灵魂石——出现在这个名单上(分别是第 15 和第 16 位)。令人惊讶的是，将灭霸带到地球的心灵之石不在名单上。
要访问空间中的实体，请读取Doc的属性ents，如下所示:
每个字符对的口语行之间的相似性
当我们讨论每个角色的顶级动词时，我们意识到，与名词不同，大多数动词都非常相似，并且传达相同的感觉。像*【走】【来】这样的术语给我们一种运动的印象或我们的角色想要去或到达某个特定地方的感觉，而像【杀死】【停止】*这样的动词暗示着确实有一个巨大的威胁必须被阻止。
带着这个想法，为了进一步研究相似性的概念，我计算了每个角色对的台词之间的相似性分数。
NLP 中的相似性概念描述了两段文本的合成或句法意义有多接近或相对——通常，相似性得分的范围是从 0 到 1，其中 0 表示完全不相似，1 表示完全相似(或者两个文本相同)。从技术上讲，相似度是通过测量单词向量之间的距离来计算的，即单词的多维表示。对于那些有兴趣了解这个主题的人，我推荐搜索 word2vec，这可能是用于生成这些单词嵌入的最常见的算法。下图显示了相似性矩阵。
Again, Ebony Maw is the most unique character
说实话，我没想到会是这个结果。一方面，我可以接受所有的相似性接近于 1，因为毕竟，电影有一个主要情节，它应该有跨对话的关联。然而，让我感到不安的是，分数非常非常相似。我的意思是，看看灭霸的分数；我希望看到一个恶棍的分数与那些想要阻止他的人的分数相比有很大的不同。从更积极的角度来看，有趣的是，彼得·帕克(蜘蛛侠)的得分变化最大。毕竟，他只是一个在混乱中被抓住的孩子，所以，有这个结果是意料之中的。
这是一个如何计算空间中两个Doc之间相似性的例子:
总结和结论
在电影《复仇者联盟 3：无限战争》中，我们跟随一群超级英雄踏上了阻止灭霸的旅程，这个威胁的目标是毁灭宇宙一半的生命。通过这部电影，我们了解到大多数英雄都有拯救世界的动机和动力，这反映在他们表达自己的方式上。在本文中，我们在 Python、NLP 和 spaCy 的帮助下，探索我们的英雄和反派如何通过研究他们的每一句台词来表达和交流。通过关注他们最常用的动词和名词等特征，我们学习、确认并重温了钢铁侠对地球的忠诚、奇异博士保护时间之石的宣誓职责、托尔复仇的渴望以及灭霸完成自己命运的雄心。
复仇者联盟周快乐，等残局剧本出来后再见。
在我的 GitHub 上可以找到用于生成这个实验的代码，链接如下:
在 GitHub 上创建一个帐户，为 juandes/infinity-war-spacy 开发做贡献。
github.com](https://github.com/juandes/infinity-war-spacy)
如果你有任何问题，评论，疑问，或者想聊天，请在这里或 Twitter 上留下评论，我将很乐意帮助你。
[## 胡安·德迪奥斯·桑托斯(@ jdiosantos)|推特
胡安·德迪奥斯·桑托斯的最新推文(@ jdiossantos)。机器学习/数据工程师。还有，口袋妖怪大师，还有…
twitter.com](https://twitter.com/jdiossantos)
啊哦…我变成了一只恐龙
用这些技能保持相关的金融专业。
据金融高管国际公司称，
【两年前】78%的首席财务官认为精通 Excel 是他们 FP & A 团队最重要的技能。今天只有 5%的人有同样的感觉。相反，首席财务官将对新技术的适应能力视为新员工的首要技能。
科技正在改变一切。它改变了传统的工作，创造了新的角色，并带来了对更新技能的需求——包括当前的 FP & A 专业人员。根据哈克特集团(Hackett Group)2018 年的数字技能调查，在不断变化的人才需求方面，金融业远远落后于形势，目前似乎没有一种跟踪供需趋势的既定方法。
如今，FP&A 的员工需要掌握以前根本不存在的新技能，否则就有被更年轻的新员工超越的风险。进入金融科技专业人士的时代——在金融、技术和 IT 领域拥有强大背景的个人。虽然今天这种知识受到高度重视，但在不久的将来，我们会看到它从重视转变为必需。我们应该期待在短短几年内看到财务政策和会计职位的角色和组成发生巨大变化。所需的后端技能无疑会发生变化。更多的数据科学家和有分析思维能力的个人将是必要的，我们将看到以前的金融巨星变得过时。
虽然今天这种知识受到高度重视，但在不久的将来，我们会看到它从重视转变为必需。
在 EY 对金融领导者的一项调查中，近 80%的受访者表示，他们不认为自己具备必要的技能来满足 10 年后的工作需求。这应该是一个紧急的警钟，大声喊出迫切需要的技能提升。
既然是这样，我们决定涵盖未来的 FP&A 的技术技能。如果你想让自己成为一名有吸引力的员工，并确保自己不会过时，可以考虑学习下面我们介绍的技能。为了在财务政策和会计领域保持相关性，这里有一些工具和技巧可以帮助你达到并保持这种状态。
“越来越多的公司和职位描述将要求财务职位描述中的大数据[和] 建模经验”——Geetanjali tan don，AFP FP &咨询委员会，拜耳作物科学数字& IT 转型财务主管
数据和技术敏锐度=职业生存的关键
理想情况下，FP&A 专业人员应该对大量数据感到舒适，理解数据是如何构造的，理解集成是如何发生的，并且能够有效地使用包括云应用在内的多个系统。
编程能力
能够做一些简单的编码或者理解使用一个系统是必要的，以便用不同的技术改进流程是一种优势。与 FP & A 专业人士最相关的编码语言有:
Java
Python
c#
Murex
由于金融科技员工的需求越来越大，以下是与每个角色相关的编程语言:
区块链专家和开发者: C、C++、Java
App 开发者: C#、C++、Java、Python
数据科学家: 数据科学工具包:R、Python、Weka、NumPy、MatLab
数据可视化工具:D3.js、ggplot
查询语言使用熟练程度:SQL、Hive、Pig
NoSQL 数据库:MongoDB、Cassandra、
有很多在线学习平台让学习编码语言变得容易和有趣。Coursera 就是这样一个平台，它提供的课程从、【面向所有人的 Python】、到、【使用 SQL 的现代大数据分析】你所要做的就是留出一些时间，带着开放的心态来。
数据可视化和分析技能 大数据工具执行数据分析，以便从大型数据集中提取洞察。FP &专业人士可以从使用可视化来解释他们的数据中受益匪浅。学习如何使用领先的分析工具可以帮助你发展你的可视化和分析技能。
你应该熟悉的一个可视化工具是 Power BI 。微软的业务分析工具 Power BI 提供交互式可视化功能，帮助各部门创建报告和仪表板。向管理层展示所有最新的财务报告，以便他们做出最佳、最明智的决策。
考虑看看 Udemy 的微软 Power BI 入门课，了解这个平台，熟悉它的特性。
财务专业人员将继续需要在分析技能和技术悟性方面具有适应性、灵活性和技能提升。
熟悉技术 为了增加自己的市场价值，FP &专业人士应该熟悉一系列与行业相关的技术。这类技术包括 SPSS、Excel、SQL、SAS、R、MatLab、Linux、Hadoop 等。与金融敏锐度相结合，熟悉这些技术对 FP & A 专业人士来说是一个巨大的优势。
此外，对自动化的日益依赖促使更多的招聘经理寻找具备适当 IT 技能的金融专业人士来利用新的金融系统。能够证明自己在预测分析或应付账款自动化等领域的知识和能力的候选人会发现自己很受欢迎。
根据毕马威的说法，“没有改变的是金融的基础作用。财务部门仍需要为组织的其他部门提供洞察力，确保有效的控制和风险管理，并推动自身和组织的效率……这些目标(预计)将保持不变。然而，不同的是在每一项上花费的时间和精力的比例以及交付方式。”
一个大大改善了 FP&A 专业人员工作方式的自动化平台是 DataRails 。DataRails 通过自动整合基于 Excel 的数据，消除了传统手动流程中的繁重工作。这节省了以前浪费在繁琐的手动工作上的大量时间。通过最大限度地将时间花在分析数字上，而不是浪费时间将它们放在一起，以高超的洞察力打动管理层。
“自动化将继续让财务和财务与会计专业人员真正反思他们的工作方式，减少他们在日常和重复的手动密集型工作上花费的时间，并腾出更多时间作为重要的业务合作伙伴。金融专业人士将继续需要在分析技能和技术悟性方面具有适应性、灵活性和技能提升。但有如此多的资源可以提供帮助，所以现在不应该是担心或恐惧的时候，相反，对于那些相信终身学习的专业人士来说，这应该被视为一个巨大的机会。”
- Rachele Collins 博士，APQC 首席研究员
用 Keras 估计剩余寿命
从时间序列到图像…问 CNN“下一次故障什么时候发生？”
在许多人工智能解决方案中，预测罕见事件正成为研发的重要课题。生存分析、客户流失预测、预测性维护和异常检测是处理罕见事件的最受欢迎的应用领域的一些例子。考虑到这些情况，我们可以将罕见事件想象为在特定条件下发生的特定状态，不同于正常行为，但在经济利益方面起着关键作用。
在这篇文章中，我开发了一个机器学习解决方案来预测特定发动机组件的剩余使用寿命( RUL )。这类问题在预测性维护领域起着关键作用，这里的目的是说’ 距离下一次故障还剩多少时间？’。为了实现这个目标，我在 Keras 中开发了一个卷积神经网络，以图像的形式处理时间序列。
数据集
对于数据科学家来说，在处理这种任务时，最重要的问题是缺乏可用的观察形式的罕见事件。因此，实现良好性能的第一步是尝试处理最丰富的数据集，以处理各种可能的情况。
由 NASA 提供的涡扇发动机退化模拟 数据集 ，正在成为一个同类型(共 100 台)发动机机群剩余使用寿命( RUL )估算的重要基准。数据以时间序列的形式提供:3 个操作设置、21 个传感器测量值和周期，即工作寿命的时间观察值。
发动机在每个时间序列开始时正常工作，在序列中的某个点出现故障。在训练集中，故障在数量上增长，直到系统故障。在测试集中，时间序列在系统故障前的某个时间结束。目标是预测测试组中故障前的剩余运行循环数，即发动机将继续运行的最后一个循环后的运行循环数。
为了更好地理解这个解释，我们试着看一下数据:
train_df.id.value_counts().plot.bar()
Life Cycles for each engine before a fault in train set
发动机有不同的寿命。列车数据中的平均工作时间为 206 个周期，最少 128 个周期，最多 362 个周期。
在单个发动机的列车设置中，绘制了操作设置和传感器测量值，亲爱的:
engine_id = train_df[train_df['id'] == 1]engine_id[train_df.columns[2:]].plot(subplots=True, sharex=True, figsize=(20,30))
settings 1–3, sensors 1–9 for engine1
sensors 10–21 for engine1
绘图总是一个好主意……这样，我们可以对我们所掌握的数据有一个印象深刻的总体概述。在大部分系列的结尾，我们可以观察到一个不同的行为，这预示着未来的失败。
准备数据
为了预测每台发动机的 RUL，我们采用了一种分类方法，通过以下方式自行生成标签:
从 0(故障)到 15 个剩余周期，我们标记为 2；从 16 到 45 个周期，我们标记为 1，其余的(> 46)标记为 0。很明显，在现实场景中，标记为 2 的类别是最有经济价值的。预测该类的良好性能将允许我们运行一个适当的维护程序，避免未来的故障并节省资金。
为了让我们能够处理列车的最大数量的数据，我们使用固定窗口和 1 步滑动来分割序列。例如，引擎 1 具有 192 个训练周期，窗口长度等于 50，我们提取长度为 50 的 142 个时间序列:
*窗口 1 - >从周期 0 到周期 50，窗口 2 - >从周期 1 到周期 51，…，窗口 142 - >从周期 141 到周期 50，窗口 191。*每个窗口都标有该窗口所考虑的最终循环的相应标签。
sequence_length = 50def gen_sequence(id_df, seq_len, seq_cols):data_matrix = id_df[seq_cols].values
n_elem = data_matrix.shape[0]
for a,b in zip(range(0,n_elem-seq_len), range(seq_len,n_elem)):
yield data_matrix[a:b,:]
def gen_labels(id_df, seq_len, lab):
data_matrix = id_df[lab].values
n_elem= data_matrix.shape[0]
return data_matrix[seq_len:n_elem,:]
从时间序列到图像
为了让事情变得更有趣，我决定将这个系列转化成图像。以便用它们来填充我们的分类模型。
我根据这个惊人的资源创建了这些图片。这个概念很简单…当我们试图将时间序列转换成图像时，我们总是利用声谱图。这个选择很聪明，但并不总是最好的选择(你可以在这里阅读)。在这篇文章中，作者解释了他对用声谱图表示处理音频系列的困惑。他谈到了声音，但意思可以在我们的场景中翻译。光谱图是强大的，但它们的使用可能会导致信息的丢失，特别是如果我们试图用计算机视觉的方式来解决问题。为了有效，2D CNN 需要空间不变性；这建立在一个假设上，即一幅经典图像(比如一张照片)的特征无论在哪里都有相同的含义。另一方面，声谱图意味着由两个不同的单位(频率和时间)构成的二维表示。
出于这些原因，我决定利用循环图来转换我的时间序列窗口(长度为 50 个周期)。它们很容易用 python 实现，只需几行代码，利用 Scipy。
from scipy.spatial.distance import pdist, squareformdef rec_plot(s, eps=0.10, steps=10):d = pdist(s[:,None])
d = np.floor(d/eps)
d[d>steps] = steps
Z = squareform(d)
return Z
使用这个函数，我们能够为我们处理的每个时间序列生成一个 50x50 的图像(我已经排除了方差为 0 的恒定时间序列)。因此，每一次观察都是由一组大小为 50x50x17 (17 是没有零方差的时间序列)的图像组成的，如下所示。
Example of a train observation
模型
至此，我们已经准备好构建我们的模型了。我采用了经典的 2D CNN 架构:
model = Sequential()
model.add(Conv2D(32, (3, 3), activation='relu', input_shape=(50, 50, 17)))
model.add(Conv2D(32, (3, 3), activation='relu'))
model.add(MaxPooling2D(pool_size=(2, 2)))
model.add(Dropout(0.25))
model.add(Conv2D(64, (3, 3), activation='relu'))
model.add(Conv2D(64, (3, 3), activation='relu'))
model.add(MaxPooling2D(pool_size=(2, 2)))
model.add(Dropout(0.25))
model.add(Flatten())
model.add(Dense(256, activation='relu'))
model.add(Dropout(0.5))
model.add(Dense(3, activation='softmax'))model.compile(loss='categorical_crossentropy', optimizer='adam', metrics=['accuracy'])
我只对 10 个时期进行了拟合，达到了 0.832%的精度。
从混淆矩阵中我们可以看到，我们的模型可以很好地区分发动机何时接近故障(2 个标签:<16 cycles remaining) or when it works normally (0 label: > 45 个循环)。一点点噪音出现在中级(> 15，<46 cycles). We are satisfied to achieve a great and clear result for the prediction of class 2 — i.e. near to failure.
SUMMARY
In this post, we try to solve a 预测性维护问题。估算发动机的 RUL 由于收集这类数据的难度，我们有意识地应对罕见事件。我们提出了一个有趣的解决方案，利用循环图将时间序列转换成图像。通过这种方式，我们能够区分处于工作寿命末期的油井发动机。
保持联系: Linkedin
遥感基础:归一化差异植被指数
卫星图像在生态学研究中的应用
NDVI visualization of continental US (Source: GIS Geography)
如果你还没有，请查看我之前的帖子，它总结了我在大会数据科学沉浸式课程中的顶点项目: 巴西亚马逊地区的土地使用和森林砍伐 。它很好地介绍了我在机器学习、遥感和生态学方面的一些兴趣，也是这篇文章的一个很好的序言。
那么，我的下一步是什么？
我认为这个帖子是对自己的承诺。我承诺将继续寻找方法来建立我的工程和机器学习能力，并推进对生态系统的研究和理解。
我深深地投入到学习这些工具和我们可以用它们做的一切，并且被那些既推进技术又增加科学家使用它进行保护的公司所鼓舞。他们的努力让我能够继续追求我的兴趣。
归一化植被指数
卫星图像的一个应用是归一化差异植被指数(NDVI)。这是一个基本但非常强大的指数，可以在任何包含近红外颜色通道的卫星图像中计算。该指数将反射的近红外光与反射的可见红光进行比较，并通过以下等式简单计算:
The scale for NDVI is -1 to 1, with values near -1 representing water and those near 1 representing dense and healthy vegetation.
植被指数因影像内容而异。这是因为陆地物体对不同波长的光有不同的反应(见我关于卫星图片内容的技术贴)。植被反射大量的近红外光，因此植被密度越大，图像的植被指数就越大。
**这对卫星图像分类中的目标或土地类型识别有影响。**一个想法是计算图像中每个像素的 NDVI，并使用这些值作为图像中的附加颜色通道。这可能是传递到神经网络的额外信息层，并可能增加网络检测图像边缘的能力。例如，检测图像从密集到稀疏的植被、城市景观到乡村环境、或者森林到清晰的农业的区域。
NDVI 的另一个应用是用来跟踪生态系统随时间的变化。我们可以通过计算植被指数随时间的变化来衡量林业对生态系统的影响。现在，我们可以每天访问整个地球的图像，我们可以通过植被指数等指标来跟踪气候变化的影响。这对于理解不同尺度的气候变化的影响是有价值的。植被指数的变化将因地方和地区而异。这可以将保护工作引向变化较大的地区，或者有助于管理和规划工作。
野火对加州植被指数的局部影响
我想在更高的水平上使用植被指数来构建工具，以进一步丰富我在数据工程和数据科学方面的经验。在完成我的顶点项目之前，我对使用卫星图像来了解 2018 年加州野火对周围植被的影响产生了兴趣。直觉上，被野火烧毁的地区经历了植被密度的剧烈变化。这些变化将反映在 NDVI(没有双关语的意思)。利用遥感，我们可以跟踪植被随时间的变化，并回答诸如“生态系统从野火中恢复需要多长时间？”。
Satellite Image from USGS of the Camp Fire in November 2018: Paradise, California
这可能适用于最近的加州火灾，以及前几年的火灾。感谢像 Planet 的 Open California 这样的开放访问数据程序，我们可以使用 2015 年的图像立即构建这些应用程序。
我很兴奋去探索！
用异常检测去除拉曼光谱中的尖峰:Python 中的 Whitaker-Hayes 算法
Detecting and removing spikes from Raman spectra
摘要。在这篇文章中，我想对 D. Whitaker 和 K. Hayes 在《化学计量学和智能实验室系统杂志》上提出的一种去除拉曼光谱中尖峰信号的新方法进行评论。在他们的出版物中，作者使用了一种基于改进的 Z 分数异常值检测的算法来定位这种尖峰，如果存在的话，随后使用简单的移动平均来消除它们。他们不是计算光谱强度的 Z 值，而是计算一次差分光谱的 Z 值。他们还提供了用 r 实现算法的代码。在这里，我用 Python 实现了算法，并展示了如何将它应用于不同的数据集。
拉曼光谱是一种广泛使用的分析技术，它提供来自分子和固体的结构和电子信息。它适用于实验室和大规模生产规模，并在许多不同的领域，如物理，化学，生物，医学或工业应用。
拉曼光谱学中已知的一个典型问题是拉曼光谱有时被尖峰“污染”。尖峰是出现在频谱上随机位置的正的窄带宽峰值。它们起源于高能宇宙射线撞击用于测量拉曼光谱的电荷耦合器件探测器。这些尖峰是有问题的，因为它们可能会妨碍后续分析，特别是在需要多变量数据分析的情况下。因此，处理拉曼光谱数据的第一步是清除尖峰。
首先，加载所需的 Python 包:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
图 1 显示了石墨烯的拉曼光谱。在过去的几年里，石墨烯由于其卓越的物理性能，包括优异的电子、热、光学和机械性能，已经成为非常受欢迎的材料。其特征拉曼光谱由如图所示的几个峰组成。从它们的形状和相关强度，可以获知掺杂、应变或晶界等大量信息。该光谱是被尖峰污染的光谱的明显例子。
# load the data as data frame
df = pd.read_csv(‘folder_name.../spectrum.txt’, delimiter = ‘\t’)# Transform the data to a numpy array
wavelength = np.array(df.Wavelength)
intensity = np.array(df.Intensity)# Plot the spectrum:
plt.plot(wavelength, intensity)
plt.title(‘Spectrum’, fontsize = 20)
plt.xticks(fontsize = 15)
plt.yticks(fontsize = 15)
plt.xlabel(‘Wavelength (nm)’, fontsize = 20)
plt.ylabel(‘Intensity (a.u.)’, fontsize = 20)
plt.show()
Figure 1. A spike was recorded while measuring the typical Raman spectrum of graphene (characterized by the G and 2D bands).
基于 Z 分数的拉曼光谱尖峰检测方法
尖峰强度通常高于光谱中其他拉曼峰的强度，因此使用基于 z 分数的方法可能是一个好的起点。z 得分以标准差为单位表示一个值离平均值有多远。因此，如果已知总体均值和总体标准差，则原始得分 x(i)的标准得分计算如下:
z(i) = (x(i)-μ) / σ
其中 μ 是平均值*，σ* 是总体 x 的标准偏差(x(i)代表单个拉曼光谱的值)。关于如何使用 Z 分数方法检测异常值的更多详细信息可以在参考文献[3]中找到。
让我们计算光谱中点的 z 分数:
def z_score(intensity):
mean_int = np.mean(intensity)
std_int = np.std(intensity)
z_scores = (intensity — mean_int) / std_int
return z_scoresintensity_z_score = np.array(z_score(intensity))
plt.plot(wavelength, intensity_z_score)
plt.xticks(fontsize = 15)
plt.yticks(fontsize = 15)
plt.xlabel( ‘Wavelength’ ,fontsize = 20)
plt.ylabel( ‘Z-Score’ ,fontsize = 20)
plt.show()
Figure 2. Z-scores of the spectrum plotted in figure 1.
然后需要一个阈值来判断一个值是否是异常值。该阈值的典型值为 3.5，由美国质量控制协会提出作为异常值标记规则的基础，而该出版物的作者使用 6。为了应用阈值来排除峰值，必须采用 Z 得分的绝对值:
|z(i)| = |(x(i)-μ) / σ|
计算为
def z_score(intensity):
mean_int = np.mean(intensity)
std_int = np.std(intensity)
z_scores = (intensity — mean_int) / std_int
return z_scoresthreshold = 3.5intensity_z_score = np.array(abs(z_score(intensity)))
plt.plot(wavelength, intensity_z_score)
plt.plot(wavelength, threshold*np.ones(len(wavelength)), label = ‘threshold’)
plt.xticks(fontsize = 15)
plt.yticks(fontsize = 15)
plt.xlabel( ‘Wavelength’ ,fontsize = 20)
plt.ylabel( ‘|Z-Score|’ ,fontsize = 20)
plt.show()
Figure 2b. Z-scores of the spectrum plotted in figure 1. The threshold still cuts the 2D Raman peak.
然而，z 分数非常类似于原始光谱，并且阈值仍然切断主要拉曼峰。让我们使用 3.5 的阈值来绘制检测到的峰值:
threshold = 3.5# 1 is assigned to spikes, 0 to non-spikes:
spikes = abs(np.array(z_score(intensity))) > thresholdplt.plot(wavelength, spikes, color = ‘red’)
plt.title(‘Spikes: ‘ + str(np.sum(spikes)), fontsize = 20)
plt.grid()
plt.xticks(fontsize = 15)
plt.yticks(fontsize = 15)
plt.xlabel( ‘Wavelength’ ,fontsize = 20)
plt.ylabel( 'Z-scores > ' + str(threshold) ,fontsize = 20)
plt.show()
**Figure 3.**Using modified Z-Scores, 18 spectral point are above the threshold.
用于拉曼光谱中尖峰检测的改进的基于 Z 分数的方法
第二种选择是利用稳健统计并计算光谱的修正 z 分数。这种改进的 Z 得分方法使用中位数(M)和中位数绝对偏差(MAD ),而不是平均值和标准差:
z(i) = 0.6745 (x(i)-M) / MAD
其中 MAD = median(|x-M|)和|…|代表绝对值。中位数和 MAD 分别是集中趋势和离差的稳健度量。乘数 0.6745 是标准正态分布的第 0.75 个四分位数，MAD 将 to⁴.收敛到该四分位数
def modified_z_score(intensity):
median_int = np.median(intensity)
mad_int = np.median([np.abs(intensity — median_int)])
modified_z_scores = 0.6745 * (intensity — median_int) / mad_int
return modified_z_scoresthreshold = 3.5intensity_modified_z_score = np.array(abs(modified_z_score(intensity)))
plt.plot(wavelength, intensity_modified_z_score)
plt.plot(wavelength, threshold*np.ones(len(wavelength)), label = 'threshold')
plt.xticks(fontsize = 15)
plt.yticks(fontsize = 15)
plt.xlabel('Wavelength (nm)' ,fontsize = 20)
plt.ylabel('|Modified Z-Score|' ,fontsize = 20)
plt.show()
Figure 4. Modified Z-scores of the spectrum plotted in figure 1. The threshold still cuts the 2D Raman peak.
仍然可以观察到相同的问题:
# 1 is assigned to spikes, 0 to non-spikes:
spikes = abs(np.array(modified_z_score(intensity))) > thresholdplt.plot(wavelength, spikes, color = ‘red’)
plt.title(‘Spikes: ‘ + str(np.sum(spikes)), fontsize = 20)
plt.grid()
plt.xticks(fontsize = 15)
plt.yticks(fontsize = 15)
plt.xlabel( ‘Wavelength’ ,fontsize = 20)
plt.ylabel( 'Modified Z-scores > ' + str(threshold) ,fontsize = 20)
plt.show()
Figure 5. Using modified Z-Scores, 24 spectral points are above the threshold.
2D 拉曼峰值仍然被检测为尖峰，因此需要更灵敏的方法。
Whitaker 和 Hayes 的 改进的基于 Z 分数的拉曼光谱尖峰检测方法
Whitaker 和 Hayes 建议利用尖峰的高强度和小宽度，因此使用连续光谱点之间的差异∇x(i) = x(i)-x(i-1 来计算 z 分数，其中 x(1)，…，x(n)是以等间距波数记录的单个拉曼光谱的值，i = 2，…，n。该步骤消除了线性和缓慢移动的曲线线性趋势，而尖锐的细尖峰将被保留。现在，
z(i) = 0.6745 (∇x(i)-M) / MAD
因此只包括一个额外的步骤，其中包括连续值之间的差。
# First we calculated ∇x(i):
dist = 0
delta_intensity = []
for i in np.arange(len(intensity)-1):
dist = intensity[i+1] — intensity[i]
delta_intensity.append(dist)delta_int = np.array(delta_intensity)
# Alternatively to the for loop one can use:
# delta_int = np.diff(intensity) intensity_modified_z_score = np.array(modified_z_score(delta_int))
plt.plot(wavelength[1:], intensity_modified_z_score)
plt.title('Modified Z-Score using ∇x(i)')
plt.xticks(fontsize = 15)
plt.yticks(fontsize = 15)
plt.xlabel('Wavelength (nm)', fontsize = 20)
plt.ylabel('Score', fontsize = 20)
plt.show()
Figure 6. Modified Z-scores using ∇x(i) of the spectrum plotted in figure 1.
同样，为了应用阈值来排除尖峰，必须采用修改后的 Z 得分的绝对值:
|z(i)| =|0.6745 (∇x(i)-M) / MAD|
导致
threshold = 3.5intensity_modified_z_score = np.array(np.abs(modified_z_score(delta_int)))
plt.plot(wavelength[1:], intensity_modified_z_score)
plt.plot(wavelength[1:], threshold*np.ones(len(wavelength[1:])), label = 'threshold')
plt.title('Modified Z-Score of ∇x(i)', fontsize = 20)
plt.xticks(fontsize = 15)
plt.yticks(fontsize = 15)
plt.xlabel('Wavelength (nm)', fontsize = 20)
plt.ylabel('Score', fontsize = 20)
plt.show()
Figure 7. Abosolute value of the modified Z-scores of ∇x(i) of the spectrum plotted in figure 1.
同样，检测到的尖峰数量可以计算如下
# 1 is assigned to spikes, 0 to non-spikes:
spikes = abs(np.array(modified_z_score(intensity))) > thresholdplt.plot(wavelength, spikes, color = ‘red’)
plt.title(‘Spikes: ‘ + str(np.sum(spikes)), fontsize = 20)
plt.grid()
plt.xticks(fontsize = 15)
plt.yticks(fontsize = 15)
plt.xlabel( ‘Wavelength’ ,fontsize = 20)
#plt.ylabel( ‘Spike(1) or not(0)?’ ,fontsize = 20)
plt.show()
Figure 8. Using modified Z-Scores with ∇x(i), 17 spectral points are above the threshold = 3.5.
对于 3.5 的推荐阈值，分配了许多错误的尖峰。然而，与拉曼峰相比，通过这种方法得到的值比以前高得多。
plt.plot(wavelength[1:],
np.array(abs(modified_z_score(delta_int))), color='black', label = '|Modified Z-Score using ∇x(i)|')
plt.plot(wavelength, np.array(abs(modified_z_score(intensity))), label = '|Modified Z-Score|', color = 'red')
plt.plot(wavelength, np.array(abs(z_score(intensity))), label = '|Z-Score|', color = 'blue')
plt.xticks(fontsize = 15)
plt.yticks(fontsize = 15)
plt.xlabel( 'Wavelength' ,fontsize = 20)
plt.ylabel( 'Score' ,fontsize = 20)
plt.legend()
plt.show()
Figure 9. Comparison between the three different approaches.
一般来说，必须根据数据集选择合适的阈值。对于这种情况，阈值= 7 已经足以获得清晰的选择。
Figure 10. Absolute values of the modified Z-scores with ∇x(i) of the spectrum plotted in figure 1.
Figure 11. Using modified Z-Scores with ∇x(i), 3 spectral points are above the threshold = 7.
固定拉曼光谱
一旦检测到尖峰信号，下一步就是移除它们并修正光谱。为此，通过计算每个候选波数的近邻的平均值(在 2m+1 个值窗口内)，在每个候选波数处获得插值。
def fixer(y,m):
threshold = 7 # binarization threshold.
spikes = abs(np.array(modified_z_score(np.diff(y)))) > threshold
y_out = y.copy() # So we don’t overwrite y for i in np.arange(len(spikes)):
if spikes[i] != 0: # If we have an spike in position i
w = np.arange(i-m,i+1+m) # we select 2 m + 1 points around our spike
w2 = w[spikes[w] == 0] # From such interval, we choose the ones which are not spikes
y_out[i] = np.mean(y[w2]) # and we average their values
return y_out# Does it work?
plt.plot(wavelength, intensity, label = 'original data')
plt.plot(wavelength, fixer(intensity,m=3), label = 'fixed spectrum')
plt.title('Despiked spectrum',fontsize = 20)
plt.xticks(fontsize = 15)
plt.yticks(fontsize = 15)
plt.xlabel('Wavelength (nm)' ,fontsize = 20)
plt.ylabel('Intensity (a.u.)' ,fontsize = 20)
plt.legend()
plt.show()
Figure 12. Original and despiked Raman spectrum using modified Z-Scores with ∇x(i) with the threshold = 7.
因此，在计算了∇x 的修改的 z 分数，并且通过设置适当的阈值来进行阈值处理之后，通过应用移动平均滤波器来去除和平滑尖峰。
其他例子
最后，给出了两个具有不同强度信号的尖峰信号的例子，显示了这种方法的优点。
Figure 13. Examples of two more spectra in the first row of graphs, their Z-scores (blue lines) , modified Z-scores (red lines) and modified Z-scores using ∇x(i) (black lines) in the second row of graphs and thei despiked spectra in the third row, are presented.
**承认。**我要感谢豪尔赫·路易斯·希塔和安娜·索拉古伦-比斯科花时间校对这些笔记。
[1]惠特克，达伦 a .，和凯文海斯。"消除拉曼光谱尖峰的简单算法."化学计量学和智能实验室系统179(2018):82–84。
[2] Iglewicz 和 D. Hoaglin。" ASQC 质量控制基本参考:统计技术."如何检测和处理异常值16(1993):1–87。
[3] Colin Gorrie 的数据故事:检测离群值的三种方法。http://colingorrie.github.io/outlier-detection.html
[4]若昂·罗德里格斯。“离群值让我们发疯:单变量离群值检测”。中等。https://medium . com/James-blogs/outliers-make-us-go-mad-单变量-离群值-检测-b3a72f1ea8c7
生产中数据科学的会合体系结构
如何构建一个前沿的数据科学平台来解决数据科学中的真正挑战:生产化。集合点架构介绍。
第 1 部分:数据科学的真正挑战
不可能错过数据领域如何获得一些新的流行词汇。它曾经是“大数据”，但最近，它被数据湖、机器学习、人工智能、区块链和数据科学独角兽所主导。随着数据科学的大肆宣传，围绕它的生态系统不断扩大，越来越多的大数据平台承诺让数据科学变得简单，并为大众带来商业成功。
但是，数据科学有什么不同之处，以至于它需要新的方法和所有这些新的平台和产品？
2016 年的一项调查得出结论: 80%的数据科学都在准备和清理数据(80/20 法则)。数据科学家的这项调查很快流行起来，并发展成为数据科学领域公认的问题陈述。不幸的是。因为如果你从这篇博文中拿走一样东西
“模型的生产是数据科学中最困难的问题” (Schutt，R & O’Neill C.《从第一线直接做数据科学》, O’Reilly 出版社。加利福尼亚州，2014)
而不是准备数据。
为什么数据科学不同于任何软件项目。
Source: Kenneth Jensen, (wikimedia)
为了理解为什么数据科学是不同的，我仍然从臭名昭著的 80/20 规则和数据挖掘的跨行业标准过程(CRISP-DM)开始。尽管 CRISP-DM 早于今天的数据科学学科，但它完美地描述了数据科学家日常工作的很大一部分。自 1996 年以来，似乎发生了变化的方面是大量集成不良的数据源，以及对数据准备步骤的日益重视。因此，80/20 法则变得如此流行就不足为奇了。
但是在数据科学中，CRISP-DM 只是图片的一半！模型不是普通的软件，也不仅仅是部署数据科学(在此插入 meme！).模型具有独特的生命周期管理，以及独特的集成、自动化和操作要求。
生产中的模型是动态的，具有复杂的生命周期。不幸的是，生产中的模型管理很少被讨论或解决。重点仍然是数据准备和 80/20 法则的重要性。
然而，管理单个模型也不是全部，数据科学的复杂性仍在继续。在某个时候，将会出现一个挑战者模型:
解决这个问题的一种方法是使用传统的 A/B 测试场景和流量分流。但是你只能孤立地看到一个单一的模型响应。假设您正在尝试测试一个严重依赖于被评分用户的上下文的推荐模型。对所有模型的请求进行并行评分，然后简单地比较它们的评分，这样不是更好吗？此外，你的挑战者模型可能不会工作得很好。不幸的是，你必须让 10%的用户经历一次糟糕的体验才能找到答案。这显然需要一个更好的方法。
但是我们还没有完成。在成功的数据科学研究中，开发和生产的界限越来越模糊。一个单一的业务目标将会看到许多不同的现任和挑战者模型，以及在一个不断变化的环境中对这些模型进行再培训的要求，其中包括持续的基准测试、监控和生产移交。真实世界的场景将很快变得更像这样:
这些都是让任何 IT 部门都物有所值的要求。挑战在于解决数据科学工具包日益增长的异构性、数据基础架构所需的弹性和灵活性，以及复杂数据科学模型生命周期的端到端自动化。
第 2 部分:解决方案
认识集合建筑
总结前面对问题陈述的介绍，我们正在寻找一些架构来
- 并行评估大量现有车型和挑战者车型
- 管理模型生命周期
- 处理日益多样化的数据科学工具包
- 允许在不影响用户体验的情况下进行生产实验，并将业务目标与数据科学目标分离
- 将 SLA 和 GDPR 等企业需求从数据科学模型中分离出来
- 扩展到每分钟 50+k 页面负载的峰值，而无需雇佣大批 DevOps 工程师
这个非常真实的数据科学问题的解决方案来自于 Ted Dunning 和 Ellen Friedman 在他们的书《机器学习物流》(O’Reilly 2017)( 免费电子书)中命名为 Rendezvous Architecture 。该架构的核心是作为流的输入数据和使用发布-订阅消息模式订阅数据的模型。通过一个流接收评分请求允许我们轻松地分发一个请求并并行评估大量的模型。
但是它没有解决如何从大量模型中健壮地返回响应，同时支持企业需求。
为了返回一个响应，架构还将模型分数视为一个流，并与名称 lendingRendezvous service相结合。从概念上讲，rendezvous 服务订阅原始的评分请求，将它们与来自模型响应流的模型评分相结合。
在 rendezvous 服务中订阅评分请求将从成功返回评分输出的模型中返回响应的能力分离开来。它还允许 rendezvous 服务实现 SLA:在模型失败或被延迟并且没有评分输出到达 SLA 内的 rendezvous 服务的情况下，rendezvous 服务可以例如发送有意义的默认响应。这是可能的，因为会合服务首先知道评分请求。rendezvous 服务实现了模型评分输出与可配置 SLA 和超时策略的竞争条件。
关于 rendezvous 架构更详细的介绍，你应该读一下 Ted Dunning 和 Ellen Friedman 的伟大著作:https://mapr . com/ebooks/machine-learning-logistics/ch03 . html
采用这种概念架构并将其应用于我们提出的要求，我们得到了以下数据科学平台:
我们根据需要将架构划分为一个业务目标层次结构，该层次结构实现了企业需求和数据科学团队的建模目标。
整个平台是用容器化的无状态服务构建的，这提供了期望的解耦和可伸缩性。甚至有状态模型也可以在模型数据丰富器的帮助下在平台上实现，如下面步骤 3 中所讨论的。
- 当评分请求到达时，入口控制器将该请求分配给集合服务实例的评分评估服务，该集合服务实例实例化异步评估策略。这为模型评分输出设置了竞争条件。每个策略实例就像一个邮箱，收集它们对应的评分请求的模型评分输出。
- 评分评估服务将策略实例 ID 和集合服务 ID 附加到评分请求，并将其发布到该业务目标的输入数据流。
- 模型数据丰富器订阅带有评分请求的输入数据流。该服务附加了额外的数据点，如元数据或有状态模型的外部状态。例如，推荐模型可能会考虑用户以前的行为来创建推荐。因此，数据丰富器可以交叉引用用户 ID 和以前用户交互的缓存。模型数据富集器然后将富集的评分请求发布到建模目标信息流。
- 每个模型的消息消费者订阅建模目标输入流，并将评分请求发送到相应的模型实例。每个模型目标可以拥有一系列不同的模型，只要它们共享输入数据模式。每个模型目标还拥有一个诱饵模型和一个金丝雀模型，前者将评分请求记录到日志流中，后者用于检测输入数据中的漂移，并提供一致的基线来比较模型性能。
- 消息生产者基于集合服务 ID 将模型评分输出发送到相应的集合流，该集合服务 ID 已经由评分评估服务附加到消息。
- 分数收集服务从其相应的集合流中读取所有进入的模型评分输出。每条消息都有一个策略实例 ID，该 ID 已由评分评估服务附加到消息上。消息策略实例 ID 是分数收集服务和分数评估服务之间的共享卷上的 Unix 套接字的名称，其允许收集服务将模型评分输出发送到相应评分请求的等待策略实例。
- 等待策略实例通过 unix 套接字服务器接收模型分数，并评估其策略逻辑是否以特定分数响应发起的评分请求。
实现集合服务
该平台的大部分相当简单，因为它由无状态的容器化服务和简单的发布/订阅消息模式组成。在生产中，这可以使用舵图和水平 pod 自动缩放器以及循环负载平衡在 Kubernetes 上轻松部署。实际上，业务和建模目标中的平台图中的每个分组实体都是水平缩放的豆荚。
主要的挑战是实现有状态的集合服务。从概念上讲，最好将 rendezvous 服务描述为邮箱。每个活动评分请求都有一个等待模型评分响应的活动邮箱。这些邮箱接收其相应评分请求的模型评分，并执行附加的评分评估策略。
当模型分数被发布到流时，每个邮箱可以订阅该流，并过滤出它们相关的模型分数进行评估(扇出模式)。这将为每个订阅的邮箱复制邮件。在每分钟 50+k 个请求的规模下，大量现有和挑战者模型并行返回分数，这将是不切实际的。为了实现可伸缩性，我们必须保持每个运行的策略实例尽可能的轻量级。
我们用 Python 实现了分数评估服务，使用 asyncio 库实现了异步事件驱动的策略评估。以下代码是该服务的简化版本:
import asyncio app.route(‘/topic/test’, methods=[‘POST’])def incoming_scoring_request(): scoring_request = SomeMessageDeserialisation( request.get_json(), os.environ['THIS_RENDEZVOUS_SERVICE_ID'] ) loop = asyncio.new_event_loop() asyncio.set_event_loop(loop) scoring_response = loop.run_until_complete( async_score_evaluation_service, scoring_request, loop ) return jsonify(scoring_response)async def async_score_evaluation_service(scoring_request, loop): policy = SomePolicy(loop.create_future()) policy_coro = asyncio.start_unix_server( policy.handle_connection, policy.socket_path, loop=loop ) policy_task = asyncio.ensure_future(policy_coro) extended_scoring_request = SomeMessageSerialisation( scoring_request, policy.socket_path ) message_task = syncio.ensure_future( publish(extended_scoring_request) ) try: scoring_response = await asyncio.wait_for(
policy.result_future, timeout=policy.sla
) except asyncio.TimeoutError: scoring_response = policy.default_value() return scoring_response
如代码片段所示，模型分数评估策略启动一个异步 unix 套接字服务器。套接字的路径由策略实例的 uuid 给出。这是获得模型分数所需的全部内容。除了评估服务，sidecar 分数收集服务订阅模型响应流。
import pulsarimport socket pulsar_consumer = get_message_consumer(os.environ['THIS_RENDEZVOUS_SERVICE_ID']) while True: msg = pulsar_consumer.receive() scoring_response = SomeMessageDeserialisation( msg.data() ) policy_socket = scoring_response.policy_socket_path sock = socket.socket(socket.AF_UNIX, socket.SOCK_STREAM) sock.connect(policy_socket) sock.sendall(SomeMessageSerialisation(scoring_response)) sock.close() consumer.acknowledge(msg)
它根据策略实例 ID 将消息发送到相应分数评估策略实例的 unix 套接字。在将扩展的评分请求发布到评分评估服务中该业务目标的数据流之前，我们将该 ID 附加到扩展的评分请求中。
为了在 Score Collection sidecar 和 Score Evaluation Service 之间创建连接，我在 unix 套接字的两个容器之间使用了一个共享卷。这种实现模式比简单的扇出模式伸缩性好得多，单个 Rendezvous 服务实例可以托管大量并发的分数评估策略。
我使用Apache Pulsar(link)作为分布式发布-订阅消息系统。Pulsar 不仅仅是一个高性能的消息传递系统，它可以轻松扩展到大量的主题，并以高度可扩展的消息持久存储为主要目标。主题压缩和分层持久存储(链接)允许 rendezvous 平台在未来无限期地保存消息，并重放评分请求，例如测试和验证新的挑战者模型或预热保持内部状态的模型。
下图显示了 minikube 开发环境中的一个示例:
该请求贯穿会合平台的所有元素，总时间为 51 毫秒。测试主题通过模型数据丰富器获取一个随机数，模型目标 test-0 运行一个简单的模型，该模型返回评分请求中所有数值的总和。我用那个装置来测试会合平台的元件。实际模型是集装箱化和解耦的，这允许我们独立地测试它们。**这意味着 rendezvous 架构的开销在 40-50 毫秒之间。**相对于该架构提供的广泛优势，这是一个非常小的成本，并且留下了充足的时间来获取外部数据点和评估复杂的模型。
做好这件事的代价是 40 毫秒
minikube 开发环境的负载测试表明，每秒 100 个请求也不是什么难事:
参考实现展示了 rendezvous 架构的优势，并展示了非常有前途的性能。我目前正在研究 Kubernetes 的生产，关注舵手和操作员。一旦我们在生产 Kubernetes 集群上部署了我们的平台，我将在这篇博客文章之后跟进我们的生产负载测试结果。
信用
我们正站在巨人的肩膀上。信用到期时的信用:
- Ted Dunning 和 Ellen Friedman 以及他们的励志书籍《机器学习物流》(O’Reilly 2017)介绍了概念上的会合架构
- 来自www.advancinganalytics.co.uk的 Terry McCann 和 Christopher Conroy 他们是我努力的早期信徒、贡献者和支持者
来自数据科学节的录音:
www.datasciencefestival.com Day 2 at Zoopla
Jan 是公司数据转型方面的成功思想领袖和顾问，拥有将数据科学大规模应用于商业生产的记录。他最近被 dataIQ 评为英国 100 位最具影响力的数据和分析从业者之一。
在 LinkedIn 上连接:https://www.linkedin.com/in/janteichmann/
阅读其他文章:【https://medium.com/@jan.teichmann】
歌曲中的重复:Python 教程
一首艾德·希兰歌曲的个案研究
Credit: Unsplash
每个人都听过一首歌或者知道一首歌听起来像什么。我可以漫不经心地说，每个人都可以用自己的话来定义一首歌。为了方便起见，一首歌曲*(根据维基百科)*是一首音乐作品，通常由人声演唱，使用声音和沉默以及各种形式(通常包括部分的重复)来表现独特和固定的音高和模式。
在他名为“**歌曲的复杂性”的期刊文章中，**计算机科学家 Donald Knuth 利用了流行歌曲从冗长且内容丰富的歌谣演变为高度重复的文本的趋势。由于有些人可能会立刻同意他的观点，这确实提出了一些问题，如:重复真的有助于歌曲成为热门吗？随着时间的推移，音乐真的越来越重复吗？
为了以案例研究的形式教授一些基本的 python 代码，我将测试这个假设*(流行歌曲真的是重复的吗？)*用我最喜欢的一首歌。检验这一假设的一种方法是找出独特的词，并计算这些词占歌曲总词数的比例。
在本教程中，我们将介绍:
- 变量和数据类型
- 列表和词典
- 基本算术运算
- 内置函数和循环
必备知识
为了充分利用本教程，您可以自己运行代码。
- 我们将要演奏的音乐是艾德·希兰的《完美》。这里可以复制歌词。然而，我在这个分析中使用的歌词被清理掉了，以得到一个结论性的结果。例如，我把像*“我们会”这样的词改成了“我们会”等等。*你可以在这里得到我版本的歌词
- 使用的编辑器是木星笔记本。这里是如何安装和使用它的快速教程。
出于本案例研究的目的，我们将通过提出两个主要问题来简化我们的假设:
- 与我们的案例研究歌曲《艾德·希兰的完美》的整个歌词相比，使用了多少独特的词？
- 在整首歌中重复使用最多的词是什么，使用了多少次？
让我们开始分析吧
基本的
1。 **A String** 是字符列表。字符是你可以在键盘上一次击键输入的任何东西，比如一个字母、一个数字或一个反斜杠。但是，Python 将字符串识别为在字符或文本的开头和结尾用双引号(" ")或单引号(’ ')分隔的任何内容。例如:“你好，世界”
对于这个案例研究，一个字符串是我们的歌词，如下所示
2。 **Variables** 通常是用于赋值或存储值的描述性名称、单词或符号。换句话说，它们是任何数据类型的存储占位符。为了在任何时候引用一个值，这是非常方便的。变量总是被赋予一个等号，后面跟着
变量的值。(查看代码输出的一种方法是使用打印功能。您可能已经知道，Jupyter 笔记本没有打印功能也可以查看输出)
为了存储歌词，我们将赋予它一个名为perfect_lyrics的变量。
3。 **Lists** 只需将不同的逗号分隔值放在方括号中，即可创建。它可以有任意数量的项，并且它们可以是不同的类型(整数、浮点、字符串等。).它甚至可以将另一个列表作为一个项目。例如:
list1 = [1,'mouse', 3.5, [2,'is',5.0]]
#3.5 is a float
现在我们已经对列表有了一个大概的了解。让我们回到我们的数据。
因为我们的目标之一是计算出使用的独特单词的数量，这意味着我们需要做一些计算，即计算每个单词。为了实现这些，我们不仅要把我们的字符串放到一个列表中，还要用一个**.split()**方法来分隔每个单词。因此，我们的数据集将看起来像这样
输入
输出
从上面的输出中，您会注意到每个单词都被分成了独立的字符串。并且整个歌词组成了名为 split_lyrics *(也是一个变量)*的列表
我们还需要将其余单词和计数中的独特单词分开。为此，我们需要使用 python 函数。
功能是一组有组织的、可重用的代码，用于执行单个相关的动作。Python 有很多内置函数，比如 print()-显示输出，list()-创建列表，len( )-计算列表或字符串中的字符数，等等。Python 还允许你创建自己的函数。但是在这个案例研究中，我们不会创建自己的函数。
我们将在案例研究中使用几个函数，但我们将从 set()函数开始。为了从整个歌词中分离出独特的单词，我们需要一个 set()和 print()函数，即
输入
输出
为了统计整首歌词的字数，我们需要一个 len()函数
输入&输出
对提取的唯一单词进行同样的操作
我们的上述分析帮助回答了我们的第一个问题:*How many unique words were used as compared to the whole lyrics of the song?*简单地说，在总共超过 290 个单词中使用了 129 个独特的单词。
我们的下一个目标是解决问题的第二部分What are the most repetitive words used and how many times were they used?
为了回答这个问题，我们需要学习更多的数据结构。
中间的
4。 **Dictionaries** 在 Python 中是无序的物品集合。数据结构只有值作为元素，而字典有一个**key:value**对。每个键-值对都将键映射到其关联的值。字典经过优化，可以在已知键的情况下检索值。您可以通过用大括号({})括起逗号分隔的键值对列表来定义字典。冒号(:)将每个键与其关联的值分开。一个简单的德-英词典就是一个例子:
colours = {"red" : "rot", "green" : "grün", "blue" : "blau", "yellow":"gelb"}print en_de["red"]#output will be the value of 'red'( the key) which is 'rot'empty_dict = {} #a dictionary can also be empty waiting to be filled up with information
5。 **Loops**在试图一遍又一遍地运行同一块代码时非常有用。在 python 中，有两种类型的循环:for 和 while。对于这个分析，我们将更多地关注for loop。for 循环用于迭代序列中的元素。当你有一段代码想要重复 n 次的时候，经常会用到它。
它是这样工作的:“对于列表或字典中的所有元素，这样做”
例如
list_number = [2,3,4]#input
for item in list_number:
multiple = 2 * item
print(multiple)
回到我们的数据集
要知道重复次数最多的单词，我们需要知道每个单词在歌词中出现的次数。为此，我们需要使用字典*(存储每个单词及其对应的计数)和 for 循环(在每个单词出现时迭代计数过程)*
首先，我们需要将唯一的单词存储在字典中
输入&输出
然后我们需要再次使用 for 循环来统计每一个唯一单词在整个歌词中出现的次数。
输入&输出
现在我们已经知道了每个单词出现的次数，让我们使用 sorted ( ) 函数将它们从最高到最低排序
只出现过一次的词好像太多了。因为我们的目标是找到最流行的单词，所以我们将通过切片将列表缩小到 10 个单词。你可以在这里了解更多关于切片的
输入&输出
然后变回字典
另一个问题是:艾德·希兰的歌曲《完美》中最流行的 10 个词是什么？通过使用字典数据结构下的 key()方法，然后创建这些单词的列表，我们可以很容易地提取这些信息。
输入&输出
从上面的输出代码中，我们可以肯定地说，在艾德·希兰创作的歌曲《完美的 T9》中，出现了 24 次的单词“我”
让我们进一步改进我们的分析
先进的
6。 **Visualizing data with python.**已经开发了各种技术来可视化地呈现数据，但是在这个分析中，我们将把数据可视化严格地集中在 Python 中的一个库上，即 Matplotlib。
让我们快速直观地分析一下案例研究中最流行的 10 个词
见解摘要
我们可以从这个案例研究中获得哪些关键的见解？我可以找出三个
- 根据选择的歌曲，我们看到独特的词只占总词数的 44%
- 在使用的 129 个独特的单词中，大约有 70 个单词出现过一次，约占总使用单词的 24%。
- “我”这个词在整首歌中被使用了 24 次，也就是说，有 8 %的时间，一个“我”被使用了。
结论
这个案例研究确实有助于发现歌曲中有很大一部分是由重复的词组成的。然而，我们不能仅凭一首歌就断定唐纳德·克努特理论是正确的。我们需要分析更多的歌曲，才能得出结论说热门歌曲是一首歌中单词重复的结果。尽管指出艾德·希兰的实际上是一首热门歌曲是件好事……(至今仍在《我的世界》中)。
延伸阅读
如果您想提高您的 Python 技能，这些文章可能会有所帮助
附注:像我一样，任何人都可以学习成为一名数据分析师，如果你想了解我的下一个项目或我的学习更新，请随时注册我的 简讯
使用 Python-2 构建高级会计模型
使用正确的数据科学工具充分利用会计公式
Photo by Carlos Muza on Unsplash
背景
如果你看了我之前的文章这里，你就知道我们可以很容易地从雅虎财经刮出财务报表。但是现在呢？当您进行任何类型的业务评估和财务分析时，我们的数据框架中的大多数特征都是重要的指标。它们都是公开的和经过审计的信息，最终决定了股票价格的变动。管理层希望他们看起来不错，所以他们经常操纵他们，并确保数字符合公众的期望。
现在，这就是数据科学和您的领域知识发挥作用的地方。为了解释会计比率，你需要理解这些特征的含义。这往往需要上几堂大学会计课。你可能已经学过一些比率，如每股收益、利润率、速动比率等。今天，我将重点讨论**股本回报率、**与股本相关的企业盈利能力。
股本回报率=净收入/股本
数据清理
在进行任何数值操作之前，我们需要整理数据。以下是我们的目标:
- 将数据类型从 object 转换为 numeric。
- 将日期索引转换为 DateTime 以进行时序分析。
我不会在清理数据上做太多的细节，因为这是另一个话题，但是可以随意检查我写的这个函数。当你真正检查数据帧时，它们更有意义。你可以暂时忽略它。
def convert_to_int(df):
col = df.columns
for col in df:
temp = df[col].to_string().replace(‘,’,’’).split(‘\n’) # Replace comma with space and clean out extra spaces temp.pop(0) #Remove first erronous element df.index = pd.to_datetime(df.index) # Convert index to datetime df[col]= [i[10:].strip() for i in temp] # stripping off nth characters df[col] = pd.to_numeric(df[col], errors=’coerce’).fillna(0).astype(int) # Transform columns from object to numeric return df
简介
简单的杜邦模型将净资产收益率分解为 3 个部分
净资产收益率=(净收入/销售额)(销售额/总资产)(总资产/股东权益)
第一个组成部分侧重于经营业绩:每笔销售产生多少净收入。第二个衡量公司利用其资产创造销售收入的效率。第三个关注股权管理:你每股产生多少资产。
然而，如果你仔细观察引擎盖下发生的事情，你会注意到那些组件非常嘈杂。例如，净收入由财务费用组成，这意味着运营 KPI 与投资 KPI 相结合。第二，什么样的资产最能反映整体经营状况？我们想看看运营方面可以产生销售的资产。
现在你可能会问，我们为什么要关注运营？我们为什么不考虑投资和财务方面呢？我们确实有来自投资和贷款的现金流入。嗯，要看行业和商业模式。一般来说，经营现金流是企业的血液。其他一切都只是支撑。
高级杜邦模型
净资产收益率=净资产收益率+财务杠杆*利差
上述公式的想法是，我们希望将运营和财务部分分开。 **1。**我们想知道每股股票的回报率是多少。好了，我们现在来看经营中资产的回报，基本上就是净收入加上利息之类的财务费用(因为是从投资方来的，这是一个不稳定不可靠的指标)和少数股权(实际上不属于你的普通股权益份额的东西)。
**2。**好了，现在是财务杠杆部分。我买了你的股份，但是你用了多少来偿还债务呢？这样，在扣除成本后，有多少金融债务实际上产生了利润。
不要被不熟悉的公式吓倒，下面我来做个分解。一旦你理解了逻辑流程，一切都会变得清晰。
净营运资产回报率=净营运利润率*净营运资产周转率
净营业利润率=(净收入+净财务费用)/净销售额
净营业资产周转率=净销售额/平均资产
财务杠杆=平均净财务责任/平均普通股权益
净借款成本=净财务费用/平均净财务义务
净财务责任=总债务+少数股东权益+优先股
净财务费用=净利息费用* (1-税率)+优先股股息+盈利中的少数股权
利差=净经营资产回报率-净借款成本
数据操作
在这里，我创建了一个新的数据框架来包含所有的比率计算，以得出最终的股本回报率。为了更清楚，我用粗体突出了重要的部分。
ratio = pd.DataFrame()
接下来，我发现有效税率是多少。
ratio[‘effective tax rate’]= Income_st[‘Income Tax Expense’]/Income_st[‘Income Before Tax’]
然后我从运营部分减去财务部分。记住，我们想从我们的税收中拿回那部分。
**ratio[‘net operating margin’]= (Income_st[‘Income Before Tax’]+Income_st[‘Interest Expense’]*(1-ratio[‘effective tax rate’]))/Income_st[‘Total Revenue’]**
在这里，我使用今年和去年之间的平均时间点来更好地表示总负债和权益的潜在价值。
ratio[‘net operating asset turnover’] = Income_st[‘Total Revenue’]/((balance_st[‘Total liabilities and stockholders\’ equity’]+balance_st[‘Total liabilities and stockholders\’ equity’].shift(1))/2)**ratio[‘Return on NOA’]= ratio[‘net operating margin’]*ratio[‘net operating asset turnover’]**ratio[‘Net financial expense’] = Income_st[‘Interest Expense’]*(1-ratio[‘effective tax rate’])ratio[‘Net borrowing cost’]= ratio[‘Net financial expense’]/(balance_st[‘Total Liabilities’]+balance_st[‘Total Liabilities’].shift(1)/2)**ratio[‘spread’]= ratio[‘Return on NOA’]-ratio[‘Net borrowing cost’]****ratio[‘leverage’] =(balance_st[‘Total Liabilities’]+balance_st[‘Total Liabilities’].shift(-1)/2) / (balance_st[‘Total liabilities and stockholders\’ equity’]+balance_st[‘Total liabilities and stockholders\’ equity’].shift(1))/2**
这是最终的净资产收益率。
ratio[‘ROE(Advanced dupont)’] = ratio[‘Return on NOA’]+ ratio[‘leverage’]*ratio[‘spread’]
数据可视化
这是我们新的数据框架和一些简单的数据可视化。
你注意到有很多 NaN 值，这只是因为我在计算中使用了平均值，而我们没有 2015 年的数据。如果您有更多的数据，并按季度或月份进行分析，您会看到一个更复杂的数据框架。
用 Python 替换 Excel
在与我的初恋情人 Excel 相处了近十年后，是时候继续前进，寻找一个更好的另一半了，他在我的日常任务中与我同甘共苦，比我好得多，也比我快得多，并且可以在充满挑战的技术时代给我带来优势，在这个时代，新技术正以非常快的速度被新事物抛弃。这个想法是在 Python 中复制几乎所有的 excel 功能，无论是使用简单的过滤器，还是从行中创建数据数组并处理它们以获得有趣结果的复杂任务
这里遵循的方法是从简单的任务开始，然后转移到复杂的计算任务。我鼓励你自己重复这些步骤，以便更好地理解。
创造这种东西的灵感来自于一个免费教程的不可用性。我大量阅读并关注 Python 文档，您会从该网站中找到很多灵感。
GitHub 回购链接
https://github.com/ank0409/Ditching-Excel-for-Python
将 Excel 文件导入熊猫数据框架
第一步是将 excel 文件导入 DataFrame，这样我们就可以在它上面执行所有的任务。我将演示熊猫的 read_excel 方法，它支持 xls 和 xlsx 文件扩展名。
read_csv 与使用 read_excel 相同，我们不会深入讨论，但我会分享一个例子。
虽然 read_excel 方法包含了数百万个参数，但是我会让你熟悉那些在日常操作中非常有用的最常见的参数。
虽然 read_excel 方法包括数百万个参数，但我会让您熟悉在日常操作中会非常方便的最常见的参数
我将使用 Iris 样本数据集，该数据集在网上免费提供，用于教学目的。
请点击以下链接下载数据集，并确保将其保存在保存 python 文件的同一文件夹中
https://archive.ics.uci.edu/ml/datasets/iris
第一步是用 Python 导入必要的库
我们可以使用以下代码将电子表格数据导入 Python:
pandas.read_excel(io，sheet_name=0，header=0，names=None，index_col=None，parse_cols=None，usecols=None，squeeze=False，dtype=None，engine=None，converters=None，true_values=None，false_values=None，skiprows=None，na_values=None，keep_default_na=True，verbose=False，parse_dates=False，date_parser=None，千位=None，comment=None，skip_footer=0，skip _ footer = 0
因为有太多的参数可用，让我们看看最常用的一个。
重要熊猫 read_excel 选项
如果我们使用本地文件的路径，默认情况下它由“\”分隔，但是 python 接受“/”，所以要更改斜线或简单地将文件添加到 python 文件所在的同一文件夹中。如果你需要以上的详细解释，请参考下面的文章。https://medium . com/@ ageitgey/python-3-quick-tip-the-easy-way-to-deal-of-file-paths-on-windows-MAC-and-Linux-11a 072 b 58 d5f
我们可以使用 Python 扫描目录中的文件，并挑选出我们想要的文件。
导入特定的工作表
默认情况下，文件中的第一个工作表按原样导入到数据框中。
使用 sheet_name 参数，我们可以明确地提到我们想要导入的工作表。默认值为 0，即文件中的第一张图纸。
我们既可以提到工作表的名称，也可以传递一个整数值来引用工作表的索引
使用工作表中的列作为索引
除非明确说明，否则在数据帧中添加一个索引列，默认情况下从 0 开始。
使用 index_col 参数，我们可以操作数据帧中的索引列，如果我们将值从 none 设置为 0，它将使用第一列作为我们的索引。
跳过行和列
默认的 read_excel 参数假设第一行是列名的列表，它作为列标签自动包含在数据帧中。
使用 skiprows 和 header 等参数，我们可以操作导入数据帧的行为。
导入特定的列
使用 usecols 参数，我们可以指定是否在数据帧中导入特定的列
这并不是可用功能的终结，而是一个开始，您可以根据自己的需求来使用它们
让我们来看看从 10，000 英尺高度的数据
现在我们有了数据框架，让我们从多个角度来看这些数据，以便掌握其中的窍门/
Pandas 有很多我们可以使用的功能。我们将使用其中的一些来浏览一下我们的数据集。
“头”对“尾”:
查看第一个或最后一个的五个行。
默认为 5，然而该参数允许我们使用一个特定的数字
查看特定列的数据
获取所有列的名称
信息方法
给出了数据帧的概要
形状方法
返回数据帧的维度
查看 DataFrame 中的数据类型
切片和切块，即 Excel 过滤器
描述性报告是关于数据子集和聚合的，当我们稍微理解我们的数据时，我们开始使用过滤器来查看更小的数据集或查看特定的列，可能会有更好的理解。Python 提供了许多不同的方法来分割数据帧，我们将使用其中的几种方法来了解它是如何工作的
查看特定列
选择列有三种主要方法:
- 使用点符号:例如 data.column_name
- 使用方括号和列名，例如 data[‘column_name’]
- 使用数字索引和 iloc 选择器 data.loc[:，’ column_number’]
查看多列
查看特定行的数据
这里使用的方法是使用 loc 函数进行切片，在这里我们可以指定由冒号分隔的开始和结束行
记住，索引从 0 开始，而不是从 1 开始
将行和列切片在一起
筛选列中的数据
过滤多个值
使用列表筛选多个值
筛选不在列表中或不等于 Excel 中的值
在多列中使用多个条件进行筛选
输入应该总是列表
我们可以用这种方法复制 excel 中的高级过滤功能
使用数字条件过滤
在 Excel 中复制自定义筛选器
组合两个过滤器以获得结果
包含 Excel 中的函数
从 DataFrame 中获取唯一值
如果我们想要查看具有唯一值的整个数据帧，我们可以使用 drop_duplicates 方法
排序值
按某一列对数据进行排序，默认情况下是升序排序
数据的统计摘要
数据帧描述方法
生成描述性统计数据，总结数据集分布的集中趋势、离散度和形状，不包括 NaN 值
字符列的汇总统计信息
数据聚合
计算特定列的唯一值
产生的输出是一个序列。您可以将其称为单列透视表
细胞计数
对每列或每行的非 NA 细胞进行计数
总和
汇总数据以获得行或列的快照
复制针对每一行添加总计列的方法
向现有数据集中添加汇总列
特定列的总和，使用 loc 方法并传递列名
或者，我们可以使用下面的方法
不喜欢新列，使用 drop 方法将其删除
在每列下添加总和
上面已经做了很多，我们使用的方法是:
- Sum_Total:对列求和
- T_Sum:将串行输出转换为数据帧并转置
- 重新索引以添加缺少的列
- Row_Total:将 T_Sum 附加到现有数据帧
基于标准的总和，即 Excel 中的 Sumif
苏米夫斯
平均 if
最大
福建话
Groupby，即 Excel 中的小计
数据框架中的数据透视表，即 Excel 中的数据透视表
谁不喜欢 Excel 中的数据透视表呢？它是分析数据的最佳方式之一，可以快速浏览信息，用一个超级简单的界面帮助你分割数据，帮助你根据数据绘制图表，添加计算列等。不，我们不会有一个接口来工作，我们将不得不显式地编写代码来获得输出，不，它不会为您生成图表，但我不认为我们可以在不了解数据透视表的情况下完成一个教程。
一个简单的数据透视表向我们显示了行中的 SepalWidth 值、列中的 SepalLength 值和列标签中的 Name 值的总和
让我们看看能不能把它弄复杂一点。
通过使用 fill_value 参数，空格现在被替换为 0
我们可以使用字典方法对值进行单独计算，也可以对值进行多次计算
如果我们使用 margins 参数，我们可以添加总计行
纵向查找函数
Excel 中的 vlookup 是一个多么神奇的公式啊，我想这是每个人在学习加法之前都想学的第一件事。当有人应用 vlookup 时看起来很迷人，当我们得到输出时看起来很神奇。让生活变得简单。我可以非常自信地说，它是在电子表格上进行的所有数据争论的支柱。
很遗憾我们熊猫没有 vlookup 功能！
因为我们在 Pandas 中没有“Vlookup”功能，所以使用 Merge 作为替代，这与 SQL 相同。总共有四个合并选项可用:
- ‘left’—使用左侧数据帧中的共享列并匹配右侧数据帧。将任何 N/A 填写为 NaN
- ‘右’—使用右侧数据帧中的共享列，并与左侧数据帧匹配。将任何 N/A 填写为 NaN
- ’ inner’ —仅显示两个共享列重叠的数据。默认方法。
- ‘outer’—当左数据帧或右数据帧中存在匹配时，返回所有记录。
上面的例子可能不是支持这个概念的最好例子，但是工作原理是一样的。
正如他们所说，“没有完美的教程存在”，我的也没有:)
用 Python 从各种工作表中提取数据
或者如何学习统一 Google 工作表、Excel 和 CSV 文件——代码指南
Widespread tabular data storage file formats — CSV, Microsoft Excel, Google Sheets
Python 通常被称为粘合语言。这是因为随着时间的推移，开发了大量的接口库和特性——这是由它的广泛使用和令人惊叹的、广泛的开源社区推动的。这些库和特性提供了对不同文件格式以及数据源(数据库、网页和 API)的直接访问。
这个故事着重于数据的提取部分。下周的故事将更深入地分析综合数据，以获得有意义和令人兴奋的见解。
了解如何使用 Seaborn 快速生成漂亮而有见地的图表
towardsdatascience.com](/how-to-explore-and-visualize-a-dataset-with-python-7da5024900ef)
但是不要让这阻止你自己分析数据。
您将学到的内容:
- 从 Google Sheets 中提取数据
- 从 CSV 文件中提取数据
- 从 Excel 文件中提取数据
这篇文章是写给谁的:
- Python 初学者
- 不得不经常争论数据的人
由于本文旨在作为一篇代码文章，您应该设置您的开发环境(我推荐 Jupyter Notebook/Lab)并开始一个新的笔记本。你可以在这里找到源代码和文件。
如果你不知道如何使用 Jupyter/Python。请查看该指南:
到底是什么阻止了你？下面是如何开始！
towardsdatascience.com](/get-started-with-python-e50dc8c96589)
情况:
在今天的故事中，我将带你进入一个虚构但可能非常熟悉的场景。你要结合各种来源的数据来创建一个报告或运行一些分析。
免责声明:下面的例子和使用的数据完全是虚构的
你的任务是找出如何提高你的销售团队的业绩。在我们假设的情况下，潜在客户有相当自发的需求。当这种情况发生时，您的销售团队会在系统中输入一个订单线索。然后，您的销售代表会尝试在订单线索被发现时安排一次会议。有时在之前，有时在之后。你的销售代表有一个费用预算，并且总是把会议和他们付钱的一顿饭结合在一起。销售代表报销他们的费用，并将发票交给会计团队处理。在潜在客户决定是否接受你的报价后，勤奋的销售代表会跟踪订单线索是否转化为销售。
对于您的分析，您可以访问以下三个数据源:
- 100，000 条订单线索(谷歌表单)
- 约 50.000 份餐费发票(Excel 文件)
- 公司和负责这些公司的销售代表的列表(CVS 文件)
获取 Google Sheets 数据:
访问 Google Sheets 是三者中最复杂的，因为它要求您设置一些使用 Google Sheets API 的凭证。理论上，您可以抓取一个公开可用的 Google Sheet(即，提取 HTML 源代码)，但是您必须使用 Beautiful Soup 之类的工具进行大量的数据操作，才能将 HTML 转储转换成有用的东西。我确实试过这个，但是结果很乱，不值得努力。那就 API 吧。此外，我们将使用 gspread 更加无缝地转换成熊猫数据帧。
获取 OAuth2 凭据
前往 谷歌开发者控制台 并创建一个新项目(或者选择一个已有的)。单击“创建项目”按钮。如果你的公司使用谷歌邮件，你可能想换成你的私人账户，以避免潜在的权限冲突。
为你的项目选择一个名字(名字不重要，我称之为矿媒数据提取)
点击 API&服务并前往库
**启用 Google Sheets API。**点击结果，在下一页点击启用 API。
**创建服务账户&密钥文件。**服务账户是用于具有有限访问权限的编程访问的专用账户。服务帐户可以而且应该由 project 设置，并具有尽可能具体的权限以及手头任务所必需的权限。
创建一个 JSON(另一种文件格式)密钥文件。对于角色，选择“项目->查看器”
如果在上一步中没有设置角色，现在就设置。
Note: Setting “Viewer” is somewhat restrictive. If you want to create google sheets programmatically, you’ll have to choose a different setting
您的私有 JSON 密钥文件将可以下载或者自动下载。我建议将该文件重命名为“medium _ Data _ Extraction _ key . JSON ”,并将该文件移动到您的 Jupyter 笔记本的文件夹中，因为这将使下面的示例无缝工作。JSON 文件包含您最近创建的服务帐户的凭证。
完美。你差不多完成了。
下载数据
首先，您必须通过在笔记本中运行以下命令来下载并安装附加的软件包。
!pip install gspread
!pip install oauth2client
其次，您必须确保将之前创建的 JSON 密钥文件移动到运行 Jupyter 笔记本的文件夹中，如果它还不在那里的话。或者，您可以指定一个不同的 GOOGLE_KEY_FILE 路径。
WORKBOOK_KEY 是我为本次会议准备的 Google 工作表的工作簿 id。
WORKBOOK_KEY = '10HX66PbcGDvx6QKM8DC9_zCGp1TD_CZhovGUbtu_M6Y'
该表是公开的。如果您想下载不同的数据，您需要更改 WORKBOOK_KEY。id 通常可以在问题 URL 的 Google 工作表的最后两个反斜杠之间找到。
获取 CSV 数据
我们既可以通过传统方式从 repo 下载 CSV 数据，也可以使用下面的代码片段。同样，您可能需要像这样安装缺失的请求包(在您的笔记本中运行):
!pip install requests
CSV 数据的美妙之处在于 Python / Pandas 可以开箱即用地处理它。对于 Excel，还需要额外的库。
获取 Excel 数据
在我们开始之前，你很可能必须安装 openpyxl 和 xlrd ，这使得你的熊猫也能打开 Excel 表格。
!pip install openpyxl
!pip install xlrd
完成之后，我们以同样的方式获得 Excel 数据，并将其加载到另一个数据框架中。
**搞定！**您已经创建了三个不同的 Pandas 数据框，并且可以在同一个 Jupyter 笔记本中访问它们:
- 销售 _ 数据
- 销售 _ 团队
- 发票
请继续关注下周的文章，届时我们将探讨如何组合这些数据框架并分析数据。
在 Excel 中用 Java 替换 VBA
几乎每个工作场所都有 Excel。从顶级投资公司和大规模工程公司到个体商人，人们都使用 Excel 完成工作。
本文将探讨使用 Excel 的一些问题和优势，以及如何使用嵌入 Excel 的Java来克服这些问题。
你不用找很远就能找到对 Excel 的批评和它的错误使用导致公司重大损失的案例。在过去的几年里，有许多案例表明，Excel 电子表格中的错误至少是令人尴尬且代价高昂的错误的部分原因。
一个简单的电子表格错误让一家公司损失了 2400 万美元。这个错误导致了加拿大大型电力公司 TransAlta
www.theregister.co.uk](https://www.theregister.co.uk/2003/06/19/excel_snafu_costs_firm_24m/) [## 伦敦鲸的崩溃部分是由于使用 Excel 时的一个错误
这是人们开始在博客圈谈论的事情，应该让整个华尔街暂停下来。结束…
www.businessinsider.com](https://www.businessinsider.com/excel-partly-to-blame-for-trading-loss-2013-2?r=US&IR=T)
有这样的风险，为什么公司还在用 Excel，有什么办法可以防止类似的情况发生？
Excel 被如此大量使用的主要原因是 工人生产率 。使用 Excel 的人可以比使用任何其他工具更有效地完成大量复杂的工作。
软件开发人员经常争辩说，用他们最喜欢的编程语言也能达到同样的效果。从技术上讲，他们可能是正确的，但这是假设每个人都有时间和兴趣学习成为一名开发人员(这需要我们大多数人很多很多年)。
大多数公司没有资源让开发人员专门为每个业务用户服务，即使他们有，那么沟通需要什么和迭代以获得期望的结果将很难与个人“拼凑”Excel 电子表格相竞争。
Excel 的问题
Excel 有什么地方容易出错？这本身没有什么，但是开发人员已经习惯了各种事情来减少非 Excel 解决方案中的风险。下面列出了 Excel 电子表格开发中的一些缺点。
- 太复杂。一个电子表格可能开始时相当简单，只有几个单元格和公式。然后它一点一点地成长。范围被复制以处理越来越多的情况或多组数据，直到很难理解发生了什么。这项任务可能仍然相当简单，但是由于需要复制，而且每个单元格只能容纳一个数据单元，电子表格可能会变得过于复杂。
- 最少或无测试。Excel 电子表格很少受到测试。通常没有为 VBA 代码编写的单元测试，唯一的功能测试通常由电子表格的作者在有限的输入集上进行。当另一个用户不得不使用同一个电子表格时，很有可能因为他们不熟悉其工作原理的细微差别，他们会偶然发现一些从未测试过的错误。
- 陈旧或过时的数据。将 Excel 连接到外部数据源可能很棘手。它可以直接连接到数据库，但是如果您需要的数据是由另一个系统产生的，或者您没有直接访问所需数据的数据库，该怎么办？数据经常从其他系统的报告中复制并粘贴到 Excel 中，通常无法知道数据上次是何时复制的，甚至无法知道数据是否完整。
- 虫子失控蔓延。即使您知道一个电子表格中有错误，并且您已经知道如何修复它，您如何找到复制并粘贴了相同 VBA 代码的其他电子表格的所有实例，或者甚至是完全相同的电子表格的副本？很可能你不能。电子表格被复制并通过电子邮件发送，电子表格、数据和代码之间没有分离。
- 版本控制地狱。当你正在处理一个大的电子表格，并且它已经稳定运行了，你会怎么做？最有可能的答案是你保存一份拷贝——也许在文件名上加上日期。这就是 Excel 中的版本控制。如果一个变化有意想不到的后果呢？你怎么知道这个变化是什么时候引入的？几乎不可能！
我们如何解决这些问题？
这些问题都源于采用表面上很简单的东西(相关数字的网格)，并将其推到极难推理其行为和正确性的程度。
本质上，Excel 是一种表达事物之间关系的方式。 A1 是 B1 和 C1 之和，例如*。当这些关系变得越来越复杂时，问题就开始出现了。*
如果你想计算“ A1 是时间序列的每日收益的方差 X ”，在 Excel 中会是什么样子？如果您是一位经验丰富的 Excel 用户，您可能会想象一个表示时间序列 B 的表，其中包含用于计算回报的额外列和用于计算方差的公式。但是，如果现在我们想计算另一个 N 时间序列的回报呢？复制并粘贴每个新时间序列的公式？这就是错误开始蔓延的原因！
更好的方法是将计算时间序列日收益方差的算法封装到一个函数中。然后我们可以重复调用这个函数，次数不限，不会有中间单元格被编辑或复制错误的风险。
现在，设想用 Excel 中的单个单元格来表示时间序列，而不是数据表。如果这可以实现，那么我们又回到了两个单元之间的简单关系—“A1 = daily _ returns _ variance _ of(B1)”。突然间，我们的电子表格开始看起来不那么复杂了！
我们仍然有时间序列数据必须来自某处的问题。与其从另一个系统或数据库复制粘贴，不如我们有一个函数直接从那个系统或数据库加载时间序列？这样，每次我们计算电子表格时，我们都知道数据是最新的和完整的！继续前面的例子，我们可能有“ B1 =负载时间系列(股票，开始日期，结束日期)”。稍后，我们将讨论如何在单个单元格中存储整个数据集。
通常不仅仅是使用 Excel 的人编写他们使用的函数。通过为最终用户提供一组可靠的 Excel 函数，技术团队可以比仅仅编写仅具有浅层 Excel 集成的应用程序(如导出报告)更有效地支持业务。
Java 对我们有什么帮助？
通过思考我们的电子表格，将算法放入函数中，并通过直接获取数据而不是复制和粘贴，我们已经解决了 Excel 电子表格的几个大问题。我们还没有触及如何编写这些函数，以及围绕测试、错误修复和版本控制的问题。
如果我们决定在 VBA 编写我们所有的函数(相信我，很多人都这么做！)那么我们就不能利用过去 20 年中软件开发的任何进步！
Java 与现代软件开发并驾齐驱，在 VBA 上有很大贡献。
- 测试。Java 有许多不同的测试框架，都有不同的优点和缺点。无论您选择哪一种，能够在您的代码库上运行自动化测试套件都会让您确信它在做正确的事情。这在 VBA 是不可能的。
- 广泛的库支持。编写 VBA 代码通常是编写网上找到的相当标准的算法，并把它们转换成 VBA。想做一些琐碎的事情，比如对一组数据进行排序？在 Java 中这不成问题，但在 VBA，你要负责确保你的排序算法有效，并且不需要任何测试。现在想象一下写一个复杂的衍生品定价模型！
- 将代码放在 Excel 之外。 VBA 代码通常保存在工作簿中，这就是为什么在共享工作簿时，错误变得如此难以追踪。如果您的电子表格引用了一个编译过的 Java 库(JAR ),那么它就在所有引用它的电子表格之外，可以很容易地更新。
- 版本控制。 Java 源代码只是文本，因此很容易被签入版本控制系统。大多数 Java IDEs 对此有很好的支持，因为这是现代软件开发的标准部分。
- **发展环境。**VBA 编辑器(VBE)已经很多年没变了。它只提供了一个非常基本的文本编辑器，具有基本的调试功能。另一方面，Java 有一系列优秀的 ide 可供选择。
但是 Java 不是 Excel 的一部分！
的确如此，但是 Excel 有一个“加载项”的概念，允许开发人员扩展 Excel 的功能。一个这样的插件是 Jinx，Excel Java 插件。
使用 Jinx ，你可以完全抛弃 VBA，完全用 Java 编写工作表函数、宏和菜单。
用 Java 编写工作表函数就像在 Java 方法中添加 Jinx 的@ExcelFunction 注释一样简单:
**package** **com.mycompany.xladdin**;
**import** **com.exceljava.jinx.ExcelFunction**;
*/***
** A simple Excel function.*
**/*
**public** **class** **ExcelFunctions** {
*/****
** Multiply two numbers and return the result.*
**/*
**@ExcelFunction**
**public** **static** double multiply(double x, double y) {
**return** x * y;
}
}
您可以返回您期望的所有基本类型，以及 1d 和 2d 数组。对于更复杂的类型，您可以编写自己的类型转换器，或者您可以将 Java 对象作为对象句柄直接返回给 Excel，以便传递给另一个 Java 方法。
jinx免费下载。参见 Jinx 用户指南了解更多关于如何用 Java 替代 VBA 的信息。
将时间序列作为单个单元格返回是怎么回事？
Jinx 函数可以返回所有你期望的标准类型(整型、双精度型、数组等。)，但是也可以返回 Java 对象！当一个复杂的 Java 对象(比如一个表示从数据库加载的时间序列的类)被返回时，它将作为一个对象句柄返回到 Excel。然后，可以将该对象句柄传递给其他 Jinx 函数，并将从第一个函数返回的原始 Java 对象传递给 Java 方法。
对于简化电子表格来说，这是一种非常有用的技术，可以避免电子表格中涉及到复杂的数据。需要时，可以使用 Jinx 数组函数将对象扩展为 Excel 数组。
您可以在用户指南的 Jinx 对象缓存部分了解更多关于这些对象句柄的信息。
其他语言
这些技术不是 Java 独有的。
其他 JVM 语言也是如此，比如 Scala T1 或 T2 kot Lin T3。Jinx 适用于所有 JVM 语言，不仅仅是 Java。
另一种编写 Excel 插件的流行语言是 Python。这可以通过使用 Excel 的 PyXLL 插件来实现。
复制多伦多图书语料库数据集——一篇综述
通过从原始来源收集和预处理图书，复制不再公开的多伦多图书语料库数据集
Front page of the Smashwords website
虽然 Toronto BookCorpus (TBC)数据集不再公开提供，但它仍然经常用于现代 NLP 研究(例如，像 BERT、RoBERTa、XLNet、XLM 等的变压器)。)因此，对于我们这些无法在数据集离线前获取其副本的人来说，这篇文章提供了一种尽可能复制原始 TBC 数据集的方法。
🔍了解原始的多伦多图书语料库数据集
因此，为了尽可能好地复制 TBC 数据集，我们首先需要查阅介绍它的 原始论文 和 网站 以更好地理解它的内容。
在论文中，朱等人(2015)写道:*“我们从网络上收集了 11038 本书的语料库。[……]我们只收录了超过 2 万字的书籍，以过滤掉可能更嘈杂的短篇故事。”*接下来，作者给出了一些汇总统计数据:
从网站上，我们了解到网站 Smashwords 是数据集中收集和使用的 11038 本书的原始来源。
📚收集书籍
现在我们已经(大致)知道了要收集多少本书，从什么来源收集，我们就可以开始收集这些书了。为此，我写了一些抓取 Smashwords 网站的代码，在 这个 GitHub 资源库 中公开。代码很快(并发)，结构良好，有据可查，所以应该非常容易使用。
Inspector mode on a Smashwords book page (accessible through “Inspect Element” or F12 on Firefox)
🔗获取纯文本图书 URL
为了获得要下载的明文书籍的 URL 列表，我们首先需要从 Smashwords 的首页获取书籍页面的 URL(每本书在 Smashwords 上都有自己的页面)。接下来，我们可以从这些书籍页面中抓取明文书籍的 URL 进行下载。你可以使用我的代码 在这里 找到这样做的说明。
📥下载明文书籍
现在我们有了一个要下载的明文书籍列表(使用它们的 URL)，我们需要…下载它们！这有点棘手，因为 Smashwords(暂时)阻止任何在一定时间内下载太多书(> 500 本)的 IP 地址。然而，这可以通过在下载书籍和经常切换 IP 地址(大约 30 次)时使用 VPN 来规避。你可以使用我的代码 在这里 找到这样做的说明。
⚙️预处理书籍
为了获得 Toronto BookCorpus 数据集在大小和内容方面的真实副本，我们需要对我们刚刚下载的明文书籍进行如下预处理:1 .句子标记的书籍和 2。将所有书籍写入一个文本文件，每行一句话。你可以使用我的代码 在这里 找到这样做的说明。
**就这样，大功告成！**🙌现在，您可以使用您刚刚创建的 Toronto BookCorpus 数据集的副本，亲自尝试 NLP 中最新最棒的内容。🤗
[1] Zhu，y .，Kiros，r .，Zemel，r .，Salakhutdinov，r .，Urtasun，r .，Torralba，a .，& Fidler，s .，对齐书籍和电影:通过观看电影和阅读书籍实现类似故事的视觉解释(2015)，《IEEE 计算机视觉国际会议论文集》(第 19–27 页)。
使用 Keras 和 TensorFlow 报告时间执行预测
这篇文章的目的是用实践的方式向软件开发者解释机器学习。该模型基于企业系统中的一个常见用例—预测生成业务报告之前的等待时间。
Source: Pixabay
用例
业务应用程序中的报告生成通常需要时间，可能是几秒钟到几分钟。生成报告需要时间，因为通常会获取和处理许多记录，这个过程需要时间。用户经常感到沮丧，他们不知道要等多久才能完成报告，可能会关闭浏览器等。如果我们能够在提交报告请求之前通知用户——执行它需要多长时间，这将是一个很大的可用性改进。
我已经使用 Keras 回归实现了机器学习模型，根据训练数据(来自过去报告执行的日志信息)计算预期的报告执行时间。Keras 是一个将 TensorFlow 复杂性包装成简单且用户友好的 API 的库。
这篇文章的目标是为软件开发人员提供一步步的指导，告诉他们如何构建一个简单而有用的机器学习模型。
Python 源代码和训练数据可以在我的 GitHub repo 上获得。这段代码基于 Keras 教程。
数据
使用 Pandas 库从 CSV 文件中获取训练数据:
column_names = ['report_id','report_params','day_part','exec_time']
raw_dataset = pd.read_csv('report_exec_times.csv')
dataset = raw_dataset.copy()
数据有四个特征，其中 exec_time 是目标特征(我们将要学习预测的那个)。 Report_id (报告类型标识符)和 day_part (一天标识符的一部分——上午、中午或下午)是分类特征，我们需要对这些特征进行编码，以帮助机器学习模型从中学习。基于报表类型，不同的报表具有不同的复杂性，模型将学习计算报表执行时间。如果指定了更多的报告参数，则报告的执行速度会更快(在当前的训练数据中假设，但这可以更改)-要处理的数据会更少。假设报告在中午和下午执行速度较慢(同样，如果您的场景不同，这可以更改)。
数据处理
我们需要对分类特征进行编码( report_id 和 day_part ，并对 report_params 特征进行规范化。这将有助于机器学习模型以更少的错误从数据中学习。
分类特征编码以一种简单的方式完成-查找所有唯一值，创建与唯一值一样多的列，并为当前行中与值匹配的列赋值 1。
从数据集中移除分类要素:
report_id = dataset.pop('report_id')
day_part = dataset.pop('day_part')
当值与当前行中的值匹配时，通过向数据集添加新列并设置 1 来编码 report_id 和 day_part 功能:
dataset['report_1'] = (report_id == 1)*1.0
dataset['report_2'] = (report_id == 2)*1.0
dataset['report_3'] = (report_id == 3)*1.0
dataset['report_4'] = (report_id == 4)*1.0
dataset['report_5'] = (report_id == 5)*1.0dataset['day_morning'] = (day_part == 1)*1.0
dataset['day_midday'] = (day_part == 2)*1.0
dataset['day_afternoon'] = (day_part == 3)*1.0
拆分行的子集进行测试是一个好主意，这将有助于评估模型定型质量。通常的做法是使用大约 20%的数据进行测试:
train_dataset = dataset.sample(frac=0.8,random_state=0)
test_dataset = dataset.drop(train_dataset.index)
下一步是为训练数据集生成统计数据(我们对目标变量的统计数据不感兴趣):
train_stats = train_dataset.describe()
train_stats.pop("exec_time")
train_stats = train_stats.transpose()
train_stats
创建标签(目标变量)—我们将训练模型来预测标签(执行时间):
train_labels = train_dataset.pop('exec_time')
test_labels = test_dataset.pop('exec_time')
归一化数据-将要素转换为相似的比例(数字之间不会有大的差异)，这将有助于模型更好地学习。使用先前计算的数据集统计信息完成归一化:
def norm(x):
return (x - train_stats['mean']) / train_stats['std']normed_train_data = norm(train_dataset)
normed_test_data = norm(test_dataset)
Keras 车型
模型用 Keras API 定义。建立一个具有两层数据处理和第三层数据输出的神经网络。根据我对给定训练数据的测试，当每层有 50 个单元时，模型表现良好。使用 SGD(随机梯度下降)优化器(把它想成是一种从一步到另一步改进训练的工具)。训练质量——损失，通过均方误差度量来测量。这表明训练过程进行得有多顺利:
def build_model():
model = keras.Sequential([
layers.Dense(50, activation='sigmoid', input_shape=[len(train_dataset.keys())]),
layers.Dense(50, activation='sigmoid'),
layers.Dense(1)
])optimizer = keras.optimizers.SGD(0.001)model.compile(loss='mean_squared_error',
optimizer=optimizer,
metrics=['mean_absolute_error', 'mean_squared_error'])
return model
模型训练和评估
Keras 提供了一个选项来定义提前停止回调。这有助于防止模型过度拟合。如果在 10 次迭代中没有改进，您可以指定停止训练:
model = build_model()# The patience parameter is the amount of epochs to check for improvement
early_stop = keras.callbacks.EarlyStopping(monitor='val_loss', patience=10)history = model.fit(normed_train_data, train_labels, epochs=EPOCHS,
validation_split = 0.2, batch_size=40, verbose=0, callbacks=[early_stop, PrintDot()])plot_history(history)
使用验证分割很重要，它会自动分配一部分训练数据，以便在训练过程中不断评估模型质量。
训练结果，由对照验证数据的训练误差表示:
Training results
让我们通过针对测试数据集运行模型来评估训练模型的质量:
loss, mae, mse = model.evaluate(normed_test_data, test_labels, verbose=0)print("Testing set Mean Abs Error: {:5.2f} Report Execution Time".format(mae))
结果——误差约为 4 秒，这非常合理:
Testing set Mean Abs Error: 3.65 Report Execution Time
新数据预测
我们已经准备好针对新数据运行我们的模型，新数据是训练或测试数据集中不存在的数据。我将使用类型 1 的报告，并告诉报告在下午执行(类型 3)。报告参数的数量等于 15，这样数量的参数在训练数据中不可用。让我们看看模型将如何执行，如果指定了更多的参数，它应该会选择更快的报告执行趋势。
数据帧结构:
headers = ['report_id', 'report_params', 'day_part']
dataset_input = pd.DataFrame([[1, 15, 3]],
columns=headers,
dtype=float,
index=['input'])
标准化:
report_id = dataset_input.pop('report_id')
day_part = dataset_input.pop('day_part')dataset_input['report_1'] = (report_id == 1)*1.0
dataset_input['report_2'] = (report_id == 2)*1.0
dataset_input['report_3'] = (report_id == 3)*1.0
dataset_input['report_4'] = (report_id == 4)*1.0
dataset_input['report_5'] = (report_id == 5)*1.0dataset_input['day_morning'] = (day_part == 1)*1.0
dataset_input['day_midday'] = (day_part == 2)*1.0
dataset_input['day_afternoon'] = (day_part == 3)*1.0
编码:
normed_dataset_input = norm(dataset_input)
使用定型的 Keras 模型预测报表时间执行:
res = model.predict(normed_dataset_input)
结果:429.1053772 秒
让我们看看训练数据。具有 10 个参数的类型 1 和日期时间类型 3 的报告被记录为在 431 秒内执行。在同一时间使用更多参数的相同报告执行速度更快，这意味着模型运行良好。
模型保存和重用
了解如何将训练好的模型保存到文件中以供将来重用是很有用的:
model.save("report_exec_times_model.h5")
模型还原:
from keras.models import load_modelmodelFromFile = load_model('report_exec_times_model.h5')
modelFromFile.summary()
您可以在恢复的模型上运行预测功能，就像在刚刚训练的模型上一样:
res = modelFromFile.predict(normed_dataset_input)
基巴纳的报告和数据可视化
简介
在本文中，我想根据我的个人经验展示一下 Kibana 的功能。我的目的是对 Kibana 的功能给出一个高层次的概述，并分析定制插件对旋转和报告的需求。在接下来的文章中，我希望对 Kibana 进行一次的深入探索，并更详细地介绍数据可视化过程*。*
什么是基巴纳
凭借 GitHub 上超过 11k 颗星星， Kibana 赢得了全世界开发者的心，多年来一直是弹性搜索数据可视化的最佳平台。不是没有原因的。
那么是什么让 Kibana 成为弹性搜索的必备工具呢？
它的主要目的听起来很简单，但它确实很强大:
- Kibana 旨在通过提供一个单一的接口，使与弹性堆栈的交互变得简单而省时，从而帮助您更好地理解您的数据。 Kibana 为您做了很多繁重的工作，包括通过 REST API 查询 Elasticsearch 的数据。通过这种方式，它消除了手工编写查询的需要。
- 你会喜欢的另一个特性是，Kibana 是 100% 开源的——它的源代码可以在 GitHub 库上找到，每个人都可以为它的开发做出贡献。由于这个原因，产品被积极地更新、维护和增强。
- 此外，社区非常强大并且支持——你可以在基巴纳讨论论坛上感受到这一点。
可视化能力
让我简单介绍一下可视化功能。
当您将数据导入到 Kibana 中时，神奇的事情就开始了。
有多个核心可视化，如区域、饼图、线条和条形图、直方图、日光、热度、区域和坐标图、量规、
控件值得特别关注——使用它们，您可以添加交互式输入(下拉菜单和单选滑块)并在实时中过滤仪表板的内容。
除了基本的可视化，Kibana 还支持Vega——一种高级语法，允许快速构建定制的可视化。
如果你对分析时间序列数据感兴趣，你可以利用 Timelion 提供的功能。Timelion 是一种特殊类型的可视化工具，它可以获取原始时间序列数据，并以某种方式呈现出来，从而帮助您从数据中获得可操作的见解。
一旦定义了索引模式并选择了可视化，就可以应用度量聚合和桶聚合来对符合特定搜索标准的文档进行分类。
人们可能会对各种口味和用途的 Y 轴聚合印象深刻，尤其是统计聚合:um、平均值、最小值、最大值、计数(唯一计数)、标准偏差、中值、百分位数、百分位数排名、最高命中率和地理质心。
对于 X 轴，您可以从桶聚合中受益，如日期直方图、范围、术语、过滤器和重要术语。除了聚合之外，您还可以通过应用后续的子聚合来进一步划分数据。
除了上述聚合，您还可以定义父管道和同级管道聚合。要详细了解所有可能的聚合，我推荐阅读这篇综合文章。
为了只显示符合特定标准的文档中的数据，您可以添加字段过滤器。
过滤有两个选项:
- 通过使用基巴纳查询语言 (KBL)或 Lucene 作为语法来编写搜索查询。
- 通过点击可视化元素。
添加过滤器的界面非常舒适——来看看吧:
没有任何东西限制您编写基于多个字段和条件的复杂过滤器。
可视化样本
以下是一个堆积条形图的示例，您可以创建该图来了解承运人飞往特定目的地的航班的平均延误时间:
这种可视化建立在样本飞行数据的基础上。
保存结果
一旦您对结果感到满意，您就可以保存创建的可视化。它已准备好添加到仪表板中，并与您的队友共享。可视化是完全可重复使用的。
向仪表板添加组件
下面是一个可以用 Kibana 构建的仪表板示例:
您可以通过查看底层的原始数据或聚合数据来检查每一个可视化，如果有必要，甚至可以将这些数据导出到 CSV。
为了数据分析过程的完整性，生成的仪表板可以转换为 PDF 和 PNG 格式的定制报告。
个人印象
最酷的一个优点是，Kibana 让你可以自由控制仪表盘的几乎每个方面。从提供给可视化效果的数据源到调整特定元素的单个颜色。其他强大的特性是聚合和过滤功能。
但是我特别喜欢仪表盘的元素是完全交互式和可配置的。每个组件都可以根据您的喜好定制—您可以更改颜色、图例、标题等。
增强了基巴纳语的报告功能
作为商业智能的工具，Kibana 是完美的。它的界面允许在几分钟内创建一个仪表板，并在它的帮助下分析数据。
尽管有大量的可视化工具，如果基巴纳有一个专门为高级报告目的设计的工具，那就太好了。它应该能够通过切片和切块从不同角度查看数据。在我看来，数据透视表最适合这项任务。Kibana 有一个插件，很容易安装和使用 Elasticsearch 数据。它叫 Flexmonster 。虽然最初是作为一个可以集成到任何使用 JavaScript 的应用程序中的数据透视表组件创建的，但它也可以作为 Kibana 的一部分。你可以将它连接到 Elasticsearch index ，从中获取文档并开始探索数据。
透视表的特性
Flexmonster 为有效的报告提供了多种功能，如聚合函数、过滤器、排序以及内置的将报告导出为 PDF、Excel 和其他格式的**。**
对于聚合，有 sum、count、distinct count、average、min 和 max 函数可用。
过滤选项
您可以按成员名称或报表筛选器来筛选数据，以保持报表井井有条。
但是我想特别强调的一个特性是通过 UI 完全交互地改变切片的能力。设置好一个报告后，你可以将 拖放到 网格的右边。要知道哪些记录支持聚合值，您可以钻取单元格。
一个额外的优势是组件的性能。人们可以注意到大量的行被足够平滑地渲染。
数据透视表
我喜欢的另一件事是在网格和图表模式之间切换的能力。通过单击工具栏上的按钮，您可以立即在图表中显示报告的层次结构。图表也支持过滤。
定制组件
数据透视表提供了通过在紧凑视图和经典视图之间切换来配置布局的自由。另一种定制组件外观的方法是应用主题，改变 CSS 样式和工具栏。
个人印象
今天，我尝试介绍了 Kibana 的数据可视化功能，并展示了一个名为 Flexmonster 的 web 报告插件。
在我看来， Flexmonster 补充了 Kibana 中可用的数据可视化——它们是很好的搭配。这个工具完全符合我在数据分析方面的要求。有了它，我可以专注于报告，并揭示隐藏在数据中的见解。
我建议尝试一下。希望它们能让你的数据分析大放异彩。
感谢您的阅读！
反馈
在 Kibana 中，你最喜欢哪种可视化？你使用任何第三方插件进行数据分析吗？我很高兴听到你的经历。
扫一扫
721
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
