原文：TowardsDataScience Blog

协议：CC BY-NC-SA 4.0

原文：https://towardsdatascience.com/predicting-cancer-with-logistic-regression-in-python-7b203ace16bc?source=collection_archive---------18-----------------------

Source

在我的第一个逻辑回归分析中，我们仅仅触及了表面。讨论的只是高级概念和双变量模型示例。在这一分析中，我们将着眼于更具挑战性的数据，并学习更先进的技术和解释。

1.数据背景

2.数据探索/清理

3.数据可视化

4.构建模型

5.测试模型

测量体内某些蛋白质的水平已经被证明可以预测癌症的发展。医生可以进行测试来检查这些蛋白质的水平。我们有 255 名患者的样本，希望获得关于 4 种蛋白质及其与癌症生长的潜在关系的信息。

我们知道:

• 每位患者测得的每种蛋白质的浓度。
• 每个患者是否已经被诊断患有癌症(0 =没有癌症；1=癌症)。

我们的目标是:

通过从我们样本中的蛋白质水平和癌症之间的关系中提取信息来预测未来的患者是否患有癌症。

我们将关注的 4 种蛋白质:

甲胎蛋白

癌胚抗原

癌抗原 125

癌抗原 50

我从 UAB 的 MBA 项目 @ 那里得到了这套用于教育目的的数据。

让我们通过引入数据和导入必要的模块来开始分析。

%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
import seaborn as sns
df = pd.read_excel(r"C:UsersAndrewDesktopcea.xlsx")
df.head()

Figure 1

df.head()为我们提供了数据集的前 5 个特征。每行是一名患者，每列包含一个描述性属性。类别(Y)描述了患者是没有癌症(0)还是患有癌症(1)。接下来的 4 列是在病人血液中发现的蛋白质水平。

df.describe()

Figure 2

我们可以从描述方法中检索关于样本的一些基本信息。该数据集中有 255 行，具有不同的标准差和均值。值得一提的是，与平均值相比，这 4 种蛋白质都具有较低的 50%值。这意味着大多数蛋白质水平很低。这 4 种蛋白质还具有高范围，这意味着存在高值异常值。比如我们看 AFP，均值 4.58，中位数 3.86，最高值 82.26。

df.isnull().sum()

检查数据中是否有空值总是好的。这可以通过多种方式解决。幸运的是，这次我们不用担心空值。

让我们想象一下我们的每个变量，并根据我们所看到的假设可能会发生什么:

yhist = plt.hist('class (Y)', data = df, color='g')
plt.xlabel('Diagnosis (1) vs no diagnosis (0)')
plt.ylabel('Quantity')
plt.title('Class (Y) Distribution')

Figure 3

这个样本中有更多的无癌患者。如图 2 所示，“类别(Y)”变量的平均值为 0.44。

**#setting the axes**
axes = plt.axes()
axes.set_xlim([0,(df['AFP'].max())])**#making histogram with 20 bins**
plt.hist('AFP', data = df, bins = 20)
plt.xlabel('AFP Level')
plt.ylabel('Quantity')
plt.title('AFP Distribution')

Figure 4

如前所述，甲胎蛋白水平低的患者相对较多。

**#color palette**
pal = sns.color_palette("Set1")**#setting variable for max level of protein in dataset**
lim = df['AFP'].max()**#setting bin size to be 20**
bins = np.linspace(0,lim,(lim/(lim*(1/20))))**#creating new column in df with bin categories per feature**
df['AFP_binned'] = pd.cut(df['AFP'], bins)**#creating a crosstab stacked bar chart variable**
chart = pd.crosstab(df['AFP_binned'],df['class (Y)'])**#normalizing chart and plotting chart**
chart.div(chart.sum(1).astype(float), axis=0).plot(kind='bar', color = pal,stacked=True)
plt.xlabel('Bins')
plt.ylabel('Quantity')
plt.title('Normalized Stacked Bar Chart: AFP vs Class(Y)')

Figure 5

作为分布直方图的补充，上面的堆积条形图显示了随着 AFP 水平的增加，1 与 0 的比例。该图表与分布直方图一样，也分为 20 个区间。结合我们从上述目标变量的分布和比例中获得的知识，我们可以直观地确定从该蛋白质中可能没有获得太多的预测知识。

让我们一步一步来。大多数患者的 AFP 值低于 10，如图 4 中的前 2 条所示。因为大多数患者都在前两个柱中，所以图 5 中它们之间 Y 的变化比其他柱中 Y 的变化更重要。

从柱 1 到柱 2，癌症患者的比例略有增加。癌症患者的比例从第 2 栏下降到第 3 栏。bar 3 之后，剩下来分析的患者少之又少，对趋势影响不大。从这里我们可以看到，目标变量看起来基本上与 AFP 的变化无关。最显著的变化(条 1 到 2)非常轻微，之后的变化不在同一方向。

让我们看看其他蛋白质是什么样子的。

Figure 6

Figure 7

东航似乎有一个不同的故事。图 6 显示了与 AFP 相似的分布形状；然而，图 7 显示了癌症发病率的不同变化。就像 AFP(由于分布形状)一样，柱之间最显著的癌症变化在柱 1 和柱 2 之间。

从柱 1 到柱 2 的变化从大约 63%的非癌性变为 18%的非癌性(或者换句话说，37%的癌性变为 82%的癌性)。此外，从 bin 2 到 bin 3 的变化是相同的方向，更多的癌症。从具有 100%癌症的 bin 5 开始的异常值加强了较高 CEA 可能指示癌症的趋势。

Figure 8

Figure 9

CA125 有点棘手。柱 1-2 表明，像癌胚抗原一样，这种蛋白的高水平可能导致癌症。然而，这一次似乎有两种趋势。随着几乎所有的后一类人群变成非癌症人群，这一趋势发生逆转。我们稍后将更详细地研究这个变量。

Figure 10

Figure 11

CA50 看起来没什么前途。前 4 个条柱似乎表明了较高癌症发病率的趋势。然而，趋势看起来在第 7-9 条中发生了逆转。CA50 水平和癌症之间的关系可能很小或可以忽略。

让我们把这个模型放在一起，看看回归能告诉我们什么。

**#importing module**
import statsmodels.api as sm**#setting up X and y**
cols= [‘AFP’,’CEA’,’CA125',’CA50']
X= df[cols]
y = df[‘class (Y)’]**#filling in the statsmodels Logit method**
logit_model = sm.Logit(y,X)
result = logit_model.fit()
print(result.summary())

Figure 12

突出显示的值是本报告中重要的内容:我们的 4 个独立变量及其 p 值。AFP 和 CA50 的 p 值较高。如果我们的α值为 0.05，那么 AFP 和 CA50 的值太高，无法拒绝我们的无效假设(我们的无效假设是蛋白质水平对癌症发病率没有影响)。然而，CEA 和 CA125 通过了测试，并且被确定为显著的。AFP 和 CA50 都是基于我们在堆积条形图上看到的假设而被忽略的，所以这是有道理的。让我们去掉这些变量，再次进行回归分析:

**#deleted p values above the 0.05 alpha threshold**
cols= ['CEA','CA125']
X= df[cols]
y = df['class (Y)']
logit_model = sm.Logit(y,X)
result = logit_model.fit()
print(result.summary())

Figure 13

有了我们的最终系数，我们对每个剩余蛋白质和癌症之间的关系有了更多的了解。CEA 的阳性关系比 CA125 的阴性关系强约 3 倍。随着癌胚抗原的增加，癌症的可能性增加。随着 CA125 的增加，癌症的可能性降低。

我们将把样本数据分为训练和测试，以测试我们的回归结果。

from sklearn.linear_model import LogisticRegression
from sklearn import metrics**#shuffling df**
df = df.sample(frac=1).reset_index(drop=True)**#redefining columns** 
cols= ['CEA','CA125']
X= df[cols]
y = df['class (Y)']**#Dividing into training(70%) and testing(30%)**
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=0)**#Running new regression on training data**
logreg = LogisticRegression()
logreg.fit(X_train, y_train)**#Calculating the accuracy of the training model on the testing data**
accuracy = logreg.score(X_test, y_test)
print('The accuracy is: ' + str(accuracy *100) + '%')

形象化上面计算的精确度的一个好方法是使用混淆矩阵。下面是混淆矩阵的概念框架。

编辑:我和一个生物统计的朋友谈论我的分析，这个领域的惯例是这种疾病被认为是阳性的。我武断地将癌症设定为阴性，因为我当时并不知道这一点。

Figure 14

“迷惑”是很多人的关键词。试着一次看一行:第一行是一个好的起点。这一行告诉我们有多少实例被预测为良性的。如果我们查看这些列，我们可以看到该预测中实际值的拆分。请记住，行是预测值，列是实际值。

from sklearn.metrics import confusion_matrix
confusion_matrix = confusion_matrix(y_test, y_pred)
print(confusion_matrix)

将上面的矩阵与图 14 进行匹配，了解其含义:

• 我们模型的 39 个猜测是真阳性:模型认为病人没有癌症，他们确实没有癌症。
• 我们的模型的猜测中有 18 个是真阴性:模型认为病人患有癌症，他们确实患有癌症。
• 14 我们的模型的猜测是假阴性:模型以为病人得了癌症，但实际上他们并没有得癌症
• 6 我们的模型的猜测是假阳性:模型认为病人没有癌症，但他们实际上患了癌症。

我们总数据的 30%给了测试组，剩下 255(.3) = 77 个被测试的实例。矩阵的和是 77。将“真实”的数字除以总数，将得到我们模型的精确度:57/77 = 74.03%。

请记住，在执行这个测试之前，我们随机打乱了数据。我运行了几次回归，得到了 65%到 85%的准确率。

最后，我们将执行接收器操作特性 (ROC)分析，作为测试我们模型的另一种方式。该测试的两个目的是

1. 确定最佳“截止”点在哪里。
2. 通过另一个称为“曲线下面积”(AUC)的指标来确定模型的分类效果。

我们将从头开始创建我们的 ROC 曲线。以下是用于格式化新数据框架以计算 ROC、截止点和 AUC 的所有代码。

**#Formatting y_test and y_predicted probabilities for ROC curve**
y_pred_prob = pd.DataFrame(y_pred_prob)
y_1_prob = y_pred_prob[1]
y_test_1 = y_test.reset_index()
y_test_1 = y_test_1['class (Y)']**#Forming new df for ROC Curve and Accuracy curve**
df = pd.DataFrame({ 'y_test': y_test_1, 'model_probability': y_1_prob})
df = df.sort_values('model_probability')**#Creating 'True Positive', 'False Positive', 'True Negative' and 'False Negative' columns** 
df['tp'] = (df['y_test'] == int(0)).cumsum()
df['fp'] = (df['y_test'] == int(1)).cumsum()
total_0s = df['y_test'].sum()
total_1s = abs(total_0s - len(df))
df['total_1s'] = total_1s
df['total_0s']= total_0s
df['total_instances'] = df['total_1s'] + df['total_0s']
df['tn'] = df['total_0s'] - df['fp']
df['fn'] = df['total_1s'] - df['tp']
df['fp_rate'] = df['fp'] / df['total_0s']
df['tp_rate'] = df['tp'] / df['total_1s']**#Calculating accuracy column**
df['accuracy'] = (df['tp'] + df['tn']) / (df['total_1s'] + df['total_0s'])**#Deleting unnecessary columns**
df.reset_index(inplace = True)
del df['total_1s']
del df['total_0s']
del df['total_instances']
del df['index']
df

为了理解下面的数据框中发生了什么，让我们一行一行地分析它。

• 索引:这个数据帧是按照 model_probability 排序的，为了方便起见，我重新编制了索引。
• CA125 和 CEA :蛋白质水平原始检测数据。
• model_probability :该列来自我们的训练数据的逻辑模型，该模型基于输入的测试蛋白质水平输出其被分类为“1”(癌性)的概率预测。第一行是最不可能被归类为癌症的实例，其 CA125 水平高而 CEA 水平低。
• y_test :我们用来检查模型性能的测试数据的实际分类。

其余的专栏仅仅基于“y_test”，而不是我们的模型的预测。把这些值想象成它们自己的混淆矩阵。这些将有助于我们确定以后的最佳分界点。

• tp(真阳性):该列从 0 开始。如果 y_test 为‘0’(良性)，则该值增加 1。它是所有潜在真阳性的累积追踪器。第一行就是一个例子。
• fp(假阳性):该列从 0 开始。如果 y_test 为‘1’(癌变)，则该值增加 1。它是所有潜在误报的累积跟踪器。第四行就是一个例子。
• tn(真阴性):该列从 32 开始(测试集中 1 的总数)。如果 y_test 为‘1’(癌变)，则该值减少 1。它是所有潜在真阴性的累积追踪器。第四行就是一个例子。
• fn(假阴性):该列从 45 开始(测试集中 0 的总数)。如果 y_test 为‘0’(良性)，则该值减少 1。它是所有潜在假阴性的累积跟踪器。第四行就是一个例子。
• FP _ Rate(False Positive Rate):这是通过获取行的误报计数并除以总的阳性数(在我们的例子中是 45)计算出来的。它让我们知道我们可以通过在该行设置分界点来分类的假阳性的数量。我们希望尽可能降低成本。
• tp_rate(真阳性率):也称为敏感度，计算方法是获取行的真阳性计数，然后除以阳性总数。它让我们知道我们可以通过在那一行设置分界点来分类的真阳性的数量。我们希望尽可能地保持高水平。
• 准确性:真阳性和真阴性的总和除以总实例数(在我们的例子中是 77)。我们一行一行地计算潜在的精确度，基于我们混淆矩阵的可能性。

Figure 15

我粘贴了整个数据帧，因为它值得研究一段时间，并从所有移动的片段中找出意义。看完之后，试着找出最高的准确率。如果您可以找到它，您可以将其与相应的 model_probability 进行匹配，以发现我们数据的最佳分界点。

**#Plotting**
plt.plot(df['model_probability'],df['accuracy'], color = 'c')
plt.xlabel('Model Probability')
plt.ylabel('Accuracy')
plt.title('Optimal Cutoff')**#Arrow and Star**
plt.plot(0.535612, 0.753247, 'r*')
ax = plt.axes()
ax.arrow(0.41, 0.625, 0.1, 0.1, head_width=0.01, head_length=0.01, fc='k', ec='k')
plt.show()

Figure 16

模型概率为 54%,其中精确度最高为 75%。这可能看起来违背直觉，但这意味着如果我们在将患者归类为癌症患者时使用 54%而不是 50%,这实际上会更准确。如果我们想最大限度地提高准确性，我们应该将阈值设置为 54%，但是，由于癌症的极端性质，将阈值降低到 50%以下可能是明智的，以确保无论如何都可以检查出可能患有癌症的患者。换句话说，当涉及到癌症时，假阳性比假阴性更重要！

最后，让我们绘制 ROC 曲线并找出 AUC:

**#Calculating AUC**
AUC = 1-(np.trapz(df[‘fp_rate’], df[‘tp_rate’]))**#Plotting ROC/AUC graph**
plt.plot(df[‘fp_rate’], df[‘tp_rate’], color = ‘k’, label=’ROC Curve (AUC = %0.2f)’ % AUC)**#Plotting AUC=0.5 red line**
plt.plot([0, 1], [0, 1],’r — ‘)
plt.xlabel(‘False Positive Rate’)
plt.ylabel(‘True Positive Rate (Sensitivity)’)
plt.title(‘Receiver operating characteristic’)
plt.legend(loc=”lower right”)

Figure 17

黑色的 ROC 曲线显示了我们的测试数据的真阳性率和假阳性率之间的权衡。穿过图表中心的红色虚线是为了提供一种最坏可能模型看起来像 ROC 曲线的感觉。

ROC 线越靠近左上方，我们的模型越有预测性。它越像红色虚线，预测性就越低。

这就是曲线下面积(AUC)的由来。AUC，顾名思义，是 ROC 曲线下的空间面积。直观上，这个值越接近 1，我们的分类模型就越好。虚线的 AUC 是 0.5。完美模型的 AUC 应该是 1。我们的 AUC 为 0.82，相当不错。

如果您觉得这很有帮助，请订阅。

我的其他文章，如果你想了解更多:

原文：https://towardsdatascience.com/predicting-carcinogens-with-logistic-regression-knn-gradient-boosting-and-molecular-e7952294a08c?source=collection_archive---------42-----------------------

T21 世纪毒理学计划(Tox21)是几个联邦机构之间的独特合作，旨在开发快速测试物质是否对人类健康产生不利影响的新方法。Tox21 中检测的物质包括各种产品，如:商业化学品、杀虫剂、食品添加剂/污染物和医用化合物。

p53 基因编码一种同名的蛋白质，被称为肿瘤抑制蛋白。当细胞经历 DNA 损伤时，p53 蛋白在细胞中表达，这可以将正常细胞转化为癌细胞。为了抵消这种影响，p53 可以导致生长停滞、修复 DNA 或开始细胞死亡过程。因此，当 DNA 损伤发生时，p53 表达显著增加。这种蛋白质表达的增加是不规则细胞健康的良好指标。Tox21 数据是通过测试在 p53 细胞机制控制下产生荧光报告基因产物的细胞系产生的。通过测量针对各种化合物的报告基因产物的水平，研究人员能够确定化合物是否是 p53 途径的激动剂(激活剂)。

指纹是一种将分子结构和性质表示为二进制位串(0 和 1)的方法。这种表示最初被开发并应用于搜索具有特定子结构的分子的数据库，但它也可以应用于机器学习。哈希函数是一个随机数生成器，应用于分子的每个特征，例如存在的键和分子的类型，这意味着它们充当函数的种子。

下图说明了如何通过生成分子指纹来评估两个分子的相似性(使用 Tanimoto 指数)。每个分子首先应用哈希函数，然后根据特征生成指纹。下图中的指纹生成器查看某个焊接距离半径以及该距离内的特征。

[1]

指纹的种类

这个项目将使用四种不同类型的指纹，并比较它们的预测能力。

• 摩根圆形指纹— 该指纹属于扩展连接指纹(ECFPs)家族，使用摩根算法生成[2，3]。这些指纹通过圆形原子邻域(键半径)来表示分子结构和亚结构的存在。另一个重要的特征是它们能够确定分子功能性的存在与否，这可以在对分子进行分类时进一步帮助区分。

• 类似日光的指纹 —该指纹生成器(使用 RDKit)生成的指纹与使用日光指纹算法生成的指纹相似。一个过于简单的解释:该算法沿着分子内的路径(拓扑路径)散列键[4][5]。下图显示了如何识别和描述焊接路径。

[6]

• 原子对指纹 —顾名思义，原子对指纹是使用原子对(作为特征)及其拓扑距离构建的。通过首先识别重原子和它们之间的最短距离，生成如下所示的原子对指纹。一对原子中的单个原子的特征(如键的数量)被编码。然后，这些编码特征被转换成位串，并表示为一个数字。然后，这个串联的整数串被传递给哈希函数，哈希函数将其赋值为 1 或 0。

[7]

• 拓扑扭转指纹— 是通过识别四个非氢原子的键路径而产生的 2D 结构指纹。然后这些四路碎片有了它们的特征，比如计算出的每个原子的 π 电子数、原子类型和非氢分支数。每个原子由 8 位来表征，这 4 个原子将它们的位存储在 32 位数组中。

Tox21 数据已经标记了活动(1)和非活动(0)状态，因此除了加载、格式化列名和为每个分子生成 mol 文件(这些文件通常被归类为包含分子数据信息、原子、键、坐标和连接信息的数据文件)之外，在此步骤中实际上没有太多事情要做。

要继续操作，请编辑数据在机器上的路径。我们加载一个 sdf 文件——它是一个结构数据文件,包含化合物的相关数据。我们实际上不需要这些细节，因为它们不会增加 p53 激动剂的分类，因此我们放弃了它们。例如,“公式列”没有给我们任何有用的信息，因为它是经验公式，不包含 SMILES 字符串提供的有价值的结构信息。我们确实需要分子的摩尔表示，所以我们生成了它们。

我们还应该确保数据中没有任何重复。特别是这个数据集，总共有 8629 条记录，只有 6826 条唯一记录，这意味着大约 20%的数据是重复的，需要丢弃。

指纹生成

这是一个使用摩根算法生成指纹的例子。为“mol”列中的每种化合物生成指纹，半径为 2，位长为 2048。还为每个 mol 生成具有相同参数的列表。然后，for 循环遍历该列表，创建一个 numpy 数组，并在名称中添加一个由“np”表示的新列表。当运行分类器时，这些列表将成为我们的 x 变量。

设置 x 和 y 变量:

当分类标签(在我们的例子中:活动= 1；inactive = 0)是成比例偏斜的，偏差被引入到模型中。例如，如果数据集包含 95%的活性标记，那么模型可能将非活性分子分类为活性分子，简单地说，准确度将反映类别分布。确定阶级不平衡是否是直接的——计算标签的分布。一旦我们确定了少数民族在我们的数据集中有多普遍，我们就可以继续平衡它。

上面的计数图显示了我们的数据中存在不成比例的类分布，比例为 15:1(非活动与活动),这意味着大约 93.84%的类是非活动的。

有两种方法可以平衡数据:过采样和欠采样。过采样通过使用多种算法之一生成合成数据点来实现平衡。我们今天要用的是 SMOTE(合成少数过采样技术)算法的一个变种。SMOTE 的工作原理是搜索给定点数据的 k-最近邻，然后沿着邻点之间的线生成数据点。由于这些点落在真实数据点之间的线上，所以存在一些相关性，这就是为什么我们将使用 ADASYN 算法(自适应合成)。ADASYN 的工作方式基本上与 SMOTE 相同，但是它向点添加随机生成的小值，以增加方差并模拟更真实的数据。

• 注意 —我选择不使用欠采样，因为它通过从多数类中移除数据点来平衡类，因此有数据损失，可能影响分类器的区分能力。

[9]

在这里，我们将 ADASYN 算法应用于每个 x 变量，这将最终导致每个指纹的不同少数类计数，但范围很窄，计数的变化很小，所以这不会造成问题。

我们可以绘制一个重新采样的分布图来直观地表示类别平衡:

为了训练、测试和再次测试，我们需要创建三个独立的数据分区。训练数据将用于训练我们的模型——本质上，模型将使用该数据来学习哪些指纹(2048 位模式)可能属于 p53 激动剂。测试集将用于评估我们模型的准确性，最后一个验证集将用作无偏数据集来测试我们模型的预测能力。

我们的数据将被分成 85/15 训练/测试集，然后训练数据将被进一步分成 85/15 训练/验证集。

逻辑回归算法根据分类反应(结果)变量与预测特征的关系，将其分类为 1 到 0。相反，线性回归输出的响应变量是连续的，可以是任何实数。最重要的是，线性回归不输出概率，而是拟合最佳超平面。因此，逻辑回归是像我们这样的二元分类问题的自然选择。

交叉验证是一种统计技术，用于降低数据过度拟合的风险。当分类器记忆数据并拟合噪声和误差而不是变量之间的潜在关系时，发生过拟合。如果模型过于复杂，并且用于训练的数据有限，就会出现这种情况。我们的模型似乎没有这两个限制，但无论如何，执行交叉验证是一个好的实践。

上面的图片展示了一个模型是如何“记忆”和拟合数据的。下图显示了交叉验证是如何创建 10 个折叠的，每个折叠都被分成测试集和训练集以适应模型，然后对每个折叠的误差进行平均。

Over-fitting (source) | cross validation (source)

以下代码行显示了如何执行逻辑回归，获得每个折叠的交叉验证分数(只需对训练数据进行一次迭代，以使用选择的分类器估计分数)。然后，我们在执行 k-fold (10)交叉验证的同时，对 1000 次迭代的逻辑回归进行建模。

对于其余的指纹，我们可以简单地编写一个函数来模拟数据:

测试集预测和混淆矩阵

混淆矩阵是描述训练模型对测试数据的性能的表格。矩阵由四个象限组成:

• 真阳性(TP): 这些是预测为阳性(活性 p53 激动剂)并且实际上是阳性的化合物。
• 真阴性(TN): 这些化合物被预测是无活性的，实际上是无活性的。
• 假阳性(FP): 预测化合物有活性但实际上没有活性的情况(I 型错误)。
• 假阴性(FN): 预测化合物无活性但实际上有活性的情况(II 型错误)

我们可以绘制测试集和验证集预测的混淆矩阵以及准确性，准确性通过(TP+TN)/总案例来确定。左边的混淆矩阵是测试集预测的混淆矩阵，右边的是验证集的混淆矩阵。

k-最近邻算法(knn)取一个数据点并查看 k 个最近的数据点(k =7 意味着 7 个最近的)，然后根据 k 个最近点的多数标签将该数据点标记为 0 或 1。例如，如果给定一个 k=7 的数据点，并且最接近的点恰好是五个 0 和两个 1，那么该点将被分类为 0(无效)。

Knn 在处理较少的特征时表现最佳，从而降低了过度拟合的可能性。正常情况下，我们必须执行 PCA 或其他形式的降维，然而，由于我们使用的是单一的预测特征，我们不需要这样做。

knn

通过交叉验证执行网格搜索，为我们的数据确定最佳邻居数量。我们将测试 9 个候选数字(3、5、7、9、11、13、15、17、19)，并执行 5 重交叉验证。然而，这可能需要一段时间，因为 9 个邻居 x 每个邻居的 5 倍验证等于 45 次拟合。k 倍可以减少，这样我们就可以保持一个体面的邻居范围来尝试。

注意:将“n_jobs”设置为-1——这允许在最大 CPU 使用率下进行并行处理

网格搜索建议使用的最佳邻居数量为 3。

为了测试 k 个邻居数的范围，我们可以用不同的 k 值来拟合不同的指纹数据。使用原子对训练数据，我测试了 k = 3，5，7，9，这按顺序产生了以下准确度分数:

0.846
0.846
0.815
0.785

看起来 3 和 5 有相同的分数，但是 3 可能是因为计算时间较少而被选中的。最后，用 k=3 训练模型，首先计算测试数据的准确度，然后计算验证数据的准确度

到目前为止，我们已经为每个分类器的训练数据拟合了单个模型。梯度增强将许多单独的模型结合起来，创建一个精确的模型，这就是所谓的集合。增强是用于创建合奏的方法。它通过将数据与初始模型拟合来工作，然后建立后续模型，该模型用于准确预测初始模型未正确分类的标签。本质上，每个后续模型都致力于最小化前一个模型的预测误差，这导致预测误差的总体降低。因此，最终模型做出的预测是所有先前模型计算的预测的加权和的结果。

By Alex Rogozhnikov (source)

像我们之前做的那样，执行网格搜索来调整参数:

检查测试和验证数据的准确性

曲线下面积(AUC)值是评估模型在二元分类任务上的性能的另一种度量，并且是最广泛使用的评估指标之一。计算 AUC 需要给定模型的两个属性——敏感性和特异性。敏感度是正确分类的阳性数据点与所有分类为阳性的数据点的比例。敏感度是被正确分类为阴性的数据点与被分类为阴性的所有数据点的比例。

1. 一个模型的灵敏度:真阳性/(假阴性+真阳性)。
2. 模型的特异性:真阴性/(假阳性+真阴性)

当真阳性率相对于假阳性率作图时，可以计算该曲线的 AUC，并从本质上对模型区分两类的能力进行评分。较高的 AUC 意味着归类为 p53 激动剂的分子确实更有可能是激动剂。

模型的 f1 分数是使用回忆(像敏感度一样计算)和精确度计算的，它是以下各项的比率:

• 真阳性/(假阳性+真阳性)

F1 分数可通过以下方式获得:

F1 分数在模型召回率(较少的预测点被证明是假阴性)和精确度(较少的预测点被证明是假阳性)之间找到了平衡。

下表显示了使用四个指纹时每个分类器的预测能力。包括的指标有训练、测试和验证数据集的准确度分数、验证数据的 AUC 分数以及验证数据的 f1 分数。对于每个指纹和分类器，最高分以绿色突出显示。

就验证准确性而言，逻辑回归得分一直最高。接下来，当评估 AUC(可以说是最重要的指标)时，梯度增强产生了最高的分数，但是逻辑回归具有最高的平均分数。最后，最高的个体和平均 f1 分数属于逻辑回归，表明它产生的模型在精度和稳健性之间达到了最大的平衡。

赢家:逻辑回归

文献表明，尽管拓扑扭曲可能不如日光样指纹和摩根指纹那样受欢迎，但与其他指纹相比，它们表现出更高的性能[8]。

结果显示，日光样指纹在逻辑回归和 knn 中都具有最高的 AUC 得分——当使用梯度增强时，原子对指纹得分仅略微高一些。然而，当考虑梯度增强时，原子对指纹在验证准确性、AUC 和 F1 分数方面超过了日光样指纹。

获胜者:日光般的指纹

由于相似的分子亚结构，指纹中有时会发生图案之间的意外碰撞。当相同的位由多个模式设置时，会发生这种情况(将两个相似结构的特定特征散列为相同的位)。例如，拓扑扭转指纹产生了八个相同的指纹，这可以通过增加位长(从 2048 到 4096)来避免。

[1] Lo，Y.-C .，Rensi，S. E .，Torng，w .，& Altman，R. B. (2018 年)。化学信息学和药物发现中的机器学习。今日药物发现，23(8)，1538–1546。https://doi.org/10.1016/j.drudis.2018.05.010

[2] Morgan，H. L.《化学结构的独特机器描述的产生——化学文摘服务处开发的一种技术》。化学博士。医生。1965, 5: 107–112.

[3]罗杰斯，d；扩展连接指纹。化学博士。Inf。模型。2010, 50(5): 742–754.

[4]https://www.rdkit.org/docs/GettingStartedInPython.html

[5]https://www . daylight . com/day html/doc/theory/theory . finger . html

[6]https://docs . eyes open . com/toolkits/python/graphsimtk/fingerprint . html

[7] Jelínek，j .，koda，p .，& Hoksza，D. (2017 年)。利用氨基酸结构邻域知识库预测蛋白质相互作用位点。BMC 生物信息学，18(S15)。https://doi.org/10.1186/s12859-017-1921-4

[8]斯柯达公司和霍克斯扎公司(2015 年)。拓扑挠指纹探索。2015 年 IEEE 生物信息学和生物医学国际会议(BIBM)。doi:10.1109/bibm。58606.88868688666

[9] Lee，h .，Kim，j .，& Kim，S. (2017 年)。基于高斯的 SMOTE 算法求解偏斜类分布。国际模糊逻辑和智能系统杂志，17(4)，229–234。https://doi.org/10.5391/ijfis.2017.17.4.229

原文：https://towardsdatascience.com/predicting-click-through-rate-for-a-website-7cd2a892d26e?source=collection_archive---------4-----------------------

想象一个用户访问一个网站，并执行工作搜索。从显示的结果集合中，用户点击他/她感兴趣的某些结果，并且在检查工作描述之后，她进一步点击其中的应用按钮以进入应用页面。申请率被定义为申请的分数(在访问职位描述页面后)，目标是使用下一节中描述的数据集来预测此指标。

这篇文章将为任何机器学习新手提供完整的指南。我的目标是为您提供一个应用机器学习的端到端蓝图，同时尽可能保持它的可操作性和简洁。

1. 数据收集
2. 探索性数据分析 —首先，“了解”数据。这一步应该快速、高效、果断。
3. 数据清理——然后，清理你的数据以避免许多常见的陷阱。更好的数据胜过更好的算法。
4. 功能工程 —接下来，通过创建新功能来帮助您的算法“专注于”重要的东西。
5. 算法选择 —选择最佳、最合适的算法，不浪费您的时间
6. 模型训练和调优 —最后，训练你的模型。一旦你完成了前 5 步，这一步就相当公式化了
7. 洞察力 —最终将推动业务发展的最终结果

你可以从这里下载数据。

Dataset

数据集中的每一行都对应于一个职务列表的用户视图。它有 10 列，如下所述。
1。职位接近度 tf idf:衡量查询和职位的接近度。
2。描述接近度 tf idf:衡量查询和工作描述的接近度。
3。主查询 tf idf:与用户查询与职位名称和职位描述的接近程度相关的分数。
4。查询 jl 得分:衡量查询和职务列表对的流行程度。
5。查询职称得分:衡量查询和职称对的受欢迎程度。
6。城市匹配:指明职务列表是否与用户(或用户指定的)位置匹配。7。职务年龄天数:指明已发布职务列表的年龄。8。申请:指明用户是否已申请此职务列表。9。太平洋搜索日期:活动日期。10。类别 id:单击的职位的类别 ID。

这个问题有两个部分。

1. 我们必须只关注前 7 列，并使用这些功能来预测有多少用户申请该网站。
2. 我们必须考虑将最后一列添加到特征集(“类 id”)，并检查分类性能是否提高。

探索性分析的目的是**“了解”**数据集。预先这样做将使项目的其余部分更加顺利，主要有三个方面:

1. 您将获得关于数据清理的有价值的提示(这将决定您的模型的成败)。
2. 你会想到特征工程的想法(可以让你的模型从好变得更好)。
3. 您将获得对数据集的“感觉”，这将有助于您交流结果并产生更大的影响。

但是，机器学习的探索性分析应该是快速、高效、果断 …而不是冗长！

不要跳过这一步，但也不要在这一步上卡住。

你看，有无限可能的图、图表和表格，但是你只需要一把的来“了解”足够好的数据来使用它。

在这一步中，我们将向您展示为您的投资带来最大收益的可视化效果。

首先，您需要回答一组关于数据集的基本问题:

• 我有多少观察值？
• 有多少功能？
• 我的要素的数据类型是什么？它们是数字吗？绝对的？
• 我有目标变量吗？

该数据集总共包含 1200890 个观察值，其中 789586 个观察值是在 2018 年 1 月 21 日至 2018 年 1 月 26 日之间收集的，而其余的观察值是在 2018 年 1 月 27 日收集的。总共有 9 个特性和 1 个目标变量(“应用”)。没有明确的特征。

Dataset info

接下来，我们可以检查数据集中有多少缺失值。从上表可以清楚地看出，标题邻近 tfidf、描述邻近 tfidf 和城市匹配包含空值。

由于该数据集是关于有多少访问网站的客户点击应用，它可能不是很多，我们可以假设这将是一个不平衡的数据。让我们看看我们是否正确。

Data Distribution (1 = Apply, 0 = Not Apply)

我们可以看到，这种分布确实是不平衡的。为了处理这种情况，有各种技术，如欠采样、过采样、SMOTE，使用 AUC 等指标。但这些将在以后讨论。

Data Description

以下是该表中的观察结果:

1. 大多数预测值的 75%瓦片和最大值之间存在显著差异。
2. ’ title_proximity_tfidf ‘，’ description_proximity_tfidf ‘，’ main_query_tfidf ‘，’ query_jl_score ‘，’ query_title_score ‘，’ job_age_days '的中值低于平均值
3. 因此，观察值 1 和 2 表明数据中存在大量异常值

相关性允许您查看数字特征和其他数字特征之间的关系。

相关性是一个介于-1 和 1 之间的值，表示两个要素协调移动的程度。你不需要记住数学来计算它们。只要知道以下直觉:

1. 正相关意味着随着一个特征的增加，另一个特征也增加。例如孩子的年龄和身高。
2. 负相关性意味着一个特征增加，另一个特征减少。例如花在学习上的时间和参加的聚会。
3. 接近-1 或 1 的相关性表示强关系。
4. 那些接近 0 的表示弱关系。
5. 0 表示没有关系

Correlation between variables

大多数特征彼此之间并不高度相关。仅有的两个具有某种相关性的特征是主查询 tfidf 和标题接近度 tfidf(大约 0.8)。

Outliers

探索性数据分析中最重要的步骤之一是异常值检测和处理。机器学习算法对数据点的范围和分布非常敏感。数据异常值会欺骗训练过程，导致训练时间更长，模型更不准确。离群值被定义为与剩余数据显著不同的样本。这些点位于分布的整体模式之外。均值、方差和相关性等统计指标很容易受到异常值的影响。

异常值的性质:

由于以下原因之一，数据集中可能会出现异常值，

1. 数据集中真正的极高值和极低值
2. 由于人为或机械错误而引入
3. 通过替换缺失值引入

异常值检测

• 极值分析
• z 分数法
• k 表示基于聚类的方法
• 可视化数据
• 箱线图

异常值处理

• 均值/中值或随机插补
• 整理
• 顶部、底部和零编码
• [数]离散化

然而，在本文中，我将使用箱线图方法检测异常值。如果您想深入了解如何检测和处理异常值，请参考这篇文章。箱线图也称为晶须图，是一种图形方法，通常用四分位数和四分位数间距来描述，有助于定义上限和下限，超出上限和下限的任何数据都将被视为异常值。

简而言之，分位数是分布中与该分布中值的等级顺序相关的点。对于给定的样本，您可以通过对样本进行排序来找到任何分位数。排序样本的中间值是中间分位数或第 50 个百分位数(也称为样本的中位数)。

下面是我们数据集的箱线图。

如您所见，数据集中有许多异常值。删除它们可能会导致重要信息的丢失。

数据清理是每个人都在做但没有人真正谈论的事情之一。当然，这不是机器学习中“最性感”的部分。不，没有隐藏的技巧和秘密要揭开。

然而，适当的数据清理可以成就或破坏您的项目。专业数据科学家通常会在这一步花费大量时间。

为什么？因为机器学习中一个简单的道理:

更好的数据胜过更好的算法。

换句话说……垃圾进来，垃圾出去。即使你忘记了本课程的其他内容，请记住这一点。

事实上，如果您有一个适当清理的数据集，即使简单的算法也可以从数据中获得令人印象深刻的见解！

显然，不同类型的数据需要不同类型的清理。然而，本文中介绍的系统化方法总是可以作为一个很好的起点。

数据清理的第一步是从数据集中移除不需要的观测值。这包括重复的或不相关的观察值。我们的数据集包含相当多的重复条目，这些条目将被删除。

在应用机器学习中，缺失数据是一个看似棘手的问题。

首先，明确一点， 你不能简单地忽略数据集中的缺失值。你必须以某种方式处理它们，因为大多数算法不接受缺失值。

“常识”在这里是不理智的

以下是处理缺失数据的最常用方法:

1. 删除有缺失值的观察值
2. 输入基于其他观察的缺失值
3. 插值和外推
4. 使用 KNN
5. 均值/中位数插补
6. 回归插补
7. 随机回归插补
8. 热卡插补

如果你想更详细地了解他们，请参考这篇文章。在我们的例子中，两个特征(标题邻近 tfidf 和描述邻近 tfidf)主要包含 0，因此我将用 0 替换缺少的值。对于城市匹配特征，1 和 0 分布几乎相等。在这里，我可以选择删除包含 null 的值，或者使用 mean 替换。我已经删除了这种情况下的值。

特征工程是关于从现有的输入特征中创建新的输入特征。

一般来说，你可以把数据清理看成是一个减法的过程，把特征工程看成是一个加法的过程。

这通常是数据科学家为提高模型性能所能做的最有价值的任务之一，原因有三:

1. 您可以隔离并突出显示关键信息，这有助于您的算法“专注于”重要的内容。
2. 可以自带 领域专长 。
3. 最重要的是，一旦你理解了特性工程的“词汇表”，你就可以引入其他人的领域专长！

以下是我们可以执行特征工程的一些方法，但是请注意，这并不是所有特征工程的详尽概要，因为这个步骤有无限的可能性。好消息是，随着你获得更多的经验，这项技能自然会提高。

1. 注入领域知识
2. 创建交互式功能
3. 组合稀疏类
4. 添加虚拟变量
5. 移除未使用的功能

在我们的例子中，由于没有太多关于数据集的领域知识，我们在特征工程的应用中受到限制。我应用的唯一特征工程是将两个相关的特征(title_proximity_tfid 和 main_query_tfidf)相乘，以创建一个名为 main title tfidf 的新列。

影响模型选择的一些因素有:

1. 模型是否满足业务目标
2. 模型需要多少预处理
3. 模型有多精确
4. 这个模型的解释力有多强
5. 模型有多快:建立模型需要多长时间，模型做预测需要多长时间。
6. 模型的可伸缩性如何

影响算法选择的一个重要标准是模型复杂度。一般来说，较复杂的模型是:

1. 它依赖于更多的特征来学习和预测(例如，使用两个特征对十个特征来预测目标)
2. 它依赖于更复杂的特征工程(例如，使用多项式、相互作用或主成分)
3. 它有更多的计算开销(例如，单个决策树与 100 棵树的随机森林相比)。

除此之外，基于参数的数量或一些超参数的选择，相同的机器学习算法可以变得更加复杂。举个例子，

1. 一个回归模型可以有更多的特征，或者多项式项和交互项。
2. 决策树可以有或多或少的深度。

使相同的算法更复杂会增加过度拟合的机会。

逻辑回归

逻辑回归模型符合一条“直线”。实际上，他们很少表现良好。对于大多数机器学习问题，我们实际上建议跳过它们。

它们的主要优点是易于解释和理解。但是，我们的目标不是研究数据，写研究报告。我们的目标是建立一个可以做出准确预测的模型。

在这方面，逻辑回归有两个主要缺陷:

1. 输入特征多了容易过拟合。
2. 它不能轻易表达非线性关系。

正规化

如上所述，逻辑回归遭受过拟合和处理非线性关系的困难。正则化是一种通过人为惩罚模型系数来防止过度拟合的技术。

• 它可以阻止大系数(通过抑制它们)。
• 它还可以完全移除要素(通过将它们的系数设置为 0)。
• 惩罚的“力度”是 可调 。(明天会有更多相关内容……)

正则化的类型有套索(L1)、脊(L2)和弹性网(脊和套索的折衷)

决策树

决策树将数据建模为分层分支的“树”。它们分支，直到到达代表预测的“叶子”。由于它们的分支结构，决策树可以很容易地模拟非线性关系。

不幸的是，决策树也有一个重大缺陷。如果你允许他们无限制地增长，他们可以完全“记住”训练数据，只是从创建越来越多的分支开始。因此，单个无约束决策树非常容易过度拟合。​

那么，我们如何利用决策树的灵活性，同时防止它们过度拟合训练数据呢？

树系综

集成是用于组合来自多个独立模型的预测的机器学习方法。有几种不同的组装方法，但最常用的两种是:

1. 装袋:试图减少复杂模型过度拟合的机会。

Bagging

• 它并行培养了大量的“强”学习者。
• 一个 强学习者 是一个相对无约束的模型。
• Bagging 然后将所有强学习者结合在一起，以便“平滑”他们的预测。
• 常用的技术是随机森林

2.助推:试图提高简单模型的预测灵活性。

Boosting

• 它按顺序训练出大量的“弱”学习者。
• 一个 弱学习器 是一个约束模型(即你可以限制每个决策树的最大深度)。
• 序列中的每一个都着重于从之前的错误中学习。
• 然后，Boosting 将所有弱学习者组合成一个强学习者。
• 常用的技术是 XGBoost 和 LightGBM

3. LightGBM: Light GBM 是一个使用基于树的学习算法的梯度推进框架。轻型 GBM 垂直生长树而其他算法水平生长树意味着轻型 GBM 逐叶生长而其他算法逐级生长。它会选择 delta 损失最大的叶子来生长。当生长相同的叶子时，逐叶算法可以比逐层算法减少更多的损失。

How LightGBM works

How other boosting algorithm works

还有许多其他算法，如支持向量机、神经网络等。但我们不会在这里接受。

对于我们的例子，我将使用 XGBoost、Random Forest 和 LightGBM。

培训模式:2018 年 1 月 21 日至 2018 年 1 月 26 日

测试型号:2018 年 1 月 27 日

我们将使用的指标是 AUC。我们得到的初始 AUC 值

1. XGBoost: 0.5803
2. LightGBM: 0.5807
3. 随机森林:0.5806

在几乎所有实际情况下，都需要搜索机器学习模型的参数以获得最佳交叉验证性能，从而获得具有最佳泛化估计的模型。scikit-learn 中的一个标准方法是使用 GridSearchCV 类，它为每个要尝试的参数取一组值，并简单地枚举参数值的所有组合。随着新参数的增加，这种搜索的复杂性呈指数增长。一种更可扩展的方法是使用 RandomizedSearchCV，但是它没有利用搜索空间的结构。

Scikit-optimize 为 GridSearchCV 提供了一个替代方案，它利用贝叶斯优化，其中一个称为“代理”的预测模型用于对搜索空间建模，并用于尽快获得良好的参数值组合。

贝叶斯优化，一种基于模型的寻找函数最小值的方法，最近被应用于机器学习超参数调整，结果表明这种方法可以在测试集上实现更好的性能，同时比随机搜索需要更少的迭代。

在应用贝叶斯优化和交叉验证后，AUC 值:

1. XGBoost: 0.5819
2. LightGBM: 0.5819
3. 随机森林:0.5810

尽管改进并不显著，但贝叶斯优化器能够以更快的速度执行调优操作。

堆叠是一种集成学习技术，通过元分类器或元回归器组合多个分类或回归模型。基于完整的训练集来训练基础级模型，然后在基础级模型的输出上训练元模型作为特征。

基础层通常由不同的学习算法组成，因此堆叠集成通常是异构的。在我们的例子中，我将组合 XGBoost 和 LGBoost 输出，并将使用投票分类器，这是 Scikit-learn 中提供的一个包。

如下图所示，叠加两个输出后，AUC 得分从 0.5819 提高到 0.5848。

ROC Graph for XGB, LGB, and Stacking

1. 如此低的 AUC 分数 0.5848 可能是由于我们在数据集中没有太多的特征，这使得算法很难正确地对目标变量进行分类。

1. 包含最后一列(class_id)以改善结果

如果你对代码感兴趣，你可以在这里找到我的笔记本。

[## 利用机器学习检测信用卡欺诈

用数据科学抓坏人

towardsdatascience.com](/detecting-credit-card-fraud-using-machine-learning-a3d83423d3b8) [## 第五章:机器学习中的算法选择-数据科学入门

在本指南中，我们将向您展示如何在众多选项中选择最有效的机器学习算法…

elitedatascience.com](https://elitedatascience.com/algorithm-selection) [## Python 中的自动机器学习超参数调整

Python 中使用贝叶斯优化进行自动超参数调优的完整演练

towardsdatascience.com](/automated-machine-learning-hyperparameter-tuning-in-python-dfda59b72f8a)

原文：https://towardsdatascience.com/predicting-customer-churn-using-pyspark-6a78a78a8412?source=collection_archive---------25-----------------------

Photo by Kaboompics from Pexels

客户流失是当今商业世界的一个主要商业问题，数据科学家正在快速采用工具和技术来有效预测客户流失并实时避免客户流失。

Apache Spark，尤其是 PySpark，是执行探索性分析和机器学习来解决这一分类问题的最快工具之一。

在本文中，主要目标是探索用户活动数据集，并使用 Spark 建立一个预测用户流失的机器学习模型。

**喜欢你读的书吗？**关注我上LinkedIn 或 中

点击获取数据服务产品

数据已经出来了。json 格式由别名为 Sparkify 的音乐公司的日志组成。这是用户与音乐平台交互时捕获的日志。

加载数据

path = “s3n://xxxxxxx-xxxx/sparkify/mini_sparkify_event_data.json”
df = spark.read.json(path)

数据架构:

df.printSchema()

在本节中，我们将执行一些数据预处理，以回答一些业务问题。

删除任何空行或缺少 userId 的行

#Drop null rows
df = df.dropna(how = “any”, subset = [“userId”, “sessionId”,”ts”])#Drop rows with missing userId
df = df.filter(df["userId"] != "")

删除任何重复的行

df.select(“userId”).dropDuplicates()

添加流失标签

churn_indicator = udf(lambda c: 1 if c == 'Cancellation Confirmation' else 0, IntegerType())
df = df.withColumn('churn_indication', churn_indicator('page'))#Add churn columns to indicate users who have churned
windowval = Window.partitionBy('userId')
df = df.withColumn('churn', max('churn_indication').over(windowval))

那么我们的客户流失数据框架将会是:

df_churn = df.filter('churn == 1')

下一步，我们将回答以下业务问题:

性别对流失有影响吗？

#Group by users and gender to aggregate count of usersdf_churn_by_gender = df.select(["userId", "gender","churn"]).groupby(["churn", "gender"]).count().sort("churn").toPandas()#Plot a barplotsb.barplot(x='churn', y='count', hue='gender', data=df_churn_by_gender)
plt.title('What is the churn comparison by gender', fontsize= 16);
plt.xlabel('Churn');
plt.ylabel('Number of Users');

从上图可以看出，男性比女性有更多的麻烦，因为女性是这项服务的主要用户。

用户流失多久？

为了回答这个问题，我们添加了一个以天为单位的总时间列，并计算每个用户总天数的最大值。

df = df.withColumn('total_time_days', (df.ts-df.registration)/1000/3600/24)total_time_df = df.select('UserId','churn','total_time_days').groupBy('userId','churn').agg(max('total_time_days').alias('total_time')).toPandas()sb.boxplot(data=total_time_df, x='churn', y='total_time', orient='v');
plt.ylabel('Total days');
plt.xlabel('Churned');
plt.title('After how long do Users churn?');

大多数用户在使用音乐平台的第二个月就会流失

哪一天和哪一小时的客户流失率高？

我们从添加小时和工作日列开始。为了对上述问题有所了解，我们将汇总每天每小时的用户数量，然后以工作日为中心，如下图所示。

#Create an hour and weekday column
calc_hour = udf(lambda t: dt.datetime.fromtimestamp(t / 1000.0).hour)
df_churn = df_churn.withColumn(“ts_hour”, calc_hour(df.ts))calc_weekday = udf(lambda t: dt.datetime.fromtimestamp(t / 1000.0).strftime(“%A”))
df_churn = df_churn.withColumn(“ts_weekday”, calc_weekday(df.ts))df_churn_by_time = df_churn.select(['userId', 'ts_weekday','ts_hour','churn']).groupby(["userId","ts_weekday","ts_hour"]).agg(count(df.userId).alias('count'))#.toPandas()df_churn_by_hr_week = df_churn_by_time.groupBy('ts_hour').pivot('ts_weekday').sum('count').sort('ts_hour')df_churn_by_hr_week = df_churn_by_hr_week.withColumn('ts_hour',df_churn_by_hr_week.ts_hour.cast('int'))
df_churn_by_hr_week = df_churn_by_hr_week.toPandas().set_index('ts_hour').sort_index(axis=0,ascending=True)#Plot a Heat Mapplt.figure(figsize=(16,10))
sb.heatmap(df_churn_by_hr_week, fmt='d',  cmap='viridis_r', annot_kws={"size": 12},  cbar_kws={'label': 'Number of Churns'})
plt.title("Which Hour and Day has High Customer Churn?", y=1, fontsize=18)
plt.xlabel('Day of the week', labelpad=8)
plt.ylabel('Hour (24hr)', labelpad=8)
plt.yticks(rotation=360);

让我们开始构建一些有希望的特性，我们将使用它们来训练我们的模型

1. 性别特征。男性或女性

gender_ft = df.select("userId", "gender").dropDuplicates().replace(['M', 'F'], ['0', '1'], 'gender').select('userId', col('gender').cast('int'))

2.订阅级别—付费或免费

level_ft = df.select("userId", "level").dropDuplicates().replace(['free', 'paid'], ['0', '1'], 'level').select('userId', col('level').cast('int'))

3.用户在平台上的总时间(生命周期)

total_time_ft = df.select(‘UserId’,’total_time_days’).groupBy(‘userId’).agg(max(‘total_time_days’).alias(‘total_time’))

4.用户听过的歌曲总数

total_songs_ft = df.select(‘UserId’,’song’).groupBy(‘userId’).agg(count(‘UserId’).alias(‘total_songs’))

5.平台上每个用户会话的歌曲数量

songs_per_session_ft = df.filter(df.page=="NextSong").groupBy('UserId','sessionId').count().groupBy('userId').agg(avg('count').alias('songs_per_session'))

6.艺术家总数

total_artists_ft = df.filter(df.page==”NextSong”).select(‘UserId’,’artist’).dropDuplicates().groupBy(‘userId’).agg(count(‘UserId’).alias(‘total_artists’))

7.总收听时间

total_listening_ft = df.select(‘UserId’,’length’).groupBy(‘userId’).sum().withColumnRenamed(‘sum(length)’, ‘total_listen_time’)

基于页面导航日志的更多功能

8.添加的好友数量

friends_added_ft = df.filter(df.page==”Add Friend”).groupBy(‘userId’).agg(count(‘UserId’).alias(‘friends_added’))

9.帮助页面访问

help_ft = df.filter(df.page==”Help”).groupBy(‘UserId’).agg(count(‘UserId’).alias(‘help_access’))

10.添加到播放列表的歌曲数量

songs_playlist_ft = df.filter(df.page==”Add to Playlist”).groupBy(‘UserId’).agg(count(‘UserId’).alias(‘songs_on_playlist’))

11.赞数—竖起大拇指页面

likes_ft = df.filter(df.page==”Thumbs Up”).groupBy(‘UserId’).agg(count(‘UserId’).alias(‘likes’))

12.不喜欢的数量—拇指向下页面

dislikes_ft = df.filter(df.page==”Thumbs Down”).groupBy(‘UserId’).agg(count(‘UserId’).alias(‘dislikes’))

最后，

13.流失标签

churn_label_ft = df.select(‘UserId’, col(‘churn’).alias(‘label’)).dropDuplicates()

最后，我们将使用外部连接将我们的特性连接到一个数据帧中，清除一些空值并删除不必要的 UserId 列

#combine all datasources into a single data frame
feature_df = gender_ft.join(level_ft,’userID’,’outer’) 
 .join(total_time_ft,’userID’,’outer’) 
 .join(total_songs_ft,’userID’,’outer’) 
 .join(total_artists_ft,’userID’,’outer’) 
 .join(songs_per_session_ft,’userID’,’outer’) 
 .join(total_listening_ft,’userID’,’outer’) 
 .join(friends_added_ft,’userID’,’outer’) 
 .join(help_ft,’userID’,’outer’) 
 .join(songs_playlist_ft,’userID’,’outer’) 
 .join(likes_ft,’userID’,’outer’) 
 .join(dislikes_ft,’userID’,’outer’) 
 .join(churn_label_ft,’userID’,’outer’)#Drop unnecessary userid column and fill null valuesfeature_df = feature_df.drop(‘userID’).fillna(0)feature_df.show(5)

我们将首先使用 VestorAssembler(一种将多个列合并为一个向量列的要素转换器)和 StandardScaler(通过移除平均值并缩放至单位方差来标准化要素)对要素数据集进行矢量化和标准化

#Vectorize features
cols = [‘gender’, ‘level’, ‘total_time’, ‘total_songs’, ‘total_artists’, ‘songs_per_session’, ‘total_listen_time’, ‘friends_added’, ‘help_access’, ‘songs_on_playlist’, ‘likes’, ‘dislikes’]
assembler = VectorAssembler(inputCols=cols, outputCol=”vfeatures”)
data = assembler.transform(feature_df)# standardize features
scaler = StandardScaler(inputCol="vfeatures", outputCol="features", withStd=True)
scalerModel = scaler.fit(data)
data = scalerModel.transform(data)

将数据集分为训练、测试、验证和测试，如下所示:

# train test split
train, rest = data.randomSplit([0.6, 0.4], seed=42)
validation, test = rest.randomSplit([0.5, 0.5], seed=42)

下一步包括比较模型的性能并选择性能最佳的模型。我们将使用默认参数训练模型，并对最佳模型进行改进。最佳模型将用于扩展到更大的数据集。

1. 逻辑回归

#estimator
lr = LogisticRegression(maxIter=10)
#evaluator
mce_f1_evaluator = MulticlassClassificationEvaluator(metricName=’f1')# build an empty paramGrid for now
paramGrid = ParamGridBuilder().build()lr_cv = CrossValidator(estimator=lr,
 evaluator=mce_f1_evaluator, 
 estimatorParamMaps=paramGrid,
 numFolds=3)#Train the model
lr_cv_model = lr_cv.fit(train)
#fit the model
svm_results = svm_model.transform(validation)#Evaluate using f1 score
evaluator = MulticlassClassificationEvaluator(predictionCol = "prediction")
print('F1 Score:{}'.format(evaluator.evaluate(svm_results, {evaluator.metricName: "f1"})))

2.支持向量机

#estimator
svm = LinearSVC(maxIter=10)
#evluator
mce_f1_evaluator = MulticlassClassificationEvaluator(metricName=’f1')
# build an empty paramGrid for now
paramGrid = ParamGridBuilder().build()svm_cv = CrossValidator(estimator=svm,
 estimatorParamMaps=paramGrid,
 evaluator=mce_f1_evaluator,
 numFolds=3)#train
svm_model = svm_cv.fit(train)
#fit
svm_results = svm_model.transform(validation)
#evaluate
evaluator = MulticlassClassificationEvaluator(predictionCol = "prediction")
print('F1 Score:{}'.format(evaluator.evaluate(svm_results, {evaluator.metricName: "f1"})))

3.随机森林分类器

#classifier
rf = RandomForestClassifier()#evaluator
mce_f1_evaluator = MulticlassClassificationEvaluator(metricName=’f1')# build an empty paramGrid for now
paramGrid = ParamGridBuilder().build()rf_cv = CrossValidator(estimator=rf,
 estimatorParamMaps=paramGrid,
 evaluator=mce_f1_evaluator,
 numFolds=3)#train
rf_model = rf_cv.fit(train)#validate
rf_results = rf_model.transform(validation)#evaluate
evaluator = MulticlassClassificationEvaluator(predictionCol = "prediction")
print('F1 Score:{}'.format(evaluator.evaluate(rf_results, {evaluator.metricName: "f1"})))

4.梯度推进树

#classifier
gbt = GBTClassifier(maxIter=10,seed=42)#evaluator
mce_f1_evaluator = MulticlassClassificationEvaluator(metricName=’f1')# build an empty paramGrid for now
paramGrid = ParamGridBuilder().build()gbt_cv = CrossValidator(estimator=gbt,
 estimatorParamMaps=paramGrid,
 evaluator=mce_f1_evaluator,
 numFolds=3)#train
gbt_model = gbt_cv.fit(train)#validate
gbt_results = gb.model.transform(validation)#evaluate
evaluator = MulticlassClassificationEvaluator(predictionCol = "prediction")
print('F1 Score:{}'.format(evaluator.evaluate(gbt_results, {evaluator.metricName: "f1"})))

4 个分类器的结果如下:

Classifier Results

梯度提升树是迄今为止最好的评分模型，F1 得分为 0.9，准确度得分为 0.9，其次是随机森林分类器。我们还可以观察其他指标，比如模型训练所需的时间。可以说，我们希望提供最好的结果，因此基于准确性和 f1 分数，我们可以选择我们在这种情况下的最佳模型。

在这最后一部分，我们希望通过将超参数添加到梯度推进树和随机森林分类器模型来调整我们的 2 个最佳模型

渐变提升树

以下是默认参数

用参数建立模型

gbt = GBTClassifier(maxIter=10,seed=42)# build param Grid
paramGrid_gbt = ParamGridBuilder() 
    .addGrid(gbt.maxIter,[5, 10,15]) 
    .addGrid(gbt.maxDepth,[2,3,4,5,6,7,8]) 
    .build()
# set evaluator
mce_f1_evaluator = MulticlassClassificationEvaluator(metricName='f1')gbt_hpt_cv = CrossValidator(estimator=gbt,
                          estimatorParamMaps=paramGrid_gbt,
                          evaluator=mce_f1_evaluator,
                          numFolds=3)gbt_hpt_model = gbt_hpt_cv.fit(train)#Test the model using test data
gbt_model_results = gbt_hpt_model.transform(test)

随机森林分类器

#Default parameters

用参数建立模型

#classifier
rf = RandomForestClassifier()#evaluator
mce_f1_evaluator = MulticlassClassificationEvaluator(metricName=’f1')# build an empty paramGrid for now
paramGrid = ParamGridBuilder() 
 .addGrid(rf.impurity,[‘entropy’, ‘gini’]) 
 .addGrid(rf.maxDepth,[2,3,4,5,6,7,8]) 
 .build()rf_cv = CrossValidator(estimator=rf,
 estimatorParamMaps=paramGrid,
 evaluator=mce_f1_evaluator,
 numFolds=3)
#train model
rf_model = rf_cv.fit(train)
#Test using test data set
rf_model_results = rf_model.transform(test)

在测试数据上测试两个调整的模型之后，随机森林分类器具有体面的 0.85 准确度分数和 0.84 f1 分数。

梯度提升以 0.78 的准确度分数和 0.78 的 f1 分数位居第二

与梯度相比，随机森林分类器在参数调整后的测试数据集上得分最高。当把时间作为一个衡量标准时，训练也需要更少的时间。因此，在 4 种算法中，随机森林是最健壮的

我首先加载数据，删除丢失的值，删除重复的行，为数据探索创建额外的列，并使用可视化工具回答了一些业务问题。然后开始研究特征生成和建模。该项目的最有趣和最困难的部分是 spark ML 的实现，使用 PySpark 预处理和特征工程数据；考虑到这是第一次使用这个工具。能够在云上部署 ML 解决方案(在本例中是 AWS EMR 服务)是一个很好的技能补充。

另一个具有挑战性的部分是特征生成。虽然我只能想出 12 个功能，但他们有可能创造更多。因此，有了更大的数据集和更多的特征，该模型在预测客户流失方面会更加准确和稳健

该模型也可以作为 API 部署，在业务环境中使用。

更多细节请查看我的 Github 库

原文：https://towardsdatascience.com/predicting-customer-churn-with-neural-networks-in-keras-f904c7113fcd?source=collection_archive---------12-----------------------

这对于任何地方的组织都是一个大问题，也是我们看到机器学习高采用率的主要领域之一，这可能是因为我们正在预测客户行为。

“流失”一词用于描述客户停止使用某个组织的服务。对于订阅服务来说，这是一个非常强大的 KPI，订阅服务的大部分收入来自重复支付(通常是每月支付)。

网飞是订阅公司的一个很好的例子，他们也是一个接近技术前沿的组织。网飞使用深度学习技术来预测客户是否会在实际离开之前离开，这意味着他们可以采取预防措施来确保他们留下来。

不用深入兔子洞，开始吧…

网飞收集了很多关于个人的数据，你看了什么，你什么时候看的，你喜欢和不喜欢的一切等等。他们可以将这些数据与深度学习分类技术结合使用，计算出他们认为客户何时会离开。简单的解释可能是这样的“如果一个客户 N 天没看任何东西，那么他们很快就会流失”。

使用神经网络分析所有收集到的数据将使组织能够根据他们现有的数据建立他们客户的档案。一旦一群用户被分类，网飞可以决定采取什么行动，如果一个客户是一个可疑的流失。例如，他们可以了解他们想要留住哪些客户，并为这些人提供折扣和促销，然后还可以确定哪些客户是“注定失败的”，可以让他们离开。

IBM 电信数据集已经在互联网上流传了一年多，所以我想现在是尝试使用它来预测客户流失的好时机。

我需要从别处借用一些代码。我直接从这篇博文中摘录了数据准备代码。这让我可以花更多的时间调整模型，并(试图)获得更高的精确度。

这是非常简单的，它几乎是一个直接的复制/粘贴到您自己的笔记本上的工作。在网上搜索时，我发现一大堆人在分发导入了错误库的代码。这里要注意的主要事情是我导入 Keras 模块的第二个块。

很多帖子(我假设都来自同一个帖子)混淆了导入 Keras 库和 TensorFlow 库。一个好的经验法则是确保你是:

• 只进口 Keras 的东西:他们会从“Keras”开始。例如，从 keras.models 导入…
• 仅导入 TensorFlow Keras 库，这些库以 TensorFlow.Keras 开头，例如来自 TensorFlow。Keras.models 导入…

import numpy as np 
import pandas as pd 
from matplotlib import pyplot as plt 
from sklearn.model_selection import train_test_split 
import keras 
from keras.models import Sequential 
from keras.layers import InputLayer 
from keras.layers import Dense 
from keras.layers import Dropout 
from keras.constraints import maxnorm

读取 pandas 数据框架中的测试数据(通过单击顶部的 IBM 数据集链接来获取 CSV)。

data = pd.read_csv('churn.csv')

一旦我们在 Pandas DF 中有了数据，我们就可以使用窃取的预处理代码将数据转换成针对神经网络优化的格式。本质上，我们正在做以下事情:

• 将分类变量转换为数值变量(例如，是/否转换为 0/1)。
• 正确设置列格式，使它们都是数字格式
• 填充任何空值。

data.SeniorCitizen.replace([0, 1], ["No", "Yes"], inplace= True) data.TotalCharges.replace([" "], ["0"], inplace= True) data.TotalCharges = data.TotalCharges.astype(float) data.drop("customerID", axis= 1, inplace= True) data.Churn.replace(["Yes", "No"], [1, 0], inplace= True)

一旦我们有了一个干净的数据集，我们就可以使用 pandas get_dummies 功能将所有分类列替换为‘dummy’数字‘indicator’列。类似于上面的代码所做的，但是这将去掉整个数据集。

data = pd.get_dummies(data)

像任何模型一样，我们应该将数据分成训练和验证(或测试集)。

首先我们把数据集分成 X 和 y。

• x 包含了我们用来做预测的所有变量。
• y 只包含结果(无论客户是否翻炒)。

X = data.drop("Churn", axis= 1) y = data.Churn

接下来，我们使用标准的 train_test_split 将数据分成训练和测试(验证)集。

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size= 0.2, random_state= 1234)

现在我们已经准备好数据集，让我们建立一个模型。

Keras 的一个伟大之处在于，我们可以非常简单地建立一个基于层的神经网络。看起来有点像这个图。

首先，我们告诉 Keras 我们想要使用什么类型的模型。大多数情况下，你会使用序列模型。

model = Sequential()

接下来让我们首先构建输入层。图中的红色层。这里有几点需要注意。

我们在整个神经网络中使用密集层，我不会深入不同层和神经网络之间的差异，但大多数时候你会使用密集层或 LSTM 层。点击此处了解更多信息。

第一个数字— 16 是节点(或上图中的圆圈)的数量。我们从 16 开始，以后我们可以改变它，看看它如何影响我们的精度。

设置 input_dim 参数很重要，它必须与 X_train 数据集中的列数相匹配。我们可以简单地计算数据集中的列数，然后像我下面所做的那样键入它——但实际上，您会希望使用 X_train.shape[1]来自动计算值。

激活函数——你真的应该在某个时候仔细阅读这些——我可能会专门就这些写另一篇文章，但目前只知道我们的层都使用“Relu ”,除了输出层(我们将在后面讨论)。作为一般的经验法则*——“如果你不确定，使用 relu”*

model.add(Dense(16, input_dim=46, activation='relu', kernel_constraint=maxnorm(3)))

我们增加了一个辍学层。dropout 层确保我们在每次迭代神经网络时删除设定百分比(在本例中为 0.2%或 20%)的数据。这是可选的，但值得包含在您的代码中，以防止它过度适应。你可以改变不同的辍学率来进行实验，以获得更高的精确度。

model.add(Dropout(rate=0.2))

现在，我们添加了所谓的隐藏层，我们可以拥有尽可能多的隐藏层，但值得考虑的是，使具有大量隐藏层的神经网络工作所需的额外计算能力。

注意到“内核约束”参数了吗？这涉及使“下降”层有效工作所需的权重的缩放。文档中还有更多的信息，但我想你需要知道的是，从本质上讲，神经网络会自动试错数据集中变量的所有不同权重(这就是为什么我们首先使用神经网络)，kernel_constraint 参数增加了对这一过程的控制。

model.add(Dense(8, activation='relu', kernel_constraint=maxnorm(3)))

添加另一个脱落层，以避免过度拟合。

model.add(Dropout(rate=0.2))

最后，我们添加一个输出层:这定义了神经网络的最终输出。这里有几件事你需要记住:

• 我们的第一个参数是数字 1。这是因为我们的 NN 的 out 是一个包含指示符的列，该指示符将指定我们的客户是否会流失。
• 激活功能不同。对于单个(是/否)分类模型，我们使用“sigmoid ”,有许多不同的功能可用，例如，如果我们正在构建一个对多个结果进行分类的网络(例如将客户分组到特定的组),我们可以使用“softmax”激活功能。

model.add(Dense(1, activation='sigmoid'))

现在我已经描述了每一个单独的层，我将展示挤在一起的整体。

model = Sequential() 
model.add(Dense(16, input_dim=46, activation='relu', kernel_constraint=maxnorm(3)))
model.add(Dropout(rate=0.2)) 
model.add(Dense(8, activation='relu', kernel_constraint=maxnorm(3))) model.add(Dropout(rate=0.2)) model.add(Dense(1, activation='sigmoid'))

接下来，我们编译我们的模型(把它们粘在一起，告诉它应该如何工作)。

我们使用编译方法来实现这一点。它需要三个参数:

• 优化器。这可能是现有优化器的字符串标识符(我们使用“Adam”)，可以使用这些标识符，看看哪个是最好的，您也可以将其拆分出来，并手动调整学习速率以受益于更高的准确性，但我们将使用默认值。参见:优化师。
• 损失函数。这是模型试图最小化的目标。由于这是一个单一的分类问题，我们将使用“二元交叉熵”。参见:损失。
• 指标列表。对于任何分类问题，我们都将其设置为“准确性”。

model.compile(loss = "binary_crossentropy", optimizer = 'adam', metrics=['accuracy'])

快速拟合模型，请注意，我们在这里将测试和验证数据集都拟合到了我们的模型中，这样 Keras 将立即告诉我们在这两个数据集上的表现。

这里我们需要注意几个论点。

• Epochs =这是我们在神经网络中来回移动的次数——更多的 Epochs 意味着更高的准确性，但也意味着更多的处理时间，而且太多的 epochs 会导致过度拟合。您可以随时对此进行调整，所以尝试一些不同的值，看看结果会发生什么变化。
• batch_size =这是一次通过神经网络的记录批数。批量越小，精确度越低，但速度越快。同样，您可以通过在代码中直接更改这个数字来进行试验。

history = model.fit(X_train, y_train, validation_data=(X_test, y_test), epochs=40, batch_size=10)

一旦我们运行上面的代码，我们将得到模型执行情况的指示——您将看到它在各个时期运行，并为训练(acc)和测试(val_acc)集提供准确性分数。

一旦它完成了它的任务，剩下要做的就是观想结果。

plt.plot(history.history['acc']) plt.plot(history.history['val_acc']) 
plt.title('model accuracy') plt.ylabel('accuracy') plt.xlabel('epoch') plt.legend(['train', 'test'], loc='upper left') plt.show()

我们需要做的最后一件事是保存模型，这样我们就可以在以后将它部署到生产环境中——AWS SageMaker/Azure/Google 都有不同的方法来做到这一点，但是通常你需要一个 JSON 和一个权重文件。

# serialize model to JSON 
model_json = model.to_json() with open("model.json", "w") as json_file: json_file.write(model_json) # serialize weights to HDF5 model.save_weights("model.h5") print("Saved model to disk")

所以我们有它。一个神经网络产生了大约 79%的客户流失准确率。当然，我们可以花更多的时间研究学习率、激活函数、节点数量、时期数量等，以使其更加准确，但希望这是开始研究的坚实基础。

接下来，我将研究如何将这些代码部署到生产环境中——这变得很棘手。

原载于 2019 年 5 月 27 日【http://drunkendatascience.com。

原文：https://towardsdatascience.com/predicting-customer-churn-with-pyspark-95cd352d393?source=collection_archive---------27-----------------------

你怎么称呼一群分散的恐龙繁殖岛？

侏罗纪火花…Ba Dum Tss

O 对于采用订阅式商业模式的公司来说，最重要的问题之一是客户流失。客户因为各种原因降级或停止服务，服务提供商通常直到客户离开才知道他们何时或为什么离开！

如果公司能够在客户离开之前预测到他们何时会流失，会怎么样？

这将是强大的！如果我们能够可靠地预测客户是否会流失，我们就有机会通过促销、宣传新功能等方式留住这些客户。这是一种主动的方法来留住客户，而不是一种被动的方法来挽回失去的客户。

Udacity 提供了两个类似格式的用户活动数据集(128MB 和 12GB)，来自一家虚构的音乐流媒体公司 Sparkify，我使用这些数据来更好地了解 Sparkify 的客户，然后预测客户是否会流失，准确率超过 80%。

我无法在本地对我的机器上的 12GB 数据集进行分析和建模，所以我在 Jupyter 笔记本中使用 PySpark 的本地模式 在 128MB 数据集上探索并原型化了我的工作流，这让您 模拟在一台机器上的分布式集群上工作的 。

然后，我把工作流变成了一个 Python 脚本 ，旋出了一个 4 节点 AWS EC2 集群 ，在集群上通过 AWS EMR 执行脚本。

以下是我在这篇文章中将要涉及的内容:

1. 特征工程&在本地探索数据
2. 本地数据预处理
3. 本地机器学习
4. 在 AWS 分布式集群上运行步骤 1–3

这篇文章是写给谁的？

本文旨在向您介绍我所采取的步骤和一些有用的代码，这些步骤和代码带我从一个小型 Jupyter 笔记本分析到一个 4 节点 AWS EMR 集群上的 12GB 数据集分析。所以，我假设你对Python&机器学习 。

此外，在这篇文章中，我展示了许多我使用**【py spark】的代码片段，这些代码片段与通过 Pandas 和 Scikit-Learn 库处理数据的非常不同。如果你不熟悉 PySpark，那就略读一下吧！没什么大不了的。

最后，如果你想对分析和机器学习如何通过我的代码进行更详细的分析和演练，请查看我的完整 Jupyter 笔记本(托管在我的网站)或 GitHub Repo 。

我希望你喜欢这个！

现在让我们开始吧…

首先，我将 128MB 的数据集从 JSON 格式转换成 Spark 数据帧。以下是模式，以便您了解数据的结构:

Original imported features.

日期特性工程| 经过一些简单的探索，我使用“ts”(时间戳)属性通过 PySpark UDF(用户定义函数)创建了一些日期特性。

我为新的日期特性创建了一个计数图，以便更好地理解用户在整个数据集期间的行为。

Note this smaller dataset (128MB) only contains data for about two months.

事件特性工程| 我使用了“页面”特性来标记用户访问的特定页面或执行的特定操作。这些标志特性将在以后帮助我在用户级别而不是事件级别聚合用户活动。

This is identical to one-hot encoding based on this feature. In hindsight, this would have been faster with PySpark’s native OneHotEncoder. Oh well! It works.

运行计数特征工程| 我还利用了数据的事务性质(例如，某个用户在某个时间点执行的每个事件对应一行)来创建一些窗口计算特征(例如，一个月内的运行收听计数等)。)这里有一个例子，说明这样的功能是如何创建的。

Creating a window for each user per month, counting rows whenever a user listens, and joining back to DF.

这将创建一个类似于图中第四列的列:

Last Column: Note the null values

由于这是一个连续计数，空值应该在前面填充。然而，鉴于 Spark 的弹性分布式数据集(rdd)的分布式本质，没有前置填充空值的原生函数。为了解释，如果一个数据集在云中的几个机器之间随机分区，那么机器如何知道前一个值何时是数据集中的真直接前一个值？

我最好的解决方法是创建一个 running listens 列的 lag，然后用这个 lag 迭代地替换空值。连续的空值需要多次执行这种练习，因为我们实际上是用先前的值替换当前的空值(对于连续的空值，先前的值为空)。我是这样实现的:

Sharing this snippet in case anyone finds it useful. This is the only way I found to front-fill null values that was computationally cheaper than using a complicated cross-join.

上面的代码片段用其先前的值填充空值 6 次，然后用 0 填充剩余的值以节省计算时间。这可能不是前置填充的最佳方式，但如果要用相对较少的空值填充大型数据集，这可能会很有用。这是可行的，因为该操作需要 Spark 连接分区来计算 lag(即紧接在前面的值)。

定义流失| 最后，我将客户流失定义为每当用户访问“取消确认”或“降级”页面时。

This labels the “Churn” column of the DataFrame as 1 if a user visits the aforementioned pages.

现在让我们退一步，想想我们拥有什么，我们想要实现什么。

我们所拥有的: 一个“事件”数据集，指定用户在给定的时间点执行什么活动。

我们想要的: 预测客户流失。

凭直觉，我不认为在我们现有的数据集中预测变动是一个好主意。基于用户事件进行训练和预测将是非常昂贵的，这可能是非常非常昂贵的*因为你每小时可以有数千个事件(或者更多！)*****

因此，我的方法是简单地汇总每个用户的数据。具体来说:

我使用了 、事件级数据 到 来创建 一个基于 指标的用户级矩阵 ，它本质上是在一个矩阵中总结每个用户的活动，该矩阵的行数与唯一用户的行数一样多。

这可能无法从文本中很好地翻译出来，所以让我给你看看我是怎么做的。下面是一个代码聚合示例，用于查看每个用户执行的全部活动:

这种类型的聚合帮助我创建一个整洁的数据框架和可视化，如下所示…

Truncated DF | One row per user with user-specific summary statistics.

Using the summary statistics to compare our churned users vs. our non-churned users.

有趣的是，看起来我们被激怒的用户比我们没有被激怒的用户更活跃。

现在我有了一个很好的、紧密的矩阵中的数据，我用它作为机器学习模型的基础。

我的建模方法简单明了:尝试一些算法，选择一个看起来最有希望的，然后调整这个模型的超参数。

这里有一个用户矩阵的模式，这样您就知道正在建模什么了:

One user and several metrics/aggregations for this user per row.

PySpark ML 要求数据采用非常特殊的数据帧格式。它需要它的 特征 在一个 列的向量 中，其中向量的每个元素代表它的每个特征的值。它还要求其 标签 在其 自有列 中。

下面的代码将原始矩阵转换成这种 ML 友好的格式，并标准化这些值，使它们具有相同的比例。

现在我们已经准备好了数据，下面的代码显示了我最初是如何评估模型的。我任意选择了逻辑回归、随机森林和梯度提升树作为最终流失模型的候选。

在初始模型评估后，我发现 gbt 分类器(梯度增强树)表现最好，准确率为 79.1%，F-1 得分为 0.799。

因此，我使用 PySpark 的原生 CrossValidator 和 ParamGridBuilder 将该模型调优为网格搜索，网格搜索使用 K-Fold 验证选择最佳超参数。这里，由于计算时间昂贵，我使用了超参数的小网格和 3 重验证。

经过调优后，重新评估模型，准确率和 F-1 分别提高到 82.35%和 0.831！最佳模型的最大深度为 3，最大箱数为 16。如果在我的笔记本电脑上运行不需要这么长时间(大约需要 1 个小时)，我会尝试更广泛的网格搜索。

PySpark 自动计算这些基于树的模型的特征重要性。

注意:出于好奇， 下面是 一篇有趣的文章，解释了特征重要性是如何计算的，以及 为什么它实际上不是那么准确 。这真的不在这个项目的范围内，所以我只是敷衍一下。

Feature Importances calculated for the tuned GBTClassifier.

看来前三个 预测流失最重要的特征 是:

• 每次会话的平均歌曲播放次数
• 总赞数
• 最连续几天不播放歌曲

看起来，与错误和歌曲播放总量相关的指标与预测客户流失完全无关。

到目前为止，我一直在描述我在笔记本电脑上使用 Spark 的本地模式对小型 128MB 数据集执行的分析。

为了对 12GB 的数据集进行同样的分析，我需要更多的处理能力。这就是 AWS 弹性地图简化(EMR)和弹性云计算(EC2)的用武之地。

注意:EC2 可以让你在云中的机器上工作。EMR 用于轻松管理已经安装了 Spark 的 EC2 集群。

在设置好 EMR 并通过我的终端访问 AWS 之后，我将我的脚本上传到我的集群的 Hadoop 分布式文件系统(HDFS)，从 Udacity 的 S3 存储桶加载包含 12GB 数据集的数据，并从我的终端在主节点上运行该脚本。

注意: 这里是我用来设置 EMR 的 和 这里是我用来设置我在终端上访问 AWS 的 。

Don’t even try to read this. Just FYI: this is what it looks like when you run a Spark app through EMR from your CLI/terminal.

这个脚本需要很长时间才能完成(我在 4 台机器上花了大约 5 个小时)。为了不浪费所有这些辛苦的工作，我在脚本的最后保存了模型和预处理过的用户矩阵。

如果您感到好奇，在完整的 12GB 数据集上运行此分析和工作流会产生非常高的准确性和 F-1 分数！

This was a HUGE improvement from the model trained on 128MB dataset.

更多的功能可以从用户活动中设计出来，比如每天/每周/每月的竖起大拇指数，竖起大拇指与不竖起大拇指的比率，等等。特征工程可以比简单地优化一种算法更好地改善结果。

因此，可以做进一步的工作，从我们的交易用户数据中提取更多的特征来改进我们的预测！

一旦一个模型被创建，也许它可以被部署到生产中，并且每隔 x 天或 x 小时运行一次。一旦我们预测到用户可能会流失，我们就有机会进行干预！

为了评估这个假设部署的模型做得有多好，我们可以运行一些概念验证分析，并且在给定的测试期间不干预它的预测。如果它预测的用户将以比普通用户更高的速度流失，这可以表明我们的模型工作正常！

学习 PySpark 和 AWS 好像是一场噩梦(我从这次经历中知道)。然而，如果你已经熟悉 Python 和机器学习，你不需要知道太多的来开始使用 PySpark 和 AWS。

你知道如何处理数据。你知道机器学习是如何工作的。战斗的另一半是知道如何在 PySpark 和 AWS 集群上执行这些相同的任务，这可以通过 Google 搜索和教程找到！我建议慢慢浏览我的(或其他人的) Jupyter 笔记本以便你对 PySpark 有所了解，然后查找如何在 AWS EMR 上运行 Spark 应用程序。

这看起来工作量很大，但是我根本不可能用笔记本电脑处理 12GB 的数据集。每天都有令人难以置信的大量数据被创建，如果你想处理这么多数据，Spark 是一个很好的学习工具！

原文：https://towardsdatascience.com/predicting-customer-churn-with-spark-4d093907b2dc?source=collection_archive---------16-----------------------

对于许多公司来说，客户流失是一个主要问题。一些人停止使用这项服务是很自然的，但如果这个比例变得太大，就会阻碍增长，不管收入来源如何(广告销售、订阅或两者兼而有之)。考虑到这一点，企业通过识别处于风险中的客户来预测客户流失的能力至关重要，因为这使他们能够采取某些行动，如个性化的优惠或折扣，以尝试和减少客户的流失。基于历史数据建立的机器学习模型可以让我们洞察客户流失的信号，并帮助我们在它发生之前预测它。

对于这个例子，我们使用一个虚构的音乐应用公司 Sparkify 的日志数据。这是一个玩具数据集，相对较小，因此可以由一台计算机处理。尽管如此，为了模拟真实世界，我们使用 Spark(在本地模式下)来处理数据和构建模型。Spark 是大数据处理和建模的领先解决方案之一，通过 DAG 的内存处理和计算延迟评估来提高速度。你可以点击这里了解更多。通过在这个项目中使用 Spark，尽管这不是绝对必要的，但我们建立了一个可扩展的框架来分析任何大小的数据的变动，因为代码和分析可以很容易地扩展，只要它们部署在能够处理所需计算的集群(如 AWS、IBM Cloud 或 GCP)上。

在这篇博文中，我将总结包含在这个 GitHub 库中的分析。在简要概述手头的数据后，我们将介绍该模型、其结果以及它对 Sparkify 的意义。

数据

在本地模式下启动 Spark 会话后，我们可以加载数据集。它包含 2018 年 10 月 1 日至 2018 年 12 月 3 日之间 226 个不同用户的信息。它是 JSON 格式的(关于 JSON 格式的更多信息在这里)，可以很容易地用以下命令加载:

path = “mini_sparkify_event_data.json”df = spark.read.json(path)

这些数据捕捉了各种各样的行为，如听歌、竖起大拇指、点击主页、更改帐户设置或向播放列表添加歌曲。因此，尽管它是一个很小的用户子集，但数据集仍然包括 278，251 行。

在数据集中出现的 226 名用户中，有 52 名用户最终出现了混乱。为了正确地训练和评估我们的模型，我们解决了这种不均衡，以确保我们的预测可以准确地预测这两个类别，而不是过度倾向于预测不存在流失。我们通过一种称为上采样的技术来做到这一点，即从搅动的用户群体中进行替换采样，直到我们得到两个大小相当的群体。

型号

为了构建这个模型的功能，我让数据探索决定了我的方法，以及我在一家非常相似的(真实的)音乐应用公司工作时获得的领域知识。特别是，我在寻找那些在翻炒用户和不翻炒用户之间价值差异很大的特性。寻找这样的功能突出了帐户类型(免费与付费)以及其他帐户相关信息的重要性，如状态和注册日期，让我们对用户有所了解。

Users who churned are more likely to have a free account

另一组特性集中在人们在平台上的行为。这些元素，如会话长度、每次会话的歌曲数量以及拇指向上/拇指向下、添加到播放列表或添加朋友，为我们提供了额外的洞察力。建模之前的直观解释是，这些特征捕捉到了与用户参与度相关的潜在变量，较低的参与度与较高的搅动可能性相关联。例如，停止搅动的用户每月来平台 9 次，而留下来的用户每月来 14 次。我们将在后面看到这种前概念是否被数据所证实。

The distribution in terms of session length and number of items in session differs between both groups

那些对创建这些特性的技术细节感兴趣的人可以参考介绍中链接的代码，但是由于 Spark 的管道，我们可以高效地处理所有这些特性，并将其转换为适合分析的形式。您可以在下面看到数据集的前几行及其所有功能。

First few rows of features dataset

检查这些特性的分布，我们可以看到 itemInSession、thumbsUp、addFriend 和 addToPlaylist 非常分散地分布在平均值周围。这一点很重要，因为模型依赖于可变性来学习和做出预测。另一方面，每日会话或长度等特征的变化较小，因此我预计这些特征在预测中的权重较小。

完成后，我们测试三种不同的分类模型(随机森林、逻辑回归和梯度推进),并评估它们在测试集上的准确性和 F1 得分。考虑这两者是很重要的(不仅仅是准确性)，因为后一个度量允许我们调整测试集中存在的类不平衡，并且通过扩展，在真实世界中也是如此。由于三个模型之间的结果没有显著差异，我们选择进一步调整逻辑回归模型，因为它具有更好的可解释性。我们通过交叉验证，利用网格搜索算法找到最佳的参数组合。特别是，我们测试了以下值:

• minInfoGain(在树节点处被考虑的分割的最小信息增益):0，1
• maxDepth(树的最大深度):5，10
• numTrees(树的数量):20，50

我特别选择了这些参数，因为它们与防止过度拟合有关。

结果

优化后，最优超参数是 50 棵树，0 最小信息增益和最大深度 10。我们有一个模型，在测试集上达到 73%的准确率，F1 值为 0.72。这两个指标合在一起非常令人鼓舞:只有关于 191 个用户(在我们的训练集中)的数据，我们就能够有效地将用户分类到这两个类别中，而在预测一个或另一个方面的性能没有显著差异。有趣的是，网格搜索后我们模型的性能并没有提高，很可能是因为我们的数据集很小。我们甚至通过训练和预测不同的随机状态来评估模型的稳健性，并发现模型的准确性在它们之间非常一致。

看看特征重要性，我们之前的直觉得到了证实:静态变量(注册月份、地理位置)和行为(添加好友)在我们的预测中都有很大的分量。这应该鼓励 Sparkify 记录尽可能多的信息，因为在试图预测客户流失时，所有信号都很重要。

结论

Spark 为我们提供了一个预测客户流失的通用框架。它可以为任何公司处理大数据，只要它部署在能够处理所需计算的集群上。如果这种分析应用于具有更多可用计算能力的更大数据集，我认为甚至会达到更好的精确度/F1，因为我们将能够在更大的超参数空间中为更多用户进行搜索。我们甚至可以在一个非常大的超参数空间上组合随机搜索，以产生一个子集，网格搜索将在这个子集上寻找最佳组合，以便进一步加快计算速度和提高性能。最后，为了更深入地了解模型，我们可以利用 SHAP 值或排列重要性来了解各个特征如何影响模型预测。

根据一小部分客户的历史数据，我们建立了一个模型，能够以 73%的准确率识别出有流失风险的用户。它可以定期(每天/每周，取决于现有的计算基础设施)应用于用户群，并标记可能很快离开服务的用户。记住这些信息，Sparkify 可以采取缓解措施，例如发送个性化信息或提供每月折扣。所有这些都可以自动化，并将对收入和增长产生巨大影响。应该通过 A/B 测试来确定要采取的具体行动。

最后，随着缓解措施的实施以及用户群的增长和演变，必须定期对该模型进行重新训练，以使该模型适应不断变化的条件。

原文：https://towardsdatascience.com/predicting-customer-churns-without-any-labelling-data-in-e-commerce-e54cb70944fe?source=collection_archive---------19-----------------------

Beautiful Cali

当你花了很大的力气去吸引你的客户，却不知道谁会去吸引时，你会感到困惑甚至沮丧吗？当你打算花 1 万美元，却只能获得 5 万美元的收入时，一切都感觉崩溃了(可能会让情况变得更糟)。我们发现在帮助你限制营销成本方面非常有用的一件事是定义和预测你的流失客户。有了一个精确的模型，你将能够削减高达 80%的成本，并减少你的竞选战略的不确定性的一大部分。让我们直接进入今天的话题吧！

首先，我们想谈谈为什么很难预测搅动。但在此之前，我们需要了解一下两种不同类型的 2C(对客户)商业模式。№1 是人们订阅和预付费的地方，这称为订阅模式。№2 是人们来买东西时没有留下任何关于他们何时会回来下一个订单的暗示，这被称为交易模式。

订阅模式有哪些例子？我会暂停 30 秒让你思考至少 2 秒。我的答案是苹果音乐，Spotify，亚马逊 Prime，纽约时报，甚至你的电，煤气，垃圾，下水道，有线电视都算在订阅里。订阅模式与众不同的是一个明确的契约。这意味着双方都知道钱和时间。因此，在客户流失建模中，这是一个相对简单的案例。

事务模型的一些例子是什么？每一个销售有形产品而没有签约重复的品牌都是一种交易模式。您现在将立即理解为什么在这种情况下更难对搅动进行建模，因为搅动实际上是不可见和不透明的。

我们将在事务模型中讨论搅动，因为我们喜欢挑战。但是为了克服挑战，我们需要变得更聪明。我们知道在监督学习的框架下解决这个问题是不可能的。为什么？因为作为一名员工，99%的时候你会发现自己在为一个项目争取多一天的时间。当你有一个监督学习的问题，但你没有足够的标记数据时，通常正确的解决方案是找到另一种方法，而不是为自己自我标记 1000，2000 个数据。时间很重要。(参见我的另一篇帖子关于作为初级数据科学家你肯定会陷入的一些陷阱。)

那么我们如何解决这个问题呢？我们需要转向营销中另一个更广为人知和使用的问题，顾客终身价值(CLV)。

CLV 是做什么的？CLV 是一个相当简单的模型，它可以根据顾客的历史行为来估计你能从他/她那里得到多少。让它简单有效的是 CLV 分析中只有三个组成部分，频率、最近和货币价值。频率是到目前为止的订单数量。最近是第一个订单和最后一个订单之间的天数。几乎每个人都很难记住这个定义，因为 recency 与“最近”有关，后者本质上代表你最后一次做某事是什么时候？这与新近的定义相反。你需要翻转它来记住它。货币价值是所有订单的平均价格。在我们的客户流失模型中，我们不太在乎钱，所以我们暂时不考虑它。

在继续之前，我想请读者暂停 30 秒，思考一下为什么这种模式在营销中如此有用和受欢迎。我的答案是 a)它对每个人都有意义 b)因为它直觉上有意义，它一定抓住了人类行为的本质。

现在是时候真正了解这个模型了。我会把所有的数学都交给最初的论文，因为不是这里的每个人都同意这种复杂性。我们更关心可用性而不是理论。该模型的基本思想是这样的:我们通过某个时间单位(天)观察每个客户，每个客户在任何给定的一天都有一定的订购概率，每个客户在任何给定的一天都有一定的流失概率。我们假设这两个概率是独立的，并且随时间保持不变。然后，我们可以建立一些方程来描述交易行为，并使用数学来估计概率。使用 Python，我们可以充分利用这个包。

首先，您需要使用以下命令安装这个包:

pip install lifetimes

该产品包包括 4 个可用于我们客户流失模型的模型，如果您感兴趣，可以查看gamma gamma fitter以全面评估客户生命周期价值建模之旅。****

我个人建议在最初的尝试中使用 BetaGeoFitter 或 ModifiedBetaGeoFitter 。原因是它们相对更容易收敛，训练速度更快(这也是原作者开发它们而不是经典的paretonbfilter的重要原因之一)。我使用betageobitabinomfitter的经验告诉我，这个模型也比其他两个模型需要更多的时间来训练。

我们现在需要一些训练数据，您可以通过在 python 中键入以下命令来获得

from lifetimes.datasets import load_cdnow_summary
data = load_cdnow_summary(index_col=[0])

我在这篇文章中使用了一些其他数据源，但是您可以预期使用这个示例数据集生成类似的曲线图、图表和数字。

查看数据集，我们有以下结构:

print(data.head())
*"""*
 *frequency   recency      T*
*ID*
*1    2           30.43       38.86
"""*

频率是某个客户的历史订单数量。最近度是第一次发现和最后一次购买之间的天数。 T 是客户被发现的天数。如果你使用自己的数据集，你需要以同样的方式组织你的数据。

我将在我的帖子中使用 ModifiedBetaGeoFitter 。

mbgf = ModifiedBetaGeoFitter(penalizer_coef=0.001)# penalizer_coef is helpful in converging the trainingmbgf.fit(summary_train['frequency'], summary_train['recency'], summary_train['T'])# display some parameters of the model
display(mbgf)# <lifetimes.ModifiedBetaGeoFitter: fitted with xxx subjects, a: 1.10, alpha: 3.27, b: 0.05, r: 0.85>

这个模型已经训练成功，我们希望得到它的一些视觉效果:

plot_probability_alive_matrix(mbgf)

The brighter, the more likely a customer has not churned

该图描述了新近度和频率与流失可能性之间的关系。给定一定的频率，较高的新近度增加了不流失的可能性。这是有道理的，因为高新近性意味着客户最近的活动表明他们坚持使用该品牌。给定一定程度的新近性，更高的频率增加了流失的可能性。这可能看起来不直观，因为我们通常认为顾客买得越多，他们就越喜欢这个品牌，离开的可能性就越小。但是如果我解释给你听的话，它就不一样了:

假设你经营一家超市，你有一些忠实的顾客，你至少每隔一天就会见到他们。我们把其中一个叫做贝拉吧。她一年前就开始来了，她从来没有错过一天，所以你真的认为她很忠诚。有一天你没有见到她，第二天你又没有见到她，你开始自言自语“她一切都好吗？我会永远失去她吗？”。当你从一两天前开始不再见到一些忠实的客户时，你的担心就会立即产生。因为你非常关心你的敬业的人，所以你不会太关心一年只来一次的人。现在，根据你对他们是否会回来的担心程度，想想他们回来的可能性。一个忠诚的人不再拜访你，这比你 7 个月没见一个人，但你只希望一年见一次更强烈的信号。正是这种行为的改变引起了我们的担忧，我们真的感觉到了。

现在我们完全理解了这个模型能给我们什么样的预测，接下来我们真正的预测是什么:

summary_train['mbgf_lh'] = mbgf.conditional_probability_alive(summary_train['frequency'], summary_train['recency'], summary_train['T'])summary_train['mbgf_lh'].plot.hist()

Most customers have low probability to return, but you can manage it before it happens

这张图显示了未来返回可能性的分布。我们注意到，大多数人不太可能回来(在电子商务中，大多数人都是买了就走，再也不回来)。现在我遵守了我的诺言，带你走到了这段旅程的终点。

我想给你留一个作业，你可以在下面评论:你如何利用这个情节来管理你的电子商务业务？你能从中赚钱吗？如果有，如何实现？

一些有用的资源:

彼得·法德尔教授在沃顿商学院谈论顾客终身价值。

他的论文:

*盘点你的客户:他们是谁，接下来会做什么？*T3，

原文：https://towardsdatascience.com/predicting-customer-lifetime-value-with-buy-til-you-die-probabilistic-models-in-python-f5cac78758d9?source=collection_archive---------1-----------------------

Tokyo, Japan - photo credit: Pexels

客户的价值是什么？在搅动之前，客户还会购买多少次？他在未来 3 个月内流失的可能性有多大？最重要的是，我们应该期望客户“存活”多久？

虽然这些问题在营销、产品、风险投资和公司财务专业人士中很常见，但总是很难用准确的数字来正确回答。

在的非契约性商业环境，顾客可以随时终止与零售商的关系，无需事先通知，这可能更加棘手。
亚马逊的图书(或任何其他没有订阅的产品类别)，Zalando 的服装，Booking.com 的酒店都是非合同商业设置的例子。对于所有这三种电子商务，我们无法通过查看客户合同的结束日期来了解他是“活着”(将来会购买)还是“死了”(永远不会再购买)。我们只能依靠客户过去的购买和其他不太典型的事件(网站访问、评论等)。).
但在这种情况下，我们如何决定客户是会回来还是永远离开呢？

“买到死”概率模型通过评估客户未来交易的预期次数及其“活着”的概率，帮助我们量化客户的终身价值。

为了理解“买直到你死”模型是如何工作的，我们把重点放在预测现实生活数据的最佳选择上:BG/NBD 模型。
贝塔几何/负二项分布模型于 2004 年由 P. Fader 的论文提出，是对 Schmittlein 等人于 1987 年开发的帕累托/NBD 模型(第一个 BTYD)的改进。

特别是，为了预测未来的交易，该模型将客户的购买行为视为掷硬币游戏。
每位顾客有 2 枚硬币:一枚购买硬币控制顾客购买的概率，一枚骰子硬币控制顾客退出且不再购买的概率。

让我们通过模型假设来理解一切是如何进行的。

假设 1: 活跃时，客户的交易笔数遵循一个 泊松过程 ，交易率λ(=一个时间间隔内的预期交易笔数)。

A customer’s purchasing behavior observed over a period of 12 months, where the number of transactions is distributed as a Poisson Process with unobserved transaction rate c

在特定时间间隔(12 个月)的每个子周期(1 个月)，每个顾客投掷他的购买硬币，根据结果，他购买或不购买。
我们在周期内观察到的交易数量(头数)取决于每个客户在λ附近的概率分布。
让我们在客户的泊松概率分布下方绘图，以形象化我们刚才所说的内容。

Poisson Probability Mass Function of a customer with λ = 4.3

这里我们假设我们的随机客户的交易率λ = 4.3。
因此，他有 19%的概率在随机的 12 个月内购买 4 次，有 4%的概率购买 8 次，等等。

假设 2: 客户间交易率的异质性遵循一个 伽玛分布 。

这相当于说每个顾客都有自己的购买硬币(有自己的头尾概率)。 为了更好地理解这一假设，我们模拟了 100 个客户的泊松分布，其中每个λ都用伽马分布建模，参数为:形状=9，比例=0.5。

Simulation of 100 customers Poisson Probability distributions where each customer’s λ depends on a Gamma distribution with shape = 9 and scale = 0.5

如假设中所述，在给定的时间间隔内，每个顾客都有自己购买 x 次的概率。

假设 3: 在任何一笔交易之后， 一个客户以概率 p 变得不活跃 因此，客户“退出”的点按照一个 (移位)几何分布 跨交易分布。

每次交易结束后，每位顾客都会投掷第二枚硬币，即骰子硬币。 鉴于 p 是“死亡”的概率，那么我们可以定义 P(活着)= 1-p. 再一次，我们来绘制一个随机的客户概率分布，更好地把握这个假设的含义。

Shifted Geometric Probability Mass Function for a customer with p = 0.52

假设我们的客户变得不活跃的概率 p = 0.52，那么他在第二次交易后变得不活跃的概率是 25%，他在第三次交易后变得不活跃的概率是 12%。

正如你所看到的，顾客买的越多，他活着的概率就越高。

***假设 4:***p 中的异质性遵循一个 贝塔分布 。

至于购买硬币，每个顾客都有自己的硬币，在特定数量的交易后，它有自己的存活概率。
我们可以在下面看到如何寻找 10 个客户的模拟，其中 p 遵循 Beta 分布，其中 α = 2， β = 3。

Simulation of 10 random customers Geometric Probability distributions where p is follows a beta distribution with parameters α = 2 and β = 3

假设 5: t 交易率 λ 和退出概率 p 在客户之间独立变化。

最后，通过对历史客户数据拟合前面提到的分布，我们能够推导出一个模型，该模型为每个客户提供:

• P(X(t) = x | λ，p) - 在长度为 t 的时间段内观察到 X 笔交易的概率
• E(X(t) | λ，p) - 长度为 t 的时间段内的预期交易数
• P(τ > t) - 客户在 τ期间变得不活跃的概率

然后，拟合的分布参数被用于前瞻性的基于客户的分析中，以找到具有由 x，t ₓ ，T — 定义的过去观察到的行为的个人在长度为 t 的未来时间段中的预期交易次数，其中 x =历史交易次数， t ₓ =最后购买时间， T =客户年龄。
下面是数学爱好者的最终公式(详细推导见 P. Fader 论文的附录):

The expected number of transactions in a future period of length t for an individual with past observed behavior (X = x, tₓ, T; where x = n. historical transactions, tₓ = time of last purchase and T = Age of a customer) given the fitted model parameters r, α, a, b

既然我们已经了解了“买到死”模型是如何工作的，我们终于准备好从理论到实践，将 BG/NBD 模型应用于真实的客户数据。

在可用于实现模型的各种备选方案中，我强烈推荐 Python 中的生存期包(在此使用)和 r 中的 BTYD 库，这些包在将模型的方程包装成方便的函数方面做得非常好，使我们的生活变得非常容易。

如前所述，BG/NBD 模型将几种分布拟合到历史客户购买数据中。为此，我们需要建立一个数据集，为每个客户提供以下三个信息:

最近(来源于 t ₓ) :顾客最后一次购买时的年龄，等于顾客第一次购买和最后一次购买之间的时间。

频率(x) :客户重复购买的次数。

客户年龄(T) :客户在研究期间结束时的年龄，等于客户第一次购买和数据集中最后一天之间的持续时间。

下面是您的数据应该是什么样子的示例:

您是选择几天、几个月还是几年，很大程度上取决于您典型的客户购买周期。例如，送餐业务往往会在同一周内遇到重复的顾客，所以他们可能会去几天。在上面的例子中，我用了月，因为它更合适。

一旦我们创建了数据集，我们就可以将它传递给模型并打印出摘要。

如果您需要一些真实的客户数据，请随意使用 UCI ML 存储库提供的在线零售数据集(我将使用一个虚构的数据集，不提供任何有关我分析的公司的敏感信息)。

酷！这是什么？我们将假设的分布拟合到历史数据中，并推导出模型参数: alpha 和 r 用于伽马分布(假设 2)，而 a 和 b 用于贝塔分布(假设 4)。
在总结中，我们还为每个参数提供了一个置信区间，我们可以用它来计算每个客户的预期未来交易的置信区间。

现在我们已经建立了一个模型，我们可以检查它是否真的有意义。第一种方法是根据拟合的模型参数人工生成具有预期购买行为的客户，并将其与真实数据进行比较。

就我们所见，人为的客户分布与真实数据非常相似。
在这一步，我还建议计算总体百分比误差(=预测交易/实际交易-1)和 在校准周期 中完成的每笔交易的百分比误差。
通过这种方式，您可以量化模型与现实的接近程度，以及它是否更适合某些客户。例如，该模型在 5、6 和 7 个校准交易桶中放置的客户可能比实际少，这可能最终导致整体严重低估。

现在我们有了一个拟合的模型，我们可以查看它的频率/最近矩阵，以检查客户最近(上次购买时的年龄)、频率(重复交易的次数)和下一个时间段的预期交易次数之间的预期关系(基于我们的拟合模型参数)(下图左侧)。
我们还可以根据客户的近期情况和频率来想象客户存活的预期概率(下图)。

事实上，我们看到，如果一个客户购买了超过 25 次，并且他们最近一次购买是在他们超过 25 个月大的时候(右下角)，那么他们是你最热门的客户，最有可能活着并购买。相反，你最冷的客户在右上角:他们很快买了很多，我们几个月没见了。

一旦您验证了模型足够接近实际数据，我们就可以看到它在预测未来购买方面有多好。

得益于 Lifetimes '*calibration _ and _ holdout _ data()*函数，我们可以将简单的事务数据集快速拆分为校准期和维持期。我们首先将模型拟合到 2 年的校准期，然后预测下一年的事务，最后比较预测的事务与维持的事务。

下面是 cal_hold 数据帧的样子:

通过比较下图中的平均实际购买量和预测购买量，我们可以注意到，对于在校准期内有 0 到 3 次重复交易的客户，预测购买量和实际购买量非常接近，而对于重复交易次数较多的客户，预测购买量和实际购买量会越来越大。

同样在这种情况下，由于图表本身可能会产生误导—对于有很多客户的时段，小的差异可能会导致大的误差—为了正确地评估预测，我们应该查看总体百分比误差(图表中的预测误差)、，在我的场景中它占-6.3%(模型预测的交易比实际少 6.3%)。
通过查看校准期内重复测试的误差百分比，我们注意到，我们低估了重复测试 3 次或更多次的客户(-7.3%至-30.3%)，以及没有重复测试的客户(-12.2%)。我们还强烈高估了 1 个回头客(+19.7%)。

尽管我们在更细的层次上预测失误，但这是一个相当好的结果，因为我们感兴趣的是交易的总体预测量。尽管如此，仅仅在一个时期内的交叉验证并不能让我们理解未来会发生什么。未来预期的-6.3%的误差合理吗，或者这是一个难以置信的幸运机会？

为了更好地评估模型，我们将在几个期间进行交叉验证，然后检查每个期间的错误。为此，我们简单地构建一个 For 循环，在几个周期内迭代交叉验证。
具体来说，我们在具有 6 年交易的数据集上运行它，其中对于每次迭代，我们在选定的 2 年校准期内对客户子集进行采样，并预测未来 1 年的交易。因此，我们最终有 4 个交叉验证的期间。
然后我们标绘结果。

如下图所示，每年的预测误差始终在 4.1%到-7.9%之间，2018 年是预测最好的一年(-2.3%预测误差)。这是非常好的，尤其是与绝对预测误差大于 10%的普通队列模型相比。

一旦你建立了模型并验证了它的有效性，你就可以很容易地查看单个客户的预测以及他们存活的概率。这非常有价值，因为你可以将 CLV 预测用于营销活动、预测或更普遍的防止流失。

例如，在下面，我们绘制了客户存活的历史概率。

从图中我们可以观察到顾客每增加一次购买，他活着的概率就增加，然后又开始下降；但是速度较慢，因为每个新的交易增加了他的频率和他的新近度。

总之，预测 CLV 总是一项棘手的任务，通常历史频率模型无法区分过去购买次数接近的客户。
购买直到你买概率模型来拯救我们，它允许我们仅使用 3 个客户的信息(客户的频率、最近和年龄)来建立相当准确的预测。

在本文中，我们有意没有提及一些与 CLV 有关的重要话题，但在让您了解之前，让我对其中的三个话题做一个简短的评论:

• 我们预测了未来的交易，但是我们忽略了 CLV 等式中的“价值”部分。通常伽马-伽马子模型用于 BG/NBD 模型之上，以估计交易的货币价值。
• 我们使用 BG/NBD 模型来预测未来交易总量，但如果您打算在用户层面采取行动，您应该正确衡量单个客户预测的准确性。
• 如果模型不太适合您的客户(并且这些假设对您的业务来说是合理的)，那么考虑让它适合客户群体(例如。按用户国家划分)，和/或 将其与具有附加功能的线性模型(例如。客户的网站访问量、自上次访问以来的时间、产品评论、获取渠道等。).

感谢阅读，用数据不断改造世界！

建设性的反馈和激励性的谈话总是受欢迎的。随意连接在Linkedin上打招呼！

原文：https://towardsdatascience.com/predicting-dengue-in-singapore-f9be0a761ce4?source=collection_archive---------19-----------------------

你正在公园里享受一个愉快的夜晚，野餐垫已经摆好，手里拿着一杯可口的冷饮。

“啪！”

你刚刚杀死了一只埃及伊蚊。你欣赏它手上独特的黑白条纹身体。你认为只是无害的咬一口。

接下来的一周，你躺在床上发高烧，肌肉疼痛，出皮疹。究竟发生了什么事？你从两周前孵化的伊蚊身上感染了登革热。

登革热是一种可怕的疾病。如果没有及时或有效的干预，散发的登革热病例很容易演变成流行病。伊蚊通过吸食受感染的人类或其母亲感染登革热，后者在其两周的寿命内可以产下数百个卵。

如果流行病发生，早期预警系统可以帮助医疗保健提供者更好地管理登革热患者的增加。

天气信息:雌性伊蚊喜欢炎热多雨的天气。它们在死水中产卵，温暖的温度促使它们的后代生长得更快。

**谷歌趋势:**包含“登革热”一词的搜索词已被成功用于预测多个国家的登革热病例。这可能是因为那些对病毒更敏感的人，也更有可能搜索它。

从网络上收集了 2014 年至 2018 年新加坡特定的登革热病例数据、天气信息和搜索兴趣。

使用 Selenium 和 BeautifulSoup 从政府网站上搜集登革热和天气数据:

1. 每周登革热病例数。报废时只有 2014 年到 2018 年的数据。此限制为所有其他数据源设置时间段。
2. 每日天气信息。只有樟宜气象站的数据被废弃；这个车站是政府的主要历史参考点。

Google Trends search interest data for search term ‘dengue’ from 2014 to 2018.

以下词语的每周搜索兴趣数据来自 Google Trends:

• 疼痛或疼痛
• 登革热
• 发热
• 头痛
• 呕吐
• 恶心
• 发疹
• 目痛

数据集相对干净，除了重复和缺失的值，它们都被删除了。登革出血热病例也从登革热数据集中删除，因为这类病例很少。

主要挑战是创建一个跨数据集的公共时间要素来合并它们。使用流行病学周(e 周，例如 2018 年第 52 周)对登革热数据集进行了日期确定。 Epiweeks 软件包用于将数据集的日期特征转换为 e-weeks。

对于天气数据集，使用 groupby 获得每个 e 周的平均温度和降雨量。

在将数据拟合到模型之前，删除了所有日期时间功能。

时滞被用来解释感染者出现症状和伊蚊孵化所需的时间。

Assumed timeline of events

对于每个每周的登革热观察，其对应的搜索兴趣数据将来自一周前，而天气数据将来自两周前。原因如下:

• 登革热的症状通常在被感染后的 4 到 7 天内出现。我认为在此期间，潜在的患者可能会在谷歌上搜索登革热及其症状。
• 天气变化被认为会影响蚊子的繁殖速度。假设是较高的温度促使伊蚊产卵；伊蚊大约需要一周的时间孵化和生长。

合并的数据框架被分成 20%的测试数据和 80%的训练验证数据(用于 k 倍交叉验证)。

通过取最高和最低温度之间的差值创建了温度范围特征。我预计较小的温度范围预示着登革热病例，尤其是在较温暖的几周。

对登革热病例数(目标)和对登革热的搜索兴趣(特征)进行对数转换，以校正偏差并提高其相关性。

在模型拟合之前，使用标准定标器对所有变量进行标准化。

在模型构建过程的不同阶段进行特征选择。

1. **在构建模型之前:**为了避免多重共线性，高度相关的要素被移除。
2. 交叉验证确定线性模型可能适用于数据后: Lasso 回归用于特征选择。使用 LassoCV 获得 alpha 值，并应用于模型。移除了具有零值系数的要素。
3. 选择合适的线性模型后:使用向后逐步方法选择特征。具有最高 p 值的要素被移除，直到模型的校正 r 减小。

以下是这些变量及其被删除原因的摘要:

Table of variables, and their status (removed/included in model)

由于观察值数量较少(248，在去除了因对数据应用时间滞后而产生的缺失值后)，因此模型选择使用了 5 k 倍交叉验证。

以下是所用线性模型的总结，以及它们在 k 倍范围内的 r 均值和标准差。对于脊和套索，alpha 值设置为默认值 1。

Simple linear model r^2: 0.755 +- 0.041 
Polynomial model r^2: 0.441 +- 0.153
Ridge model r^2: 0.755 +- 0.043
Lasso model r^2: -0.010 +- 0.012

为了在简单模型和脊模型之间进行选择，在使用 Lasso 回归进行特征选择并应用从 RidgeCV 获得的 alpha 值之后，我重新运行了交叉验证。

Simple linear model r^2:: 0.757 +- 0.043 
Polynomial model r^2: 0.665 +- 0.053
Ridge model r^2: 0.757 +- 0.046
Lasso model r^2: -0.010 +- 0.012

两个模型的 r 均值和标准差相似，因此选择了一个简单的线性回归。

使用 StatsModel 对整个 80%的训练验证数据重新训练模型。为了进一步简化模型，进行了向后逐步特征选择。

以下是最终模型的结果:

StatsModel Output for Training-Validation Set

将测试数据拟合到模型后，得到的调整后的 r 为:0.760。这个模型似乎概括得很好。

Left to Right: Scatterplot of y-predicted and y-residuals; Q-Q plot of y-residuals; Scatterplot of e-week and y-residuals.

散点图中没有可辨别的模式，y 残差似乎呈正态分布，除了尾部的一些偏差。

从 StatsModel 输出和上述图表判断，线性回归模型似乎适用于该数据。

该模型可以解释新加坡每周登革热病例记录中 77%的差异，并具有良好的泛化能力。

在预测未来登革热病例时，对“登革热”的搜索兴趣似乎是最重要的。这一专题可能已经包含了一些关于蚊子繁殖和传播的信息，因为新加坡新闻媒体倾向于报道繁殖集群、登革热感染和相关死亡。

因此，医疗保健提供者可以使用登革热搜索兴趣 作为粗略的指标，而不是监控大范围的变量。

我对头痛和疼痛的显著负系数感到惊讶。新加坡对登革热的认识很高。我怀疑，当蚊子叮咬后出现这些症状或登革热病例增加时，意识到登革热的人可能会直接搜索“登革热”。因此，更好地理解搜索模式有助于改进模型。

*** GitHub ***

** * 这是一个 METIS 数据科学训练营项目。线性回归和网页抓取是项目需求***

原文：https://towardsdatascience.com/predicting-earthquakes-with-xgboost-5feca64d28f6?source=collection_archive---------36-----------------------

我最近创建了一个 Plot.ly Dash web 应用程序，它处理地震的震级和深度，以预测它们在美国和美国领土内的位置。

很多时候，当我开始一个新的机器学习项目，并开始查看我的数据时，我会先发制人地考虑模型和验证。尽管我喜欢筛选 NaNs 和脏数据，但我更喜欢拟合模型和预测结果。尽管如此，我倾向于最初选择一个最小可行的产品，然后向上滚动开始数据清理过程。数据清理和模型选择一样重要，这一点在这个项目中得到了很好的证明，在这个项目中，我从一个随机的森林分类器开始，得到了低得多的分数。最后，在最后一秒，我将模型切换到 XGBoost 分类器，这使我的分数提高了约 0.01%，但对复制的 gap-per 预测大幅下降。

一开始，我很快意识到我的模型将围绕着根据地震的属性对州名进行分类。预期的结果将是明确的，并使用地震的属性进行预测。我从一个快速、极简的数据清理开始，为了避免一些严重的泄漏，我做了一些删减，包括经度和纬度列，这会使算法完全不相关。接下来，我开始了一些经典的数据探索。我决定用直方图开始我的分析，以确定州数据中某些结果的频率。这个直方图对于确定我的过程中的下一步，以及在开始拟合模型之前，数据下一步需要去哪里是必不可少的。

California is separated marginally from the rest of the data.

我对这一结果的最初反应是恐惧，这是理所当然的；如果加州比其他地方的数据要对更多的地震负责，那么真的有必要应用一个模型吗？这个问题很好地延续到任何数据科学家将对他们的数据做的下一件事。所以我用 iter-tools 写了一个函数，它能给我一个很好的主意，如果没有一个真实的模型，预测会有多准确，这当然也被称为多数类基线。这将让我知道，如果我们提供的任何数据(x)只猜测加利福尼亚，这个模型会有多准确。此外，您可能已经注意到直方图中有一个“CA ”,所以我很快将其替换为“California ”,以使命名更加一致，并且对该州有一个单独的预测。有趣的是，这是这个特定问题暴露出来的唯一地点。在将查询“CA”替换为“California”之后，我继续获取我的多数类基线。

def most_common(L):
  SL = sorted((x, i) for i, x in enumerate(L))
  groups = itertools.groupby(SL, key=operator.itemgetter(0))
  def _auxfun(g):
    item, iterable = g
    count = 0
    min_index = len(L)
    for _, where in iterable:
      count += 1
      min_index = min(min_index, where)
    return count, -min_index
  return max(groups, key=_auxfun)[0]
def baseliner(y,ypr):
  r = []
  m = most_common(y)
  for i in ypr:
    r.append(m)
  return(r)

我快速编写了两个函数，并使用测试集中的目标和特性运行该函数，产生了以下结果:

from sklearn.metrics import accuracy_score# Setting our majority class accuracy
acc = accuracy_score(ytest,y_pred)
# Printing out the results:
print("Baseline Accuracy:",acc)

结果是:

Baseline Accuracy: 0.5108583247156153

至少可以说，我非常震惊！令人高兴的是，这表明数据对于一点点机器学习仍然是可行的，生活中没有比这更令人兴奋的事情了。考虑到这一点，我用一个顺序编码器、迭代估算器和最后的 XGBoost 分类器设置并安装了我的管道。当内核计算结果时，我积极地预期结果，希望模型最初会获得很高的分数，因为仍然有很大的改进空间，使我能够获得非常高的训练准确度分数，这使我对模型的验证非常兴奋。

pipeline.fit(Xtrain,ytrain)

我安装了管道，设置了模型，并准备获得一些验证准确性，看看这个项目的未来会发生什么。这听起来好像我是在杞人忧天，但这真的是“数据科学家的事情”，你不会明白的。我在从我的 test-x 预测我的 test-y 时初始化了管道。

ypr = pipeline.predict(Xtest)
accuracy_rf = accuracy_score(ytest,ypr)

我终于得到了我的准确度分数，为了制造悬念，我不得不为我最初的准确度做一个夸张的展示:

print("===================================================================")
print("___Unscaled, continous encoded Xgb Classifier Accuracy___")
print("_________________________",accuracy_rf,"______________________")
print("===================================================================")

这产生了:

===================================================================
___Unscaled, continous encoded Random Forest Classifier Accuracy___
_________________________ 0.8749038445433443 ______________________
===================================================================

在这之后，我回去清理了数据——稍微深入一点，我最终回去删除了“时间”一栏，因为害怕泄露。此外，我优化了一些超参数，以小幅提高精度，最高可达:

===================================================================
_______________________Xgb Classifier Accuracy______________________
_________________________ 0.9118407445708376 ______________________
===================================================================

这表明，当数据猴子不断提醒你不断清理数据的重要性时，它们可能终究是正确的。

在获得一个相当准确的模型后，我决定在将所有这些都转移到 Dash-app 之前，用 Plotly 对数据进行更多的探索，从一些幅度的依赖关系图开始。

虽然这肯定不是该项目的必要获得，也没有必要说明任何不可能假设的事情，但它确实显示了一些关于这个特性在实际预测中有多少作用的事情，并且看到它与其他特性相比的位置是非常令人惊讶的。接下来，a 用各州的震级标绘了一些点:

我认为这可能是一个更重要的可视化，因为现在我们可以看到为什么多数类不是一个非常准确的预测器，通常某一级别的地震在位置上是一致的，至少在某种程度上是如此。此外，值得注意的是，气泡的大小是地震的间隙，(这在悬停数据中更明显，在停滞的截图中不太明显)。在 3.0-4.0 的震级范围内，加利福尼亚州不太可能发生地震，被多米尼加共和国忽视了，该国也没有发生-1.0-2.5 的地震。

虽然这看起来像是旧闻，特别是对那些有领域知识的人来说，但了解一下为什么准确度可以如此容易地非常高是非常有趣的，尽管该模型需要的工作比清理数据要少。

当然，下一步是将这个模型放入 Dash 应用程序，提供一个用户友好的体验，使用这个模型预测地震的准确率达到 91%。但是我认为现在是提供一个笔记本链接的好时机，它也将出现在页脚，因为我们正在从笔记本转移到 flask/dash。

对我来说，这是这样一个项目的另一部分，我觉得比清理数据更有趣。很多时候，整天埋头于笔记本会变得乏味，感觉有点做作。当然，我不得不从矢量图形开始，但是我是数据科学家，不是图形设计师，所以我不太热衷于提供更多关于 Adobe Illustrator 设计的信息。

为了创建 Dash 应用程序，我需要知道更多关于数据的事情。由于数据主要是振动扰动(seismec)，它主要以 double 数据类型存储，所以我的范围主要由一些非常精确的小数组成。此外，我需要分类特征观测值中的所有唯一值。我决定快速创建一个新的笔记本，我可以用它来填充我的半个屏幕，以便快速简单地将最小值和最大值输入到我的 IDE 中。这让我真的希望有一种简单的方法在 Dash 中将列表添加到组合框中。

这里有一个笔记本的链接。我使用了 Numpy idx-min 和 max，以及 Numpy unique，这两个都是熊猫内置的，这里没有什么真正严重或复杂的事情发生。

至少对于 slider 的来说，更自动化的方法是可能的，通过在 DCC 中使用 min()和 max()作为 slider 构造函数中的 max 和 min 值。

由于这不是一个情节性的 Dash 教程，(也许有一天)，我不会解释如何使用 Dash 进行机器学习。但如果你愿意，你可以看看这个应用程序，它将于今天下午(9 月 27 日)推出，这里。而至于来源，这里有。

我真诚地希望这次旅行至少能读起来有趣，或者甚至有一点教育意义。构建 Dash 应用程序是一种简单、快速、有趣且富有创造性的方式，可以展示您的数据科学实力，同时还可以产生一个伟大的项目，并在此过程中学习一些东西。将来我肯定会在空闲时间赚更多的钱。

原文：https://towardsdatascience.com/predicting-electricity-demand-in-la-outperforming-the-government-a0921463fde8?source=collection_archive---------13-----------------------

这是一篇小博文，描述了我为 Springboard 的数据科学职业轨道所做的第一个顶点项目。我比较了一些根据天气特征预测电力需求的机器学习模型。我发现大多数模型对电力需求的预测比政府预测要准确 5%。查看我的 GitHub 库，获取该项目的代码、图表和数据集。

为了训练和测试我的模型，我使用了 2015 年 7 月至 2018 年 9 月洛杉矶约 27，000 个小时的电力需求测量值。该数据由洛杉矶水电部记录，并通过能源信息管理局的公共 API 发布。

对于天气数据，我从美国国家海洋和大气管理局(NOAA)网站检索每小时的当地气候数据(LCD) 。在 LCD 数据集中，我们有熟悉的测量方法，如压力和温度，以及不太直接相关的量，如风向。增加了一个额外的功能:每小时的冷却和加热程度。加热/冷却度天数分别是低于/高于 65 华氏度(负值设置为零)的度数，代表该温度所需的加热/冷却总量。我们特别预计 CDD 将成为洛杉矶电力需求的强劲代表，因为洛杉矶的主要能源消耗在夏季会降温。图 1 显示了 2015 年 7 月至 2018 年 9 月 LA 的电力需求时间流程。请注意夏季的需求高峰。

Figure 1. Electricity demand versus time in LA. Note the sinusoidal pattern, with the summer season experiencing more demand because of air conditioning.

经过一天的初步清理和准备，我进行了一些探索性的数据分析，以查看独立变量(天气)和非独立变量(电力)之间的初始关系。图 2 显示了各特性与电力需求之间的皮尔逊相关系数和图 3r值和值。我们发现，日冷却度天数与电力需求有很强的正相关关系——绝大多数电力使用是在洛杉矶的夏季用于冷却。

# Code to produce figure 2
import pandas# df is our main dataset with weather + electricity data
print('DEMAND CORRELATIONS (PEARSON) FOR LA')
print(df.corr()['demand'].sort_values(ascending=False)[1:])

Figure 2. Pearson correlation coefficients for all weather features with electricity demand. Cooling degrees is a strong predictor of electricity demand.

# Code to produce figure 3
import scipy
import pandas as pd# get r^2 values per column and print
demand_r = {}
for col in df.columns:
 if col != 'demand':
  slope, intercept, r_value, p_value, std_err = scipy.stats.stats.linregress(df['demand'], df[col])
  demand_r[col] = r_value**2print('DEMAND CORRELATIONS (r^2) FOR LA')
demand_r_df = pd.DataFrame({'col': demand_r.keys(), 'r^2': demand_r.values()})
print(demand_r_df.sort_values(by='r^2', ascending=False))

Figure 3. r² values for all weather features with electricity demand.

我们的一些天气特征有共线性；有三个压力特征(气象站、海平面、高度计)和另外三个温度特征(湿球、干球、露点)。共线性增加了预测建模的难度，因此我们删除了具有最低 r 值的两个温度和压力特征。

由于人们在白天和晚上的行为是不同的，不管那天有多热或多冷，我添加了一个“hourlytimeofday”列，代表白天或晚上(早上 6 点到下午 6 点之间为 0，其他时间为 1)。温度数据对于识别加热和冷却峰值很重要，但是回归可以开始依赖温度作为白天的代理。这个特征大大增加了我们的多元回归 r，表明了它的重要性。

在使用其他机器学习技术之前，我们对所有特征与电力需求进行多重普通最小二乘(OLS)回归，以隔离不重要的特征。图 4 显示了结果。高 p 值特征证实了我们的直觉——加热对洛杉矶并不重要。天空条件，即云覆盖和风速，直观上不是预测电力需求的重要特征。这些不太有用的要素由高 p 值表示，p 值大于 0.1 的任何要素都将从数据集中删除。

Figure 4. Multivariate OLS regression for all the features + constant. Shown are a list of coefficients with errors, t-statistics, confidence intervals.

现在我们为机器学习建模准备数据集。我们的第一步是将时间序列数据转换成适合监督学习的形式。从这篇博文中大量提取，我转换了我们的数据集，这样我使用了前一时间步的所有特征(即电力需求加上所有天气特征)来预测当前时间步的电力需求。由于前一小时的电力需求似乎是当前电力需求的相关代理，因此该过程似乎是直观的。在以这种格式重新构建我们的数据后，我们将数据分成训练和测试集，然后开始评估我们的机器学习模型。为了使分析与我们后来使用的递归神经网络和长短期记忆一致，我们将第一年的数据指定为训练集，其余数据指定为测试集(约 32%训练和约 68%测试)。这种分割避免了过度拟合，并且在训练和测试上更传统的 80/20 分割产生了几乎相同的模型性能。建模的一般过程是:1 .)在所有功能上使用最小-最大缩放器以使它们更具可比性，2。)适合训练集，以及 3 .)计算测试集上的相关度量，并在最后比较所有模型。在我们的分析中，我们没有对模型进行超参数调整(sklearn 的默认设置),但仍然取得了非常显著的结果。

我使用 sklearn 模块训练和测试了大量的标准机器学习模型，如上所述。建模的具体过程和结果可以在本 IPython 笔记本中找到。对于直接开箱即用的模型，我获得了令人印象深刻的准确性，在我尝试的一半以上的模型中， r 值超过 0.95。此外，训练和测试的 r 值相差不超过百分之几，这表明我没有过度适应。

除了开箱即用的模型，我们还使用了循环神经网络(RNN)和长短期记忆(LSTM)。本 IPython 笔记本中的详细介绍了 RNNs 和 LSTM 的基础知识以及建模过程。使用这个模型，我获得了与我之前测试的最好的机器学习模型相当的性能。图 5 显示了使用 LSTM 建模的训练集和测试集的损失(平均绝对误差)与历元的关系。测试损失仍然略高于训练损失，表明模型没有过度拟合。

*# Code to produce figure 5
from keras.models import Sequential
from keras.layers import Dense, LSTM
import matplotlib.pyplot as plt# design network
model = Sequential()
model.add(LSTM(50, input_shape=(train_X.shape[1], train_X.shape[2])))
model.add(Dense(1))
model.compile(loss='mae', optimizer='adam')
# fit network
history = model.fit(train_X, train_y, epochs=50, batch_size=72, validation_data=(test_X, test_y), verbose=2, shuffle=False)
# plot history
plt.plot(history.history['loss'], label='train')
plt.plot(history.history['val_loss'], label='test')
plt.xlabel('Epoch')
plt.ylabel('Loss (MAE)')
plt.legend()
plt.show()*

Figure 5. Loss measured in mean-absolute error (MAE) versus epoch for both training and testing of our LSTM network. Testing error remains slightly above training, indicating that we are not overfitting. The LSTM network was optimized using the Adam implementation of stochastic gradient descent with MAE loss.

EIA 还提供了一天前的电力需求预测，我们希望将其与我们的模型进行比较。EIA 用各种方法预测需求。天气预报很大程度上取决于天气预报和一周中的第天。使用相同的测试集，我们发现 EIA 的前一天需求预测的准确性比我们的所有模型都差，只有一个模型除外，这表明了简单、开箱即用的机器学习模型在回归问题上的能力。

结果汇总在按 r 排序的表 1 中。总的来说，集成机器学习模型比其他模型表现得更好，尽管简单的线性回归本身表现得非常好。所有模型的训练和测试误差保持在彼此的百分比水平之内，这表明我们没有过度拟合数据。由于交叉验证方法在处理相关的时间流数据时变得更加棘手，我们在分析中省略了交叉验证，以保持所有模型的可比性和一致性。我们发现梯度增强模型具有最好的性能，尽管前五个模型的性能都相似(它们之间的性能差异约为 1%)。

Table 1. Table of results of machine learning models. The best performing model (gradient boosting) is highlighted.

选择了性能最佳的模型(梯度推进)后，我们现在想看看哪些特性在预测电力需求时最为重要。这显示在图 6 中。“特征重要性被计算为通过到达该节点的概率加权的节点杂质的减少。节点概率可以通过到达该节点的样本数除以样本总数来计算前一小时的需求是预测当前小时需求的最重要特征，但如我们所料，其他替代因素如降温天数也是一个强有力的预测因素。

Figure 6. Feature importance for each feature in the gradient boosting model. Demand at the previous time step is four times more important than cooling degree days. Reframing our supervised learning problem in this way greatly improved the predictive power of the models.

我们表明，通过机器学习，超越 EIA 的日前电力需求预测实际上相当容易。值得注意的是，除了一个模型(k-NN)之外，我们的所有模型都比 EIA 的模型表现得更好。来自公共数据源的相对简单的数据采集、清理和机器学习建模导致比传统方法更精确的建模。我们报告说，与 EIA 的模型相比，总体性能提高了约 5%;能源公司可以利用这些模型更好地应对高需求电力高峰，并最终增加利润。既然已经确定了最佳模型，那么可以执行额外的超参数调整和交叉验证技术来提高性能，尽管这种提高很可能是微不足道的，并且由于我们的模型已经优于 EIA 的模型，所以我们在没有调整的情况下得出结论。请查看这个项目的 GitHub以获得更详细的报告、代码、情节和文档。