移动广告点击率预测
Photo credit: Pixabay
在线广告，Google PPC，AdWords 活动，移动广告
在互联网营销中，点击率(CTR)是一种衡量广告客户在每一次展示中获得的点击量的指标。
移动已经与所有渠道无缝连接，移动是推动所有商业的驱动力。预计今年移动广告收入将达到 10.8 亿美元，比去年增长 122%。
在这项研究分析中，Criteo Labs 分享了 10 天的 Avazu 数据，供我们开发预测广告点击率(CTR)的模型。给定用户和他(或她)正在访问的页面。他(或她)点击给定广告的概率是多少？此分析的目标是为 CTR 估计的最精确 ML 算法设定基准。我们开始吧！
数据
数据集可以在这里找到。
数据字段
- id:广告标识符
- 单击:0/1 表示非单击/单击
- 小时:格式为 YYMMDDHH，因此 14091123 表示 UTC 2014 年 9 月 11 日 23:00。
- C1——匿名分类变量
- 横幅 _ 位置
- 站点 id
- 站点 _ 域
- 站点 _ 类别
- app_id
- app_domain
- app _ 类别
- 设备 id
- 设备 _ip
- 设备型号
- 设备类型
- 设备连接类型
- C14-C21-匿名分类变量
EDA 和特征工程
训练集包含超过 4000 万条记录，为了能够在本地处理，我们将随机抽取其中的 100 万条。
import numpy as n
import random
import pandas as pd
import gzipn = 40428967 #total number of records in the clickstream data
sample_size = 1000000
skip_values = sorted(random.sample(range(1,n), n-sample_size))parse_date = lambda val : pd.datetime.strptime(val, '%y%m%d%H')with gzip.open('train.gz') as f:
train = pd.read_csv(f, parse_dates = ['hour'], date_parser = parse_date, dtype=types_train, skiprows = skip_values)
Figure 1
由于匿名化，我们不知道每个特性中的每个值意味着什么。此外，大多数特征是分类的，并且大多数分类特征具有很多值。这使得 EDA 不太直观更容易混淆，但我们会尽力而为。
特性
我们可以将数据中的所有特征分为以下几类:
- 目标特征:单击
- 站点特征:站点 id，站点域，站点类别
- 应用功能:应用 id、应用域、应用类别
- 设备特征:设备标识、设备 ip、设备型号、设备类型、设备连接类型
- 匿名分类特征:C14-C21
import seaborn as sns
import matplotlib.pyplot as pltsns.countplot(x='click',data=train, palette='hls')
plt.show();
Figure 2
train['click'].value_counts()/len(train)
Figure 3
整体点击率约为。17%，而大约。83%没有被点击。
train.hour.describe()
Figure 4
该数据涵盖了从 2014 年 10 月 21 日到 2014 年 10 月 30 日的 10 天点击流数据，即 240 个小时。
train.groupby('hour').agg({'click':'sum'}).plot(figsize=(12,6))
plt.ylabel('Number of clicks')
plt.title('Number of clicks by hour');
Figure 5
每小时的点击模式每天看起来都很相似。然而，有几个高峰时间，一个是在 10 月 22 日中午，另一个是在 10 月 28 日中午。10 月 24 日午夜时分，点击量非常低。
日期时间特征的特征工程
小时
从日期时间特征中提取小时。
train['hour_of_day'] = train.hour.apply(lambda x: x.hour)
train.groupby('hour_of_day').agg({'click':'sum'}).plot(figsize=(12,6))
plt.ylabel('Number of clicks')
plt.title('click trends by hour of day');
Figure 6
一般来说，点击数最高的是 13 和 14 小时(下午 1 点和 2 点)，点击数最低的是 0 小时(午夜)。对于粗略估计，这似乎是一个有用的特性。
让我们考虑一下印象。
train.groupby(['hour_of_day', 'click']).size().unstack().plot(kind='bar', title="Hour of Day", figsize=(12,6))
plt.ylabel('count')
plt.title('Hourly impressions vs. clicks');
Figure 7
这里没有什么令人震惊的。
既然我们已经看了点击和印象。我们可以计算点击率(CTR)。CTR 是广告点击与印象的比率。它衡量每个广告的点击率。
每小时一次
import seaborn as snsdf_click = train[train['click'] == 1]
df_hour = train[['hour_of_day','click']].groupby(['hour_of_day']).count().reset_index()
df_hour = df_hour.rename(columns={'click': 'impressions'})
df_hour['clicks'] = df_click[['hour_of_day','click']].groupby(['hour_of_day']).count().reset_index()['click']
df_hour['CTR'] = df_hour['clicks']/df_hour['impressions']*100plt.figure(figsize=(12,6))
sns.barplot(y='CTR', x='hour_of_day', data=df_hour)
plt.title('Hourly CTR');
Figure 8
这里一个有趣的观察是最高的点击率发生在午夜，1，7 和 15。如果你记得的话，午夜前后的浏览量和点击量最少。
星期几
train['day_of_week'] = train['hour'].apply(lambda val: val.weekday_name)
cats = ['Monday', 'Tuesday', 'Wednesday', 'Thursday', 'Friday', 'Saturday', 'Sunday']
train.groupby('day_of_week').agg({'click':'sum'}).reindex(cats).plot(figsize=(12,6))
ticks = list(range(0, 7, 1)) # points on the x axis where you want the label to appear
labels = "Mon Tues Weds Thurs Fri Sat Sun".split()
plt.xticks(ticks, labels)
plt.title('click trends by day of week');
Figure 9
train.groupby(['day_of_week','click']).size().unstack().reindex(cats).plot(kind='bar', title="Day of the Week", figsize=(12,6))
ticks = list(range(0, 7, 1)) # points on the x axis where you want the label to appear
labels = "Mon Tues Weds Thurs Fri Sat Sun".split()
plt.xticks(ticks, labels)
plt.title('Impressions vs. clicks by day of week');
Figure 10
周二的浏览量和点击数最多，周三次之，周四次之。周一和周五的浏览量和点击量最少。
星期几中心
df_click = train[train['click'] == 1]
df_dayofweek = train[['day_of_week','click']].groupby(['day_of_week']).count().reset_index()
df_dayofweek = df_dayofweek.rename(columns={'click': 'impressions'})
df_dayofweek['clicks'] = df_click[['day_of_week','click']].groupby(['day_of_week']).count().reset_index()['click']
df_dayofweek['CTR'] = df_dayofweek['clicks']/df_dayofweek['impressions']*100plt.figure(figsize=(12,6))
sns.barplot(y='CTR', x='day_of_week', data=df_dayofweek, order=['Monday', 'Tuesday', 'Wednesday', 'Thursday', 'Friday', 'Saturday', 'Sunday'])
plt.title('Day of week CTR');
Figure 11
虽然周二和周三的点击数和点击数最高，但它们的点击率却是最低的。周六和周日享受最高的 CTR。显然，人们在周末有更多的时间点击。
C1 专题
C1 是匿名分类特征之一。虽然我们不知道它的意义，但我们还是想看看它的分布。
print(train.C1.value_counts()/len(train))
Figure 12
C1 值= 1005 拥有最多的数据，几乎是我们正在使用的所有数据的 92%。让我们看看 C1 的值是否能说明一些关于 CTR 的东西。
C1_values = train.C1.unique()
C1_values.sort()
ctr_avg_list=[]
for i in C1_values:
ctr_avg=train.loc[np.where((train.C1 == i))].click.mean()
ctr_avg_list.append(ctr_avg)
print("for C1 value: {}, click through rate: {}".format(i,ctr_avg))
Figure 13
train.groupby(['C1', 'click']).size().unstack().plot(kind='bar', figsize=(12,6), title='C1 histogram');
Figure 14
df_c1 = train[['C1','click']].groupby(['C1']).count().reset_index()
df_c1 = df_c1.rename(columns={'click': 'impressions'})
df_c1['clicks'] = df_click[['C1','click']].groupby(['C1']).count().reset_index()['click']
df_c1['CTR'] = df_c1['clicks']/df_c1['impressions']*100plt.figure(figsize=(12,6))
sns.barplot(y='CTR', x='C1', data=df_c1)
plt.title('CTR by C1');
Figure 15
重要的 C1 值和 CTR 对有:
C1=1005: 92%的数据和 0.17%的点击率
C1=1002: 5.5%的数据和 0.21%的点击率
C1=1010:数据的 2.2%和 0.095 的点击率
C1 = 1002 比平均 CTR 高得多，C1=1010 比平均 CTR 低得多，这两个 C1 值似乎对预测 CTR 很重要。
横幅位置
听说影响你横幅广告效果的因素有很多，但最有影响的还是横幅位置。让我们看看这是不是真的。
print(train.banner_pos.value_counts()/len(train))
Figure 16
banner_pos = train.banner_pos.unique()
banner_pos.sort()
ctr_avg_list=[]
for i in banner_pos:
ctr_avg=train.loc[np.where((train.banner_pos == i))].click.mean()
ctr_avg_list.append(ctr_avg)
print("for banner position: {}, click through rate: {}".format(i,ctr_avg))
Figure 17
重要的横幅位置有:
位置 0:数据的 72%和 0.16 CTR
位置 1: 28%的数据和 0.18%的 CTR
train.groupby(['banner_pos', 'click']).size().unstack().plot(kind='bar', figsize=(12,6), title='banner position histogram');
Figure 18
df_banner = train[['banner_pos','click']].groupby(['banner_pos']).count().reset_index()
df_banner = df_banner.rename(columns={'click': 'impressions'})
df_banner['clicks'] = df_click[['banner_pos','click']].groupby(['banner_pos']).count().reset_index()['click']
df_banner['CTR'] = df_banner['clicks']/df_banner['impressions']*100
sort_banners = df_banner.sort_values(by='CTR',ascending=False)['banner_pos'].tolist()
plt.figure(figsize=(12,6))
sns.barplot(y='CTR', x='banner_pos', data=df_banner, order=sort_banners)
plt.title('CTR by banner position');
Figure 19
虽然横幅位置 0 具有最高的印象和点击量，但横幅位置 7 享有最高的点击率。增加横幅位置 7 的广告数量似乎是一个好主意。
设备类型
print('The impressions by device types')
print((train.device_type.value_counts()/len(train)))
Figure 20
train[['device_type','click']].groupby(['device_type','click']).size().unstack().plot(kind='bar', title='device types');
Figure 21
设备类型 1 获得最多的展示和点击，其他设备类型仅获得最少的展示和点击。我们可能想要查看关于设备类型 1 的更多细节。
df_click[df_click['device_type']==1].groupby(['hour_of_day', 'click']).size().unstack().plot(kind='bar', title="Clicks from device type 1 by hour of day", figsize=(12,6));
Figure 22
正如所料，大多数点击发生在营业时间，来自设备类型 1。
device_type_click = df_click.groupby('device_type').agg({'click':'sum'}).reset_index()
device_type_impression = train.groupby('device_type').agg({'click':'count'}).reset_index().rename(columns={'click': 'impressions'})
merged_device_type = pd.merge(left = device_type_click , right = device_type_impression, how = 'inner', on = 'device_type')
merged_device_type['CTR'] = merged_device_type['click'] / merged_device_type['impressions']*100merged_device_type
Figure 23
最高的 CTR 来自设备类型 0。
用同样的方法，我探索了所有其他分类功能，如网站功能，应用程序功能和 C14 C21 功能。探索的方式大同小异，详细内容可以在 Github 上找到，在此不再赘述。
建筑模型
哈希简介
一个散列函数是一个将一组对象映射到一组整数的函数。当使用散列函数时，这种映射将任意长度的键作为输入，并输出特定范围内的整数。
我们的精简数据集仍然包含 100 万个样本和～2M 特征值。哈希的目的是最小化功能的内存消耗。
如果你想了解更多，这里有一篇由卢卡斯·伯纳蒂撰写的关于杂凑技巧的优秀文章。
Python 有一个内置函数，它执行一个名为hash()的散列。对于我们数据中的对象，散列并不奇怪。
def convert_obj_to_int(self):
object_list_columns = self.columns
object_list_dtypes = self.dtypes
new_col_suffix = '_int'
for index in range(0,len(object_list_columns)):
if object_list_dtypes[index] == object :
self[object_list_columns[index]+new_col_suffix] = self[object_list_columns[index]].map( lambda x: hash(x))
self.drop([object_list_columns[index]],inplace=True,axis=1)
return self
train = convert_obj_to_int(train)
LightGBM 模型
lightGBM_CTR.py
训练后的最终输出:
Figure 24
Xgboost 模型
Xgboost_CTR.py
它将训练到 eval-logloss 在 20 轮中没有改善。最后的输出是:
Figure 25
Jupyter 笔记本可以在 Github 上找到。周末愉快！
MobileNetV2:反向残差和线性瓶颈
2017 年 4 月，谷歌的一组研究人员发表了一篇论文，介绍了一种针对移动设备优化的神经网络架构。他们努力寻求一种模型，在保持参数和数学运算尽可能低的同时提供高精度。为了给智能手机带来深度神经网络，这是非常必要的。
被称为 MobileNet 的架构围绕着使用深度方向可分离卷积的思想，它由一个深度方向卷积和一个点方向卷积组成。如果你对这个操作的细节有点模糊，请随意查看我的另一篇文章，它详细解释了这个概念。 MobileNetV2 用两个主要想法扩展了它的前身。
反向残差
残余块用跳过连接来连接卷积块的开始和结束。通过添加这两种状态，网络有机会访问在卷积块中未被修改的早期激活。事实证明，这种方法对于构建深度网络至关重要。
A residual block connects wide layers with a skip connection while layers in between are narrow
当稍微靠近观察跳过连接时，我们注意到原始残差块遵循关于通道数量的宽- >窄- >宽方法。输入具有大量通道，这些通道通过廉价的 1x1 卷积进行压缩。这样，下面的 3×3 卷积的参数就少得多。为了最终增加输入和输出，使用另一个 1x1 卷积再次增加通道的数量。在喀拉斯，它看起来像这样:
def residual_block(x, squeeze=16, expand=64):
m = Conv2D(squeeze, (1,1), activation='relu')(x)
m = Conv2D(squeeze, (3,3), activation='relu')(m)
m = Conv2D(expand, (1,1), activation='relu')(m)
return Add()([m, x])
另一方面，MobileNetV2 遵循窄->宽->窄的方法。第一步使用 1×1 卷积来加宽网络，因为接下来的 3×3 深度方向卷积已经大大减少了参数的数量。之后，另一个 1x1 卷积挤压网络，以便匹配初始的信道数量。
An inverted residual block connects narrow layers with a skip connection while layers in between are wide
在喀拉斯，它看起来像这样:
def inverted_residual_block(x, expand=64, squeeze=16):
m = Conv2D(expand, (1,1), activation='relu')(x)
m = DepthwiseConv2D((3,3), activation='relu')(m)
m = Conv2D(squeeze, (1,1), activation='relu')(m)
return Add()([m, x])
作者将这种想法描述为反向剩余块，因为在网络的狭窄部分之间存在跳跃连接，这与原始剩余连接的工作方式相反。当您运行上面的两个代码片段时，您会注意到反向块的参数要少得多。
线性瓶颈
我们在神经网络中使用非线性激活函数的原因是多个矩阵乘法不能简化为单个数值运算。它允许我们建立多层神经网络。同时，通常在神经网络中使用的激活函数 ReLU 丢弃小于 0 的值。这种信息丢失可以通过增加信道数量来解决，以便增加网络的容量。
对于反向残差块，我们做相反的事情，并挤压跳跃连接所链接的层。这损害了网络的性能。作者引入了线性瓶颈的概念，其中剩余块的最后卷积在被添加到初始激活之前具有线性输出。将它写入代码非常简单，因为我们只需丢弃卷积模块的最后一个激活函数:
def inverted_linear_residual_block(x, expand=64, squeeze=16):
m = Conv2D(expand, (1,1), activation='relu')(x)
m = DepthwiseConv2D((3,3), activation='relu')(m)
m = Conv2D(squeeze, (1,1))(m)
return Add()([m, x])
ReLU6
上面的片段显示了一个卷积块的结构，它包含了反向残差和线性瓶颈。如果您想尽可能地匹配 MobileNetV2，您还需要另外两个部分。第一个方面只是在每个卷积层后面增加了批量归一化，这是你现在可能已经习惯的。
第二个附加项不太常见。作者使用 ReLU6 而不是 ReLU，这将激活值限制在最大值…嗯…6。只要在 0 到 6 之间，激活就是线性的。
def relu(x):
return max(0, x)def relu6(x):
return min(max(0, x), 6)
这在你处理定点推理时很有帮助。它将小数点左边的信息限制为 3 位，这意味着我们可以保证小数点右边的精度。这也在最初的 MobileNet 论文中使用。最后一个构建块看起来像这样:
def bottleneck_block(x, expand=64, squeeze=16):
m = Conv2D(expand, (1,1))(x)
m = BatchNormalization()(m)
m = Activation('relu6')(m)
m = DepthwiseConv2D((3,3))(m)
m = BatchNormalization()(m)
m = Activation('relu6')(m)
m = Conv2D(squeeze, (1,1))(m)
m = BatchNormalization()(m)
return Add()([m, x])
建筑
The MobileNetV2 architecture
现在我们已经了解了 MobileNetV2 的构建块，我们可以看一下整个架构。在表格中，您可以看到瓶颈块是如何排列的。 t 代表通道的膨胀率。如您所见，他们使用了系数 6，而不是我们示例中的系数 4。 c 表示输入通道的数量， n 表示该块重复的频率。最后, s 告诉我们块的第一次重复是否为下采样过程使用了步幅 2。总而言之，这是一个非常简单和常见的卷积块集合。
Performance of MobileNetV2 and other architectures on ImageNet
总结想法
我特别高兴的是 MobileNetV2 提供了与 NASNet 类似的参数效率。NASNet 是几项图像识别任务的最新技术。它的构建模块相当复杂，这使得它为什么工作得这么好变得相当不直观。NASNet 的构建模块不是由人类设计的，而是由另一个神经网络设计的。引入一个简单的架构，比如 MobileNetV2，它表现出了相当的效率，这让我更加相信下一个大的架构可能也是由人类设计的。
使用 Mockito 模仿同一个测试类中的方法
在测试驱动开发(TDD)中，单元测试是隐含实现质量的子部分。在使用 junit 进行单元测试时，您会遇到想要模拟类的地方。当您调用具有外部通信(如数据库调用或 rest 调用)的类的方法时，就完成了模仿。通过模仿，您可以显式地定义方法的返回值，而无需实际执行方法的步骤。在这篇文章中，我将讨论在你编写测试用例的同一个测试类中模仿方法。
假设您有一个 Person 类，它有外部通信并相应地返回值。
public class Person {
private String name;
private int age; public Person(String name, int age) {
this.name = name;
this.age = age;
}
public String getName() {
return name;
}
public void setName(String name) {
this.name = name;
}
public int getAge() {
return age;
}
public void setAge(int age) {
this.age = age;
}
public boolean runInGround(String location) {
if(location.equals("ground")) {
System.*out*.println("The person runs in the " + location);
return true;
} else {
System.*out*.println("The person doesn't run in the " + location);
return false;
}
}
public boolean isPlay() {
if(this.runInGround("ground")) {
System.*out*.println("The person plays.");
return true;
}
else {
System.*out*.println("The person doesn't play");
return false;
}
}
}
出于解释的目的，让我们假设 runInGround(字符串位置)方法通过检查数据库值来返回值，这是一个外部通信。这个特殊的方法是在同一个 Person 类的 isPlay()方法中调用的，它影响 isPlay()方法的返回值。
那么如何在不让 runInGround(字符串位置)执行并与数据库对话的情况下测试 isPlay()。
我们可以模仿 PersonTest 类中的 *runInGround(字符串位置)*方法，如下所示。
import org.junit.Assert;
import org.junit.Test;
import org.mockito.Mockito;
public class PersonTest{
@Test
public void playTest() {
Person person = new Person("name", 15, "23435678V");
Person person1 = Mockito.*spy*(person);
Mockito.*doReturn*(true).when(person1).runInGround("ground");
Assert.*assertEquals*(true, person1.isPlay());
}
}
这里我们需要使用 Mockito.spy()来模拟我们正在测试的同一个类，而不是使用 mock(class)。然后我们可以如下模拟我们想要的方法。
Mockito.*doReturn*(true).when(person1).runInGround("ground");
希望这将有所帮助。用单元测试快乐编码:d。
延伸阅读
[1]http://www.vogella.com/tutorials/Mockito/article.html
https://dzone.com/articles/introduction-to-java-tdd
模态测试和核密度估计
python 中的多模态测试
当处理大量可能具有不同数据分布的数据集时，我们面临以下考虑因素:
- 数据分布是单峰的吗？如果是，哪种模型最接近它(均匀分布、T 分布、卡方分布、柯西分布等)？
- 如果数据分布是多模态的，我们能否自动识别模态的数量并提供更细粒度的描述性统计？
- 如何估计一个新数据集的概率密度函数？
本笔记本处理以下主题:
- 直方图与概率密度函数近似
- 核密度估计
- 最佳带宽的选择:Silverman/ Scott/ Grid 搜索交叉验证
- 单峰分布的统计检验
- 单峰性倾斜试验
- 基于核密度估计的数据分布模式数识别
直方图和 pdf
正如这篇博文中所解释的，https://mgl Lerner . github . io/posts/histograms-and-kernel-density-estimation-kde-2 . html直方图的缺点是在不合适大小的容器中隐藏了实际数据分布的一些细节。
def plotHistogramAndPdf(data, x, pdf):
ax = plt.gca()
plt.hist(data, bins = 4, alpha = 0.4, label = 'histogram of input values');
plt.ylabel('Frequency')
plt.xlabel('x values')
ax2 = ax.twinx()
plt.plot(x, pdf, c = 'red', label = 'probability density function');
plt.ylabel('PDF')
[tl.set_color('r') for tl in ax2.get_yticklabels()]
ax.legend(bbox_to_anchor=(0.4, 1.15))
ax2.legend(bbox_to_anchor=(1.15,1.15))
plt.savefig('figures/hist.jpg', bbox_inches='tight')
plotHistogramAndPdf(data, x, true_pdf)
核密度估计
核密度估计依赖于任意带宽，该带宽决定了返回的近似的平滑程度。以下示例说明了各种带宽值的影响:
def getKernelDensityEstimation(values, x, bandwidth = 0.2, kernel = 'gaussian'):
model = KernelDensity(kernel = kernel, bandwidth=bandwidth)
model.fit(values[:, np.newaxis])
log_density = model.score_samples(x[:, np.newaxis])
return np.exp(log_density)for bandwidth in np.linspace(0.2, 3, 3):
kde = getKernelDensityEstimation(data, x, bandwidth=bandwidth)
plt.plot(x, kde, alpha = 0.8, label = f'bandwidth = {round(bandwidth, 2)}')
plt.plot(x, true_pdf, label = 'True PDF')
plt.legend()
plt.title('Effect of various bandwidth values \nThe larger the bandwidth, the smoother the approximation becomes');
plt.savefig('figures/bw.jpg', bbox_inches='tight')
核密度估计最佳带宽的选择方法
为了确定最佳带宽，有几种方法:
- Silverman 的经验法则:假设未知密度为高斯分布。它不是最佳带宽选择器，但可以用作非常快速、相当好的估计器，或者用作多级带宽选择器中的第一估计器。更精确的求解方程插件规则使用积分平方密度导数泛函的估计来估计最佳带宽。它们需要大量计算来使用迭代方法求解非线性方程。他们用腐烂作为第一估计
- Scott 的经验法则:对于正态分布数据的随机样本是最佳的，在某种意义上，它最小化了密度估计的积分均方误差。
这两种方法具有计算速度快的优点，但是它们通常给出太少的面元，并且很可能对底层数据分布进行欠拟合。这两种方法都已经在 statsmodels 包中实现，如下图所示。
- 基于交叉验证的方法:statsmodels 带有一个 cv 带宽参数。或者，我们可以实现网格搜索交叉验证。与前两种方法不同，执行网格搜索可能需要更多的计算，尤其是对于较大的数据集
from statsmodels.nonparametric.bandwidths import bw_silverman, bw_scott, select_bandwidthsilverman_bandwidth = bw_silverman(data)# select bandwidth allows to set a different kernel
silverman_bandwidth_gauss = select_bandwidth(data, bw = 'silverman', kernel = 'gauss')scott_bandwidth = bw_scott(data)def bestBandwidth(data, minBandwidth = 0.1, maxBandwidth = 2, nb_bandwidths = 30, cv = 30):
"""
Run a cross validation grid search to identify the optimal bandwidth for the kernel density
estimation.
"""
from sklearn.model_selection import GridSearchCV
model = GridSearchCV(KernelDensity(),
{'bandwidth': np.linspace(minBandwidth, maxBandwidth, nb_bandwidths)}, cv=cv)
model.fit(data[:, None])
return model.best_params_['bandwidth']cv_bandwidth = bestBandwidth(data)print(f"Silverman bandwidth = {silverman_bandwidth}")
print(f"Scott bandwidth = {scott_bandwidth}")
print(f"CV bandwidth = {cv_bandwidth}")
正如预期的那样，第一个 Silverman 和 Scott 返回了更大的带宽值，这导致了更大的箱，从而丢失了关于数据分布的信息。
Statsmodels 允许基于交叉验证和最大似然运算符自动搜索最佳带宽:
from statsmodels.nonparametric.kernel_density import KDEMultivariate
stats_models_cv = KDEMultivariate(data, 'c', bw = 'cv_ml').pdf(x)
画出不同的近似值
plt.figure(figsize= (14, 6))
plt.plot(x, true_pdf, label = 'True PDF')kde = getKernelDensityEstimation(data, x, bandwidth=silverman_bandwidth)
plt.plot(x, kde, alpha = 0.8, label = f'Silverman bandwidth')kde = getKernelDensityEstimation(data, x, bandwidth=scott_bandwidth)
plt.plot(x, kde, alpha = 0.8, label = f'Scott bandwidth')kde = getKernelDensityEstimation(data, x, bandwidth=cv_bandwidth)
plt.plot(x, kde, alpha = 0.8, label = f'CV bandwidth')plt.plot(x, stats_models_cv, alpha = 0.8, label = f'Statsmodels CV maximum likelihood')plt.legend()
plt.title('Comparative of various bandwidth estimations for KDE');
plt.savefig('figures/comp_bw.jpg', bbox_inches='tight')
单峰分布的统计检验
有许多统计测试可以解决数据形态问题:
- 倾斜试验
- 过量质量测试
- 地图测试
- 模式存在测试
- 矮小测试
- 跨度测试
- 马鞍试验
不幸的是，在 python 开源库中实现的并不多。
倾斜试验
以下 python 包https://github.com/BenjaminDoran/unidip提供了倾角测试的实现，以及利用单峰性的哈迪根倾角测试在数据中广泛提取密度峰值的功能。
from unidip import UniDip
import unidip.dip as dipdata = np.msort(data)
print(dip.diptst(data))
intervals = UniDip(data).run()
print(intervals)
确定并绘制 KDE 的局部最大值
一旦我们有了核密度函数的估计，我们就可以确定该分布是否是多峰的，并识别对应于这些模式的最大值或峰值。
这可以通过识别一阶导数改变符号的点来实现。默认情况下，getinflexinpoints 方法可以返回所有拐点(最小值+最大值)，或者只返回一个选择(typeOfInflexion = ‘max’/ 'min ')。
下图描绘了可能对应于多种数据分布模式的最大值。可以通过基于峰的高度设置阈值来继续分析，以便过滤掉一些不太重要的值。
def getExtremePoints(data, typeOfExtreme = None, maxPoints = None):
"""
This method returns the indeces where there is a change in the trend of the input series.
typeOfExtreme = None returns all extreme points, max only maximum values and min
only min,
"""
a = np.diff(data)
asign = np.sign(a)
signchange = ((np.roll(asign, 1) - asign) != 0).astype(int)
idx = np.where(signchange ==1)[0]if typeOfInflexion == 'max' and data[idx[0]] < data[idx[1]]:
idx = idx[1:][::2]
elif typeOfInflexion == 'min' and data[idx[0]] > data[idx[1]]:
idx = idx[1:][::2]
elif typeOfInflexion is not None:
idx = idx[::2]
# sort ids by min value
if 0 in idx:
idx = np.delete(idx, 0)
if (len(data)-1) in idx:
idx = np.delete(idx, len(data)-1)
idx = idx[np.argsort(data[idx])]
# If we have maxpoints we want to make sure the timeseries has a cutpoint
# in each segment, not all on a small interval
if maxPoints is not None:
idx= idx[:maxPoints]
if len(idx) < maxPoints:
return (np.arange(maxPoints) + 1) * (len(data)//(maxPoints + 1))
return idx
我们注意到获得的值对应于生成的分布的初始锚。
资源
- https://jakevdp . github . io/blog/2013/12/01/kernel-density-estimation/
- https://mglerner . github . io/posts/histograms-and-kernel-density-estimation-kde-2 . html
- https://en.wikipedia.org/wiki/Multimodal_distribution
Titanic 数据集上的模型比较
在为期 12 周的训练营中，我刚刚接触了数据科学的服务，我最喜欢的练习之一是查看模型比较。开始模型选择和 EDA 需要深思熟虑，并为我提供任何合理的见解执行。随着我成为一名数据科学家，我认为划分哪些模型最适合哪些类型的数据非常重要。这需要我自己进行大量的分析和探索。在这一点上，让我们看看泰坦尼克号数据集作为一个例子。
我应该提到，我正试图根据上面列出的几个因素来确定谁幸存了下来。在我深入研究之前，我想通过上面的表格指出，我们可以把性变成一个虚拟变量。接下来，我们应该执行和 EDA 定义我们想要的 X 和 y。
现在，我们应该导入并准备一个训练测试分割，并创建一个名为模型评估的功能，该功能查看准确性分数、混淆矩阵和分类报告。
现在让我们从 K 个最近的邻居开始。
现在 KNN 使用网格搜索。
这给了我们一个 0.636363636363635 的最终分数。现在让我们试试 KNN 装袋！
现在，逻辑回归
现在我们将运行一个决策树，但首先我们要运行一个网格搜索。
在这之后，我们应该在 DT 上使用装袋。
最后是随机森林和多余的树。
现在我们可以开始比较我们的模型了！首先，我们可以使用训练/测试分割来确定哪个模型执行得最好:
代替训练测试分割，我们可以看一个分层的 K 折叠来看模型在那里如何排列。
在这两方面，我认为决策树上的网格搜索包做得最好。使用训练测试分割，它排名很高，虽然它接近分层 k-fold 的下端，但它的误差比一些更高的预成型模型低得多！
通过师生知识蒸馏介绍 PyTorch 模型压缩
通过知识蒸馏的模型压缩可以节省推理时间、功率效率和模型大小。
Knowledge River Delta
在资源受限的移动和实时系统中服务 ML 模型可能是一个真正的问题。ML 社区一直在开发解决方案来压缩由较大的服务器集群生成的模型的大小。模型压缩承诺节省推理时间、功率效率和模型大小。所有这些都可以让飞行救援无人机在一次充电后覆盖更多的土地表面，同时不会耗尽移动应用用户的电池。
模型知识提炼是一种在不损失太多预测能力的情况下减少模型规模的方法。
Geoffrey Hinton 在 2018 年深度学习峰会上关于使用 【知识蒸馏】 (KD)的演讲让我去寻找另一类问题的当前技术状态:推荐系统(RecSys)。
这让我想到了唐嘉熙在 2018 年 KDD 发表的关于 排名蒸馏 (RD)的优秀作品，他在其中讨论了他将知识蒸馏应用于排名任务的相关方法。
在这篇博客中，我在 Movielens 100K 数据集 上复制了这个排名提炼工作的一小部分。在这方面工作是一种领悟。即使 KD 是一个从一个模型到一个小模型提取知识的坚实的概念框架，将它应用于推荐系统的排序任务也不是一个简单的任务。
第一个挑战是我们在一个比常见的 fit/predict API 更低的抽象层次上工作，这些 API 存在于 Scikit-learn 和 Keras 等更高层次的库中。这是因为实现这个 KD 所需的改变是在损失函数公式本身。为了解决这个问题，我跟随第三篇论文的脚步，使用优雅的 PyTorch API 在 RecSys 中构建这个 KD。
第二个挑战是，即使 PyTorch 是一个优雅的库，我们也需要一个更高层次的框架，专门研究带有 PyTorch 的 RecSys。这些天的选择框架似乎是来自库拉的聚焦。我强烈推荐它，API 设计易于使用，它让用户自定义我们这个实验需要的大多数方面。
我们走吧！
定制分级蒸馏反向传播流
目标是从 Movielens 100K 数据集生成 3 个模型:学生模型、带蒸馏的学生模型和教师模型，并比较它们的 MAP@K 指标以及物理磁盘大小。
我们需要解释我们将要使用的策略，从教师模型到学生模型，用蒸馏来教授一些黑暗知识。以下是对培训过程中将要发生的事情的解释:
Flow of data and forward/backward propagation during the knowledge distillation method
在上图中，我们展示了培训流程:
- 对于学生模型，我们使用传统的方法，使用带有数据标签和单个排名损失的训练数据。
- 对于教师模型，我们与学生模型类似地对其进行预训练，但我们使用更大的网络规模来实现更高的 K (MAP@K)平均精度。在完成较大模型的训练之后，我们存储预先训练的教师模型。
- 对于带有蒸馏的学生模型，我们使用带有标签和排名损失的训练数据。然而，在这个例子中，我们使用了教师模型对我们提供给学生模型的数据的预测。更准确地说，除了学生的损失之外，我们还使用教师的损失来计算和反向传播学生模型网络中的梯度。这些额外的信息应该可以提高学生模型的预测能力，同时保持模型大小与没有经过提取的学生模型相同。
结果比较
首先，我们需要一些培训数据，我们用这些数据来建立一个预培训教师模型。我们使用 movielens 100K 数据集，并且只使用电影/用户交互。我们将尝试预测用户最有可能评价的前 5 部电影。
为此，我们将使用 Spotlight 库提供的implicit factorization model。该模型使用基于嵌入的模型结构:
Implicit Factorization Model with a Bi-Linear model structure
对于损失，我们使用类似于下面负对数似然函数的方法。我们对正对和负对进行采样，并要求优化器提高正对中的排序项目(d+)并减少负对中的项目(d-):
Loss function related to the negative logarithmic of the likelihood function.
使用 200 作为 movielens 数据集上每个嵌入层的大小来训练“大型”教学模型，这为我们提供了以下指标:
让我们用一个小得多的模型做同样的尝试，用 2 作为每个嵌入层的大小:
2 次观察:
- 第一，学生模型本身序列化后的大小更小(0.10 mb vs 6.34)。这与网络的大小一致，因为嵌入的大小要小 100 倍。
- 第三，学生模型的 MAP@5 低于教师模型(0.050 vs 0.073)。较小的网络可以远离较大的网络。挑战是:我们能在保持模型尺寸不变的情况下做得更好吗？
这是我们接下来要尝试的。我们培训第三个模型，这是一个学生模型，由预先培训的教师模型推动。
为此，我们需要在损失函数中混合从两个模型中获得的两个损失。这就是 PyTorch 闪耀的地方。我们所要做的就是定义一个修正的损失函数，它总结了学生和老师的损失，让梯度下降发挥它的魔力。其核心是，如果您对使用对数 sigmoid 损失的正负损失有所了解，我们通过教师网络传递当前一批数据，获得候选预测，并使用它们来生成教师损失值。我们用于优化的最终损失是 pos/neg/teacher 这三个损失的总和。下面是组合损失函数的一个片段:
我们表现如何？
- 首先，提取模型的 MAP@5 值更接近教师模型的值，仅使用 2 作为嵌入层的大小(0.070 比 0.073)
- 第二，大小仍然是 0.10mb，类似于非蒸馏的学生模型。
这里有一个表格，列出了所有这些值以供比较
我们从这次冒险中学到了什么？
- 我很高兴看到 PyTorch 能够灵活地复制 KDD2018 论文的一小部分。
- 知识蒸馏真的很酷，也为推荐系统工作。
- 总的来说，每当两个或更多的人工智能模型之间有互动时，我对它们的结果非常感兴趣。
如果您对这种类型的跨模型交互感兴趣，我邀请您深入研究 KDD2018 论文。我们没有讨论如何通过加权教师的模型损失或仅考虑教师模型的 top-k 建议来改进这种设置。大概是以后的帖子吧。
直到那时！
谢了。
参考资料:
唐佳夕，还有王柯。排名提取:为推荐系统学习高性能的紧凑排名模型。第 24 届 ACM SIGKDD 知识发现国际会议论文集&数据挖掘。ACM，2018。
辛顿·杰弗里、奥里奥尔·维尼亚尔斯和杰夫·迪恩。**在神经网络中提取知识。**arXiv 预印本 arXiv:1503.02531 (2015)。
2017 年【https://github.com/maciejkula/spotlight】库拉聚光灯T2
2018 年 PyTorch ，https://pytorch.org/
模型评估 I:精度和召回率
source - pexels.com
为了测试像支持向量机这样的分类系统的质量，需要执行一些评估指标。支持向量机是我在内核中简单解释过的分类算法。
一点背景
精确和回忆在信息抽取中都有广泛的应用。Precision 是检索到的*的 文档的编号，Recall 是检索到的的 相关文档的编号。***
相关性是指所呈现的信息对此刻正在讨论的主题有多有用。
让我们来看一个约翰和乔希之间对话的例子:
如果约翰说“我喜欢冰淇淋”，而乔希回答“我有一个朋友叫约翰·多伊”，那么乔希刚才说的与约翰暗示的毫无关系，但是如果乔希说“我有一个朋友叫约翰·多伊，他也喜欢冰淇淋”，那么乔希的陈述就变得相关，因为它现在与约翰的陈述相关。
精确度和召回率都非常有助于理解呈现了哪组文档或信息，以及这些文档中有多少对于所提问题是实际上有用的。
尽管精确度和召回率有时会被混淆为同义词，但它们不是。
精确度和召回率彼此成反比，因此理解它们的差异对于建立有效的分类系统是重要的。
让我们考虑另一个例子
假设我在谷歌上搜索“什么是精确和召回？”不到一分钟，我就有大约 1560 万个结果。
假设在这 1560 万个结果中，与我的问题相关的链接大约有 200 万个。假设还有大约 600 多万个相关但没有被谷歌返回的结果，对于这样的系统，我们可以说它的精度是 2M/1560 万，召回率是 2M/8 万。
这意味着谷歌的算法检索* e 所有 相关 链接 的概率是 0.25 ( 召回)，所有 检索到的链接 是 相关的概率是 0.13 ( 精度)。*
思考精度和回忆的另一种方式是这样的:
如果有人让你列出你去年圣诞节收到的 5 件礼物的名字，但你记不清这 5 件礼物的名字，所以你随机猜了 7 次。在你记得的 7 个名字中，有 5 个被正确回忆起来，而 2 个是你生日时收到的礼物。即使你有 100%的召回率(5/5)，你的准确率也只有 71.4% (5/7)。
现在，我已经让您对精确和回忆有了初步的了解，让我们再深入一点。
第一类和第二类错误
谈精度和回忆而不提I 型& II 型错误 就像叙述穆罕默德·阿里的历史而跳过“拳击”——我
ouch! that punch was precise 😃 — source here
第一类错误
这是一个真 零假设(Ho)的不正确拒绝。
零假设是一个默认为真的陈述，除非被证明是相反的。
*I 型错误导致假阳性 (FP) *。例如，当一个文件被搜索引擎返回为“相关”时，结果却是“不相关”。
一个例子是，实际上并没有发生火灾，但火警却响了。这种错误等同于“相信谎言”或“虚惊一场”。
第二类错误
这是一个假 零假设(Ho)的不正确保留。
这相当于系统忽略了可能没有检索到相关文档的可能性。
这种错误会导致假阴性 (FN)。即没有检索到本应检索到的相关文档。
这种情况的一个例子是发生了火灾，但火警没有响。这种错误等同于“不相信真理”或“失误”。
考虑下表:
Table 1.0 (confusion matrix)
误报和漏报分别是精确度和召回率的两个独特特征。
Source here
减少一个意味着增加另一个
P α 1/R
一些数学公式
在分类任务中，
*Precision P = TP/(TP+ FP)*
即{真阳性(正确检索的文档)数)}/{(检索的文档总数)}
*Recall R = TP/(TP + FN)*
即{真阳性(正确检索的文档)数)}/{(检索的相关文档总数)}
从 Google 搜索的例子来看， perfect Precision 得分为 1.0 意味着搜索引擎检索到的每个结果都是相关的(但没有说明是否检索到了所有相关的文档)
而完美回忆得分为 1.0 意味着所有相关文档都从搜索引擎中检索到了(但没有说明检索到的结果中有多少是不相关的)。
这不是一个非此即彼的问题，因为就一个问题作出的决定会自动影响另一个问题。因此，对于每个系统，通常都有一个基准，它被认为是“可接受的”,而不会失去太多的特性。
对于能够实现最大精确度(无假阳性)和最大召回率(无假阴性)的任何系统，需要没有类型 I 和类型 II 错误。
精确度和召回分数不是孤立讨论的。取而代之的是，将一个度量的值与另一个度量的固定水平进行比较(例如，召回水平为 0.75 的精度)，或者将两者合并为一个度量。
结合了精确度和召回率的度量的例子是 F-measure
*F = 2 * ((Precision * Recall)/(Precision + Recall))*
这个指标通常被称为 F1 得分，它是精确度和召回率接近时的平均值。这是调和平均数。
其他相关指标有准确度、混淆矩阵和列联表。
继续探索吧！
相关链接
- https://en.wikipedia.org/wiki/Null_hypothesis
- https://en . Wikipedia . org/wiki/Precision _ and _ recall # F-measure
- https://en . Wikipedia . org/wiki/Accuracy _ and _ precision # In _ binary _ class ification
- https://en.wikipedia.org/wiki/Confusion_matrix
- https://en.wikipedia.org/wiki/Contingency_table
从头开始创建模型概述和代码
我最近在一次聚会上与一位数据科学家交谈，他提到他发现从头构建模型在智力上是有益的。注意到这让他对模型的行为有了很好的直觉，也是一个很好的 OOP 项目。我发现，通过在这个项目上花几天时间，我了解了每个模型的新特征，在堆叠、追加和整形 numpy 数组方面变得更加流畅，并且提高了我对优化算法以及使用 numpy 与标准 python 操作的效果的认识。
此外，我最近在深度学习中阅读和学习了大量新概念，我想后退一步，刷新我对常见分类算法的理解。
我决定采用 3 种难度递增的分类模型和一种额外的聚类算法:
K-最近邻
kNN 是机器学习中使用最广泛的分类模型之一，也是最容易理解的分类模型之一，因为它的决策过程是透明的。它被认为是一个懒惰的学习者，因为它根本不建立模型——记忆的训练数据集就是模型！为了进行样本外预测，kNN 找到 k 个最相似的训练样本，并对与它们相关联的结果类进行计数。并且通过多数投票或其他投票方案，算法预测未知样本的类别。
正如我提到的，kNN 是一个懒惰的学习者，因为它直到最后一秒——在预测时间——才做任何事情。这有几个影响。首先，这意味着 kNN 需要在预测时间内计算所有训练样本与每个未知点之间的相似度。这通常会导致更大的内存需求，并且可能比构建表示的其他模型更慢，计算量也更大。然而，如果结构适当，这可能是有利的。例如，如果新的训练样本不断生成，旧的样本被删除，那么 kNN 的预测将来自手头“最新”或“最相关”的数据。一个重要的注意事项是，kNN 假设特征被缩放。这是确保在确定相似性时所有特征被平均加权的要求。
超参数是 k 的值、相似性度量和加权度量。k 的值是预测时用于投票的邻居数量。它介于 1 和 n 之间，其中 n 是训练样本的总数。
**def** find_neighbors(self, new_sample):
*'''List the k neighbors closest to the new sample.*
*'''*
distances = []
**for** i **in** range(len(self.X)):
**if** self.metric == 'euclidean':
distance = self.euclidean(self.X[i], new_sample)
**if** self.metric == 'manhattan':
distance = self.manhattan(self.X[i], new_sample)
distances.append((self.y[i],distance))
distances = sorted(distances,key=operator.itemgetter(1))
neighbors = []
**for** i **in** range(self.k):
neighbors.append(distances[i][0])
**return** neighbors
相似性度量是根据距离量化未知样本和训练样本之间的关系的数学函数。您可以想象，相似点的排序以及反过来的分类严重依赖于所选择的度量。有许多相似性度量，每一个都有自己的用例。通常，决策基于数据集中的数据类型和一些经验法则。
默认的数字度量是欧几里德距离。它是两点之间的“直线距离”,由平方差之和定义。另一种方法是曼哈顿距离，这是从西 22 街和第三大道走到西 27 街和第五大道所需的距离，它受轴的限制，在数学上由绝对差之和定义。此外，分类数据和二进制数据通常使用汉明距离进行评估，这基本上类似于逻辑 and。它计算点之间相同属性的数量。每个度量标准都会产生不同的点排序，基于这些点它认为与未知点最相似。
def **euclidean**(self, a, b):
return np.sqrt(((a-b)**2).sum(axis=0))
恰好有一个距离度量来统治它们。它叫做闵可夫斯基距离。我的意思是，当 p = 2 时，它有能力表示欧几里得，当 p = 1 时，它有能力表示曼哈顿，当 p 接近 0 时，它有能力表示近汉明尔。因此，这是您在现成的 kNNs 中看到的最常见的距离度量。
def **minkowski**(self, a, b, p):
return ((np.abs(a-b)**p).sum(axis=0))**(1/p)
加权超参数旨在基于 k 个邻居在相似性线中的位置来修改对每个 k 个邻居投票的影响。它通常设置为多数投票，这是我在实现中硬编码的，但倒数或等级距离也是其他典型的选项。
**def** majority_vote(self, neighbors):
*'''Determine majority class from the set of neighbors.*
*'''*
class_votes = {}
**for** i **in** range(len(neighbors)):
sample_class = neighbors[i]
**if** sample_class **in** class_votes:
class_votes[sample_class] += 1
**else**:
class_votes[sample_class] = 1
sorted_votes = sorted(class_votes.items())
**if** self.ties:
sorted_votes = self.tie(sorted_votes)
**return** sorted_votes[0][0]
kNN 培训:
1.记忆数据
kNN 预测:
- 接受新样本
- 使用指定的度量计算所有训练点和新样本之间的距离
- 基于相似性对所有训练点进行排序。距离越小，相似度越高
- 将 k 个最近点定义为排序列表中的前 k 个点
- 基于 k 个最近点类的加权集合计算预测值
**def** predict(self, X_test):
*'''Predict class for each value in array of new samples.*
*'''*
self.tie_count = 0
y_pred = []
**for** i **in** range(len(X_test)):
neighbors = self.find_neighbors(X_test[i])
pred_class = self.majority_vote(neighbors)
y_pred.append(pred_class)
**if** self.ties:
print('**{}** ties'.format(self.tie_count))
**return** y_pred
虽然很少，但我也添加了一个警告:当多数投票产生平局时，平局将出现在训练数据中出现频率最高的类中。由于 Iris 数据集很小，这种情况经常发生，但在较大的数据集中就不太常见了。
**def** tie(self,sorted_votes):
*'''Determine when ties occur in the the neighbors. Of the tied classes,*
*choose the class most frequent in the training data.*
*Print out number of ties.*
*'''*
tie = {}
**for** pair **in** sorted_votes:
count = pair[1]
**if** count **in** tie:
self.tie_count += 1
*#print('tie')*
tie[count].append(pair[0])
**else**:
tie[count] = [pair[0]]
*#print(tie)*
tie_class_frequency = {}
**if** len(tie[count]) > 1:
*#print('tie')*
**for** tie_class **in** tie[count]:
tie_class_frequency[tie_class] = np.count_nonzero(self.y == tie_class)
max_class = max(tie_class_frequency, key=tie_class_frequency.get)
*#print(max_class)*
sorted_votes = [(max_class,1)]
**return** sorted_votes
这个练习对我来说是一个很好的方式，可以深入到 kNN 的各个层次。这是从零开始构建更多模型的一个简单的开始，给了我一个解决问题的良好基础。了解从相似性度量到加权选项的众多选项有助于我构建类，以便在我决定做更多事情时允许增量添加。
决策树
决策树是一种流行的机器学习模型，因为它们易于解释并且具有预测能力。然而，它们本身并不常用。相反，它们被用作被认为是最先进模型的集成方法的基础。决策树是通过递归地将数据集一分为二来构建的，直到分支是纯的或者满足停止标准。最后的节点称为叶节点，对应于一个输出类。其余节点表示任意输入要素和分割数据集的要素的值。为了预测新样本的类别，样本简单地沿着分裂逻辑一个节点一个节点地向下遍历树，直到到达叶节点。
**def** split(self, feat, val, Xy):
Xi_left = np.array([]).reshape(0,self.Xy.shape[1])
Xi_right = np.array([]).reshape(0,self.Xy.shape[1])
**for** i **in** Xy:
*#print(i.shape)*
**if** i[feat] <= val:
Xi_left = np.vstack((Xi_left,i))
**if** i[feat] > val:
Xi_right = np.vstack((Xi_right,i))
**return** Xi_left, Xi_right
树是贪婪的算法，这意味着在每次分裂时，做出最佳决策，希望得到的树将导致全局最大值。然而，这很少偶然发生，通常他们需要许多正则化超参数来减少过拟合。上面提到的停止标准就是这样的例子。最大深度是控制树可以向下生长的总步数的参数，最小样本数是控制发生分裂所需的样本数的参数。通过减少最大深度或增加最小样本数，可以控制树的大小。
任何给定特征的最佳分离由定义纯度的成本函数确定。最常见的是用于连续特征的基尼指数和用于分类特征的信息增益或熵。基尼指数通过所得两组的阶级构成来量化分裂，其中低分意味着接近完美的分裂，高分意味着两组的构成几乎相同。
**def** gini_score(groups,classes):
n_samples = sum([len(group) **for** group **in** groups])
gini = 0
**for** group **in** groups:
size = float(len(group))
**if** size == 0:
**continue**
score = 0.0
*#print(size)*
**for** class_val **in** classes:
*#print(group.shape)*
p = (group[:,-1] == class_val).sum() / size
score += p * p
gini += (1.0 - score) * (size / n_samples)
*#print(gini)*
**return** gini
为了计算最佳分割，模型必须以各种可能的方式计算每个特征分割的基尼系数。例如，以各种可能的方式分割特征[x]上的样本[a，b，c，d，e]将导致[a]和[b，c，d，e]，[a，b]和[c，d，e]，[a，b，c]和[d，e]以及[a，b，c，d]和[e]。将对每个分割进行评分，并保存最低的基尼系数。那么对特征[y]和[z]也是如此。然后，所有特征中最低的基尼系数将被选为最佳分割，并创建一个节点。对每个节点的所得两组中的每一组重复该过程，直到遇到点点。
**def** best_split(self, Xy):
classes = np.unique(Xy[:,-1])
best_feat = 999
best_val = 999
best_score = 999
best_groups = **None**
**for** feat **in** range(Xy.shape[1]-1):
**for** i **in** Xy:
groups = self.split(feat, i[feat], Xy)
gini = self.gini_score(groups, classes)
*#print('feat {}, valued < {}, scored {}'.format(feat,i[feat], gini))*
**if** gini < best_score:
best_feat = feat
best_val = i[feat]
best_score = gini
best_groups = groups
output = {}
output['feat'] = best_feat
output['val'] = best_val
output['groups'] = best_groups
**return** output
这个结构是通过使用递归函数快速创建的。只需评估基尼系数，在最佳分割点创建节点，保存结果数据集，然后重复。
**def** terminal_node(self, group):
classes, counts = np.unique(group[:,-1],return_counts=**True**)
**return** classes[np.argmax(counts)]
**def** split_branch(self, node, depth):
left_node, right_node = node['groups']
**del**(node['groups'])
**if** **not** isinstance(left_node,np.ndarray) **or** **not** isinstance(right_node,np.ndarray):
node['left'] = node['right'] = self.terminal_node(left_node + right_node)
**return**
**if** depth >= self.max_depth:
node['left'] = self.terminal_node(left_node)
node['right'] = self.terminal_node(right_node)
**return**
**if** len(left_node) <= self.min_num_sample:
node['left'] = self.terminal_node(left_node)
**else**:
node['left'] = self.best_split(left_node)
self.split_branch(node['left'], depth+1)
**if** len(right_node) <= self.min_num_sample:
node['right'] = self.terminal_node(right_node)
**else**:
node['right'] = self.best_split(right_node)
self.split_branch(node['right'], depth+1)
决策树训练:
- 基于所有要素上所有可能分割点的最低基尼系数分割输入数据
- 基于上述分割创建了两个数据集
- 同样，每个新数据集都是根据所有可能的最低基尼系数进行分割的
- 在所有特征上分割点
- 同样，从每个节点创建两个数据集
- 依此类推，直到根据最大深度或最小样本数标准创建叶节点
决策树预测:
- 接受新样本
- 基于根节点的拆分标准评估样本
- 前进到下一个节点
- 基于该节点的拆分标准评估样本
- 前进到下一个节点
- 依此类推，直到到达叶节点
- 基于叶节点的类确定预测
决策树要求每个新样本遍历训练数据的学习表示。然后将 n 个样本中的每一个附加到输出中。
**def** predict_sample(self, node, sample):
*#print(node)*
**if** sample[node['feat']] < node['val']:
**if** isinstance(node['left'],dict):
**return** self.predict_sample(node['left'],sample)
**else**:
**return** node['left']
**else**:
**if** isinstance(node['right'],dict):
**return** self.predict_sample(node['right'],sample)
**else**:
**return** node['right'] **def** predict(self, X_test):
self.y_pred = np.array([])
**for** i **in** X_test:
*#print(i)*
self.y_pred = np.append(self.y_pred,self.predict_sample(self.root,i))
**return** self.y_pred
从头开始创建决策树算法教会了我很多东西。首先，我从来没有意识到分割搜索是多么详尽，直到我被迫编程。此外，分裂纯洁的概念现在对我来说也更清楚了。显而易见，当你构建决策树时，很容易过度拟合。只有通过创建像最大深度和最小样本数这样的超参数，树才有机会泛化。
在未来，我会跟进一个随机森林和 Kmeans 的解释，敬请期待！
分类模型的模型性能和成本函数
分类模型是预测分类的 Y 变量的机器学习模型:
- 雇主会离开还是留下来？
- 病人到底有没有癌症？
- 该客户属于高风险、中等风险还是低风险？
- 客户会支付贷款还是违约？
Y 变量只能取 2 个值的分类模型称为二元分类器。
分类模型的模型性能通常在哪个模型性能最相关方面存在争议，尤其是当数据集不平衡时。用于评估分类模型的通常模型性能测量是准确度、灵敏度或召回率、特异性、精确度、KS 统计和曲线下面积(AUC)。
让我们基于一个预测贷款违约的例子来理解一些模型性能测量。贷款违约数据集是不平衡数据集的一个典型例子，其中两个类是贷款违约 Y 和贷款违约 n。贷款违约者的数量通常只占整个数据集的很小一部分，不超过 7–8%。这提供了一个经典的不平衡数据集来理解为什么成本函数在决定使用哪个模型时至关重要。
在我们深入研究如何制定成本函数之前，让我们看看混淆矩阵、假阳性、假阴性和各种模型性能测量的定义的基本概念。
什么是混淆矩阵？
混淆矩阵是包含二进制分类器输出的表格。让我们看看预测贷款违约的二元分类器的混淆矩阵——0 表示客户将支付贷款，1 表示客户将违约。我们进一步讨论的积极类是 1(将违约的客户)。
矩阵的行表示观察到的或实际的类别，列表示预测的类别。
Confusion matrix
Definition of TP, FP, TN, FN
现在，让我们试着理解每个模型性能度量在混淆矩阵的组成部分方面转化为什么。
准确(性)
准确性是模型通过记录总数做出的正确预测的数量。最好的准确度是 100%,表明所有的预测都是正确的。
对于不平衡的数据集，精确度不是模型性能的有效衡量标准。对于默认比率为 5%的数据集，即使所有记录都被预测为 0，该模型仍将具有 95%的准确性。但是这种模式会忽略所有的默认值，对业务非常不利。因此，在这种情况下，准确性不是衡量模型性能的正确标准。
敏感性或回忆
灵敏度(回忆或真阳性率)的计算方法是正确的阳性预测数除以阳性总数。它也被称为回忆(REC)或真阳性率(TPR)。
Sensitivity
特征
特异性(真阴性率)计算为正确阴性预测数除以阴性总数。
Specificity
精确
精度(阳性预测值)的计算方法是正确的阳性预测数除以阳性预测总数。
Precision
KS 统计
KS 统计量是正分布和负分布之间分离程度的度量。KS 值为 100 表示分数对记录进行了精确的划分，一个组包含所有的正数，另一个组包含所有的负数。在实际情况下，高于 50%的 KS 值是理想的。
ROC 图和曲线下面积(AUC)
ROC 图是 X 轴表示特异性，Y 轴表示灵敏度的图。ROC 曲线下的面积是模型性能的度量。随机分类器的 AUC 是 50%,而完美分类器的 AUC 是 100%。对于实际情况，超过 70%的 AUC 是理想的。
精确度与召回率
回忆或敏感度为我们提供了关于模型在假阴性(对将违约的客户的不正确预测)上的表现的信息，而精确度为我们提供了模型在假阳性上的表现的信息。基于预测的结果，精确度或召回率对模型来说可能更重要。
成本函数在决定哪个不正确的预测可能更有害(假阳性或假阴性)时发挥作用(换句话说，哪个性能测量更重要——精确度或召回率)。
净收入函数
净收入或成本函数是通过为每个假阳性和假阴性分配成本并基于正确和不正确的预测得出总收入而得到的。让我们假设此贷款默认数据集的成本和收入如下:
Cost components of a loan default prediction model
Allocation of FP, FN & TP cost
Net revenue
由于假阴性成本最高，最 的最优模型将是假阴性 最小的模型。换句话说，与其他模型相比，具有更高敏感度的模型将获得更高的净收入。
现在我们有了计算净收入的方法，让我们根据混淆矩阵来比较两个模型:
混淆矩阵 A
混淆矩阵 B
模型 B 由于其较低的假阴性被证明是一个更好的模型，因此可以被选择用于预测违约。
总之，在各种成本成分的相当近似的估计可用的商业场景中，基于成本函数的模型性能度量与传统的模型性能度量(如灵敏度、特异性等)相比，将给出对模型选择的更好的洞察。与给出统计术语相比，如果用成本和收入来解释模型，对客户来说也更有意义。
型号选择 101，使用 R
使用 R 对简单模型选择进行快速标记
我们在做什么？
由于这是一个非常入门的模型选择，我们假设你已经获得的数据已经被清理，擦洗和准备好了。数据清理本身就是一个完整的主题，实际上是任何数据科学家的主要时间接收器。如果您想自己下载数据并跟随学习，请阅读本文末尾！
编辑:我为这篇文章做了一个续集，关于可视化和绘制我们找到的模型，如果你想在读完这篇文章后看看的话！：
你已经有了一个模型，现在呢？
medium.com](https://medium.com/@peter.nistrup/visualizing-models-101-using-r-c7c937fc5f04)
如果你喜欢这篇文章并想看更多，请务必关注我的简介。
让我们看看管道:
这是我用来创建一个简单的 LM 或 GLM 的框架:
- 使用所有可用变量和数据创建一个基础模型
- 如果 R 不工作，分解分类变量
- 增加相关的电力变换
- 添加相关变量交互
- 删除无关紧要的变量和相关的测试标准
——重复第 3-5 步，直到你用尽所有选择 - 删除任何异常数据点
- 评估你的模型
- 将你的发现可视化
1.创建基础模型
让我们从建立一个工作空间和加载我们的数据开始。在本例中，我们正在处理一个描述女性就业状况的数据集，该数据集基于您是否是外国人、政府补贴金额(对数转换)、年龄、受教育年限和子女数量(分布在两个分类变量*【young . children】*和 school.children ):
rm(list=ls()) # “Clear current R environment”setwd(“C:/R/Workspace”) # Setting up your workspace
dat <- read.table(“employment_data.txt”, header=T) # Load that data
str(dat)
summary(dat)
哪些输出:
我们注意到的第一件事是，我们的响应变量是二项式的(显然),这表明我们有一个二项式分布，这意味着我们必须拟合一个 GLM 而不是传统的 LM :
fit.1 <- glm(employed == "yes" ~ ., data = dat, family = binomial)
summary(fit.1)
默认情况下，这种拟合是我们可以使用的最通用的拟合，它使用数据集中的每个变量( ~)来拟合二项式模型( family = binomial )关于具有值“yes”的响应变量“employed”。)给出如下输出:
好的，那么马上就有几个问题，我们不喜欢看到 p 值高于 0.05 ，更不喜欢高于 0.1 ，但是在我们不顾一切地删除它们之前，让我们先检查变量交互和幂变换！
2.分解分类变量
让我们考虑没有任何孩子和实际上有大于零的任何数量的孩子之间存在分类差异的可能性，因此我们添加分类变量，用于有 0 个孩子 : ‘因子(young.children == 0)’ ，‘因子(school.children == 0)’ 和根本没有任何孩子的组合因子*‘因子(young . children+school . children = = 0)’*
我们可以用新的变量更新我们的拟合:
tempfit <- update(fit.1, .~. + factor(young.children == 0)
+ factor(school.children == 0)
+ factor(young.children +
school.children == 0))
summary(tempfit)
所以我们已经在 AIC 方面对我们的模型进行了一点改进，从 1066.8 到 1050.2 ！让我们看看连续变量，观察可能的幂变换:
3.添加相关的电源转换
寻找二项分布的潜在幂变换的一个简单方法是使用这个自定义函数:
logit.plot <- function(x, y, h, link=’logit’, xlab=deparse(substitute(x)), yname=deparse(substitute(y)), ylab, rug=T, data, …){
if(!missing(data)){
call <- match.call()
dataPos <- match(“data”,names(call))
return(invisible(with(data, eval(call[-dataPos]))))
}
if (length(levels(factor(y)))!=2) stop(‘y must be binary’)
succes <- levels(factor(y))[2]
if (missing(ylab)) ylab <- paste(binomial(link)$link,’ P{‘,yname,’=’,succes,’|’,xlab,’}’, sep=’’, collapse=’’)
if (is.factor(y)) y <- as.numeric(y==levels(y)[2])
x.seq <- seq(min(x),max(x),length=101)
smooth <- binomial(link)$linkfun(ksmooth(x, y, ’normal’, h, x.points=x.seq)$y)
plot(smooth~x.seq, ylab=ylab, xlab=xlab, type=’l’,…)
if (rug) rug(x)
invisible(xy.coords(x = x.seq, y = smooth, xlab = xlab, ylab = ylab))
}
这实际上相当简单，它绘制了链接(E[y|x]) 与 x 的关系，其中 E[y|x]使用 Nadaraya-Watson 内核回归估计进行估计:
这需要以下参数:
x —您的解释变量
y —二元结果变量
h—Nadaraya-Watson 核回归估计的带宽
数据 —自我解释
… —您想要传递给 plot()的任何附加参数
使用该函数在不同带宽上迭代，我们得到以下类型的图:
for(i in seq(2.5,10,length.out = 6))
logit.plot(x = age, y = employed == ‘yes’, h = i, data = dat, main = paste0(“bandwidth = “,i))
在这个带有“年龄”的例子中，我们可以看到，该函数在大约 7 的带宽附近开始变得平滑，该图可以近似为 2。或者 3 年。次数多项式，对于“education”和“gov.support”也是如此，但为了简单起见，我们将考虑所有三者的形状都为 2 的情况。次数多项式:
tempfit <- update(tempfit, .~. + I(age^2)
+ I(education^2)
+ I(gov.support^2))
summary(tempfit)
就 AIC 而言，改进相当大，从 1050.2 到 1017.7 ！尽管有许多无关紧要的变量！
这是我们对“完整”模型的第一次尝试，所以让我们将其定义为“fit.2 ”,然后继续。
fit.2 <- tempfit
4。添加变量交互
检查变量相互作用的最简单的方法是使用 R 函数’add 1【T1]'，这只是定义一个测试范围以及相对于原始模型测试时使用哪个测试的情况。f 检验通常只与 LM 和 AOV 车型相关，因此我们可以“安全地”忽略该检验标准，我们将使用χ-检验** (Chisq 或 Chi)😗*
add1(fit.2, scope = .~. + .^2, test=”Chisq”)
这个范围仅仅是要求测试当前的模型(。~.)加上现有变量之间的相互作用(+。^2)，这将输出许多交互，一些具有统计上显著的 p 值，但是手动排序可能会很烦人，所以让对列表进行排序，这样我们就可以在顶部得到最低的 P 值:
add1.test <- add1(fit.2, scope = .~. + .^2, test=”Chisq”)
add1.test[order(add1.test$’Pr(>Chi)’),]
是的，所以看起来外国人和年龄变量之间可能会有交互作用。在简单地添加具有最低 P 值的交互之前，要考虑的一件事是，这在我们当前的模型中是否有意义，现在年龄实际上是我们模型中最重要的变量，因此我们可能会认为添加外国人和年龄之间的交互更直观，为了简单起见，我们将坚持使用外国人:年龄交互。
在添加变量交互外国人:年龄之后，让我们测试更多的交互**😗*
add1.test <- add1(update(fit.2, .~. + foreigner:age), scope = .~. + .^2, test=”Chisq”)
add1.test[order(add1.test$’Pr(>Chi)’),]
现在看来在**外国人:因素(young . children+school . children = = 0)**中有一个重要的互动
经过几轮后，我们最终没有看到新的具有统计意义的交互，到最后我们添加了以下交互:
- 外国人:年龄
- 外国人:因素(young . children+school . children = = 0)
+年龄:school.children - 政府支持:因素(young.children == 0)
因此，让我们用新的变量交互来更新我们的拟合，如下所示:
fit.3 <- update(fit.2, .~.
+ foreigner:age
+ foreigner:factor(young.children + school.children == 0)
+ age:school.children
+ gov.support:factor(young.children == 0))summary(fit.3)
再次对 AIC 进行小改进，从 1017.7 变为 1006.9 。
5。移除无关紧要的变量
这个过程与步骤 4 中的最后一个非常相似。我们现在将简单地使用 R 中的 drop1 函数，而不是 add1、，并且由于我们寻求移除而不是追加变量，我们寻求最高的 P 值而不是最低的(我们仍将使用 χ -test 作为我们的标准):
drop1.test <- drop1(fit.3, test=”Chisq”)
drop1.test[rev(order(drop1.test$’Pr(>Chi)’)),]
这和我们的模型告诉我们的大部分是一样的——总结，政府支持肯定没有统计上的显著性，所以我们将首先移除它，以此类推，我们最终移除以下变量:
- 政府支持
- 青少年
- 教育
- 教育
移除这些变量后，我们看到所有剩余变量都具有统计显著性，因此让我们通过移除上面列出的变量来更新我们的拟合:
fit.4 <- update(fit.3, .~.
— I(gov.support^2)
— young.children
— education
— I(education^2))summary(fit.4)
不错， AIC 的另一个改进虽然微不足道，但这个模型相对于我们之前的模型的主要优势是增加了解释变量数量减少所固有的简单性**！**
有人可能会问*“为什么我们不删除gov . support变量？在看我们模型的总结时，这显然是微不足道的！”*这是因为边际原则禁止我们移除一个无关紧要的变量，如果该变量与另一个变量有显著的相互作用，如gov . support:factor(young . children = = 0)。
鉴于明显不重要，并且与 young.children == 0 变量的交互仅略微重要( p = 0.435) ，您可能会认为移除 gov.support 将有利于模型的简单性，然而，当进一步检查从模型中移除 gov.support 时，变量交互分裂为两个变量 因此，没有给我们任何额外的简单性， AIC 、所有其他系数以及零** -和剩余偏差保持完全相同，因此，我关闭它，将其留在模型中。**
对我们新改进的模型进行 add1 和 drop1 测试表明，没有新的显著交互，所有当前变量都是显著的，所以我们完成了！最终拟合为:
glm(employed == “yes” ~ foreigner
+ gov.support
+ age
+ school.children
+ factor(young.children == 0)
+ factor(school.children == 0)
+ factor(young.children + school.children == 0)
+ I(age^2)
+ foreigner:age
+ foreigner:factor(young.children + school.children == 0)
+ age:school.children
+ gov.support:factor(young.children == 0), family = binomial, data = dat)
6。移除异常值
现在我们有了一个满意的模型，我们可以寻找对模型有负面影响的异常值。
使用“汽车”包https://www.rdocumentation.org/packages/car/versions/1.0-2我们可以使用 influencePlot()和 outlierTest()函数来查找潜在的异常值:
我们看到数据点 416 在两个测试中都被归类为异常值，我们可以查看一些图中的点，以判断是否要将其移除:
par(mfrow = c(2, 2)) # 2 row(s), 2 column(s)plot(fit.4, col=ifelse(as.integer(rownames(dat))==416, “red”, “black”), pch=20)
这似乎很可能会对我们的模型造成一点影响，请注意，我们实际上应该使用皮尔逊残差来衡量我们的模型拟合度，因此，左上角的图中没有任何接近直线的东西是好的， Q-Q 图也与这种模型无关。
让我们试着去掉这一点，看看新的合身程度:
final.fit <- update(fit.4, data = dat[-c(416),])
summary(final.fit)
从最初的 1067，下降到几乎10001000AIC当比较两个不同数据集的 AIC(因为我们移除了点 416)时，我们实际上必须得出结论，416 在初始模型中也是异常值，移除它，然后将没有点 416 的初始模型的 AIC 值与没有点 416 的最终拟合值进行比较
查看另一轮 influencePlot() 和 outlierTest() 我们发现数据点 329 也在发挥作用，但是查看实际的曲线图我们会发现，我们无法像使用 416 那样真正证明删除数据的合理性。这是我们最后的合身衣服。
7.模型评估
现在我们有了一个最终拟合，我们不能自信地添加或删除任何交互变量和其他转换是时候评估我们的模型是否真正符合我们的数据，以及我们的最终拟合和第一次“天真”拟合之间是否有统计学上的显著差异。
让我们先来看看我们的拟合值与皮尔逊残差的对比:
par(mfrow = c(1, 1)) # 2 row(s), 2 column(s)
plot(p.resid ~ fit.val)
lines(ksmooth(y = p.resid, x = fit.val, bandwidth=.1, kernel=’normal’), col=’blue’)
abline(h=0)
这是一个相当不错的拟合，让我们看看主要解释变量的拟合情况:
除了’ gov.support ‘之外，一切看起来都很好，而且’ gov.support '弯曲的原因似乎是一个孤立点，与我们所有其他数据点相比，某人获得的支持量要低得多。
检查潜在的过度拟合
过度拟合是所有统计建模的祸根，你如何确保你的模型不仅仅符合你输入的精确数据？我们的目标是建立一个通用的模型，而不仅仅是迎合现有的数据。那么，我们如何测试我们的模型是否过度符合我们的数据呢？
一个流行的测试指标是通过交叉验证生成的 delta 值，我们可以通过使用来自’ boot '包的 cv.glm 函数来计算这些值，并将我们的最终拟合与我们的第一次拟合进行比较！
cv.glm(dat[-c(416),], final.fit, K = 13)$delta
cv.glm(dat[-c(416),], update(fit.1, data = dat[-c(416),]), K = 13)$delta
在上面的代码中，我们使用了 k 倍交叉验证和 k = 13 (因为 13 是 871 的因数，871 是我们在去除异常值后的数据长度)这意味着我们将数据分成 13 个“块”。
delta 值是一个向量，其中第一个分量是预测误差的原始交叉验证估计值**，第二个分量是调整后的交叉验证估计值(设计用于补偿不使用穷举测试方法(如*)引入的偏差***
运行上面的代码会产生以下增量值，请注意，这些值受一些随机变化的影响，因此您可能不会获得完全相同的值:
对于最终拟合，预测误差低于，即使在交叉验证测试时也是如此。因此，我们可以假设我们的模型没有过度拟合我们的数据！
拼图的最后一块
因此，现在我们已经得出结论，我们的模型实际上非常适合我们的数据，但它与没有任何转换和变量交互的“天真”模型相比，在统计上有显著差异吗？为此，我们可以使用 ANOVA 测试，我们只需在两次拟合中删除相同的数据点:
这两种拟合之间肯定有显著的差异，我们可以高兴地得出结论，我们的努力得到了回报！
8.模型可视化
现在我们已经有了一个模型，我们如何真正地可视化和解释它所说的关于我们的数据中的关系？
看看下面的走线通道，它使用了与本文相同的数据和型号。：
你已经有了一个模型，现在呢？
medium.com](https://medium.com/@peter.nistrup/visualizing-models-101-using-r-c7c937fc5f04)*
结束语
请记住，这纯粹是介绍性的并且这不是详尽的分析或结论！如果我们在追求中更加严格，我们会在我们模型的每个新迭代中加入交叉验证测试和方差分析测试。每当我们增加一个新的变量，交互作用或能量转换。
如果您有任何问题，请随时给我发消息，如果您觉得我错过了什么或做错了什么，请纠正我，请记住，这是对 R 建模的介绍，我很清楚，与更高级的方法相比，这个过程是高度简化的！
如果你想试试这个数据集的运气，试试这里:https://github . com/pela 15 AE/statmod/blob/master/employment _ data . txt
数据是丹麦语，因此要将标题和分类值转换为英语，请运行这段代码:
*names(dat) <- c(“employed”, “foreigner”, “gov.support”, “age”, “education”, “young.children”, “school.children”)
levels(dat$employed)[1] <- “yes”
levels(dat$employed)[2] <- “no”
levels(dat$foreigner)[1] <- “yes”
levels(dat$foreigner)[2] <- “no”*
结束语
如果你想看和了解更多，一定要关注我的 媒体 🔍*碎碎念 🐦***
阅读彼得·尼斯特拉普在媒介上的作品。数据科学、统计和人工智能…推特:@PeterNistrup，LinkedIn…
medium.com](https://medium.com/@peter.nistrup)***
数据科学的模型跟踪工具(mlflow)
Source: https://www.recordedfuture.com/machine-learning-cybersecurity-applications/
在数据科学工作中，Jupyter 笔记本是一个众所周知的工具。除此之外，我们可能会使用 databricks 的笔记本或 Colab(由谷歌)。产品化怎么样？
如何将我们的模型部署到生产中？如果你使用 TensorFlow 库来建立模型，我们可以使用 tensorflow 服务。使用 Azure 平台的话可以用 Data Science Workbench。如果我们真的使用上面的解决方案呢？一些公司建立了内部系统，如米开朗基罗(优步)、 FBLearner Flow(脸书)。然而，那些平台是为它们的内部系统设计的。
我们有开源吗？在此之前，我们有 modeldb 。我们有了新的选择，mlflow，它刚刚在上周(2018 年 6 月 5 日)发布了 alpha 版本。让我们参观一下这个产品。
首先，您必须安装 mlflow，然后在命令行中执行以下命令来启动 UI 服务
mlflow ui
并执行以下命令将数据传递给服务器
lr = LogisticRegression(max_iter=param['max_iter'], penalty=param['penalty'], solver=param['solver'])
lr.fit(x_train, y_train)
y_pred = lr.predict(x_test)(rmse, mae, r2) = eval_metrics(y_test, y_pred)mlflow.log_metric("rmse", rmse)
mlflow.log_metric("r2", r2)
mlflow.log_metric("mae", mae)
mlflow.log_param("max_iter", param['max_iter'])
mlflow.log_param("penalty", param['penalty'])
mlflow.log_param("solver", param['solver'])mlflow.sklearn.log_model(lr, "model")
稍后，转到 http://localhost:5000/
默认情况下，它会显示“日期”、“用户”、“源”、“版本”、“参数”和“指标”。
如需源代码，请访问 github 获取完整版本
估价
用户
系统将使用 pwd (linux 命令)来获取用户 ID。如果操作系统是 Window，它将使用“未知”的默认用户 ID。我相信他们对审计目的有所顾虑，这就是为什么他们不允许将用户 ID 作为参数传递。
# Copy from mlflow source code_DEFAULT_USER_ID = "unknown"def _get_user_id():
*""" Gets the ID of the user for the current run. """* try:
import pwd
return **pwd.getpwuid(os.getuid())[0]**
except ImportError:
return _DEFAULT_USER_ID
那么我们能在 Window 中做些什么来“劫持”它呢？:)我们可以覆盖全局值(在 alpha 版本中有效)。但是，不建议这样做，因为这会破坏审计/安全考虑。
mlflow.tracking._DEFAULT_USER_ID = '[https://github.com/makcedward/](https://github.com/makcedward/)'
来源
默认情况下，带有路径的文件名将被用作源文件。强烈建议提供源代码，以便更好地使用 UI 搜索。
# Copy from mlflow source codedef _get_main_file():
if len(sys.argv) > 0:
return **sys.argv[0]**
return None
版本
一般来说，版本号可能是“1.0.0”就可以了。但是，如果您定义的版本号是“123.456.789”或“测试”。你会有麻烦，因为用户界面显示多达 6 个字符。从源代码来看，它将版本简化为“shortVersion ”,允许最多 6 个字符。
# Copy from mlflow source codestatic renderVersion(run) {
if (run.source_version) {
**const shortVersion = run.source_version.substring(0, 6);**
if (run.source_type === "PROJECT") {
const GITHUB_RE = /[:@]github.com[:/]([^/.]+)\/([^/.]+)/;
const match = run.source_name.match(GITHUB_RE);
if (match) {
const url = "https://github.com/" + match[1] + "/" + match[2] + "/tree/" + run.source_version;
return <a href={url}>{shortVersion}</a>;
}
} else {
return shortVersion;
}
}
return null;
}
储存
记录将保存在文件系统中。换句话说，它不将数据持久化到 DB 中，并且它可能在 UI 搜索和其他系统问题中具有潜在的查询性能。相信将记录保存到数据库应该在路线图中。
API
API 接口由单线程的 Flask 构建。因此，它不能同时处理大量的请求。似乎在投入生产之前，必须先解决这个问题。
关于我
我是湾区的数据科学家。专注于数据科学、人工智能，尤其是 NLP 和平台相关领域的最新发展。
中:【http://medium.com/@makcedward/】T2
领英:【https://www.linkedin.com/in/edwardma1026
github:https://github.com/makcedward
https://www.kaggle.com/makcedward
使用 Python 和 Pyomo 对每周劳动力进行建模和优化
Credit: https://unsplash.com/@rawpixel
使用 Python 和 Pyomo 逐步建模和优化劳动力设计和分配问题。
在本帖中，我们将对一个调度问题进行建模和求解，在这个问题中，工人必须被分配到不同的班次，以优化给定的标准，满足不同的工作条件约束。
当开始一个新项目，计划开一家新商店，甚至准备大学的课程表时，负责任务的人有两种选择:手工求解或建模并作为优化问题求解。
即使有一个优秀的劳动力规划者，使用优化框架来处理这个问题也有很多好处:
- 扩展:问题可能会变得和适当的环境一样大，这可能会超出一个人的能力。
- 不平衡:在很多情况下，员工之间的不平衡是不可避免的。在这种情况下，决定将不是个人的。
Pyomo 作为优化建模环境
为了解决这个问题，我们将使用 Pyomo，用他们自己的话说，这是一种“基于 Python 的开源优化建模语言，具有多种优化功能”。对于那些已经使用 Python 的人来说，用 Pyomo 建模和解决问题会很简单。
Pyomo 的主页，这是一种用于线性编程的基于 Python 的可扩展开源优化建模语言…
www.pyomo.org](http://www.pyomo.org/)
Pyomo 允许在各种解决方案中进行选择，包括开源和商业解决方案。此外，用户可以选择在 Neos 服务器上解决问题，这是一个免费的基于互联网的解决方案，可以直接从 Pyomo 使用。这让我们不用在机器上下载和安装不同的解算器，也允许免费使用商业解算器。更多信息在https://neos-server.org/neos。
我已经在我的机器上安装了带有 conda 的 Pyomo，但也可以用 pip 来完成。可以在 http://www.pyomo.org/installation 的阅读安装说明。
在这篇文章中，我们将利用硬币或项目 Cbc，【https://projects.coin-or.org/Cbc】T4，来解决一个混合整数规划问题。
问题描述
大学校园里新开了一家食品店，一周 7 天，一天 24 小时营业。每天有三次八小时轮班。早班为 6 时至 14 时，晚班为 14 时至 22 时，夜班为 22 时至次日 6 时。
晚上只有一个工人，而白天有两个，除了星期天，每班只有一个。每个工人每周最多不超过 40 小时，并且在两班之间必须休息 12 小时。
至于每周的休息日，一个周日休息的员工也会喜欢在周六也休息。
原则上，有 10 名员工可用，这显然超出了规模。需要的工人越少，其他商店的资源就越多。
模型组成
为了使用 Pyomo 并解决这个问题，需要导入一些包。
第一步是输入数据，这意味着提供模型，在这种情况下，我们正在考虑的日子，工人，班次…我们不能明确地定义它们，稍后通过使用 AbstractModel() 输入数据，但对于这个帖子，继续使用 ConcreteModel() ，因此所有数据必须在建模步骤中可用。
我们已经添加了模型所需的所有数据。我们将能够在定义约束时调用这些元素，而不需要手动插入每个元素或考虑 if、else 子句。
为了构建模型，我们需要初始化模型并创建决策变量。
将模型初始化为对象后，将变量、约束和目标函数等元素作为属性添加。正如前面所解释的，我们创建了一个混凝土模型()，因为此时数据正在被提供。我们用 Var() 变量添加到模型中，用列表索引。
- works:由工人、天数和班次索引的二元变量。1 如果工人必须在那天轮班工作。
- 需要:由工人索引的二进制变量。1 员工是否必须包括在员工队伍中。
- no_pref:由 workers 索引的二进制变量。1 如果它在周日不工作，但在周六工作。
定义完所有变量后，我们就可以添加目标函数了。
目标是找到一个最大限度地减少工人数量的时间表，一旦实现了这个目标，也是最大限度地减少在星期天工作但在星期六不工作的工人数量。我们将工人数量部分乘以一个足够大的常数，这样只有在决定了最佳工人数量后，才考虑最小化周末偏好。
我们还需要添加约束，首先创建一个调用 ConstraintList() 的约束容器，然后用函数 add 将我们想要的任何约束添加到容器中。约束的代码如下，每个约束的解释都是行内注释:
现在，创建的模型可以求解了。一个模型可以用 model.pprint() 来可视化。这个问题是由二进制变量组成的，所以一个混合整数规划求解器适合我们的要求，CBC 从 COIN-OR 将被选中，【https://projects.coin-or.org/Cbc】T2。
这需要在本地安装求解器(不要忘记将它添加到路径中，以便 Pyomo 识别它)。然而，在 neos-server 中运行它将我们从这一过程中解放出来。
提取溶液
我们这样做是为了得到一个问题的解决方案，所以这是下一步。问题的解决方案是稀疏的，从这个意义上说，我们只是想知道它们中哪些是非零的，为了可视化，我们必须处理给定的解决方案。
我们现在可以检查所获得的解，并验证所有提出的约束条件是否得到满足。以 JSON 格式保存最佳时间表的字典，解决方案如下:
最佳解决方案需要激活 7 个工人。然后，有 7 名工人，他们不满意的人数是 2 (W6，W9)，因为他们在星期天不工作，但在星期六必须工作。观察问题公式，由于周日有 3 个班次，周六有 5 个班次，人们不能期望少于 2 个周六工人在周日不工作。
结论
在这篇文章中，我们已经完成了以下目标:
- 解释 Pyomo 作为 Python 中优化建模框架的基础知识。
- 学习如何建立一个特定问题的模型，在这个问题中，劳动力计划者必须在各种限制下优化每周时间表。
这篇文章旨在友好地介绍 Pyomo 的使用和优化问题的建模，因此，解决由于问题的规模或难以制定一个好的模型(是的，建模可以被认为是一门艺术)而导致的更困难的问题可能需要更深入地了解这些主题。
感谢您的阅读。玩得开心！
为 MLB 2018 年季后赛球队建模
随着十月的到来，季后赛终于开始了。在 162 场棒球赛季的整个过程中，球队会生成数据，而且是大量的数据，因为他们努力进入季后赛，有机会赢得世界职业棒球大赛。越来越多的这种数据被逐场使用，以优化防守排列、场地选择和情况调整。
对于这个项目，我想收集更多的宏观赛季数据，看看我们是否可以预测——以合理的准确度——哪些球队将根据他们的年终统计数据进入季后赛。在建立了截至 2017 年的数据模型后，我们可以引入未经测试的 2018 年数据，看看模型的表现如何。
我们的高级流程是:
- 数据收集和清理
- 数据准备
- 模型测试和调整
- 2018 年季后赛球队预测
所有网页抓取和建模的项目代码都可以在 我的 GitHub 这里 找到。
数据收集和清理
我想要一个项目，将工作在网络抓取，棒球统计提供了充分的机会，与来自网络的数据工作。经过一番搜索，我发现 baseball-reference.com 有一套强大的可废弃的数据，并使用该网站收集信息。
为了从击球和投球的角度来处理这个问题，我决定收集每一个方面的代表性数据。Baseball-reference 有非常有用的汇总详细统计数据的表格:
Example of annual batting statistics
Example of annual pitching statistics
在检查他们的 robots.txt 页面时，我发现我想访问的页面被允许抓取，并且要求 3 秒钟的延迟，以免服务器过载。
下一步是确定分析中包括多少年。事实证明，1969 年是在网站上包含更详细的团队统计数据的转折点，因此这成为选择的自然年份，以便使用的数据既全面又一致。
击球和投球数据的 URL 格式也被证明是有益的，年份是唯一可以改变的变量。这允许我们用每个目标年份编写一个简单的for循环来捕获每个页面上的 HTML。
for i in range(1969, 2019):
site = requests.get('[https://www.baseball-reference.com/leagues/MLB/{}.shtml'.format(i))](https://www.baseball-reference.com/leagues/MLB/{}.shtml'.format(i))) soup = BeautifulSoup(site.text, 'html.parser')
.
.
.
# Rest of code
实际的数据抓取需要检查 HTML 代码，以找到在我们的抓取中调用的正确标签。同样，数据结构很有帮助，因为它在我们的目标特性中是标准化的:
很明显，我们想要捕获的大部分数据都包含在 data-stat 属性下的一个标记中。因此，我们可以定义一个函数，它接收 class 和 data-stat 并返回该标签的内容。对于任何不遵循这种模式的，我们可以只写一次性代码。
def get_attribute(class_side, data_stat):
all_instances = soup.find_all('td', attrs = {'class':
class_side , 'data-stat' : data_stat})
all_instances = [i.text for i in all_instances][:-2]
# [:-2] eliminates league total and average return all_instances# Example call to the function:batters_used = get_attribute('right', 'batters_used')
一旦获得各个数据点的正确代码就绪，每年的数据就会迭代地添加到 for 循环之外的数据框中，这样，到刮擦结束时，数据框中就有了 1969 年至 2018 年的所有联盟信息。
在这一点上，数据没有丢失值，并且相当干净。最后一步是简单地合并击球和投球数据集，为每支球队每年创建一个单独的观察。
df_all = df_pitching.merge(df_batting, on=['year','team'], how='left', suffixes=('_pitching', '_batting'))
最后，球队进入季后赛的年份取自维基百科。
数据准备
作为准备数据的第一步，我们将通过查看季后赛球队与非季后赛球队在数字列平均值方面的不同来给出变量的上下文:
我们看到季后赛球队在三振出局和跑垒(以及与全垒打和三垒打相关的变量)方面的出色表现，更少的得分/更低的 ERA，以及更稳定的投手阵容(使用更少的投手)——这不是革命性的，但对上下文有帮助。
接下来，数据准备需要两个主要步骤:标准化和缩放数据。
标准化是指变量中的某些幅度差异可能是由于所玩游戏数量的轻微变化(即 163 比 161)造成的。同样地，对于投球统计，有些是每九局报告一次，而有些则是全部。
因此，任何不在预定义范围内的变量(例如，每九局)，我们将标准化为每场比赛的统计数据。
Are differences in variables a result of different number of games played? Or actual differences?
Per-game standardization allows for valid comparisons across teams.
在这个项目中，缩放数据有一个有趣的变化。因为每个赛季都会重置一个球队的统计数据，而季后赛是由单个赛季的表现决定的，所以简单地从整体上缩放数据是行不通的。相反，我们必须在每个季节内单独扩展。
我选择的定标器是最小最大定标器。这将获取数据，并根据每列的最大值和最小值将其放在 0 到 1 的范围内。选择这一比例的原因是这种方法很容易解释——可以看出谁在某一年表现最好，并比较不同年份的相对表现。
为了实现逐季缩放的变量，我们可以创建一个空的数据框，遍历每一年，将 MinMaxScaler 应用于该年的数据，并将其附加到新创建的数据框。
df_scaled = pd.DataFrame()for yr in list(set(df_all.year)):
df = df_all[df_all.year == yr]
df = df.select_dtypes(include=['float','int'])
scaler = MinMaxScaler()
df_transformed = scaler.fit_transform(df)
df_scaled = df_scaled.append(pd.DataFrame(df_transformed))
A sample of the MinMaxScaler transformed data. 0 indicates bottom of the league for a given year while 1 indicates the best performing team in that category.
型号选择和测试
这个项目很好地解决了分类问题——如果球队在某一年进入了季后赛，那么季后赛一栏标记为 1，否则标记为 0。因此，我们可以测试一系列分类模型，看看结果如何比较。
为了评估我们模型的成功，我们将查看几个不同的指标:精度、回忆、 F-1 分数，以及曲线下接收器操作特征面积 (ROC AUC)。要更深入地阅读分类指标，这是一篇由威廉·科尔森撰写的好文章。
我们的数据是不平衡的，因为我们的样本中只有约 22%被标记为进入了季后赛，这意味着(1)当我们分成训练集和测试集时，我们应该对样本进行分层，以及(2) 22%是我们希望改进的原始基线值，如果随机选择，预计会被猜到。
分层是一种确保我们的测试集有足够的目标变量表示的方法，通过在 sklearn train_test_split函数中指定分层选项来实现。
X_train, X_test, y_train, y_test = train_test_split(df_features, df_targets, random_state = 60, stratify=df_targets, train_size = 0.7, test_size=0.3)
车型
我们将测试五种不同的分类模型，每种模型都采用不同的方法来预测数据(我还尝试了一种带有投票分类器的元估计器，但这并没有改善结果)。
我们将经历两次模型测试迭代——首先是开箱即用的模型，然后是网格搜索以优化超参数。
我们将测试的模型是逻辑回归、随机森林、KNeighbors 分类器、支持向量分类器(SVC)和 XGBoost 分类器(XGB)。
逻辑回归使用误差最小化来创建分隔数据的线性边界，随机森林构建一系列决策树来进行预测，KNeighbors 计算一个点与其“n”个邻居之间的距离，然后汇总邻居的预测。
SVC 使用超平面划分数据，并根据点相对于超平面的位置进行预测，XGBoost 使用迭代残差调整来改善预测。
既然已经收集、清理和准备好了数据(包括分成训练集和测试集)，运行和测试开箱即用模型的实际操作就只有几行代码:
from sklearn.neighbors import KNeighborsClassifierknc = KNeighborsClassifier()knc.fit(X_train, y_train)print(classification_report(y_test, knc.predict(X_test)))
print(roc_auc_score(y_test, knc.predict(X_test)))
解释精确度的方法是，对于这个模型，在预测进入季后赛的球队中，61%是准确的(精确度)。在应该被标记为进入季后赛的球队中，42%被正确识别(回忆)。F-1 分数是两者的混合平均值。ROC 绘制假阳性率对真阳性率，然后 AUC 使用微积分和积分来确定曲线下的面积(越高=越好)。
对于这五个模型中的每一个，我首先拟合和预测现成模型的概率。SVC 和逻辑回归模型实现了测试数据的最高精度(0.65)，最高的召回率和 ROC AUC 分数是 XGBoost 分类器(分别为 0.57 和 0.74)。
提高车型性能
接下来，我们将看看是否可以通过调整超参数来改进模型。因为总的数据集不是太大，我们可以用一系列不同的超参数选项进行全面的网格搜索。
网格搜索的基本前提是，它用参数的每种排列来拟合模型，然后返回得分最高的那个。我们通过传入超参数选项的字典来初始化和拟合模型，然后使用标准的 sklearn fit方法。
params_knc = {'n_neighbors':[3, 5, 7, 9, 11], 'weights':['uniform', 'distance'], 'metric':['euclidean','minkowski']}gs_knc = GridSearchCV(KNeighborsClassifier(), params_knc, cv = 5)gs_knc.fit(X_train, y_train)
为了得到最佳参数，我们可以简单地在我们拟合的模型上调用.best_estimator_，我们可以用与前面模型相同的调用得到分类报告和 ROC AUC 得分。
print(gs_knc.best_params_)
print(classification_report(y_test, gs_knc.predict(X_test)))
print(gs_knc.best_params_)print(classification_report(y_test, gs_knc.best_estimator_.predict(X_test)))print('ROC AUC Score: %.2f' % roc_auc_score(y_test, gs_knc.best_estimator_.predict(X_test)))
KNeighbors Classifier Grid Search Results
最终，超参数调整略微改进了模型，最高精度提高了 0.05 到 0.70(随机森林)，而召回没有看到任何改进。
调整后的 SVC 和逻辑回归在曲线下面积上得分最高:
预测 2018 年 MLB 季后赛球队
该项目的目标是根据年终统计数据预测哪些球队将进入 2018 年季后赛。由于季后赛正在进行，我们可以在这个看不见的数据上测试数据的表现(2018 年的数据没有在任何模型测试中使用)。
我们可以从两个角度来处理这个问题。首先，我们将让模型预测现状，看看哪些球队通过了季后赛预测的门槛。接下来，使用我们的领域知识，即 10 支球队将进入季后赛(包括外卡)，我们可以选择预测概率最高的 10 支球队，看看其中有多少是正确的。
模型预测
使用我们优化的模型，最好的结果来自 XGBoost Classsifier，其精度为 1.00(所有预测都进入了季后赛)，召回率为 0.80(在进入季后赛的 10 支球队中，模型正确识别了 8 支)。
XGB 分类器错误识别的球队是密尔沃基队和科罗拉多队，其中一支球队——科罗拉多队——是 NL 的外卡参与者，这意味着他们不会在模型适用的一些过去季后赛球队的结构下进入季后赛。
用概率预测
在模型上使用predict_proba方法，我们可以利用我们对 10 支季后赛球队的领域知识，来看看预测的前 10 支球队与进入季后赛的 10 支球队相比如何。
事实上，这种方法为优化的 SVC 模型产生了 10 个正确预测球队中的 9 个的结果，唯一的错误是模型预测华盛顿将进入季后赛，而科罗拉多不会。
最终想法
一个简单的模型在大约 22%的情况下随机做出了正确的预测，而我们最好的模型在这个数字上提高了近 160%,预测准确率为大约 57%。
此外，当将模型应用于看不见的数据时，在领域知识的帮助下，我们将正确预测的百分比提高到了 90%。
虽然调整并没有显著改善模型，但使用不同的模型确实改变了预测能力和最终结果的准确性。在性能最差和最好的分类模型之间，召回率提高了 40%以上。
改善结果的后续步骤可能包括:
- 将数据分成不同的时间段，并优化每个独特时间段的超参数，以查看季后赛因素是否变化
- 将预测概率方法应用于每年，并观察它捕捉季后赛球队的准确性
- 整合其他数据源(成功率可能是一个很高的预测值，但我想在本次分析中排除它)
- 测试神经网络或更复杂的模型
- 当前数据集中更深层次的特征工程
对我来说，这一点的一个有趣的未来应用将在明年投入生产，在整个赛季中实时跟踪球队的预测季后赛概率，并查看从开始到结束的赛季如何利用这一模型改变季后赛预测。
如果您有任何反馈或想法要分享，可以通过 jordan@jordanbean.com 联系我或在 LinkedIn 上给我发消息。
建模:教授机器学习算法以交付商业价值
如何训练、调整和验证机器学习模型
这是关于我们如何在功能实验室进行机器学习的四部分系列中的第四部分。完整的文章集可以在下面找到:
- 概述:机器学习的通用框架
- 预测工程:如何设置你的机器学习问题
- 特征工程:机器学习的动力
- 建模:教授算法(本文)
这些文章涵盖了应用于预测客户流失的概念和完整实现。project Jupyter 笔记本在 GitHub 上都有。(充分披露:我在 Feature Labs 工作，这是一家初创公司开发工具，包括 Featuretools ，用于解决机器学习的问题。这里记录的所有工作都是用开源工具和数据完成的。)
机器学习建模过程
预测和特征工程的输出是一组标签时间，我们想要预测的历史实例，以及特征，用于训练模型以预测标签的预测变量。建模的过程意味着训练机器学习算法来根据特征预测标签，根据业务需求对其进行调整，并根据维持数据对其进行验证。
Inputs and outputs of the modeling process.
建模的输出是经过训练的模型，可用于推理，对新数据点进行预测。
机器学习的目标不是一个在训练数据上做得很好的模型，而是一个证明它满足业务需求并且可以在实时数据上部署的模型。
与特征工程类似，建模与机器学习过程中的先前步骤无关，并且具有标准化输入，这意味着我们可以改变预测问题，而无需重写我们的所有代码。如果业务需求发生变化，我们可以生成新的标签时间，构建相应的特征，并将其输入到模型中。
客户流失建模的实现
在这个系列中，我们使用机器学习来解决客户流失问题。有几种方法来制定任务，但我们的定义是:
在每个月的第一天预测这个月哪些客户会流失。使用一个月的提前期，客户流失率为 31 天，没有订阅。提前期为 1 个月，这意味着我们提前 1 个月进行预测:1 月 1 日，我们预测 2 月份的客户流失。
尽管机器学习算法听起来在技术上可能很复杂，但用 Python 实现它们很简单，这要感谢像 Scikit-Learn 这样的标准机器学习库。作为一点实用的建议，经验结果已经表明机器学习模型和超参数的选择很重要，但是不如特征工程重要。
因此，理性的决定是将大部分精力放在预测和特征工程上，并插入一个预先构建的机器学习解决方案。
在这个项目中，我和 Scikit-Learn 一起快速实现了一些模型。为了让数据为机器学习做好准备，我们必须采取一些基本步骤:缺失值插补、分类变量编码，如果输入维度过大，还可以选择特征选择(详细信息参见笔记本)。然后，我们可以用标准建模语法创建一个模型:
指标和基线结果
在应用机器学习之前，最好建立一个幼稚基线来确定机器学习是否真的有帮助。对于分类问题，这可以简单到猜测保留测试数据中所有示例的训练数据中的多数标签。对于客户流失数据，猜测每个测试标签不是流失产生 96.5%的**。**
这种高精度听起来可能令人印象深刻，但是对于不平衡的分类问题——一个类比另一个类表示得多——精度不是一个足够的度量。相反，我们希望使用召回率、精确度、或 F1 分数。
回忆表示数据中实际搅动的百分比，我们的模型用天真的猜测记录 3.5%来识别。Precision 测量由我们的模型预测的和实际的的百分比，简单的分数是 1.0%。F1 分数是这些测量值的调和平均值。
由于这是一个分类问题，对于机器学习基线，我尝试了一个表现不佳的逻辑回归*。这表明问题可能是非线性的，所以我的第二次尝试使用了随机森林分类器，结果更好。随机森林训练速度快，相对可解释，高度准确，通常是可靠的模型选择。*
无机器学习、逻辑回归和具有默认超参数的随机森林的指标如下所示:
Metrics recorded by baseline models
基于时间序列分割*，使用约 30%的数据对每个模型进行了评估。(在评估时序问题中的模型时，这一点至关重要，因为它可以防止训练数据泄漏，并可以对新数据的实际模型性能提供良好的估计。)*
使模型与业务需求保持一致
尽管 ml 模型的指标比没有机器学习的要好，但我们希望根据业务需求优化给定指标的模型。在这个例子中，我们将关注召回率和精确度。我们将调整模型，通过调整阈值来实现一定的召回率，阈值是一个观察结果被分类为正面的概率——客户流失。
精度和召回调整
在机器学习中，在召回率和精确度之间有一个基本的权衡，这意味着我们只能以减少另一个为代价来增加一个。例如，如果我们想要找到客户流失的每一个实例——100%的召回率——那么我们将不得不接受低精度——许多误报。相反，如果我们通过提高精确度来限制假阳性，那么我们将识别出更少的降低召回率的实际搅动。
通过调整模型的阈值来改变这两者之间的平衡。我们可以在模型的 精确回忆曲线 中形象化这一点。
Precision-recall curve tuned for 75% recall.
这显示了不同阈值下的精度与召回率。Scikit-Learn 中的默认阈值是 0.5，但是根据业务需求，我们可以调整这个阈值以实现所需的性能。**
对于客户流失，我们将调整阈值以实现 75%的召回率。通过检查预测的概率(实际值)，我们确定阈值应该是 0.39，以达到这个标准。在阈值为 0.39 时，我们的召回率为 75%，准确率为 8.31%。
选择召回还是精确在于商业领域。这需要确定哪种成本更高，假阳性——预测客户会流失，而事实上他们不会——或者假阴性——预测客户不会流失，而事实上他们会——并进行适当的调整。
选择 75%的召回率作为优化示例，但这是可以改变的。在这个值上，与原始基线相比，我们在召回率上实现了20 倍的提升，在精确度上实现了8 倍的提升。
模型验证
一旦我们选择了对流失进行分类的阈值，我们就可以从维持测试集中绘制出混淆矩阵来检查预测。
Confusion Matrix for Tuned Random Forest
在这个阈值，我们识别了超过一半的搅动(75%)，尽管有大量的假阳性(右上)。根据假阴性和假阳性的相对成本，我们的模型实际上可能不是一个改进！
为了确保我们的模型解决了问题，我们需要使用维持结果来计算实现模型的回报。
验证业务价值
使用模型在保留测试集上的度量作为对新数据的性能评估，我们可以在部署这个模型之前计算部署它的价值。使用历史数据，我们首先计算流失造成的典型收入损失*，然后使用实现 75%召回率和 8%精确度的模型计算流失造成的减少的收入损失。*
对客户转化做一些假设(详见笔记本)我们得出以下结论:
机器学习增加了每月活跃订户的数量，并从客户流失中挽回了 13.5%的每月损失。
考虑到订阅成本，这相当于每月 130，000 美元。
通过这些数字，我们得出结论，机器学习已经解决了每月用户增加的业务需求，并提供了一个积极的解决方案。
作为模型解释的最后一部分，我们可以查看最重要的特征，以了解与问题最相关的变量。随机森林模型中 10 个最重要的变量如下所示:
Most important features from random forest model.
最重要的变量符合我们对问题的直觉。例如，最重要的特征是截止时间前一个月的总支出。因为我们使用 1 个月的提前期，这代表了预测月份前两个月的支出。客户在此期间消费越多，他们流失的可能性就越小。我们还看到了一些顶级功能，如平均交易间隔时间或支付方式 id ，这对我们的业务监控非常重要。**
进行预测和部署
随着我们的机器学习管道的完成和模型的验证，我们已经准备好预测未来的客户流失。我们没有这个项目的实时数据，但如果有，我们可以做出如下预测:
Predictions for new data based on threshold.
这些预测和功能重要性可以交给客户参与团队，在那里他们将努力留住成员。
除了在每次获得新数据时进行预测，我们还希望继续* 验证 我们的解决方案是否已经部署。这意味着将模型预测与实际结果进行比较，并查看数据以检查概念漂移。如果性能下降到提供价值的水平以下，我们可以收集更多的数据并对其进行训练，改变预测问题，优化模型设置，或者调整调整的阈值。*
建模注释
**与预测和特征工程一样，建模阶段适用于新的预测问题，并使用数据科学中的常用工具。我们使用的机器学习框架中的每一步都是分段的，这意味着我们能够实现众多问题的解决方案，而无需重写所有代码。此外， API s — 熊猫、功能工具和 Scikit-Learn —都是用户友好的，有很棒的文档，并且抽象掉了繁琐的细节。
机器学习过程的结论
机器学习的未来不在于一次性解决方案*，而在于 通用框架 允许数据科学家快速开发解决他们面临的所有问题的解决方案。这种脚手架的功能与网站模板非常相似:每次我们建立一个网站时，我们不会从头开始，我们使用现有的模板并填充细节。*
同样的方法应该适用于解决机器学习的问题:不是为每个问题建立一个新的解决方案，而是调整现有的脚手架，并用用户友好的工具填充细节。
在这一系列文章中，我们介绍了用于解决现实世界机器学习问题的通用框架的概念和用法。
该过程总结为三个步骤:
- ***预测工程:*定义业务需求，将需求转化为有监督的机器学习问题，并创建带标签的示例
- ***特征工程:*使用标签时间和原始历史数据为每个标签构建预测变量
- ***建模:*训练、调整业务需求、验证解决方案的价值，并使用机器学习算法进行预测
A general purpose framework for solving problems with machine learning.
*虽然机器学习不是只对少数人开放的神圣艺术，但由于缺乏标准化的流程，它仍然远离许多组织。该框架的目标是使机器学习解决方案更容易**开发和部署，*这将让更多的组织看到利用这一强大技术的好处。
如果构建有意义的高性能预测模型是您关心的事情，那么请联系我们的功能实验室。虽然这个项目是用开源的 Featuretools 完成的，但是商业产品为创建机器学习解决方案提供了额外的工具和支持。
模拟纽约芬格湖国家森林的树高和断面积
我尝试使用 R 包 randomForest 来创建两个树高和断面积的回归模型，这两个模型是基于一些激光雷达和在纽约手指湖国家森林中实地收集的数据。
声明:这个项目是我对 R 的第一次真正体验。在这个学期的早些时候，我对 R 如何工作，如何上传数据做了一些简单的学习，试图熟悉 R 的行话，等等。我还有很多东西要学，有时也非常依赖我那位出色的 GIS 教授，但我认为我对使用 R 的感觉肯定比几个月前要好。
下面我将试着概述一下我是如何完成这个项目的。关于背景和上下文，我将简单介绍一下我在 arcGIS 和 FUSION 中所做的工作，但这篇文章的重点将是使用 R 包 randomForest 来测试变量的重要性并建立回归模型。
开始了。
我开始用的数据(这部分和我在 R 里做的无关):
因此，我最初用的是一个 Excel 文件，其中包含实地收集的数据变量，其中包括每个地块的测量平均树高和断面积(这两个变量我决定用 R 建模)。对于激光雷达数据，我有一堆 LAS 格式的激光雷达文件，覆盖了整个手指湖国家森林。为了从激光雷达数据中获取我需要的度量值以用于在 randomForest 中创建模型，我必须使用激光雷达查看/处理软件工具 FUSION。通过使用 FUSION 的 grid metrics 功能，我获得了整个激光雷达数据集输出的一系列大约 60 个描述性统计数据，这些数据以栅格(grid)的形式存储在一个巨大的 CSV 文件中。我不打算尝试在我的模型中使用所有这些统计值。对于这个项目，我正在研究一些已经发表的关于该主题的科学论文的方法，这些论文的共同主题是只选择几个网格度量变量用于随机森林模型生成过程。所以，我最终保留了 Brubaker 等人论文中的五个变量的值。al 2018 使用:每个单元的平均值(每个单元是一个 15 x 15m 米的区域)，每个单元的第 30、90 和 99 百分位值，以及每个单元的覆盖百分比。这被保存为 CSV 文件，下面谈到的“grid.csv”。
然后，我将“grid.csv”加载到 arcGIS 中，将其转换为所有像元的 shapefile，然后在空间上将它连接到我收集的关于地块的野外数据，这就生成了下面将要讨论的文件“FieldandLidarData.csv”。
我实际上最终在 R 中使用的数据(R 代码中出现的文件的名称在引号中):
- “grid . csv”-一个 CSV 文件，其值超过 900，000(！！)15 x 15m 米激光雷达像元。关键值包括:每个 15 x 15m 米激光雷达像元的 UTM 坐标，以及我选择用于我的模型的 5 个格网度量统计数据:每个像元的平均值、每个像元的第 30、90 和 99 百分位值，以及每个像元的覆盖百分比。来源:我在 FUSION 中基于最初为 Brubaker 等人 2018 年的研究收集的激光雷达数据创建了这个文件。
- “fieldandlidardata . csv”-一个 CSV 文件，其中包含手指湖国家森林中约 100 个场地地块的值。关键值包括:田间收集的每块地的平均树高和断面积，以及与“grid.csv”文件中相同的 5 个网格度量统计的每块地的值。来源:这是我通过将“grid.csv”中的信息与 arcGIS 中的绘图数据进行空间连接而生成的文件。
因此，基本上这两个文件包含几乎完全相同的属性，除了覆盖不同区域的值:“grid.csv”具有整个国家森林的值，而“FieldandLidarData.csv”仅覆盖整个森林中随机放置的 100 个小研究地块的区域。这两个文件之间的另一个关键区别是,“FieldandLidarData.csv”文件包含野外采集的每个地块的平均树高和断面积值，而“grid.csv”文件不包含整个研究区域的这些值……这些是我试图使用 randomForest 获得的值！
我试图回答的关键问题是:
整个手指湖国家森林的树高和断面积是什么样的？
R 中的角色进展如何:
我把上面提到的两个文件保存到我电脑上的一个文件夹中，所以我只是在那里设置了目录，安装了 randomForest 包，然后去了镇上。注意:我将只显示断面积或树高中的一个的代码和图表，因为我实际上对它们两个都做了完全相同的步骤。
我做的第一件事是查看我用来计算断面积和树高的 5 个变量的重要性。我一次运行一个模型。我首先用我的第 21 行代码种植森林，见下文。数据“thetrees”是“FieldandLidarData.csv”文件的列，其中包含我的 5 个网格度量变量的值。然后我打印了结果，显示我种植了 500 棵树的回归森林，模型解释的变异百分比是树高 48.02%，断面积 26.61%，百分比不是很大，但也不错。然后，我用我的第 23-25 行代码生成了一个变量重要性图，以查看 5 个变量中哪个比其他变量更重要，下面显示了树高变量的图片。
Code for the variable importance bit for both basal area and height.
The variable importance plot I got for tree height. Clearly the 30th percentile variable didn’t matter…
接下来的议程是查看一些模拟值与每个地块的基础面积和高度的实际值的快速图。肯定有很多分散的情况，但是直观地看到模型中的东西与实际收集到的东西相比有多好是很好的。
The plot I made for tree height (it looks better than the basal area one did so I’m using it).
The code I used to generate my modeled vs. actual height plot.
最后，我进去把我的树高和断面积模型应用到整个手指湖国家森林！因此，我调用了一个巨大的“grid.csv”文件，其中包含整个森林的 UTM 坐标和网格度量数据，然后使用我之前生成的模型来预测研究区域中每个 15 x 15 m 地块的高度和断面积值。然后我将这些数据保存在一个 CSV 文件中，耶！我说完了。
My code for making my height model for the entire study area. I used the head function a few times just to see how things were going before I generated another massive CSV file.
在能够查看我的模型后，我的一些模式和想法:
如果你好奇的话，下面是一些在 arcGIS 中制作的小面积森林的地图。我想调查手指湖国家森林周围的树高和断面积值，通过显示我在 R 中得到的预测值，我能够做到这一点，看起来真的很酷。
a) is satellite imagery, b) shows how old the forest is in the area, c) shows the predicted tree height values per cell and d) shows basal area.
一般来说，树高不会超过 25 米，而断面积在 2 米/公顷以下。总的来说，模型似乎与我们在卫星图像中看到的高大/古老的树木非常匹配，而且看起来高度和断面积彼此之间以及与林龄呈正相关——我最终也研究了断面积和高度如何随林龄变化，因为手指湖国家森林的一些部分比其他部分的森林覆盖时间更长。
如果我有更多的时间，我想我会更深入地研究树高和断面积是否真的与林龄有关。我还会更多地考虑我使用了哪些变量以及使用了多少变量，以便尝试在数据中获得更多的解释变量，特别是对于基底面积，因为我们看到的变量中只有四分之一多一点被模型解释。
参考资料:
布鲁贝克，K. M .，Q. K .约翰逊和 M.W .凯耶。(2018)."用离叶和开叶激光雷达研究落叶林中乔木和灌木生物量的空间格局."加拿大森林研究杂志。48(9): 1020–1033.
墨菲，P. W. (2018)。未公布数据。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
