支持向量机
分离，但支持
历史
这篇文章是关于支持向量机的介绍，理解数学直觉，正则化，在代码中实现概念，然后了解它的应用领域。好吧，那就系紧安全带，让我们进入这个美丽的概念支持向量机(SVM)的仙境。机器学习领域的先驱之一 Vladimir Vapnik 在 20 世纪 60 年代初发现了支持向量机。但他无法从计算上证明这种算法的魔力，因为他在 20 世纪 60 年代缺乏所需的资源。正是在 20 世纪 90 年代初，贝尔实验室邀请他前来美国进行相关的研究工作。早在 1993 年，在字符数字识别领域就有一项快速跟踪研究。Vapnik 打赌，他的支持向量机工作在数字识别方面可以比神经网络做得更好。因此，我的一位同事尝试应用 SVM，结果出乎意料地创造了奇迹。最后，这项在 20 世纪 60 年代发现的工作在 20 世纪 90 年代得到了认可，当时它实际上被实施以给出更好的计算结果。
简介
支持向量机是一种用于分类和回归建模的监督机器学习算法。支持向量机的目的是识别空间中的样本，并根据类别对其进行分离。例如，如果我给你一个水果图片列表，那么算法会将输入的图片分类为苹果或草莓。记住，在上面的例子中，分类是基于类很容易分离的假设。如果图像的分布方式使得很难通过画直线来划分类别，那会怎么样呢？在我们朝这个方向分析之前，我们需要熟悉一些基础知识。
让我们进一步了解几个技术术语:
**超平面:**分隔在给定空间中分布的正负类的平面。
Margin: 在超平面的两边平行画出的一条分隔线，将正负类分开。
**支持向量:**距离位于超平面任一侧的边缘最近的正负点。
作者图片
棘手的部分！！
我确信你的脑海中可能会出现一个问题，为什么我们不能在任何其他方向上找到一个分离的超平面和边缘。是的，这是可能的，但是我们必须确保我们选择超平面的方式使得超出边界的点的距离是最大的。为什么必须是最大值？？只是为了避免该区域中的点的错误分类。分离边缘越宽，将获得更准确的分类。下图 1 显示了边距之间有足够的间隔，以避免新点的错误分类，并正确地将它们分配到各自的类中。图 2 示出了在这种方向上的超平面可能导致边缘之间的间隙不足，并最终导致错误分类的点。
作者图片
作者图片
最大化利润—数学直觉
作者图片
考虑一个向量 ŵ 在垂直于超平面的方向，如图所示。假设我们有一个未知点 u 位于超平面中。**。**由于我们有一个垂直于超平面的向量，我们可以考虑 ŵ 和 c. 任何未知点的内积，该点的存在是未知的，即，对于我们的数学假设，它是正点还是负点，我们认为它的点积至少大于常数 c.
ŵ∙u>= c
假设常数 — c = b
ŵ∙u—c>= 0
ŵ∙u+b>= 0—(1)
等式 1 是我们的决策规则，我们将在后面的阶段利用它。
假设平面右边的点是正的，表示为+1。
平面左侧的点为负，表示为-1。
因此，使用等式 1 中的判定规则，我们有以下两个等式，分别表示正负点。
ŵ∙(x+)+b>= 1—(2)
ŵ∙(x-)+b<=-1—(3)
现在为了数学上的方便，我们引入一个特征变量 y ,它满足上面的方程
ŵ∙(x+时 y = 1)+b>=+1—(4)
ŵ∙(x-时 y =-1)+b<=-1—(5)
从上面两个等式我们可以得出结论:将 Y 乘以等式(2)和(3) 任意点被正确分类的必要条件是下面给出的。
*y (ŵ∙x+b)>= 1—(6)
重新排列上面的等式，并考虑到该点仍然未知，我们可以将其等于零，得到如下
*y (ŵ∙x+b)—1 = 0—(7)
现在我们想要最大化两个边缘之间的距离，因此我们画出正点和负点边缘的正交向量 ŵ 。我们在边际上的支持向量可以认为是 X+ 和 X- 。我们的目标是计算边距上两个支持向量之间的最大距离。然而，我们绘制了一个范数为 w 的范数，这样距离( X+) — (X -)可以通过在范数 w 上描绘距离来计算。因此，将 w 的单位向量与距离相乘实际上不会改变值。
边距宽度=(x+)—(x-)=[(x+)—(😆]*(ŵ/| |ŵ| |)
作者图片
现在从方程(7)可以计算出 X+ 和 X- 的值。
通过将 X+的 Y = 1，我们得到 X+的值如下
*1 (ŵ∙(x+)+b)—1 = 0
(x+)=(1—b)/ŵ—(8)
对于 X，Y = 1，我们得到 X 的值如下
(-1)*(ŵ∙(x-)+b)—1 = 0
(x-)=(-1—b)/ŵ—(9)
将等式(8)和(9)代入页边距宽度计算:
边距宽度= (X+) — (😆
=[(x+)—(😆]*(ŵ/| |ŵ| |)
=(1-b+1+b)(1/ŵ)(ŵ/| |ŵ|)
*= 2 (1/| |ŵ|)
*我们的目标是最大化边距宽度= max(2 (1/| |ŵ|)
*以 y (ŵ∙x+b)—1 = 0 为条件
好了，我们终于得出了页边距宽度的公式。
结束了吗？？不，不是真的！！！
我们在这里优化公式，使其一般化。为了优化公式，我们必须确保它达到凸性！！如果你不知道最优化的凸性，那也没关系。现在，如果没有数学上的便利，我们就无法实现最优化。
因此，我们将页边距宽度的等式重新定义为
min()*(| |ŵ| | ^2)
*以 y (ŵ∙x+b)—1 = 0 为条件
使用拉格朗日乘数优化
现在我们知道，当我们有一个目标函数和约束，然后我们可以利用拉格朗日乘数法，找到变量的拉格朗日依赖关系。所以利用目标函数和约束条件，我们可以构造拉格朗日函数。
l=()* | |ŵ| | 2—∑αI [易(ŵ∙x+b)—1]—–(1)
现在根据拉格朗日函数
∂L / ∂ŵ = 0
∂L / ∂b = 0
在我们的主函数(1)中应用同样的方法，我们得到
∂L / ∂ŵ = 0
ŵ — ∑ αi 易 xi = 0
ŵ = ∑ αi 易 Xi—–(2)
∂L / ∂b = 0
∑αI * yi = 0——-(3)
应用(1)中的(2)和(3)
*L = (∑αI * yi * Xi)(∑αj * yj * XJ)—(∑αI * yi * Xi)*(∑αj * yj * XJ)—b ∑αI * yi+∑αI
*L =∑αI—∑αI αj * yj * yj * Xi * XJ
所以我们最后得出结论，边距宽度的优化因子取决于任意两个样本或点 x.
线性可分与否？
回到我们的第一幅图，可以看到，通过画一条直线，积极的和消极的类别被如此干净利落地分开了。这种类型的分离被称为线性可分离变量。有没有可能一直通过画一条直线把两个类分开？如果我们在试图画一条直线时有大量错误分类的点，该怎么办？考虑下图，我们可以看到这些点是不容易分开的。在这种情况下，我们使用内核。核的作用是将任何给定维度空间中的点转换到更高维度的空间。例如，在下图中，我们有 2-D 中的点，通过应用内核，我们看到它们被转换为 3-D。令我们惊讶的是，我们看到现在通过绘制一条线，这些点可以线性分离。
图片由阿拉什·赛义德·普尔在研究之门上拍摄
正规化
想象一下这样一种情况，边距宽度变得最大，但同时出现了点的两个错误分类。在这种情况下我们会怎么做？我们会尝试重新排列超平面以避免错误分类吗？如果你仔细考虑，答案必须是否定的，因为我们在这里优化和概括的情况。虽然我们可能会错误地分类一些点，但我们不能每次都通过改变超平面来折衷边缘宽度。在这种情况下，我们将应用正则化的概念。这里，我们将在超平面的对齐和不允许错误或错误分类之间找到惊人的平衡，同时尝试最大化余量。我们包括两个变量误差, c 给我们错误分类的总点数和ɛ给我们错误分类的每个误差的值。这有助于我们实现规则化效果，即告诉 SVM 不要调整或包含异常值，也有助于我们计算 c 的值，从而在使模型免受某些错误分类的影响和同时确保边距宽度最大化之间找到一个惊人的平衡。
正则化= min()*(| |ŵ| | ^2)+∑ɛi
代码实现
使用 SVM 的虹膜分类是理解 SVM 分类器概念的最好例子之一。鸢尾是一种大约有 250-300 种的花。在我们的例子中，我们有一个由花的细节组成的数据集，例如萼片宽度、萼片长度、花瓣宽度、花瓣长度。我们的目标是使用 SVM 分类器，基于数据集中可用的细节，将花分类为正确的种类。
作者图片
我们使用 sklearn 库中可用的方法将数据分成训练和测试数据。此外，我们从 sklearn 库中导入 SVC 分类器模型，并在训练数据上训练一个模型。
在集合上训练模型之后，我们然后将预测方法应用于测试数据集。使用各种准确性度量标准，我们然后使用混淆矩阵和分类报告检查并验证模型输出。
混乱矩阵作者的图片
作者图片
SVM 的应用
字符电子识别
生物信息学
图像分类
人脸检测
感谢您花时间阅读这篇文章！！
支持向量机(SVM)清楚地解释:一个 python 教程的分类问题与三维绘图
在这篇文章中，我解释了支持向量机的核心，为什么以及如何使用它们。此外，我还展示了如何在 2D 和 3D 中绘制支持向量和决策边界。
手工制作草图作者。SVM 插图。
介绍
每个人都听说过著名且广泛使用的支持向量机 (SVMs)。最初的 SVM 算法是由 Vladimir N. Vapnik 和 Alexey Ya 发明的。1963 年的切尔沃嫩基斯。
支持向量机是监督的机器学习模型，通常用于分类 ( SVC —支持向量分类)或回归 ( SVR —支持向量回归)问题。根据目标变量(我们希望预测的)的特征，如果我们有一个离散目标变量(例如分类标签)，我们的问题将是一个分类任务，或者如果我们有一个连续目标变量(例如房价)，我们的问题将是一个回归任务。
支持向量机更常用于分类问题，因此，在本文中，我将只关注支持向量机模型。
方法的核心
在本文中，我不打算详细介绍算法的每一步(因为有大量的在线资源),而是解释关于支持向量机的最重要的概念和术语。
1.决策边界(分隔超平面)
SVCs 的目标是找到最佳 超平面(也称为决策 边界)使得最佳 将(分裂)一个数据集分成两个类/组(二元分类问题)。
根据输入特征/变量的数量，决策 边界可以是一条线(如果我们只有 2 个特征)或一个超平面(如果我们的数据集中有 2 个以上的特征)。
为了获得主要思想，考虑如下:每个观察(或样本/数据点)被绘制在一个 N 维空间中，其中**N**是我们数据集中特征/变量的数量。在那个空间中，分离 超平面的是一个 (N-1)维子空间。
超平面是 N - 维-空间的一个( N-1 )- 维-子空间。
所以，如前所述，对于一个二维空间，决策，边界将只是一条线，如下图所示。
手工制作草图由作者制作。使用只有两个特征(即 x1 和 x2)的数据集来说明 SVM 分类模型(SVC)的决策边界。决策边界是一条线。
数学上，我们可以定义决定边界为:
作者写的渲染 latex 代码。
如果你想在交互式路线图和活跃的学习社区的支持下自学数据科学，看看这个资源:https://aigents.co/learn
2.支持向量
支持 向量就是距离分离 超平面最近的样本(数据点)。如果这些样本被移除，它们将改变分离超平面的位置。由此可见，这些就是最重要的** 最重要的 样本即定义****位置和方位的最佳** 决定 边界。****
手工制作草图作者。在这个玩具二维 SVM 问题中，用紫色圈起来的点代表支持向量。
3.硬边缘:SVM 如何找到最佳超平面？
几条不同的线(或者一般来说，不同的决策边界)可以分隔我们的类。但是哪一个是最好的呢？
****手工制作草图由作者制作。此图显示了分隔两个类别的 3 个候选决策边界。
超平面与最近的数据点(样本)之间的距离称为余量。目标是在超平面和任何支持** 向量之间选择一个具有最大可能** 余量的超平面。** SVM 算法寻找最佳决策边界，如边际最大化。这里最好的线是黄线，如下图所示。******
手工制作草图作者。最佳分隔线是最大化边距(绿色距离)的黄色分隔线。
在摘要中，支持向量机选择最大化到支持向量的距离的决策边界。以最大化到支持向量的距离的方式绘制决策** 边界。如果决策边界太靠近支持向量，那么它将对噪声敏感并且不能很好地概括。**
4.关于软边界和 C 参数的一个注记
有时，我们可能想要(故意)允许一些误差**(错误分类)。这就是“软 边距”背后的主旨。软边界实现允许一些样本被错误分类或者位于决策边界的错误一侧，从而允许高度一般化的模型。**
一个软余量SVM解决了下面的优化问题:
- 增加决策边界的距离到支持向量 s(即余量)并且
- 最大化训练集中正确 分类的点数。
很明显，这两个优化目标之间有一个的权衡。这种权衡由著名的 C 参数控制。简而言之，如果 C 为小，则误分类数据点的罚分为低，因此选择具有大 裕量的决策 边界，代价是更大数量的误分类。如果 C 是大， SVM 试图最小化误分类样本的数量，并导致判定** 边界具有较小 余量。******
5.当没有清晰的分离超平面(核 SVM)时会发生什么？
如果我们有一个数据集是线性可分的，那么支持向量机的工作通常很容易。然而，在现实生活中，在大多数情况下，我们手头都有一个线性不可分的数据集，这正是内核技巧提供一些魔力的时候。
内核 技巧将原始数据点投影到一个更高维的空间中，以便使它们可以线性分离(在那个更高维的空间中)。
因此，通过使用核技巧，我们可以使我们的非线性可分的数据，在一个更高维的空间中线性可分。
内核技巧是基于一些内核 函数来测量样本的相似性。这个技巧实际上并没有将数据点转换到一个新的高维特征空间，显式地**。核 SVM 根据高维特征空间中的相似性度量来计算决策边界，而不实际进行投影。一些著名的核函数包括线性、多项式、径向基函数(RBF)和 sigmoid 核。**
手工制作草图作者。内核技巧。在原始空间中，数据不是线性可分的，但是在投影到更高维空间后，它们是线性可分的。
在 scikit-learn 中使用 Iris 数据集和线性 SVC 模型的 Python 工作示例
提醒:虹膜数据集由 150 个花样本组成，每个样本具有 4 个特征/变量(即萼片宽度/长度和花瓣宽度/长度)。
2D
让我们在 2D 绘制决策边界(我们将只使用数据集的 2 个特征):
from sklearn.svm import SVC
import numpy as np
import matplotlib.pyplot as plt
from sklearn import svm, datasetsiris = datasets.load_iris()# Select 2 features / variables
X = iris.data[:, :2] # we only take the first two features.
y = iris.target
feature_names = iris.feature_names[:2]
classes = iris.target_namesdef make_meshgrid(x, y, h=.02):
x_min, x_max = x.min() — 1, x.max() + 1
y_min, y_max = y.min() — 1, y.max() + 1
xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))
return xx, yydef plot_contours(ax, clf, xx, yy, **params):
Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])
Z = Z.reshape(xx.shape)
out = ax.contourf(xx, yy, Z, **params)
return out# The classification SVC model
model = svm.SVC(kernel="linear")
clf = model.fit(X, y)fig, ax = plt.subplots()# title for the plots
title = (‘Decision surface of linear SVC ‘)
# Set-up grid for plotting.
X0, X1 = X[:, 0], X[:, 1]
xx, yy = make_meshgrid(X0, X1)plot_contours(ax, clf, xx, yy, cmap=plt.cm.coolwarm, alpha=0.8)
ax.scatter(X0, X1, c=y, cmap=plt.cm.coolwarm, s=20, edgecolors="k")
ax.set_ylabel("{}".format(feature_names[0]))
ax.set_xlabel("{}".format(feature_names[1]))
ax.set_xticks(())
ax.set_yticks(())
ax.set_title(title)
plt.show()
以上 python 代码的输出。作者创作的赋格曲。
在虹膜数据集中，我们有 3 类花和 4 个特征。这里，我们仅使用 2 个特征(因此我们有一个 2 维特征空间)T1，并且我们绘制了线性 SVC 模型的决策边界。点的颜色对应于类别/组。
三维（three dimension 的缩写）
让我们在 3D 中绘制决策边界(我们将只使用数据集的 3 个特征):
from sklearn.svm import SVC
import numpy as np
import matplotlib.pyplot as plt
from sklearn import svm, datasets
from mpl_toolkits.mplot3d import Axes3Diris = datasets.load_iris()
X = iris.data[:, :3] # we only take the first three features.
Y = iris.target#make it binary classification problem
X = X[np.logical_or(Y==0,Y==1)]
Y = Y[np.logical_or(Y==0,Y==1)]model = svm.SVC(kernel='linear')
clf = model.fit(X, Y)# The equation of the separating plane is given by all x so that np.dot(svc.coef_[0], x) + b = 0.# Solve for w3 (z)
z = lambda x,y: (-clf.intercept_[0]-clf.coef_[0][0]*x -clf.coef_[0][1]*y) / clf.coef_[0][2]
tmp = np.linspace(-5,5,30)
x,y = np.meshgrid(tmp,tmp)fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot3D(X[Y==0,0], X[Y==0,1], X[Y==0,2],'ob')
ax.plot3D(X[Y==1,0], X[Y==1,1], X[Y==1,2],'sr')
ax.plot_surface(x, y, z(x,y))
ax.view_init(30, 60)
plt.show()
以上 python 代码的输出。图由作者生成。
在虹膜数据集中，我们有 3 类花和 4 个特征。这里我们只使用了 3 个特征(所以我们有一个三维特征空间**)和只有 2 个类(二进制分类问题)。然后，我们绘制了线性 SVC 模型的决策边界。点的颜色对应于 2 个类别/组。**
绘制支持向量
import numpy as np
import matplotlib.pyplot as plt
from sklearn import svm
np.random.seed(2)# we create 40 linearly separable points
X = np.r_[np.random.randn(20, 2) — [2, 2], np.random.randn(20, 2) + [2, 2]]
Y = [0] * 20 + [1] * 20# fit the model
clf = svm.SVC(kernel=’linear’, C=1)
clf.fit(X, Y)# get the separating hyperplane
w = clf.coef_[0]
a = -w[0] / w[1]
xx = np.linspace(-5, 5)
yy = a * xx — (clf.intercept_[0]) / w[1]margin = 1 / np.sqrt(np.sum(clf.coef_ ** 2))
yy_down = yy — np.sqrt(1 + a ** 2) * margin
yy_up = yy + np.sqrt(1 + a ** 2) * marginplt.figure(1, figsize=(4, 3))
plt.clf()
plt.plot(xx, yy, "k-")
plt.plot(xx, yy_down, "k-")
plt.plot(xx, yy_up, "k-")plt.scatter(clf.support_vectors_[:, 0], clf.support_vectors_[:, 1], s=80,
facecolors="none", zorder=10, edgecolors="k")
plt.scatter(X[:, 0], X[:, 1], c=Y, zorder=10, cmap=plt.cm.Paired,
edgecolors="k")
plt.xlabel("x1")
plt.ylabel("x2")
plt.show()
以上 python 代码的输出。图由作者生成。
双 - 圆圈 点代表支持 向量。
最新帖子
使用可从《先知脸书》公开获得的预测模型预测股票价格
towardsdatascience.com](/time-series-forecasting-predicting-stock-prices-using-facebooks-prophet-model-9ee1657132b5) [## 用新冠肺炎假设的例子解释 ROC 曲线:二分类和多分类…
在这篇文章中，我清楚地解释了什么是 ROC 曲线以及如何阅读它。我用一个新冠肺炎的例子来说明我的观点，我…
towardsdatascience.com](/roc-curve-explained-using-a-covid-19-hypothetical-example-binary-multi-class-classification-bab188ea869c) [## PCA 清楚地解释了——如何、何时、为什么使用它以及特性的重要性:Python 指南
在这篇文章中，我解释了什么是 PCA，何时以及为什么使用它，以及如何使用 scikit-learn 在 Python 中实现它。还有…
towardsdatascience.com](/pca-clearly-explained-how-when-why-to-use-it-and-feature-importance-a-guide-in-python-7c274582c37e) [## 关于 Python 中的最小-最大规范化，您需要知道的一切
在这篇文章中，我将解释什么是最小-最大缩放，什么时候使用它，以及如何使用 scikit 在 Python 中实现它
towardsdatascience.com](/everything-you-need-to-know-about-min-max-normalization-in-python-b79592732b79) [## Scikit-Learn 的标准定标器如何工作
在这篇文章中，我将解释为什么以及如何使用 scikit-learn 应用标准化
towardsdatascience.com](/how-and-why-to-standardize-your-data-996926c2c832)
敬请关注并支持这一努力
如果你喜欢并发现这篇文章有用，请关注我！
有问题吗？把它们作为评论贴出来，我会尽快回复。
参考
[1]https://www.nature.com/articles/nbt1206-1565
https://en.wikipedia.org/wiki/Support_vector_machine
[2]https://sci kit-learn . org/stable/modules/generated/sk learn . SVM . SVC . html
和我联系
- 领英:【https://www.linkedin.com/in/serafeim-loukas/】T4
- 研究之门:https://www.researchgate.net/profile/Serafeim_Loukas
- https://people.epfl.ch/serafeim.loukasEPFL简介 : 美国
- 堆栈 溢出:https://stackoverflow.com/users/5025009/seralouk**
惊奇和熵的故事
今天，我从一个有趣的问题开始我的小文章。假设有一天，外面正下着雨，你走出家门，发现你的露天庭院完全干涸了。你的第一反应会是什么？
惊喜！
为什么？
因为，当外面下雨时，你已经知道你的院子会是湿的。如果你已经看到了这一点，那么你就不会“收到”任何新的信息，因此就不会有惊讶的成分。但是既然你发现它没有湿，你就“接收”了一些信息。这让人感到意外。
现在，我们可以试着在这种情况下更正式一点。假设你已经知道庭院有 0.1%的可能是干的，0.9%的可能是湿的。*(你注意到价值观的选择了吗？我选择它们的方式是，如果你把它们加起来，值就是 1.0。你能解释一下吗？)*所以，在这个设置中，我们有两个变量。
- x(庭院干燥的可能性)
- y(庭院潮湿的可能性)
现在，让我问这个问题(即使你没有完全理解它的意思，请耐心听我说)变量 X 的“熵”是什么？
在我们做数学之前，让我告诉你熵是变量不确定性的度量。越不确定，熵越高。那么基于这种直觉，你认为变量的熵会是多少？在这里停下来想一想。
好的。现在，让我们做数学。(如果对 Claude Shannon 的原论文感兴趣，可以查看一下这里我们说一个变量的总熵等于——
其中 p(x)代表事件的概率。
下面来一点 Python:)
我们知道 X 的概率是 0.1 所以，
In [1]: import mathIn [2]: (0.1) * round(math.log2(0.1), 2)
Out[2]: -0.332
根据同样的定律，Y 的熵是—
In [3]: (0.9) * round(math.log2(0.9), 2)
Out[3]: -0.135
正如所料，X 的熵大于 Y 的熵(如果我们暂时忽略负号的话)。所以，我们可以从这个例子中看到，一个事件发生的可能性越小，这个事件的熵就越大。
如果我们愿意，我们也可以从这个角度来考虑——我们对一个事件的信息越少，也就是说，它越令人惊讶，这个事件的熵就越高。
现在，我们要问另一个问题。系统的总熵是多少？其实很简单。我们是这样衡量的—
系统的总熵
所以，基本上我们要把 X 和 Y 的值加起来，然后取反。
In [4]: - (-0.332 + (-0.135))
Out[4]: 0.467
我们用称为“比特”的单位来表示这个量。所以我们可以说这个系统有 0.467 比特的熵。因此，你可以从这个系统接收那么多的信息。
为了完成这篇文章，我们要测量一次公平抛硬币的熵。它有多少比特的熵？
在公平抛硬币的情况下，我们可以像下面这样列出正面和反面的概率
- P(H) = 0.5
- P(T) = 0.5
所以，应用上面的方法，我们可以计算熵。它将是 1 比特。
In [5]: - ( (0.5 * math.log2(0.5)) + (0.5 * math.log2(0.5)))
Out[5]: 1.0
这并不奇怪！
希望你喜欢这个关于熵和信息论的小介绍。这个小等式支撑着我们现在的整个通信系统。有个想法是好事。
如果你喜欢，请按拍手图标多次。这会鼓励我写更多。此外，你可以在这里关注我的未来文章，撰写关于信息论、复杂性理论、机器学习、算法、离散数学等方面的文章。
有关更多信息(！)和进一步的阅读—
基于成本的路线优化结果令人惊讶
基于真实世界成本的路线优化。图片由 Free-Photos 发自 Pixabay 。
基于成本的路线优化不同于仅考虑距离或时间的路线优化。
基于成本的路线优化在现实世界中很有用，因为它符合任何试图降低成本和提高效率的公司的业务目标。
相比之下，仅考虑距离的路线优化提供了一个有吸引力的假设结果-地图上看起来不错的路线。但仅此而已。
这很容易用几个例子来演示。
仅基于距离成本的路线优化
考虑以下 3 辆车和 50 个地点的交付时间表，该时间表需要针对伦敦周边的当天交付进行优化。
3 辆车，50 站距离成本仅路线优化由 Optergon 提供
车辆成本用 F/D/T ( 固定，距离&时间)表示。每辆车的基于时间的成本被设置为 0，这意味着时间不会影响系统的总成本(这显然不适用于现实世界中的送货公司)。距离成本为 1 美元/公里。
优化该交货计划(可能包括任何类型的交货和提货)会产生以下结果。
3 辆车，50 站距离成本唯一路线优化解决方案由 Optergon 提供
这看起来是一个相当不错的结果。的确如此。因为它只有基于距离的成本数据可以处理。
基于仅距离路线优化来优化其车队和交付时间表的公司永远不会怀疑他们没有利用最有效的路线优化解决方案。
从本质上说，他们把钱留在桌子上而不知道，因为没有办法通过简单地看结果来判断是否有更低成本的解决方案。
基于距离和时间成本的路线优化
让我们稍微改变一下场景，为每辆车在路上花费的时间赋值。毕竟，司机的工资是按小时支付的，在运输、物流和配送公司的总成本中占很大一部分。
3 辆车，50 站距离和时间成本路线优化由 Optergon 提供
每辆车都有相同的基于时间的成本，每小时 10 美元。在优化该时间表之后(除了基于时间的成本的变化之外是相同的)，我们获得了以下结果。
由 Optergon 提供的 3 辆车、50 站距离和时间成本路线优化解决方案
长得很像。事实上，鉴于每辆车基于时间的成本是相同的，这正是我们所期望的，对吗？毕竟，如果每辆车的成本完全相同，就不应该有变化。
除了，这里有一个微妙的区别。
返回到仅距离基于成本的路线优化解决方案，并查看总时间(显示在地图上方的*)。现在将它与基于距离和时间的优化解决方案的总时间进行比较。*
花的时间更少。大约少了 4 分钟。
最优路线被稍微改动了一下，以减少总的时间为代价，以少量的距离为代价，使优化的总成本最小化。
这里或那里的几分钟可能看起来不多，但这是一个非常小的优化，只有微小的差异。即使是现在，随着时间的推移，成本的降低也会累积起来，从而带来显著的节约。
如果我们远离这个简单的例子，即我们小型车队中的所有车辆都产生相同的成本，那么考虑基于时间和距离的成本的重要性就变得更加明显。
基于可变时间和距离成本的路径优化
让我们假设我们的一辆车需要一名司机和一名合作伙伴来完成一些交付——基本上是那辆车基于时间的成本的两倍。
由 Optergon 提供的 3 辆车、50 个停车距离和可变时间成本路线优化
蓝色皮卡 III 现在的时间成本是每小时 20 美元，是其他两辆车的两倍。在此交货计划的所有其他方面保持不变的情况下，这是优化的结果。
由 Optergon 提供的 3 辆车、50 站距离和可变时间成本路线优化解决方案
同样，这看起来是合理的。仔细观察，您会注意到蓝色皮卡 III ，这款基于时间的运营成本翻了一番的车辆最终只使用了很短的时间。
这听起来可能很合理，因为我们凭直觉理解业务目标是降低车队的总体成本。
对于任何只考虑基于距离的成本的系统来说，这种类型的解决方案是完全不可能的。只有距离成本的系统仍然会返回类似于我们的初始路径优化解决方案的东西，该解决方案是时间成本不可知的。
基于真实世界时间和距离成本的路线优化
现代公司运营着拥有许多不同车辆类型的多样化车队。这些不仅有不同的运载能力，他们也有不同的相关成本每距离。
让我们假设蓝色皮卡 I 是一辆稍微旧一点的车，不太省油。
3 辆车，50 站真实世界成本路线优化由 Optergon 提供
蓝色皮卡 I 的行驶距离成本为 1.20 美元——相比之下，两款较新的车型每公里成本仅为 1 美元。
结果非常有趣，可能不是你第一眼看到的那样。
3 辆车，50 站现实世界成本路线优化解决方案由 Optergon 提供
蓝色皮卡 I 实际上花了很长的时间送货。事实上，只比日程限制规定的 8 小时少了 10 分钟——几乎与最便宜的汽车蓝色皮卡 II 花费的时间相同。
然而，他们的优化路线有很大的不同。
仔细观察每个人走过的距离。蓝色皮卡 I 行驶了 68 公里多一点，而蓝色皮卡 II 行驶了两倍多，才刚刚超过 162 公里。
仅使用时间(或仅使用距离)进行交付规划的路线优化没有必要的数据来确保每辆车的行为方式降低车队的总体成本。
在这种情况下，前两辆皮卡不得不花很长时间停车，以防止昂贵的时基车(蓝色皮卡 III )抬高成本。
与此同时，第一辆皮卡需要尽可能少地行驶，以防止其膨胀的距离成本推高车队的整体成本。
这突出了不考虑基于时间和距离的成本之间复杂相互作用的路线优化的假设交付计划与降低现实世界中发生的总成本的优化解决方案之间的差异。
你可以在 Optergon 使用一个月的免费路线优化来尝试为自己的车队制定基于成本的交付计划。
最初发表于https://optergon.com。
用 Python 命名匹配的惊人有效的方法
数据匹配、模糊匹配、重复数据删除
姓名匹配问题先睹为快，图片作者。
R 最近我偶然发现了这个数据集，在这里我需要分析数码产品的销售记录。我得到的数据集差不多有 572000 行和 12 列。我非常兴奋能够处理如此大的数据。我兴致勃勃地快速查看数据，发现同一个名字重复走不同的行。(啊！数据清理时间太长了！).
在一些产品栏上，它包含 iphone，而在另一个栏上，它写着 Iphone，或 iphone 7 +和 iphone 7 Plus，对于一些常规客户列表，有些有 Pushpa Yadav，而在其他 Pushpa Yadav(名称是伪的)。
先睹为快的问题视图(不是来自原始数据集)，作者的图片。
这些是相同的产品名称和客户名称，但采取了不同的形式，即处理相同名称的不同版本。这类问题是数据科学家在数据分析过程中要处理的常见场景。这个场景有一个名字叫做数据匹配或模糊匹配(概率数据匹配)或简称为重复数据删除或字符串/名称匹配。
为什么会有“不同但相似的数据”？
常见的原因可能是:
- 数据输入时出现打字错误。
- 缩写用法。
- 系统数据输入没有被很好地验证以检查这种错误。
- 其他人。
无论如何，作为数据科学家或分析师，我们有责任匹配这些数据，以创建主记录供进一步分析。
所以，我草草记下了解决问题的行动:
1.手动检查并解决。
2.寻找有用的库/资源，这些库/资源是由社区的精英们分享的。
对于这样一个扩展的数据集(572000* 12)来说，第一个选择确实很麻烦，所以我开始查看不同的库，并用模糊匹配器(conda 现在没有)尝试了这些方法，fuzzywuzzy 在核心上使用了 Levenshtein distance 和 difflib 。
然而，在使用它们的过程中，我发现对于这样的大数据，这是非常耗时的。因此，我需要寻找一种更快更有效的方法。
经历了许多天的挫折后，我终于知道了克里斯·范登贝格分享的解决方案。
我们将介绍
- ngram
- …向量化…
- TF-IDF
- 余弦相似度与稀疏点 topn:领先果汁
- 工作演示(代码工作)
因为我们的目标不仅仅是匹配字符串，而且是以更快的方式匹配。因此，ngram 的概念，具有余弦相似性的 TF-IDF 开始发挥作用。
在进入工作演示(代码工作)之前，让我们了解一下基础知识。
N-grams
n 元语法广泛用于文本挖掘和自然语言处理，它是给定句子或(单词文件)中一组共现的单词。为了找到 n-gram，我们向前移动一个单词(可以根据需要移动任何一步)。
例如，对于房型*“标准房海景威基基塔”*
如果 N=3(称为三元组)，那么 N 元组(三元组)将是:
- 标准客房海洋
- 海景客房
- 怀基基海滩海景
- 观看怀基基塔
公式为 Ngrams，图片作者。
从研究论文中探索更多关于 ngrams 的信息
为什么 n=3？
可能会产生一个问题，为什么 n = 3 可以被看作 n= 1(一元)或 n=2(二元)。
这里的直觉是，与单个词相比，双词和三词可以捕捉上下文信息。例如，当独立观察时，“房间海景”比仅仅“房间”、“海洋”和“景色”有更多的含义。
TF-IDF
TF-IDF 代表术语频率-逆文档频率，而 TF-IDF 权重是信息检索和文本挖掘中经常使用的两个主要关注的权重。
1.用于评估一个单词对集合中的文档有多重要
2.重要性随着单词出现的次数成比例增加
TF-IDF 背后的两个想法
词频(TF)，又名。一个单词在文档中出现的次数除以该文档中的总单词数，即它衡量一个术语在文档中出现的频率。
逆文档频率(IDF)，计算为语料库中文档数量的对数除以特定术语出现的文档数量，即，它衡量术语的重要性。
TF-IDF 的公式。
在计算 TF 时，所有项都被认为是同等重要的。然而，众所周知，特定的词，如“是”、“的”和“那个”，可能会出现很多次，但没有什么重要性。
…向量化…
将每个单词拆分后，需要将这些单词(词条)转换为 SciKit 学习算法模型可以使用的向量，因为这些算法只理解数字特征的概念，而不考虑其底层类型(文本、图像和数字等)。)允许我们对不同类型的数据执行复杂的机器学习任务。
Scikit-learn 的 tfidf 矢量器旨在将一组原始文档转换为 TF-IDF 特征矩阵，它可以一次性计算字数、IDF 和 TF-IDF 值。
图片作者。
这样，我们就生成了 tf_idf_matrix，这是一个稀疏矩阵。与 Tfidfvectorizer 一样，我们将原始文本转换为单词和 n 元语法的数字向量表示。正如已经说过的，这使得直接使用变得容易；算法只理解数字特征的概念，而不管其基本类型(文本、图像、数字等)。).
余弦相似性
两个文本文档在它们的上下文(表面接近度)和意义(即，词汇相似度和语义相似度)方面如何彼此接近被称为文本相似度，并且有各种方法来计算文本相似度，例如余弦相似度、欧几里德距离、雅克卡系数和骰子。
余弦相似性度量两个文档之间的文本相似性，而不考虑它们的大小。数学上，余弦相似性度量测量在多维空间中投影的两个 n 维向量之间的角度的余弦，并且取值范围从 0 到 1，
在哪里，
- 1 表示更相似
- 0 表示相似度较低。
数学上:
余弦相似性公式。
然而， ING 批发银行高级分析团队的数据科学家发现余弦相似度有一些缺点:
- sklearn 版本做了大量的类型检查和错误处理。
- sklearn 版本一次性计算并存储所有相似度，而我们只对最相似的感兴趣。
因此两者都会导致使用更多的内存消耗和时间。
为了优化这些缺点，他们创建了他们的库，该库仅存储每行中前 N 个最高匹配，并且仅存储阈值以上的相似性。
获取更多关于库:sparse _ dot _ topn
他们自豪地声称，这种方法提高了大约 40%的速度，并减少了内存消耗。
有了这些基本的基础知识，我们开始我们的代码工作。
我假设你熟悉 Jupyter 笔记本和下载库。如果没有，请查看我之前的帖子，在那里我已经分享了探索性数据分析的初学者指南。
工作环境: 来自 Anaconda 的 Jupyter 笔记本
数据集: 房间类型
数据集快照
数据集的快照。
现在让我们导入所需的库/方法，读取我们的数据集并将其保存到“df”中，只要快速查看数据集，就可以将“df”命名为任何名称。
# Importing libraries and module and some setting for notebookimport pandas as pd
import re
from sklearn.feature_extraction.text import TfidfVectorizer
import numpy as np
from scipy.sparse import csr_matrix
import sparse_dot_topn.sparse_dot_topn as ct # Leading Juice for us
import time
pd.set_option('display.max_colwidth', -1)# reading dataset as dfdf = pd.read_csv('room_type.csv')# printing first five rowsdf.head(5)
快速查看 5 行数据集。
ngrams
这里，我们取 n = 3，因此 3-gram 或三元组，因为大多数房间类型包含两个或三个单词，并且我们将句子分成标记，去除所有特殊字符、标点符号和单个字符(，-。/).收集那 3 克。
def ngrams(string, n=3):string = re.sub(r'[,-./]|\sBD',r'', string)
ngrams = zip(*[string[i:] for i in range(n)])
return [''.join(ngram) for ngram in ngrams] # Testing ngrams work for verificationprint('All 3-grams in "Deluxroom":')
ngrams('Deluxroom')
为了验证，我们测试了“Deluxroom”的 3-gram 结果，输出如下:
“豪华客房”的 3 克结果。
看起来不错！
现在，在拆分每个单词(标记或词条(n-gram 生成的项目) )之后，我们需要对使用 Scikit-Learn 库提供的 Tfidfvectorizer 的每个文档中出现的单词进行计数，并将一组原始材料转换为 TF-IDF 特征矩阵。
room_types = df['RoomTypes']
vectorizer = TfidfVectorizer(min_df=1, analyzer=ngrams)
tf_idf_matrix = vectorizer.fit_transform(room_types)
让我们查看生成的稀疏 CSR 矩阵。
print(tf_idf_matrix[0])
稀疏 CSR 矩阵。
看起来不错！
余弦相似度与稀疏点 topn:领先果汁
我们正在用sparse _ dot _ topn库处理一个 CSR 矩阵。****
**def awesome_cossim_top(A, B, ntop, lower_bound=0):
# force A and B as a CSR matrix.
# If they have already been CSR, there is no overhead
A = A.tocsr()
B = B.tocsr()
M, _ = A.shape
_, N = B.shape
idx_dtype = np.int32
nnz_max = M*ntop
indptr = np.zeros(M+1, dtype=idx_dtype)
indices = np.zeros(nnz_max, dtype=idx_dtype)
data = np.zeros(nnz_max, dtype=A.dtype)ct.sparse_dot_topn(
M, N, np.asarray(A.indptr, dtype=idx_dtype),
np.asarray(A.indices, dtype=idx_dtype),
A.data,
np.asarray(B.indptr, dtype=idx_dtype),
np.asarray(B.indices, dtype=idx_dtype),
B.data,
ntop,
lower_bound,
indptr, indices, data)return csr_matrix((data,indices,indptr),shape=(M,N))**
我们存储前 10 个最相似的项目，并且只存储相似度高于 0.8 的项目，并显示模型所用的时间。
**# Top 10 with similarity above 0.8t1 = time.time()
matches = awesome_cossim_top(tf_idf_matrix, tf_idf_matrix.transpose(), 10, 0.8)
t = time.time()-t1
print("SELFTIMED:", t)**
Selftime 0.0019，快！
然后，我们解包得到的稀疏矩阵；
**# unpacks the resulting sparse matrixdef get_matches_df(sparse_matrix, name_vector, top=100):
non_zeros = sparse_matrix.nonzero()
sparserows = non_zeros[0]
sparsecols = non_zeros[1]
if top:
nr_matches = top
else:
nr_matches = sparsecols.size
left_side = np.empty([nr_matches], dtype=object)
right_side = np.empty([nr_matches], dtype=object)
similairity = np.zeros(nr_matches)
for index in range(0, nr_matches):
left_side[index] = name_vector[sparserows[index]]
right_side[index] = name_vector[sparsecols[index]]
similairity[index] = sparse_matrix.data[index]
return pd.DataFrame({'left_side': left_side,
'right_side': right_side,
'similairity': similairity})**
我们观看比赛。
**# store the matches into new dataframe called matched_df and
# printing 10 samplesmatches_df = get_matches_df(matches, room_types, top=200)
matches_df = matches_df[matches_df['similairity'] < 0.99999] # For removing all exact matches
matches_df.sample(10)**
火柴看起来相当令人满意！
套房、 1 居室和豪华套房，1 居室很可能不是同一个房型，我们得到的相似度为 0.81。看起来不错！****
按照排序的顺序，我们查看我们的匹配。
**# printing the matches in sorted ordermatches_df.sort_values(['similairity'], ascending=False).head(10)**
排序匹配。
比赛看起来相当令人满意！余弦相似性给出了两种房间类型之间相似性的良好指示。大转角大床房 1 特大床和大转角大床房大概是同一个房型，我们得到的相似度为 0.93。
因此，正确的视觉评估和用这种策略进行的匹配是非常令人满意的。
因此，通过使用 ngram 进行标记化，使用 TF-IDF 进行向量矩阵，使用 TfidfVectorizer 对每个文档中出现的单词进行计数，并使用与 sparse_dot_topn 的余弦相似性，我们甚至可以最快地匹配大型数据集的字符串(使用 57200012 获得了良好的结果)。*
在我的 GitHub 中获取整个工作 Jupyter 笔记本。
喜欢这篇文章吗？成为一个中等会员继续学习，没有限制。如果您使用 以下链接 ，我将收取您一部分会员费，无需您支付额外费用。
如果你对这篇文章有任何疑问，或者想在你的下一个数据科学项目中合作，请在 LinkedIn 上 ping 我。
敬请关注下一篇与数据科学相关的帖子。
推荐读物
- 使用 TF-IDF 从噪声中识别信号
- 促进大规模数据集中最相似实体的选择
- 余弦相似度——文本相似度度量
- N-gram 简介:它们是什么，我们为什么需要它们？
- 快速探索性数据分析:熊猫概况
调查美国公司的债务
费利克斯·米特迈尔在 Unsplash 上的照片
我们用数据来评估美国公司的财务健康状况
(如果你想用我的代码， 你可以在我的 GitHub 这里找到。抱歉，还没清理/评论。此外，你需要订阅 Sharadar 美国基础和一个 Quandl API 密钥来运行它)
典型的美国公司有多健康？这是一个非常相关的问题。虽然有理由对新冠肺炎局势的改善保持谨慎乐观，但这只是第一步。接下来的事情(重新开放和恢复经济)将会同样困难或者更加困难。
当我们最终走出这条经济隧道的另一边时，许多大大小小的企业都将不复存在。今天，我们将研究潜在的财务数据，以了解美国企业的债务状况，以及是否存在任何明显的问题点。
为什么企业会违约
公司破产最常见的原因是高昂的固定成本——换句话说，当收入骤降时，无法降低成本。公司有固定成本有多种原因:
- 为了正常运转，企业需要在工厂等高成本资产上进行大量前期投资。如果一家公司刚建完新工厂，销售额就暴跌，这对公司的底线来说是灾难性的——它付出了现金，却没有办法赚回来。职能型的高固定成本在金融术语中被称为经营杠杆。
- 固定成本高的另一个原因是债务。任何一家公司，无论其成本结构如何，都可以通过负债在繁荣时期增加回报。这有点像用经纪人提供的保证金购买股票。如果股票价格持续上涨，很好，你赚了更多的钱，因为债务允许你购买更多的股票。如果股市崩盘，那你就有大麻烦了。你不仅损失了大部分投资，还需要偿还经纪人你所借的原始金额。大量借贷的公司面临同样的风险——如果企业倒闭，现金流枯竭，那么他们可能无法偿还债务。
- 这两者的混合是隐性债务，如被锁定在长期租赁或租赁合同(如租赁飞机或租赁大量商业房地产)。
今天，我们将重点探讨等式中的债务部分，因为巨额债务通常也是高运营杠杆的标志(除非一家公司像苹果一样拥有大量现金流，否则它通常需要借款才能支付前期基础设施投资)。
谁持有所有的债务？
我们先来看看哪些行业负债最多。竞争动态(利润率、基础设施需求等。)不同行业的公司差异很大，因此我们应该预计债务的用途也会有所不同。我在下面按部门绘制了总资产负债表债务(包括长期和短期)。一些观察结果:
- 金融服务部门(银行等。)雇佣了最多的债务。不过，这并不一定意味着它们的违约风险最大。银行将债务作为赚钱的工具(或原始投入)。在高层次上，它们将通过存款(如支票和储蓄账户)借入的资金转化为产生现金的贷款——因此它们基本上是用来在经济中输送资金的。只要这些贷款产生的利息收入(而不是全体违约)高于银行为其存款支付的利息支出，银行就会赚钱。对银行的财务分析应该更多地关注他们的股本缓冲(总资产减去总负债)和贷款的风险。在评估一家银行时，与其说它借入了多少美元，不如说它把借来的钱投资到了哪里。由于这些原因，在接下来的一些情节和分析中，我将省略金融服务公司。
- 借入大量资金的行业通常是那些前期投资需求较高的行业，如通信服务(康卡斯特，美国电话电报公司)、能源(埃克森美孚)、公用事业和房地产。
截至 2019 年底的行业负债情况(来源:Sharadar 美国基本面)
**总债务水平很有趣，但更重要的是看相对于收入的债务水平，这是公司偿债能力的一个因素。**我们可以通过将债务除以 EBITDA 并按行业绘制平均债务/EBITDA(如下)来实现。EBITDA 代表未计利息、税项、折旧和摊销前的收益。EBITDA 是一家公司手头用于偿债的现金流的合理指标。
- 虽然通信服务行业发行了大量债务(在我们之前的情节中，仅次于金融服务)，但它也产生了大量利润，因为它包括谷歌和脸书等利润丰厚的公司(因此，从债务/EBITDA 角度来看，它看起来相当健康，这还不包括他们的巨额现金储备)。康卡斯特和美国电话电报公司也因其在家庭互联网上的寡头垄断而获得了丰厚的利润。
- 在这个情节中，真正跳出来的是房地产。该行业的平均收入并不多，但却借了很多钱(每 1 美元的 EBITDA 就有 17 美元的债务)。恕我直言，我不确定这一领域的商业地产风险是多少(远程办公趋势和向在线零售的持续转移对写字楼和购物中心的所有者来说不是好兆头)。但看看这个领域的大型公司(按市值计算)，有像 American Tower 和 Crown Castle 这样拥有和租赁蜂窝塔的公司，还有 Equinix，它在数据中心方面也做同样的事情——今天已经是必不可少的服务，也许明天在家工作会更是如此。因此，尽管随着经济形势的恶化，相对于收入的高负债将被证明是一个挑战，但这些公司并不像乍看上去那样受到 COVID 的影响。
2019 年底按行业划分的债务/EBITDA(来源:Sharadar 美国基本面)
按行业和市值分类
让我们稍微集中一下我们的镜头，把市值也包括在我们的分析中。我们可以使用热图来直观显示按行业和市值划分的平均债务/EBITDA 水平(市值是公司规模的一个不错的代表)。
按行业和市值划分的债务/EBITDA(来源:Sharadar 美国基本面)
在我们深入研究这些数字之前，请注意，对于一家公司的债务/EBITDA 比率应该是多少，并没有正确的答案。这取决于管理层管理流动性的能力、债务的利率、现金流的周期性(如果你的公司在衰退中现金流骤降，那么它应该注意它借了多少钱)，以及用这些债务融资的项目的回报。显然，这些变量因公司而异。它们还会根据经济状况或特定行业状况而变化——你可以打赌，如果衰退加深，EBITDA 将严重下降，以前看起来可控的债务将成为巨大的负担。
认识到债务不一定是坏事是很重要的。如果一家公司有大量且稳定的现金流，能够以低利率获得债务(由于今天的超低利率环境)，并能够从债务中获得超过所支付利息的回报，那么它至少应该借一点点。只要管理层不高估现金流的规模和稳定性，像谷歌、苹果和伯克希尔哈撒韦这样的公司的超低利率债务几乎等同于免费资金。
好，回到我们的热图。总的来说，我注意到以下几点:
- 相对于收入而言，小型房地产公司和大型工业公司负债累累。我们绝对应该看一看这两个桶里的公司。
- 正如我们之前观察到的，一般来说，房地产公司的债务相对于他们的支付能力来说是很高的。我们还发现，该行业中负债最重的是规模较小的公司(房地产投资组合不太多样化)。
- 中型能源公司(我想知道那里有多少压裂设备)和大型基础材料公司相对于他们的同行来说看起来负债过高。负债累累的大宗商品公司让我感到紧张——它们对所售商品的价格没有任何控制权，因此只要价格一跌，它们就可能被彻底摧毁。
我注意到的另一件事是，除了工业类股，大型公司的平均负债水平高于同行业的大型公司(见下图)。我猜想，这是因为大型公司主导着各自的行业，享有规模竞争优势，因此利润率更高。因此，尽管以绝对美元为基础，相对于收入，他们可能会借更多的钱(由于他们更高的盈利能力)，但大型公司实际上看起来负债更少。不过，我很惊讶地看到，在基础材料等大宗商品领域，债务存在如此明显的差异。
大盘股比超级大盘股负债更多(来源:Sharadar 美国基本面)
负债累累的个别公司
下面是能源、房地产和工业部门负债最多的 30 家公司的图表 (我只包括负债至少 50 亿美元的公司，因为我想把重点放在鲸鱼上) 。
看看前几名，他们主要是抵押贷款房地产投资信托基金(如卡普斯泰德和 AGNC)。与银行类似，这些公司借款是为了投资抵押债券组合(因此债务或多或少是他们业务的基本投入)。因此，高水平的债务是可以预期的，并不一定是一件坏事——只要抵押贷款房地产投资信托基金敏锐地管理其投资组合中的违约风险(在整个商业周期内)，并恰当地匹配其资产和负债的期限(又称债务)。因此，我们应该越过房地产公司，专注于能源和工业(Brookfield Property Partners 是一个例外，它管理着办公楼、酒店、购物中心等。并可能在未来几个月陷入困境)。
抵押贷款房地产投资信托基金之后，你会看到财务管理不善的公司，如波音(太多的股票回购)和通用电气，以及能源公司，如越洋公司和阿帕奇公司。越洋公司是一个有趣的公司。2014 年前，当油价仍徘徊在 100 美元左右时，Transocean 的海上钻井船和钻机船队需求巨大。每艘船(或钻机)成本很高，但赚得更多——所以只要油价居高不下，借钱投资更大的船队是明智的。显然，随着石油价格跌至 20 年来的最低点，以及陆上页岩油产量的激增，Transocean 公司昂贵的船队在过去 5 年里已经成为一个巨大的财务拖累。
30 家负债最多的能源、房地产和工业公司
以下是除金融服务和房地产 (仅指负债至少 50 亿美元的公司) 之外的所有行业。这不是一个漂亮的列表。下面的大多数公司都受到了疫情的严重影响，即使经济重新开始，他们也会继续感受到痛苦:
- 波音和拉美航空集团受到航空旅行急剧下降的沉重打击。
- Transocean、Apache、Cheniere Energy 和 Noble Energy 将发现在当前的油价下几乎不可能赚钱。
- 凯撒娱乐公司(Caesars Entertainment)和 AMC 娱乐公司(AMC Entertainment)将很难分别吸引人们回到赌场和影院(当它们重新开业时，将不得不承担执行严格的社交距离和容量规则的成本)。
- 福特和 Clear Channel Outdoor Holdings(广告牌)将在未来几年受到远程工作加速趋势的沉重打击。更少的司机意味着更少的汽车销售和更少的通勤者观看户外广告。
请记住，企业业绩往往并不等同于股票表现，尤其是在短期内。股票表现完全取决于实际发生的情况和预期发生的情况。正如即使是伟大公司的股票也会在市场低迷时大幅抛售一样，管理不善和缺乏竞争力的公司的股票也能在牛市中表现出色，因为牛市中充满了动物精神。
30 负债最多的行业，除了金融服务和房地产(来源:Sharadar 美国基本面)
结论
这只是洋葱的第一层。我计划在不久的将来深入研究债务结构、商业模式以及公司财务报表的其他领域。
最后，在我们开始之前，让我们检查一下总体债务水平在过去 5 年中是否增加了。他们有！除了基本材料，其他所有行业的债务都有所增加。当然，这是长期经济繁荣(刚刚结束)的副产品——随着公司的销售、利润和资产增长，他们能够借更多的钱，也确实借了。但这种行为也增加了经济低迷时期的违约风险——债务水平创下历史新高，而经济状况的恶化降低了利润和支付能力。
过去 5 年，几乎所有行业的债务水平都有所上升(来源:Sharadar 美国基本面)
希望这是有见地的。当我开始分析这些数据时，我脑子里并没有一个故事，只是想看看我能找到什么。我仍在分析，所以会有更多的来。干杯，保持安全，投资愉快！
涡扇发动机预测维修的生存分析
弗拉季斯拉夫·切尔卡先科在 Unsplash 上的照片
探索美国宇航局的涡轮风扇数据集
解释危险和示例实施
<免责声明:我的目的是展示模型开发过程中不同方法和选择的效果。这些影响经常使用测试集来显示，这被认为是(非常)不好的做法，但有助于教育目的。>
欢迎来到“探索 NASA 的涡轮风扇数据集”系列的另一部分。这将是对第一个数据集(FD001)的第四次也是最后一次分析，其中所有发动机在相同的操作条件下运行，并产生相同的故障。
在我的上一篇文章中，我们深入研究了时间序列分析，并探索了用于预测性维护的分布式滞后模型。最终模型表现相当好，RMSE 为 20.85。今天我们将探讨生存分析。这是一种我渴望尝试的技术，因为我已经多次听说和阅读过这种技术，它可能是一种适用于预测性维护的方法。然而，我从未遇到过满足我好奇心的示例实现。首先，生存分析到底是什么？
生存分析入门
生存分析起源于医学领域，用来回答关于特定人群寿命的问题。如果你知道某人的年龄，可以预测某人的一生，你也可以估计那个人还能活多少时间。例如，这种技术应用于流行病学或疾病治疗研究。然而，它也可以应用于数据由持续时间和基于时间的事件组成的许多其他情况，例如流失预测和预测性维护。下面我快速总结了生存分析中使用的几个关键概念[1，2]:
事件:一个有趣现象的发生，在我们的例子中是发动机故障。
持续时间:持续时间是指开始观察到事件发生的时间或者观察停止的时间
删截:删截发生在观察已经停止但感兴趣的对象还没有他们的‘事件’的时候。
生存函数:生存函数返回在/过去时间 t 的生存概率
危险函数:危险函数返回事件在时间 t 发生的概率，假设事件直到时间 t 还没有发生
对我来说，生存分析的一个吸引人的方面是在还没有事件的模型中包括主体(或者在我们的例子中是机器)的可能性。在更传统的机器学习中，你会从数据集中丢弃“不完整的”或经过审查的主题，这可能会使结果有偏差[3]。解释了一些基础知识后，是时候开始了！
加载数据
我们将读取数据并计算剩余使用寿命(RUL ),就像我们现在习惯的那样。
train.head()的结果
正如在以前的帖子中所讨论的，我们将剪裁任何高于 125 的 RUL 值，因为这将极大地提高模型性能。此外，从先前的探索性数据分析中得出的非信息特征被丢弃。
数据准备
我们稍后将使用的模型需要一个事件列。因此，让我们添加一个细分列，指示引擎是发生故障(1)还是仍在运行(0)。
接下来，我们需要指出每个观察的开始和停止时间。由于数据集在多个时间周期内具有连续的测量值，因此每个观测值只是一个周期。我们可以使用 time_cycles 列来指示观察的结束，我们将添加一个等于time_cycles — 1的 start 列来指示观察的开始。
train.tail()的结果
注意时间周期、RUL、分解和开始列值，以检查我们所做的数据准备是否符合我们的预期，看起来不错！
在列车组中，每台发动机都出现故障，因此没有任何经过审查的观察结果。我们将通过忽略 200 个时间周期后的任何记录，人为地对我们的数据集进行右审查。这使得我们可以在一个更真实的环境中处理数据，混合使用已经发生和还没有发生故障的引擎。
卡普兰迈耶曲线
所有的数据准备工作都完成了，是时候深入了解引擎的生存时间和概率了。我们可以使用卡普兰迈耶曲线来实现这一点，它所需要的是表明持续时间(时间周期)和事件(故障或功能)的最后一次观察
卡普兰迈耶曲线
这种方法已经为我们提供了一种粗略的工具，来估计来自相同群体的发动机在过去的时间 t 内存活的概率。例如，发动机有 100%的概率在第一个 128 个时间周期中存活。在这之后，第一个发动机开始出现故障，但是仍然有 46%的可能性发动机能够存活超过 200 个时间周期。请注意，这种方法仅表示超过某一点后的生存概率，但不能超出给定的数据进行推断。
Cox 比例风险模型
提供更多信息的常见模型是 Cox 比例风险模型。它计算危险比率，例如指示故障风险，例如在特定设置下运行的发动机发生故障的可能性是在不同设置下运行的发动机的 1.85 倍。
python lifelines 包的 CoxPH 实现还附带了漂亮的“predict_expectation”方法，为您提供了一种估计事件发生时间的直接方法。由于这种 predict_expectation 方法，我尽了最大努力将 CoxPH 模型应用于我们的数据集。为了利用发动机随时间的退化，我使用“cluster_col”来表示发动机的 unit_nr，试图让模型考虑每个发动机的多个观察值。不幸的是，结果相当糟糕。下定决心要成功，我联系了《生命线》的作者卡梅伦·戴维森-皮隆。在这里，我了解到“cluster_col”并不是指与时间相关的样本，而是指具有与时间无关的观察值的组。例如，指示不同的处理组，或在不同操作设置下运行的发动机组。
为了利用数据集的时间序列特性，您必须使用时变模型。然而，这个模型的缺点是它没有一个方法来估计事件发生的时间。它预测了部分危险，解释起来不太直观，你会在下面看到。
CoxTimeVaryingFitter
训练模型非常简单，您实例化模型并调用 fit 方法，传递数据集、id_col 以指示唯一的引擎、event_col 以指示是否发生故障以及 start 和 stop 列，以便模型可以解释观察的持续时间。
# returns
# Iteration 8: norm_delta = 0.00000, step_size = 1.00000, ll = -63.86010, newton_decrement = 0.00000, seconds_since_start = 0.1
# Convergence completed after 8 iterations.# <lifelines.CoxTimeVaryingFitter: fitted with 18627 periods, 100 subjects, 54 events>
接下来，我们将检查模型。
# returns
# <lifelines.CoxTimeVaryingFitter: fitted with 18627 periods, 100 subjects, 54 events>
# event col = 'breakdown'
# number of subjects = 100
# number of periods = 18627
# number of events = 54
# log-likelihood = -63.86
# time fit was run = 2020-09-02 15:26:11 UTC
CoxTimeVaryingFitter 摘要
CoxTimeVaryingFitter 系数图
查看模型摘要，我们对对数似然、p 值和指数(coef)感兴趣。对数似然给出了拟合优度的指示，但仅在与包含较少特征的其他类似模型相比较时。更接近 0 的对数似然性被认为是更好的(不要与对数似然比相混淆！).当查看 p 值时，传感器 9 和 15 的值相当大，p > 0.50。然而，移除传感器 9 和 15 得到的对数似然值为-64.20，因此没有提高拟合优度[4，5]。
exp(coef)显示了结垢危害风险。例如，传感器 11 的传感器值增加 1 个单位，击穿的风险增加 167.43 倍[6]。该图本质上显示了特征的系数和置信区间。
预测和评估
训练好模型后，就该开始评估了。故障(或危险)的风险取决于基线危险和部分危险(见下面的公式)。
来源:生命线时变生存回归[7]
因为我们的发动机来自统一的群体(例如，所有发动机都在相同的运行条件下运行)，所以它们的基准危害是相同的。因此，我们只需检查部分或对数部分危险，就可以得到故障风险的指示。部分危险只有在与同一人群的其他部分危险相关时才有意义。部分危险为 2e⁶的发动机发生故障的可能性是部分危险为 1e⁶.的发动机的两倍因为部分危险值相当大，所以显示部分危险的日志更容易。然而，测井局部风险降低了可解释性。一台发动机相对于另一台发动机的对数局部危险每增加 1 个单位，故障概率就增加 2.718(或 e )倍。
为了开始我们的评估，我们只需要还没有坏的引擎，它们的 log_partial_hazard 和计算的 RUL。
预测结果. head(10)
当将 log_partial_hazard 与计算的 RUL 进行比较时，您可以看到它通常很好地告知了即将发生的故障(此处显示了前 10 个)。对于更有故障风险的发动机，返回更高的 log_partial_hazards。然而，它并不总是准确的，例如，发动机 16 的危险比发动机 15 的危险高得多，尽管发动机 15 会更快损坏。
将所有 log_partial_hazards 与计算出的 RUL 相对照，会生成下图，其中有清晰可见的趋势。
log_partial_hazard 与 RUL 的关系图
由于我们处理的是时间序列数据，我们还可以预测 log_partial_hazard 随时间的变化，并观察它的表现。
一段时间内 log_partial_hazard 的图表
注意:这里要做的实际事情是为 log_partial_hazard 设置一个阈值，在此阈值之后应该执行维护。然而，如前所述，这并不能真正让你了解 RUL。您可以开发一个时间序列模型来预测何时达到该阈值，以获得更多的“事件发生时间”预测。
在现实环境中，我建议使用上面建议的两个选项中的一个。然而，由于之前的模型都预测了 RUL，我将尝试将对数部分风险值与计算出的 RUL 关联起来进行比较。
回归对数-RUL 的部分风险
首先，我们将预测截尾训练集中每个观察值的 log_partial_hazard，并检查它的散点图
RUL 与对数部分风险的散点图
你可以清楚地看到我们的 RUL 削波器在图表顶部附近的影响，但是如果没有削波，扩散会更大。散点图中的关系是非线性的，可能是指数关系。有相当多的传播，使得很难根据 log_partial_hazard 精确定位 RUL，但让我们看看我们会如何进展。
我们将使用 scipy 的曲线拟合[8]拟合一个指数模型来从 log_partial_hazard 推断 RUL。推断出 RUL 后，我们将针对训练集和测试集的计算 RUL 对其进行评估，以了解其准确性。
拟合的思路:对数部分风险的指数模型对单位 nr 1 的 RUL。
首先，我们定义评估函数。
接下来，使用 scipy 的曲线拟合来定义和拟合指数模型。
z 是数据，a 和 b 是模型系数。
# returns
# array([8.85954699e+01, 4.35302167e-02])
最后，准备测试集，并评估训练和测试预测。
# returns
# train set RMSE:26.302858422243872, R2:0.5487597218187095
# test set RMSE:27.135244169256758, R2:0.5736091039386049
27.13 的 RMSE 已经比我们的基准模型提高了 15%，基准模型的 RMSE 为 31.95。现在，让我们在完整的数据集上进行训练，看看模型的表现如何。
在整个数据集上重复
# returns
# Iteration 8: norm_delta = 0.00000, step_size = 1.00000, ll = -114.77106, newton_decrement = 0.00000, seconds_since_start = 0.2
# Convergence completed after 8 iterations.# <lifelines.CoxTimeVaryingFitter: fitted with 20631 periods, 100 subjects, 100 events># train set RMSE:26.226364780597272, R2:0.6039289060308352
# test set RMSE:26.580988808221285, R2:0.5908498441209542
最终的 RMSE 是 26.58，比我们的基线提高了 16.8%，但还没有接近支持向量回归(RMSE = 20.54)或时间序列分析(RMSE = 20.85)的解决方案。部分不准确性可以通过在预测的 log_partial_hazard 上拟合另一个模型来解释，这导致错误上的错误(因为没有模型是完美的)。
我坚信，当你放弃我们一直使用的 RUL 范式，为 log_partial_hazard 设置一个阈值时，这种方法将非常适合于定义何时需要维护。
此外，您不会经常遇到在数据集中有这么多故障示例的真实用例。更少的故障使得准确预测 RUL 变得更加困难。在这种情况下，预测崩溃的可能性并让企业决定什么样的崩溃风险是可接受的，可能会产生更好的结果。
总而言之，我认为这项技术非常有趣，多学一点也无妨！
我要特别感谢 lifelines 的作者 Cameron Davidson-Pilon，他花时间为我提供了一些关于如何最好地利用 lifelines 包来处理手头数据集的指导。此外，我要感谢 Wisse Smit 和 Maikel Grobbe 为我的文章提供的意见和评论。你可以在我的 github 页面这里找到完整的笔记本。
我们对 FD001 的分析到此结束。下一次我们将深入研究的第三个数据集(没错，阅读文章找出原因)，其中发动机出现了两个故障中的一个。一如既往，请在下面的评论中留下你的问题和评论。
参考资料:
【1】https://life lines . readthedocs . io/en/latest/Survival % 20 analysis % 20 intro . html
【2】https://en.wikipedia.org/wiki/Survival_analysis
【3】https://life lines . readthedocs . io/en/latest/Survival % 20 analysis % 20 intro . html #审查
【4】https://stats . idre . UCLA
Keras 中基于深度学习的生存分析
应用生存方法估计概率密度函数
叶夫根尼·切尔卡斯基在 Unsplash 上拍摄的照片
生存分析是*统计学的一个特殊分支，用于分析一个或多个事件发生前的预期持续时间。它也被称为*‘事件时间’分析**，因为目标是估计一些感兴趣的事件发生的时间。它在机器学习中的应用没有限制:在寿命预期、客户保持、预测性维护等问题中应用这些技术是很常见的。
乍一看，生存函数估计似乎是一项艰苦的工作，但经过一些经验和基本的数学概念，它不是那么令人望而却步！最好的起点是头脑中有清晰的想法，我们需要知道什么是我们的时间和事件。只有这样，我们才能继续开发一个模型，为感兴趣的实体提供可靠的生存函数预测。通过一些技巧，我们可以利用生存函数的映射属性，并利用它的其他形式的派生(如危险函数)。
我已经在之前的文章中澄清了这些直觉，在那里我也使用梯度推进方法(比如 LGBM)来估计一个不常见场景中的生存函数。在这篇文章中，我试图遵循同样的目标:估计密度函数，但采用神经网络。
数据
我们使用的数据集包含华盛顿特区首都自行车共享系统从 2011 年到 2012 年的两年历史数据，以每日和每小时的格式提供。提供了各种预测因子来生成预报，如时间或天气回归因子。现在我们考虑小时格式，因为我们需要大量数据来拟合我们的神经网络模型。
我们的范围是预测临时用户在某种时间和天气条件下注册了多少辆每小时的自行车租赁。
乍一看，这个任务可以像简单的回归问题一样轻松地执行。我们试图让事情变得更辣，采取一种生存的方法。我们必须估计每小时临时用户数量的累积密度函数(CDF)。在这个场景中，我们的时间变量由到目前为止注册的临时用户的数量表示，而事件日期被定义为小时结束时的用户总数。我们的模型将学习注册少于一定数量的用户的概率。
累积分布函数
为了处理这个问题，维持生存方法、线性或增强方法(如前后中所做的)需要尊重一些基本假设，以便提供生存函数的正确估计。这也包括，例如，对数据进行一点处理，或者对我们的目标应用一个简单的映射。对于神经网络，我们不关心这个，我们只需要选择正确的结构，并给它提供大量的数据。
一般来说，我们在构建自己的生存深度学习模型时，可以做出两种不同的选择:用回归或者分类来处理问题。这些方法能够估计我们期望的生存函数，提供不同的损失最小化过程，包括以两种不同的形式处理我们的目标。
与同一样本相关的目标示例，用于回归和分类
回归神经网络
用神经网络和回归估计生存函数意味着直接估计每小时的 CDF。当临时用户的数量小于真实计数时，CDF 为 0，而在该值之后为 1。为了使我们的目标密度函数具有相同的长度，我们将可达用户的最大数量固定为 400。在我们使用的网络结构之下。它将我们的外部预测值(天气和时间特征)作为输入，试图复制我们的密度函数作为输出。最后一层中的 sigmoid 激活有助于将最终值限制在 0 和 1 之间。计算训练以最小化均方误差损失。
def get_model_reg():
opt = Adam(lr=0.01)
inp = Input(shape=(len(columns)))
x = Dense(512, activation='relu')(inp)
x = Dropout(0.3)(x)
x = Dense(256, activation='relu')(x)
x = Dropout(0.2)(x)
x = Dense(256, activation='relu')(x)
out = Dense(max_count, activation='sigmoid')(x)
model = Model(inp, out)
model.compile(optimizer=opt, loss='mse')
return model
我们直接根据我们的模型估计的整个概率函数来计算我们的模型的优良性。为此，我们使用连续排名概率得分(CRPS)作为性能指标。
CRPS 对我们任务的提法
我们预测每个小时样本的完整生存曲线，最终，我们在我们的测试数据上达到 0.074 CRPS，这优于基于测试生存函数的虚拟估计的简单模型，其中 CDF 来自于训练(0.121 CRPS)。
使用回归进行每小时生存函数预测的示例
分类神经网络
使用分类方法，我们将目标作为一次性实例来处理。我们的目标，长度等于可达用户的最大数量，对应于注册的精确计数等于 1，否则等于 0。我们的模型经过训练，可以预测到目前为止临时用户的准确数量。为此，我们的神经网络试图最小化分类交叉熵，并具有最终的 softmax 激活函数。
def get_model_clas():
opt = Adam(lr=0.0001)
inp = Input(shape=(len(columns)))
x = Dense(512, activation='relu')(inp)
x = Dropout(0.3)(x)
x = Dense(256, activation='relu')(x)
x = Dropout(0.2)(x)
x = Dense(256, activation='relu')(x)
out = Dense(max_count, activation='softmax')(x)
model = Model(inp, out)
model.compile(optimizer=opt, loss='categorical_crossentropy')
return model
在预测期间，为了获得我们想要的生存函数，我们必须简单地对输出进行后处理，对每个预测应用累积和变换。利用我们估计的生存函数，我们可以像以前一样计算我们在测试集上的性能。利用分类方法，我们达到 0.091 CRPS，这比基于测试生存函数的虚拟估计的简单模型更好，其中 CDF 来自训练(0.121 CRPS)。
分类的每小时生存函数预测示例
摘要
在这篇文章中，我们开发了生存神经网络结构，它能够在给定一些外部预测的情况下估计生存函数。我们工作的特殊性包括在处理标准回归问题时改变观点。我们解决了这个问题，提供了对问题的生存解释，并展示了使我们能够通过深度学习来估计生存密度函数的策略。
保持联系: Linkedin
基于 LightGBM 加泊松回归的生存分析
应用生存方法估计概率密度函数
在 Unsplash 上由 Philippe D. 拍摄的照片
生存分析经常被用来解决诸如“T4”这样的问题，一个特定的事件还要多久才会发生？’。事件和时间是这种方法的关键变量。这些事件被认为是每时每刻都可能发生的极端情况:典型的例子是搅动、死亡、故障等等。时间的定义很简单，但是有点抽象。我们认为时间是持续时间；从这个意义上说，很明显，并非所有的“极端事件”都已经发生。由于一些外部环境，我们无法查看所有的生命历史，这意味着将持续时间视为相对的，即开始可能发生在不同的时间。
传统上，生存分析被用来衡量个体的寿命。该分析不仅可以进一步应用于传统的出生和死亡，还可以应用于任何持续时间等……我们如何通过生存建模来处理分类/回归问题令人惊讶。我们需要定义的只是事件和时间。
在这篇文章中，我们开发了一个解决方案来预测在某种时间和天气条件下每天有多少自行车租赁被注册。任务是回归问题的形式，我们必须预测计数，这也可以提醒我们泊松回归问题，我们必须模拟离散变量。考虑到时间是逐渐“到达目前为止”的租赁次数，而事件是一天结束时记录的最终计数，我们的附加价值是将观点转换为生存方法。为了模拟这种现象，我们将使用一个经过调整的 LightGBM 估计器。
从倒退到生存
当时间和事件被定义时(即，我们分别有观察和目标)，我们必须思考如何才能达到我们的最终范围:估计一个生存函数。
所有战略量化时间；它是一个非负的连续随机变量，其分布由概率密度函数 (PDF) 或累积分布函数 (CDF) 表征。
生存函数是时间的函数，它定义了死亡事件在时间 t 尚未发生的概率，或者等价地，给出了事件发生时间值大于 t 的人口比例。数学上，它是 1-CDF。危险函数 (H) 是事件发生的速率。给定一个固定的时间间隔，这是一个风险度量:危险越大，失败的风险越大。数学上，它是-log(S)。
生存函数估计中涉及的所有元素的直观映射。来源:生命线
有各种方法来拟合生存模型。具体来说，我们对所谓的半参数方法感兴趣，在该方法中，我们试图用一种非常优雅和简单的技巧来学习风险函数:我们假设可以将时间细分为合理的小区间，即我们在每个区间中考虑我们的基线风险常数。在这种格式中，分段比例风险模型相当于某个泊松回归模型。生存函数可以从风险函数中导出，反之亦然，这一事实特别有用，因为它允许我们根据自己的兴趣轻松切换目标。
数据
我们使用的数据集与华盛顿特区首都自行车共享系统 2011 年和 2012 年的两年历史数据相关，这些数据以每日和每小时的格式在此公开提供。目前，我们认为每日数据将我们的注意力集中在预测一天结束时临时用户的数量上。各种预测器可用于生成预报，如时间或天气回归。
请记住，我们将时间变量定义为“到达目前为止”的租金，将事件定义为登记的最终计数，因此一个好的起点是绘制时间的 CDF。
在固定的最大可达数量(4000)下，很明显我们的 train-test 分区的 CDF 非常不同。
泊松回归下的生存光 GBM
利用 LGBM 回归器，应用半参数指数法学习风险函数相当容易。这是可能的(如上所述),因为生存问题的负对数似然性与泊松回归的负对数似然性是一对一映射的，这在 LGBM 库中是默认初始化的。
我们需要做的就是为此目的整理我们的数据。特别地，数据集必须首先被扩展:每一行被复制成多行，从 0 到偶然计数登记，这是我们的死亡事件。然后生成两个新列:
- 到目前为止的计数:从 0 到随机计数取值登记；
- 停止:除了到目前为止计数的行 = 随机计数外，其他地方都等于 0
生成的新变量示例
有了扩展的数据集，我们就可以开始训练了。拟合的计算一如既往:我们的目标由分数(一个 0-1 变量)表示，我们的回归变量是所有可用的外部变量加上迄今为止的计数，之前由我们自己生成。所有的魔法都是由我们的假设和泊松损失造成的。在预测阶段，我们必须像以前一样操作我们的测试数据，通过创建迄今为止的计数变量来扩展输入数据(得分未知)。为了使我们的生存函数有一个统一的形状并便于分数评估，一个好的做法是将目前为止的计数与其他预测因子一起扩展到一个预先设定的极限，对我们来说是 4000。还记得 LGBM 预测风险函数，为了获得相应的生存函数，我们必须操作这个简单的后处理转换:1-exp(cumsum(H))。
让我们用图表来检查一些预测:
如您所见，从风险函数中导出的预测生存函数试图复制从一天结束时临时用户的总数中构建的真实 CDF。我们预测总用户数的累积概率分布。换句话说，我们预测的每个值都表明了一天结束时最终计数小于或等于该数字的概率。
对我们结果的评估是根据这类任务的标准分数进行的。连续排序概率得分(CRPS)将 MAE 推广到概率预测的情况。
CRPS 对我们任务的提法
CPRS 是概率预测中使用最广泛的精确度指标之一。它需要评估所涉及的整个概率函数，因此需要估计每个每日样本的整个生存曲线。最后，我们可以在测试上达到 0.283 CRPS，这比基于测试生存函数的哑估计的简单模型更好，其中 CDF 来自于训练(0.291 CRPS)。
摘要
在这篇文章中，我们开发了一个生存 LGBM，它能够在给定一些外部预测和一些简单假设的情况下估计生存函数。我们解决了一个非经典的生存问题，在这个问题中，我们不需要估计在过去的时间中幸存的概率，但是我们估计了在一天结束时特定事件发生的概率。这是一个简单的例子，但它显示了如何应用生存建模技术，在一个不常见的场景中，我们的职责是估计概率密度函数。
保持联系: Linkedin
参考文献
Kaggle: NFL 大数据碗
生命线:生存分析介绍
普林斯顿大学:广义线性模型，第七章，生存模型
生存功能 101
生存分析核心部分之一的介绍。
斯蒂芬·道森在 Unsplash 上拍摄的照片
在继续写这篇文章之前，我强烈建议你了解生存分析，因为这篇文章的主要焦点是生存函数，尽管已经提供了生存分析的简要介绍。
什么是生存分析？
生存分析是一种建模技术，用于模拟事件发生的时间。它用于了解一个事件发生需要多长时间。这些事件可能是心脏病发作后的死亡时间、毕业后找工作的时间、一对夫妇保持婚姻的时间或机器正常工作的时间等。
生存分析包括三个部分:
- 收集和处理生存数据
- 选择生存函数
- 通过不同的方法进行分析
什么是生存功能
这个函数可以帮助我们知道，一个感兴趣的物体，是否会在特定的时间后继续存在。这个感兴趣的对象可以是病人、机器或任何其他东西。
从数学上讲，生存函数用 S 表示，它显然是时间的函数。生存函数也可以解释为某个感兴趣的对象在超过某个时间 t 后生存的概率。
所以可以表示为:
s(T)=P(TT)
这里 T 是随机寿命，取自必须进行生存分析的群体，因此它不为零
该函数的值介于 0 和 1 之间(包括 0 和 1 ),并且是一个非递增函数。
存活函数有时也被称为存活函数或可靠性函数。
生存函数 S(t)也可以通过以下方式与累积分布函数 F(t)相关联:
S(t) = 1 - F(t)
因此它也被称为互补累积分布函数。
如何解读一个生存函数？
让我们看一个任意的例子，看看我们如何解释一个给定的生存函数。
任意生存函数
在观察了上面的生存函数之后，我们可以说，被观察的总人口中只有 25%存活超过 4 个月。
我们可以做出更多的预测，比如有百分之多少的人能活过 6 个月，大约是 17%，等等。
**注意:**这不是解释生存函数的唯一方式，根据分析师的需要，还可以进行其他解释，如中位生存。
生存函数的类型:
生存函数主要分为两类:参数生存函数和非参数生存函数。这种区分基于可以是平均值、标准偏差、中间值等的参数。
让我们来看看一些参数生存函数:
1.指数生存函数
在这种函数中，无论感兴趣对象的年龄是多少，失败的概率在每个时间区间都是相同的，所以对象的年龄不影响概率。
指数生存函数
2.威布尔生存函数；
该函数实际上扩展了指数生存函数，以允许恒定的、递增的或递减的危险率，其中危险率是对物品失效或死亡倾向的度量，取决于其达到的年龄。
威布尔生存函数
3.其他不同的生存功能
我们可以选择其他生存函数，如正态、对数正态、对数逻辑和伽玛，这取决于数据集的分布情况。对数逻辑生存函数如下所示:
对数逻辑生存函数
**注:**第二类生存函数是非参数生存函数，这里没有讨论，它通常使用著名的 Kaplan-Meier 估计量建模，这本身就是一个广泛的话题。要了解卡普兰-迈耶估计量，请访问 Pratik Kumar 的 了解卡普兰-迈耶估计量 。
结论
存活函数构成了存活分析的主要部分，并帮助分析者获得与直到时间事件相关的有价值的知识。如果生存函数得到明智的解释，许多重要的问题都可以得到解答。
这篇文章是作为学生守则的一部分写的——在 2020 年。
学生代码是一个全球性项目，旨在向大学/大学预科学生介绍开源软件开发。
参考
- https://lifelines.readthedocs.io/en/latest/Quickstart.html
- https://www.youtube.com/watch?v=MdmWdIV5k-I
- https://en . Wikipedia . org/wiki/Survival _ function #:~:text = The % 20 Survival % 20 function % 20 is % 20a，survivor % 20 function % 20 or % 20 reliability % 20 function。
- https://www.youtube.com/watch?v=XQfxndJH4UA
- 生存分析:大卫·克莱恩鲍姆和米切尔·克莱因的自学文本。
生存指南:实施大数据分析
unsplash.com
在您的组织中实施大数据分析时，您是否感到有些迷茫？
你并不孤单。
把这看作是你创造分析文化的生存指南——包括通往成功的路线图，指引你正确方向的指南针，点燃你激情的火花。
实施大数据分析的 8 大挑战(以及如何应对挑战)
1。没有定义或传达的成功基准
没有可衡量的成功定义的分析计划比那些记录了 KPI 的计划更容易失败。但是比没有对你的成功进行基准测试更糟糕的是没有将它们传达给组织。
任何旨在改变现状的改变，都会有反对者和怀疑者。那些感觉受到威胁的人会质疑该计划的价值，并试图给它泼冷水，但当你大喊用屋顶上的数据推动商业决策的指数收益时，这很难做到。
成功是否包括使用新的流程来访问以前锁定的数据？如何让更多的终端用户获得数据控制面板？用预测分析影响关键业务目标怎么样？
为了帮助取得成功，在口袋中部署分析计划。
在各种业务单位中进行有管理的推广。让最终用户立即接触到新的解决方案，或者在新的流程在企业推出之前让几个部门参与其中。越多有经验的用户能够帮助新用户，你就能为你的团队和 IT 卸下越多的负担。
2。良好的输入不是主要关注点
虽然业务用户需求应该始终推动您的分析计划，但请注意不要忽略这些需求的潜在需求。业务用户最终关心的是分析的输出，因此很容易过分强调数据的可视化而忽略输入。
构建仪表板对于向业务领导者和利益相关者传达见解至关重要，但它是分析周期的最后一层。
仪表板是分析周期的最后一层。在数据可视化之前，必须对其进行查找、验证、准备、调查和建模。
作为一名变革者，你需要向利益相关者强调好的、干净的数据和合理的分析的重要性。对最终结果有狭隘的看法，而不事先关注和支持每一步，只会助长混乱。
有效的分析解决方案允许用户在分析生命周期的每一点了解他们的数据，而不仅仅是在结束时。
3。直到最后一分钟才被循环
由于利益相关者越来越大的压力和迫在眉睫的最后期限，让看门人不参加项目会议并希望他们会签字同意是很诱人的，甚至是很方便的。但是，如果您的业务和 IT 团队之间没有强大的协作，您如此努力点燃的运动将会在着火之前熄灭。
IT 在企业级解决方案、部署和数据治理方面拥有丰富的经验和良好的记录，知道如何实施分析计划。尽早寻求他们的建议，看着他们成为你最大的盟友。
不要忘记，目标是企业的分析能力。因此，您需要推动满足 IT 和最终用户需求的计划和软件部署。这个过程的标志应该是业务和 IT 协同工作，在美妙的、篝火般的和谐中。
4。数据访问的大门紧闭
为了传播企业分析能力，你必须致力于打开数据仓库。如果您的计划有任何回报，最终用户需要访问数据。拥有数据的终端用户越多，就意味着赢得的机会越多。请记住，这并不意味着缺乏数据治理。
结合现代技术的深思熟虑的流程将使您能够将数据保存在受控的环境中，同时仍然为用户提供他们需要的访问权限。
说你信任你的终端用户并向他们提供数据很容易，但这就是问题所在。准备好相信您的最终用户知道如何处理这些数据。他们可能会犯错误，但使用过时的流程和解决方案(如电子表格)会导致更多的错误(事实上，一项研究显示 88%的电子表格包含错误)。
将您的数据民主化并让好奇的人得到它，只会产生大量有益于业务的问答问题。
5。业务与 IT 权力之争仍在继续
根据内部历史和文化，准备调解业务用户和 IT 之间的权力斗争——他们可能会很快变得激烈。当用户自助服务时，他们通常认为没有必要，但事实上，情况恰恰相反。
对于自助服务分析这样的解决方案，一个普遍的误解是不需要它。其实完全相反。IT 的一个主要作用是为业务提供信息。
有了企业分析文化，支持变得更加重要。即使企业继续“拥有”数据，IT 部门也将负责存放和交付数据。
他们与业务用户合作，以确保他们拥有支持自己所需的资源。企业分析能力需要这种强大的业务和 IT 协作才能蓬勃发展。
6。不安全导致数据混乱
有效的分析文化将为企业带来许多新的见解。找到关于如何在几乎任何领域提高绩效的可操作信息对企业来说都是一个极好的好处，但增加透明度对一些人来说可能意味着痛苦的调整。
不习惯数据驱动文化的同事可能会对项目绩效受到关注感到不舒服。
如果一个企业没有在第一时间发现改进的机会，它就无法改善结果。
将分析和任何相关计划定位为公司的积极发展非常重要。拥抱和庆祝改进的发现，而不是拖着过去的错误。持续改进是关键词。
7。头衔和传统经验优先于数据
民主化数据和分析的承诺可能会威胁到那些拥有旧的做事方式所需的技能和培训的人。这是公司需要接受的严峻现实。
目前标准普尔 500 指数中大约有一半的公司将在 10 年后不复存在，这种趋势在全世界都存在。随着在创新竞赛中保持领先的压力越来越大，很明显，接受数字化转型和数据驱动文化的公司将在生存和发展的竞赛中胜出。
8。复杂的解决方案在实现其优势之前很久就过时了
当您的组织开发新的分析流程和软件时，很容易陷入过去的模式。大数据计划曾经庞大而昂贵。公司花费了数百万在复杂的、完全治理的解决方案上，这些解决方案需要两年时间来实施。
您的分析计划有如此多的活动部分，以至于最终看起来与 2005 年的商业智能(BI)实施惊人地相似，该实施未能获得用户采用，并且在痛苦的两年后未能获得投资回报。
在实施结束时，您的公司陷入了无法扩展、缺乏通用性且无法支持业务真正需要的广泛使用的分析解决方案或流程。
因此，当你在推动新的流程或分析软件时，选择那些有常识支持的、每个人都容易使用的、实际上让你的生活更轻松的软件。
成为数据驱动型意味着知情，而不是让事情变得过于复杂。
带路
你的 CEO 从来不会问:“你写了多少行 SQL 代码？”他们不在乎。相反，高管们会就你的分析和你的建议提出问题。
组织文化的改变需要整个组织的适应，那些接受这种改变的人将会明白，虽然方法可能会改变，但重点是不变的:从数据中寻找答案。更深刻的见解将成为点燃您的分析文化的火花。
分析实施过程中要避免的 8 种危险
1.没有确定的成功基准
2.好的输入不是主要的焦点
3.直到最后一分钟它才被循环
4.数据访问的大门紧闭
5.业务与 IT 权力之争仍在继续
6.不安全助长了数据混淆
7.头衔和传统经验优先于数据
8.复杂的解决方案在实现其优势之前很久就过时了
TATIANA SERVIN，ALTERYX 的内容营销经理
原文发布于此:https://www . alter yx . com/input/survival-guide-implementing-big-data-analytics
在数据科学训练营中生存
来自编码训练营幸存者的提示、技巧和建议
关于
首先，简单介绍一下我。
我目前正在进行熨斗学校沉浸式数据科学项目的最后一周，最后一步是一个全包的顶点项目。
整个课程非常广泛，试图将两年的研究生课程压缩到仅仅五个月的强化课程中。
顶石是为了展示我新获得的数据科学家工具包而设计的，必须从头开始创建。不允许来自 Kaggle 的预制数据集。
在做最后一个项目的时候，我得到了一份全职的工作，网页开发员和数据科学家，这是我从一开始就梦想的工作。
从我最后的角度来看，整个经历在我脑海中相对清晰，我认为这将是一个与他人分享我的经历的好时机。
我将分享我在这段时间里使用的一些概念，这些概念帮助我保持领先，并在课程结束前两周获得了一份工作。
社区
数据科学训练营，或任何训练营，或实际生活中最重要的方面是社区。
成为拥有相似目标和兴趣的社区的一员是无价的。
有一句谚语说，“不要独自编码”，这乍一看似乎微不足道，但实际上是你学习编程时可以遵循的最佳建议。
在进入熨斗学校之前，我尝试过自己的编程教育。我使用了 FreeCodeCamp.org 大学的免费课程，顺便说一句，我强烈推荐这个课程。
我在一年半的时间里完成了那个项目的大部分工作，独自完成了数千小时的全栈 web 开发材料。
如果我当时知道我现在所知道的，团队合作比单干要好得多，我会立即开始在论坛上搜索一群人一起工作。
我在熨斗中学到的第一件事是，成为一个由志同道合的人组成的强大而专注的社区的一部分是多么有价值。
我们经常一起工作，互相分享我们的想法和支持，一起努力实现同一个最终目标。
我们成立了在线学习小组，持续了一天的大部分时间，轮流教对方我们刚刚学到的东西。
我是最强的程序员，所以我会教授与编程相关的话题。诸如网页抓取、http 请求和响应、面向对象编程等主题。
我们有一个在大学教数学的人，他教我们梯度下降背后的微积分和所有不同类型的数据。
每个人都有自己的位置，成为团队的一员感觉很棒。与单独工作相比，这种激励是无法相比的。
在这样一个艰难的过程中，拥有积极的支持是无价的，即使你觉得你可以自己处理一切，我保证你会通过拥有一个团队来激发你的想法，从而获得无价的知识和洞察力。
贝佐斯规则
克里斯蒂安·威迪格在 Unsplash 上拍摄的照片
帮助我成功的第二个想法来自我不久前读到的一篇关于杰夫·贝索斯和他雇佣的第一批软件工程师的文章。
当他创办亚马逊时，他发布广告寻找有才华的开发者。
他不希望任何开发人员，但他希望最好的，他可以找到。
因此，他的广告声称他正在寻找极具天赋和动力的开发人员，但更具体地说，他正在寻找能够以三倍于其他开发人员的速度交付产品的开发人员。
当我第一次读到这本书时，我被深深地打动了。
首先，我印象深刻，因为贝佐斯要求潜在员工做到最好，而且他不怕把标准设得极高。
我们大多数人会寻找比其他开发者“更好”的开发者。
事实上，贝佐斯想要的不仅仅是更好，而是三倍好的开发人员，这改变了我制定目标和期望的方式。
其次，我被激励是因为这是一个软件开发人员可以量化的标准。
我想澄清的是，我确实理解一个糟糕的开发人员可以很快推出一个糟糕的产品，但这不是贝佐斯所说的。
相反，他期望他雇佣的开发人员能够创造出最高质量的产品，此外，他还期望他们能够在三分之一的时间内完成。
自从我读到这篇文章，我就被激励成为最优秀的开发人员之一，能够申请那份工作，并且知道我在行业的顶端。
训练营是实践这一信念的绝佳机会，也是见证它在行动中的有效性的绝佳机会。
如果一个模块预计在一周内完成，项目在下一周结束时到期，我会把在两三天内完成这个模块作为我的目标。
然后，我会用第一周的另一半时间做模块的最终项目的第一次尝试，到下周末，我将完成三个独立的项目。
这意味着当其他人第一次开始做他们的项目时，我已经完成了一次，第二次也快完成了。
我知道你在想什么。质量和数量呢？
我发现在一个项目中，每一次成功的尝试都会成倍地提高质量。不仅如此，花费的时间也减少了。
我在网上看到无数的帖子证实了我的观察，当你学习的时候，数量比质量更重要。
学习 web 开发的一般建议是创建尽可能多的应用程序，因为每一个连续的应用程序都是一个学习的机会，也是一个更熟悉你将使用的流程、库和工具的机会。
这就是我如何能够展示 Flatiron 迄今为止制作的一些最好的学生项目，以及我现在如何拥有比我的同龄人中任何其他训练营毕业生多三倍的经验。
10 倍开发人员
TK Hammonds 在 Unsplash 上拍摄的照片
最后一条类似于贝佐斯规则，这是软件开发界的一个常见神话。
与其说它是一个具体的标签，不如说它是一个想法，如果不从字面上理解，如果你努力做到最好，它可以成为一个很好的目标。
这个想法是，有开发者会用 1/10 的时间，用 1/10 的代码行数，解决你的问题。
就像我说的，把编写代码的速度和效率提高 10 倍作为目标很可能会导致失望，然而，把开发人员始终致力于在最少的时间内交付最高质量的产品作为目标会更加现实，并且会产生一个了不起的开发人员。
围绕这个想法有无数的争论，其中大多数人都照字面理解这个定义，就像我警告过不要这么做一样。
相反，我把它作为一种激励工具，作为一种在我可以轻松休假让其他开发人员“赶上”的时候继续工作的方式。
我看到很多开发人员做这种“赶进度”的事情，他们在几个小时内完成一天的工作，然后决定休息一天。
有了这个 10 倍的目标，我就能在其他人完成后继续发展。
我提醒自己，我不是想成为一名普通的开发人员，而是想成为一名 10 倍的开发人员，10 倍的开发人员不会因为完成了当天的任务就扔掉他们的笔记本电脑，而不是自我意识或来自同行的压力。
我知道这种想法并不适合每个人，如果它适合每个人，那么它将成为标准，然后我将不得不以成为 100 倍开发者为目标，所以在某种程度上，我可以感谢不是每个人都有这么高的目标。
结论
彼得罗·兰帕索在 Unsplash 上的照片
通过编码训练营绝对不容易。我注意到有几个人落在后面，退出了我的队伍。他们中的一些人放弃了，重新开始学习，一些人完全退出了，选择了另一个领域。
对我来说，我从来没有想过这些，因为这些目标和想法远远超出了最低限度，甚至失败都不是画面的一部分。
我承认，曾经有一段时间，我在个人目标上失败了。
我没有在一个模块结束时完成项目三次，而是只完成了两次。
把标准定得和我一样高，失败只意味着一种失望感，那就是我没能更快地掌握概念，尽管在课程中从来没有落后的危险。
我没有用三倍的经验构建一个项目，而是用两次重复构建一个项目。
这很令人失望，但我仍然有一个比团队中其他人完成的还要完美两倍的最终项目。
编写 bootcamp 需要大量的动力、时间、金钱，最重要的是对成功的渴望，就像许多其他事情一样，无论你付出多少努力，你都会从中获得回报。
对我来说，这是一次将我从自学成才的开发人员提升到专业人员的经历，我想向自己和他人证明，我可以在这方面做得很好。
最后，我很感激有机会参加熨斗学校，整个过程是我有过的最好的经历之一。
在项目结束前两周，我被录用了，现在我可以非常高兴地说，我是一名专业的软件开发人员。
谢谢！
感谢阅读。我希望有人发现这个信息鼓舞人心，并能够使用这里的想法，成为最好的开发人员。
请随时给我留下任何评论、问题或反馈，我真的很感谢有机会谈论我热衷的东西。
用随机搜索算法幸存僵尸启示录
2019 年波兰北部的日落。私图。
介绍
想象一个典型的好莱坞僵尸启示录场景…
一群绝望的人，在一个魅力四射的英雄的带领下，在某个建筑里设置路障，试图坚持到电影结束。大量的嗜血生物敲着每一扇门和每一扇窗户想要进去。人们试图用他们能使用的任何东西来击退他们，但是由于他们的弹药越来越少，英雄需要派人去取更多的炮弹。
必须做出决定——权衡失去一名队员的风险和弹药耗尽的风险。听起来不像是一个优化问题吗？
虽然它肯定可以用这种方式来表达，但我们会看到，并不总是可以应用标准梯度下降算法 (SGD)。相反，我们将使用一种简单的替代方法，称为随机搜索。为了说明它的简单性，我们将用纯 python 实现它(参见:代码)。
我们的生存场景
为了便于讨论，让我们假设我们的魅力英雄有一台装有文本编辑器(例如 vim，如果他是一个真正的英雄)和 python3 的笔记本电脑——足以快速运行模拟并找出该做什么。
情况呈现如下:
- 有危险的日子，在这期间他们必须挺过去。
- 第***(N+1)***天一到，军方就会来营救幸存者。
- 在那之前，团队每天消耗 一 箱弹药，完全耗尽等于不可避免的死亡。
- 此外，弹药每天都在失去抗僵尸的能力，使得它在 E 天后无法使用。
- 小组从 X 箱弹药开始。
- 在指定的日子里 D = {d_1，d_2，…，d_K} ，空军会投放额外的装有弹药箱的吊舱来支援正在挣扎的士兵，这些日子提前就知道了。
- 不幸的是，从分离舱中回收箱子是有风险的。在[0，1] 中表示为 p，代表失去一个组成员的概率。
- 此外，每从一个箱子中取出一个箱子，风险就会随着箱子数量的增加而增加*。*
- 最后，没有必要在“零日”之后积累任何弹药，因为军方会迅速将其全部带走。
- 现在的挑战是确保团队永远不会耗尽弹药，同时尽量减少潜在的人类牺牲。
建立模型
情况看起来相当复杂，所以我们的第一步是建立一个模型。我们将使用它来模拟不同的场景，这将有助于我们找到答案。这一步在某种程度上相当于 SGD 中的正向传递。
我们知道系统完全是用 N，E，X 和成对的 (d_k，p_k) 来描述的。因此，将我们的模型实现为一个类并公开方法run来执行我们的正向传递是很方便的。
让我们从构造函数开始。
E = 30 # expiry of anti-zombie serum
X = 10 # initial ammo supply
class Survival:
def __init__(self, total_days, pods, supplies_init=[E]*X):
self.total_days = total_days
self.pods = pods
self.supplies = supplies_init
pods是一个元组列表，代表具有相关回收风险的输入弹药舱(例如[(10, 0.1), (20, 0.2)])。由于供应品耗尽，我们将其表示为“剩余天数”列表。一旦列表变空——就完了！
接下来，考虑到供应减少，我们定义以下受保护的方法:
def _consume(self):
self.supplies.remove(min(self.supplies))
self.supplies = list(map(lambda x: x - 1, self.supplies))
self.supplies = list(filter(lambda x: x != 0, self.supplies))
每天消耗一箱弹药(第 3 行)，然后所有的补给离到期日又近了一天(第 4 行)，最后，过期的从列表中移除(第 5 行)。
另一件事发生在新的补给被收回的时候。
def _retrieve(self, risk, quantity):
new_ammo = [E] * quantity
cost = quantity**2 * risk
return new_ammo, cost
在每一种情况下，当我们产生更多的弹药，我们也增加了分数风险，这与我们选择从吊舱中取回多少箱子成正比。(如果成本> 1，我们可以将 is 视为损失一个以上人员的概率)。
为了便于向前传递，让我们再定义一个函数:
def _get_risk(self, day):
pod = list(filter(lambda x: x[0] == day, self.pods))
if len(pod) > 0:
return pod[0][1]
return None
该函数只是选择与为给定日期安排的 pod 相关联的风险(如果已安排)。
最后，向前传球的本质:
def run(self, retrieval_plan):
total_risk = 0
retrieval_idx = 0
for day in range(1, self.total_days + 1):
risk_involved = self._get_risk(day)
if risk_involved:
new_ammo, partial_risk = self._retrieve(
risk_involved,
retrieval_plan[retrieval_idx])
self.supplies += new_ammo
total_risk += partial_risk
retrieval_idx += 1
if len(self.supplies) == 0:
return 0, total_risk, day
else:
self._consume()
return len(self.supplies), total_risk, day
该过程从一些retrieval_plan(稍后将详细介绍)开始，总风险为零。然后，随着我们一天一天的移动，将会有“机会”获得弹药。如果在某一天有这样的机会，我们会根据回收计划尝试回收一定数量的弹药。每当涉及到风险点时，我们都会将它们添加到总体风险计数中。
然后，不管我们是否有更多的弹药，我们检查我们的供应状况。如果我们到达终点，今晚将是这个团队的最后一晚。该函数返回 0 弹药剩余，总牺牲和纪念日。如果没有，这个团队会在英勇保卫他们的要塞时消耗弹药，这个函数会返回相应的值。
这就是我们向前传球所需要的。现在，让我们看看优化。
解决办法
图一。两个 pod 的成本函数与参数空间(z)的关系示例。较低楼层区域与违反约束相关联。
在涉及标准梯度下降的典型情况下，我们现在将设计反向传播通道，并讨论如何正确实施梯度。然而，这里没有梯度！描述风险累积的函数更像是一个过程。它是不可微的。(也没有僵尸参与……)
正因为如此，也因为我们之前提到的诸多限制，我们需要想出另一种方法来降低风险。这就是随机搜索算法的由来。
基本原则如下:
- 我们定义了一个成本函数，这是一种衡量解决方案有多“差”的方法。(在我们的示例中，它是团队成员流失的总风险)。
- 我们在搜索空间中随机初始化系统( z )。
- 我们从( z )的邻域内对我们的搜索空间中的新位置(z’)进行采样。
- 我们评估成本，如果更低，我们移动到那个点(z’->z)。
我们重复 3 遍。第四。只要需要改进，我们就可以应用硬停，例如，如果成本不再下降，时间用完，等等。
当然，这种方法并不能保证我们找到最好的解决方案，但是 SGD 也不能做到这一点。然而，凭直觉，我们可以认为，如果我们继续质疑“现状”并做出小的改进，我们应该会找到一个可能更好的解决方案。
我们来看一下实现。首先，初始化:
def initialize_retrieval_plan(pods, total_days):
a = list(map(lambda x: x[0], pods))
a.append(total_days)
return [a[x + 1] - a[x] for x in range(len(a) - 1)]
该函数接受 pod 的时间表(列表***【d _ k，p_k}*** )和生存天数，并返回一个可能的计划，该计划描述了团队在每个机会可以取回多少弹药——这是一个整数列表。
也许你已经注意到这里没有涉及到随机数。没错，我们不会随机初始化它。在没有从哪里开始的偏好的情况下，随机初始化工作得很好。由于我们的僵尸危机充满了约束，我们可以从至少一个约束可以很容易满足的位置开始。因此，我们从一个保证生存的计划开始，给人们足够的子弹来完成它。
从此，我们可以优化。
def calculate_plan(total_days, pods,
epochs=300,
supplies=[E]*X,
attempts=12,
radius=4):
...
函数calculate_plan是反向传播步骤的替代方法。我们可以将epochs、attempts和radius视为超参数，其中epochs是搜索步骤的数量，attempts是我们在每个步骤下采样的点的数量，radius表示( z )的邻域的“伸展”。
因为理论上计算需要一些时间，我们要做的第一件事是检查条件是否不合适，或者计划是否不会让他们死。例如，如果第一个豆荚的日期是在团队耗尽补给之后，就没有机会再补给，死亡是不可避免的。在这种情况下，我们可以停止计算，开始祈祷。
plan = initialize_retrieval_plan(pods, total_days)
if is_illegal(plan):
print("This retrieval plan cannot make them survive.")
return
S = Survival(total_days, pods, supplies.copy())
supplies_size, cost, day = S.run(plan)
if day < total_days:
return
is_illegal函数是一个简单的检查，如果计划中没有元素是负的(def is_illegal(plan): return (True in [c < 0 for c in plan]))。
检查了可行性之后，我们终于可以开始为最小风险进行优化了。我们的参数空间的维数等于吊舱的数量(等于plan的长度)。该空间的采样意味着随机选取plan的一个元素，看看如果我们从 pod 中移除一些箱子会发生什么。然而，为了防止违反团队不得耗尽弹药的条件，我们也选择了plan的另一个元素来给这个数量。这创建了一个新的计划，我们使用S.run(plan)验证该计划以获得新的结果。
plan_range = range(len(plan))
cost_history = [9999999] # initializing cost history for early stopping
last_plan = plan.copy()
epoch = 1
while epoch < epochs:
cache = {}
for attempt in range(attempts):
i1 = random.choice(plan_range)
i2 = random.choice(plan_range)
plan = last_plan.copy()
qty = random.choice(range(1, radius + 1))
plan[i1] -= qty
plan[i2] += qty
S = Survival(total_days, pods, supplies.copy())
supplies_size, cost, day = S.run(plan)
然后一系列的if-else语句检查剩余的约束，但是如果新的计划是可行的，我们将它添加到cache中，该计划包含具有相关成本函数的步骤的所有备选计划。然后，如果存在成本函数更低的替代方案，我们选择该方案用于下一次迭代(z’->z)。
if supplies_size < 0: # ran out of ammo
continue
elif day < total_days: # ran out of days
continue
elif is_illegal(plan): # attempt to sell ammo
continue
else: # solution found
key = '-'.join([str(x) for x in plan])
cache[key] = cost
last_plan, cost = get_minimum(cache)
print (f"{epoch}, cost: {cost}, plan: {last_plan}")
epoch += 1
cost_history.append(cost)
if cost_history[epoch - 1] == cost_history[epoch - 2]:
break
在某一点上，对某一步骤将没有进一步的计划。当这种情况发生时，我们依靠最后一个可行的计划作为解决方案。
if len(cache) == 0:
print ("No solution for this case. They must die.")
return
else:
print (f"Solution found, but it'll cost = {cost}.")
return last_plan
get_minumum功能是我们定制的arg min版本:
def get_minimum(cache):
min_risk = min(cache.values())
key = list(cache.keys())[list(cache.values()).index(min_risk)]
return [int(x) for x in key.split('-')], min_risk
潜在的问题和替代方案就是这样！例如，跑步:
if __name__ == '__main__':
pods = [(10, 0.2), (15, 0.1), (35, 0.5), (50, 0.3)]
total_days = 70
plan = calculate_plan(total_days, pods)
print(plan)
我们可以看到在寻找更好的计划和成本价值下降的进展。
图二。N = 100 和 10 个 pod 的随机搜索算法的最小成本(左)和最小时间(右)，对“尝试”和“半径”的每个组合检查 10 次。获得了成本函数与回收箱子数量的线性相关性的结果。
不幸的是，如果我们多次执行这种优化，我们将很快发现最“最优”的计划会发生变化，并且成本值有时会高于其他计划。这是因为我们的算法是“盲的”。如果没有梯度，我们在每次迭代中只给它几次机会(attempts)来寻找更好的计划，而且我们是随机进行的。因此，算法到达一个地方 (z`) 就有一个很大的威胁，它不能再前进了——类似于局部最小值。
为了克服这个问题，我们可以自然地增加attempts的数量以给它更多的机会，或者增加radius以允许更多的选择来找到更好的*。*
*在极端情况下，我们可以选择循环遍历距离为radius的所有 ***(z，z’)***对。这种方法相当于强制梯度，大大增加了我们达到全局最小值的机会。另一方面，这是要付出代价的。仅用radius == 1，组合数 ***(z，z’)**与K ** 2成正比，其中K = len(plan)。计算可能需要更长时间。
结论
防御不死族的囤积的确是一个巨大的挑战。除了僵尸之外，在本文中，我们给出了一个简单的随机搜索算法的裸 python 实现。每当我们可以定义一个代价函数，但是我们不能保证它的可微性时，这个算法可以是 SGD 的一个可行的替代方案。此外，如果问题充满了约束，像这样，这里提出的方法可以应用。最后，我们还讨论了算法本身的计算复杂性，如果确实存在重大僵尸危机，这可能会影响下一个英雄的决策。
还会有更多…
我计划把文章带到下一个层次，并提供简短的视频教程。
如果您想了解关于视频和未来文章的更新，请订阅我的 简讯 。你也可以通过填写表格让我知道你的期望。回头见！
可持续的隐私合规要求严格的数据管理
照片由 xdfolio 通过 pixabay (CC0)拍摄
自从 20 世纪 80 年代数据仓库时代开始以来，数据管理，包括元数据管理、数据治理、主数据管理，一直被提倡。然而，它很难实施或执行。没有它，企业仍然可以通过他们的数据仓库项目生存。然而，项目和组织的筒仓会引入数据重复、效率低下和没有明显的事实来源，这一直是许多组织头疼的问题，尤其是那些拥有多个系统的大量数据的组织。
随着最近数据隐私问题的增加以及几项法规的出台，这给每个组织敲响了警钟，要求他们优先实施严格的数据管理。它不仅适用于现有数据资产，还能确保在获取、收集、处理、使用和存储新数据时的法规遵从性。
两个先锋隐私法规是欧盟的 GDPR ( 通用数据保护法规，于 2018 年 5 月 25 日生效)和 CCPA ( 加州消费者隐私法案，于 2020 年 1 月 1 日生效)。这两项法规都旨在增强消费者的隐私权，以保护企业收集、使用或出售的个人信息。尽管受监管公司的规模、罚款金额以及消费者行使权利的方式有所不同，但这两项法规都规定了消费者的以下权利:
- 知情权:了解企业已经收集或正在收集客户的哪些个人数据。
- 访问权:访问客户个人信息的权利。
- 被遗忘权:要求商家删除消费者的任何个人信息。
- 有权反对出于营销、销售或与服务无关的目的处理个人信息(GDPR)，或选择不出售消费者数据(CCPA)。对于后者，CCPA 要求企业在网站上提供一个按钮或链接，明确允许消费者请求“不要出售我的个人信息”
CCPA 将个人信息定义为“识别、涉及、描述、能够与特定消费者或家庭直接或间接相关联或可以合理链接的信息”GDPR 将个人数据定义为“与已识别或可识别的自然人(数据主体)直接或间接相关的任何信息，特别是通过提及标识符。”因此，CCPA 和 GDPR 之间的一个显著区别是，CCPA 的个人信息包括家庭，而 GDPR 的个人信息只包括个人。
对于一家公司来说，确保遵守一项法规已经是一项挑战，更不用说遵守两项或更多法规了。隐私法规将直接影响企业的大规模数据收集、使用和存储。对于没有工具和治理流程来精确定位包含个人信息的数据所在的位置、谁可以访问这些数据以及这些数据是如何被使用的组织来说，这可能是非常耗费资源和时间的。另一方面，如果一家公司在以下领域实施了有效的数据管理，实现法规遵从性应该是简单的。
- 精心设计的数据架构
设计良好的数据架构是成功数据管理的重要促成因素。数据架构和数据管理是不可分的。好的数据架构有清晰的文档，从业务需求、业务到数据的映射、业务术语表、命名约定和标准开始。它还应该在设计中考虑数据隐私和安全性。数据架构给出了每个数据对象到数据元素级别的清晰度、透明度和标准，包括收集、使用和存储个人信息的位置。这将带来清晰的设计、高效的实施和最佳的结果，因为 CCPA 和 GDPR 要求消费者遵守访问、遗忘和不出售其个人信息的权利。
2。元数据管理和数据治理
首先，用准确的分类和数据谱系对所有数据对象和数据元素进行编目将是实现隐私遵从性的一个极好的能力。一个很好的类比是图书目录。对于一个图书馆或多个图书馆中的大量书籍，图书馆目录在过去几十年中已经得到发展，以系统地跟踪书名、作者、位置、流派、出版日期等。不管一个图书馆有多少本书，目录允许用户搜索任何书，而不需要实际阅读它们。随着后来的数字技术的发展，用户也可以先阅读几页样本，然后再决定是借阅还是购买这本书。公司数据资产的数据目录应以同样的方式建立，以便快速搜索可以找出哪些系统中的哪些数据资产包含个人信息，它们是何时创建的，以及它们是如何使用的。
因为公司中的任何职能或部门都可以拥有数据，所以通常需要一个由这些部门的数据管理员参与的数据治理计划，以确保数据目录拥有正确和完整的信息。数据治理也是一个在整个企业中实施标准和策略的平台，以维持组织的法规遵从性。
3。清除数据保留策略
明确的数据保留政策对于公司的隐私合规性也至关重要。如果数据不再对组织有用，或者数据已经按照 GDPR 或 CCPA 等法规的要求保留了一段时间，则应将其删除。例如，GDPR 要求个人资料的保存时间必须"不超过处理个人资料的目的所必需的时间"(第 5(1)条)。随着当今数据科学和分析的快速发展，以多种方式利用数据，使用数据的用户可能会在不经意间创建许多包含敏感信息的数据对象。保留这些数据的时间超过法律要求或使用所需的时间可能会产生各种后果，例如增加数据泄露的风险以及违反与数据访问和数据隐私相关的合规性流程的风险。
4。具有安全性的集中客户信息
出于安全和隐私信息的考虑，公司的客户信息应该存储在尽可能少的地方。应该使用最佳的安全和隐私措施对其进行仔细管理。主数据管理最初旨在提供一个具有有保证的数据质量的真实来源。同时，它还实现了隐私合规性的最佳设计。有了集中的主数据即服务，组织中就不应该有其他地方来永久存储个人信息。对于主数据系统之外的任何操作，它应该到那里检索客户信息作为黄金记录，而不需要担心存储或处理它。另一方面，主数据管理将提供隐私信息的数据接收流程、处理的基本业务规则以及数据访问的仔细控制的清晰性和强制性。综上所述，主数据是公司保护其客户信息并确保完全隐私合规性的最佳设计。
5。自动检测潜在违规的强大数据控制
公司对隐私法规的遵守需要具有高度的可持续性。构建强大的数据控制流程并自动检测不合规情况是组织的重要组成部分。此流程旨在实施业务规则和程序，以便在任何个人信息被添加或填充到不应该添加或填充的位置时发出警报，或者发出任何其他潜在违规风险的警报。
结论
GDPR 和 CCPA 等隐私法规标志着隐私法规在全球的兴起。我们预计很快会有更多的规则出现，每个规则都有其独特性，同时与其他规则有一些共性。这给企业带来了挑战，尤其是那些在多个洲和州拥有客户的企业，他们需要处理各种各样的合规流程。然而，通过精心设计的数据架构、集中且维护良好的元数据管理、数据治理、清晰的数据保留政策、主数据即服务以及强大的数据控制和自动化审计流程，公司将能够快速提出解决方案，并确保以经济高效的方式遵守所有隐私法规，从而获得最佳结果。
奇异值分解——一种几何观点
对通过奇异值分解编码的视觉信息进行解包
(图片由作者提供)寻求在矩阵作用后保持垂直的垂直轴！
教科书中的奇异值分解
奇异值分解实矩阵的奇异值分解通常在线性代数的第一门课结束时讲授。通俗地说，实矩阵的奇异值分解无非是将矩阵表示为三个更好的矩阵的乘积。
说得更确切些，但仍然平淡无奇，实际结果如下:
设 A 为 n×m 实矩阵。然后是 a=uσvᵀ
- u 是一个 n×n 正交矩阵
- v 是一个 m×m 正交矩阵
- σ是 n×m 矩阵，其中除了前 r 个对角线元素之外的每个元素都是零，其中 r 是 a 的秩。
- 此外，σ的非零对角元素正是以降序排列的 A 的所谓的奇异值。
如果那样…看起来有点平淡无奇，让我们借助一个例子来尝试理解这个结果的威力和效用。
由矩阵编码的信息
这是一个非常实际的 2×2 矩阵
(图片由作者提供)A 是矩阵[[6，2]，[-7，6]]
它实际上代表了一个相对于标准基准 e₁=(1,0 和 e₂=(0,1).的线性图线性地图 t 发送通常的基本元素 e₁⇝(6,−7 和 e₂⇝(2,6).现在，我们不仅知道 t 如何作用于 e₁和 e₂，我们实际上还知道它如何作用于任何(a，b)。例如
(图片由作者提供)求 T(2，3)是什么
让我们想象一下我们的线性变换 T 是如何作用于标准坐标轴的。我们可以看到，e₁，e₂的垂直轴被扭曲了，{T(e₁),T( e₂)}的图像不再垂直。
(图片由作者提供)T 在常用轴上的作用
寻找保持垂直的垂直轴
因为使用垂直轴总是更容易，所以我们希望找到我们定义域的垂直轴(或者沿着这些轴的单位向量),即使在应用 t 之后，这些轴仍然保持垂直。显然，e₂}的{e₁不是一个好的选择。有好的选择吗？如果有，怎么找？
由彼此成对垂直的单位向量组成的向量空间的基被称为向量空间的标准正交基。因此，我们在更广泛的背景下的问题可以表述如下:
给定一个线性映射 T:R ᵐ →R ⁿ(相对于标准基用一个 n×m 矩阵 A 表示)，我们能找到rt39】ᵐ的一个正交基 {v_1，v_2，…，v_m} 使得 {T(v_1)，T(v_2
输入 SVD！
A 的奇异值分解恰恰回答了上面的问题！
设 a 的奇异值分解为 a=uσvᵀ
(图片由作者提供)寻找 T 作用后保持垂直的轴
回忆一下 A 的 SVD 中的 V 是一个 m×m 正交矩阵。因此，v 的列形成了 R ᵐ.的标准正交基{ v1，v2，…，v_m}
类似地，u 是正交 n×n 矩阵，所以 u 的列形成了 R ⁿ的正交基{u_1，u_2，…u_n}
但是 U 和 V 有什么关系呢？
将方程 a=uσvᵀ乘以 v，得到
*av =(*uσvᵀ)v =uσ(vᵀv)= uσ
所以 uσ的列正好是{t(v_1),t(v_2),…t(v_m)}⊂rⁿ
此外，由于σ是 n×m 矩阵，如果有的话，仅沿着对角线元素(σI，I)具有非零元素，所以 Uσ的列只是由σ的相应对角线元素缩放的 U 的前 m 列
(图片由作者提供)uσ长什么样
因此 T(v _ I)=(σI，i)u_i，它们相互垂直。由于 u_i 是单位向量，所以 T(v_i)的长度正好是对角线项σI，I。
摘要
给定一个线性映射 T: R ᵐ → R ⁿ(用一个 n×m 的矩阵 a 相对于标准基来表示)，设 a=uσvᵀ是 a 的奇异值分解.那么 v 的列构成了 R ᵐ.的一个正交基此外，在ⁿ.中{T(v_i)}仍然保持相互垂直
向量{T(vi)}的长度由σ的对角线项编码。沿着{T(vi)}的单位向量可以完备化为ⁿ的标准正交基，这由 u 的列精确地给出
所以，你看矩阵的奇异值分解给了我们很多非常有用的信息！
奇异值分解的唯一性
由于σ的非零对角线元素恰好是 A 的降序的奇异值，所以σ对于给定的矩阵 A 是唯一的。
然而，V(以及 U)的列不是唯一的。首先，V 的列可以变化到符号(由于 T(v_i)只是 u_i 乘以一个常数，u_i 的符号也会改变)。为了理解为什么 U，V 不一定是唯一的，也不一定只有符号，考虑一下旋转矩阵 a 的奇异值分解是什么。
奇异值分解计算代码
当然，人们可以(也应该知道如何)手动计算 SVD。但那是以后的事了。现在，我们将使用 numpy 中的 SVD 实现。
import numpy **as** np
A **=** np.array([[6,2], [**-**7,6]])
U,s,VT **=** np.linalg.svd(A)
**print**(U)
**print**(s)
**print**(VT)
上述代码给出了
>>>U
array([[-0.4472136 , 0.89442719],
[ 0.89442719, 0.4472136 ]])
>>> s
array([10., 5.])
>>> VT
array([[-0.89442719, 0.4472136 ],
[ 0.4472136 , 0.89442719]])
总结这个例子
(图片由作者提供)使用 numpy 代码的 SVD
你可以检查 100，25 是 AᵀA 的特征值，因此 10，5 是 a 的奇异值
如果你真的用手算出奇异值分解，你会发现这相当于
(图片由作者提供)手写的奇异值分解
脚注
- 方阵 x 是正交的当且仅当 XᵀX = XXᵀ是单位矩阵
- n×m 矩阵 a 的奇异值是 AᵀA 的特征值的平方根，列出了它们的代数重数。
- 根据定义，零矢量 0 垂直于每个矢量
机器学习中的奇异值分解:PCA
理解什么是主成分分析，我们如何通过 EVD 和奇异值分解实现主成分分析，以及为什么奇异值分解实现更好。
本文讨论什么是主成分分析 (PCA)，我们如何使用特征值分解 (EVD)或奇异值分解 (SVD)，以及为什么 SVD 实现更好。
主成分分析
直观地说，PCA 是一种转换过程，它将具有可能相关特征的数据矩阵转换成一组线性不相关的变量，称为主成分。每一个主成分都是原始特征的线性组合( PCᵢ = Xϕᵢ ，这里 PCᵢ 代表第 I 个主成分)，并且在与前面的成分正交时占最大可能方差(如I≦j，pcᵢ⋅pcⱼ=0)。
给定形状为 n × p 的特征矩阵 X ，即居中，即列均值已被减去，现在等于零*，*PCA 的典型用例包括:
- **降维:**找到形状 n × k (其中 k 远小于 p )的一个更低维的近似值，同时保持大部分方差，作为分类或可视化的预处理步骤。
- **特征工程:**用 p 线性不相关的特征创建一个 X 的新表示。
- **无监督学习:**提取 k 主成分(其中 k 往往比 p 小很多)。通过查看原始要素如何影响这些因素来理解数据集。
从概念上讲，重要的是要记住，主成分分析是一种多变量数据分析的方法，EVD 和奇异值分解都是数值方法。
通过特征值分解的 PCA
传统上，PCA 是基于样本协方差矩阵 C 上的 EVD。假设 X 居中:
C 形状为 p× p 。它是对称的，因此总是可对角化的。我们可以应用特征值分解:
Q 是一个正交矩阵，它的列是 C 的特征向量(q=[q₁， q₂ ，…， qₚ )。 ∧ 是对角线上特征值按降序排列的对角矩阵(λ₁≥λ₂≥…λₚ≥0):
这里 qᵢ 是 Q 的第 I 列，也称为第 I 个主方向**。 λᵢ 是∧的第 I 个对角元素。主成分是原始特征矩阵在主方向上的投影，可以通过 XQ 获得。第 I 个主成分解释的总方差比例为 λᵢ / (λ₁ + λ₂ +…+λₚ)。**
下面是通过 Python 中的 EVD 实现的 PCA:
通过奇异值分解的 PCA
对于矩阵 X，t 这里总是存在矩阵 U ， ， ， V 这样:
U 和 V 都是正交的，σ是对角的。σ的对角元素称为奇异值**(σ₁≥σ₂≥…≥σₚ≥0)。考虑协方差矩阵 C 😗*
对比上面我们知道 V 的列是主方向，第 I 个特征值 λᵢ = σᵢ /(n-1) 。主成分可以通过XV = uσ得到。
这比 EVD 实现好多少？
- ****计算效率:对于高维数据( p > > n )，用协方差矩阵 C 执行计算可能效率低下。
- ****数值精度:形成协方差矩阵 C 会导致精度损失。
- ****数值稳定性:大多数 SVD 实现采用分而治之的方法，而 EVD 的实现使用不太稳定的 QR 算法。
下面是在 Python 中通过 SVD 实现的 PCA:
无监督学习的主成分分析
我们将使用鸢尾花数据集来说明 PCA 如何作为无监督学习工具来帮助理解数据。
我们将从 scikit-learn 加载虹膜数据集。特征矩阵包含跨越 4 个属性的 150 个观察值。每行包含鸢尾花的萼片和花瓣的长度和宽度(以厘米为单位)。目标是花的类型(3 个中的 1 个)，但是我们将只使用它来进行可视化。
加载数据集，执行数据预处理，并应用 SVD。准备 PCA 的一个基本问题是是否在中心化之上标准化数据集**(即确保每列的标准偏差为 1)。从我的角度来看，如果特征是相同规模的，标准化是没有必要的。因此，我们将下面的特征矩阵居中。**
计算并可视化由四个主成分和因子载荷解释的总方差的比例。因子负载是原始特征在主成分中的权重，通过用相应特征值的平方根缩放特征向量来获得。它们帮助我们将每个主成分解释为原始特征的加权和。
我们可以看到，第一主成分解释了总方差的 90%以上，它严重依赖于萼片和花瓣的长度。这意味着我们数据中的大多数变化可以用这两个特征的线性组合来解释。
使用 PCA 作为降维或特征工程工具通常会损害结果的可解释性。当特征的数量非常多时，每个主成分通常是太多不同的原始特征的线性组合，因此难以定义。
要查看标准化的效果，请转换特征矩阵，使每一列都有一个标准差，然后重复可视化。
我们可以看到结果变化很大。现在，第一个主成分只能解释总方差的近 75%，而第二个主成分可以解释近 95%。萼片宽度变得更加重要，尤其是对 PC 2。因此，对于我们是否应该对协方差或相关矩阵执行 PCA 的问题，没有普遍的答案。“最佳”方法是基于主观选择、仔细思考和一些经验。
这篇文章是何坤宇写的。Kunyu 目前是芝加哥大学计算分析专业的硕士生。他发现理解统计建模和机器学习技术、将它们应用于真实世界的数据并帮助创建金融服务行业的端到端解决方案是一件有趣的事情。在 LinkedIn 上联系昆玉！🐷
** [## 昆玉何-即将上任的全球量化策略非周期分析师-美银美林…
芝加哥大学正在接受理学硕士计算分析项目培训的数据科学家。对…充满热情
www.linkedin.com](https://www.linkedin.com/in/kunyuhe/)**
SVM:一个最优化问题
用拉格朗日画线
1.概观
这个博客将探索支持向量机的机制。首先，让我们对这篇文章进行一个 100 英里每小时的概述(强烈建议你在阅读这篇文章之前浏览一下)。基本上，我们给定了 n 维空间中的一些点，其中每个点都有一个二进制标签，并希望用一个超平面来分隔它们。用 x 表示这个空间中的任何一点，那么，任何超平面都可以表示为:w^T x +b=0。让我们列出一些术语。
- x^i:上面提到的 d 维空间中的第 I 个点。所以，它是一个长度为 d 的向量，它的所有元素都是实数(x ∈ R^d).
- t^i:这第 I 点的二进制标号。它可以是 1 或 0。
- w:对于将空间分成两个区域的超平面，x 向量的系数。
- b:对于把空间分成两个区域的超平面，常数项。
在之前的博客中，我们导出了优化问题，如果解决了这个问题，我们就可以得到描述分离平面的 w 和 b(我们将从那里继续我们的方程编号，γ是一个虚拟变量),它最大化了“余量”或离平面最近点的距离。请注意，每个数据点有一个不等式约束。
等式(4): SVM 最优化问题
然后，我们做了一些忍术来消除γ，并简化为以下优化问题:
等式(7):简化的 SVM 优化问题
在这篇博客中，让我们来看看拉格朗日乘子的方法对于解决像这样的约束优化问题可以提供什么样的支持向量机的见解。
2.使用拉格朗日乘数的解决方案
在本系列的前一篇博客中，我们获得了两个约束优化问题(上面的等式(4)和(7 )),它们可用于获得最大化余量的平面。有一种求解带约束优化问题的通用方法(拉格朗日乘子法)。为了保持重点，我们将在这里陈述配方，并使用它来挖掘关于 SVM 问题的见解。至于这个方法为什么有效，请阅读这篇博客其中详细介绍了拉格朗日乘数。如果我们有一个一般性的优化问题，
等式(8)一般约束优化问题。
步骤 1:定义拉格朗日函数:
方程式(9)拉格朗日函数
其中α_i 和β_i 是称为“拉格朗日乘数”的附加变量。不等式α_i 对应的乘数必须≥0，而等式β_i 对应的乘数可以是任意实数。这里再次提供了一些直观的原因。那么，为了使 w 成为最佳值，必须满足的条件(称为 KKT 条件)是:
等式(10):找到约束最优化问题的解的 KKT 条件。
方程 10-e 称为互补条件，它保证如果一个不等式约束不是“紧的”(g_i(w)>0 且不等于 0)，那么对应于该约束的拉格朗日乘数必须等于零。事实证明，这是因为拉格朗日乘数可以被解释为指定不等式约束在最优点上被“推动”的力度。如果约束甚至不是严格的(有效的)，我们在解中根本没有推动它，因此，相应的拉格朗日乘数，α_i=0。等式 10b 和 10c 非常简单，因为它们简单地陈述了原始优化问题的约束应该在最优点处得到满足(几乎是一个重言式)。
2.1.关于 SVM 的观察
让我们回到支持向量机。在等式(4)和(7)中，我们根据每个点到分隔线的垂直距离(边距)为每个点指定了一个不等式约束。与线(边距)距离最小的点将是约束条件相等的点。如果有多个共享该最小距离的点，则它们的约束将根据等式(4)或(7)变成等式。这样的点被称为“支持向量”，因为它们“支持”它们之间的线(正如我们将看到的)。从问题的几何结构中，很容易看出必须有至少两个支持向量(与直线共享最小距离的点，因此具有“紧”约束)，一个具有正标签，一个具有负标签。假设不是这种情况，并且只有一个具有最小距离 d 的点。在不失一般性的情况下，我们可以假设该点具有正标签。然后，在距离为 d+δd 的线的另一侧有另一个带有负标签的点。可以沿着 w 向量向负点移动δd/2 的距离，并增加相同距离的最小边距(现在，最近的正点和最近的负点都成为支持向量)。这意味着我们开始的另一行是一个假先知；不可能是真正的最佳利润线，因为我们很容易提高利润。
但是，这完全依赖于问题的几何解释。如果我们只得到最优化问题(4)或(7)(我们假设我们知道如何从一个问题得到另一个问题)，我们会得到同样的结论吗？对于方程(4)中的问题，方程(9)中定义的拉格朗日量变成:
对γ求导，我们得到，
等式(11)
对 b 求导，
如果我们认为{I}是正标签的集合而{J}是负标签的集合，我们可以重写上面的等式:
等式(12)
等式(11)和(12)以及所有α都≥0 的事实意味着在每个正类和负类中至少有一个非零α_i。并且由于α_i 表示对应于第 I 个点的约束有多“紧”(0 表示一点也不紧)，这意味着两个类中的每一个都必须有至少两个点，约束是活动的，因此拥有最小余量(跨越点)。
2.2“支持向量”到底是什么意思？
在上一节中，我们用公式表示了方程(4)中给出的系统的拉格朗日函数，并对γ求导。现在，让我们为等式(10)给出的公式形成拉格朗日量，因为这要简单得多:
等式(13):等式(7)的拉格朗日量
根据 10-a 对 w 求导，并设为零，我们得到:
等式(14):最优 w。
与之前一样，每个点都有一个与之对应的不等式约束，拉格朗日乘数α_i 也是如此。此外，除了与分隔线具有最小可能距离的点(等式(4)或(7)中的约束是有效的)，所有其他点的α_i 都等于零(因为约束是无效的)。
从方程(14)中，我们看到这些点(α_ I = 0)对拉格朗日量没有贡献，因此对最佳线的 w 也没有贡献。同样容易看出，它们也不影响最优线的 b。因此，在定义直线时，只有最接近直线的点(因此它们的不等式约束变成了等式)才是重要的。这就是为什么这些点被称为“支持向量”。通常只有少数几个，但它们支持它们之间的分离面。
3.最简单分类问题
在这一节中，我们将考虑一个非常简单的分类问题，它能够抓住优化行为的本质。所以我们可以把正在发生的事情可视化，我们把特征空间做成二维的。让我们在上面放两个点，并像这样标记它们(绿色表示正标签，红色表示负标签):
很明显，分隔这两点的最佳位置是由 x+y=0 给出的紫色线。
现在让我们看看我们到目前为止学习的数学是如何告诉我们关于这个问题我们已经知道了什么。为什么要这么做？以便明天它能告诉我们一些我们不知道的事情。
由于我们有 t⁰=1 和 t =-1，我们从方程(12)得到，α_0 = α_1 = α。将此代入方程(14)(这是一个向量方程)，我们得到 w_0=w_1=2 α。因此，我们立即得到，对于 x 和 y，该线必须具有相等的系数。我们还知道，因为只有两个点，所以对应于它们的拉格朗日乘数必须> 0，并且对应于它们的不等式约束必须变得“紧密”(等式；见第 2.1 节最后一行)。这意味着，
减去这两个方程，我们得到 b=0。所以，分离平面，在这种情况下，是直线:x+y=0，正如所料。
3.1.添加第三个浮点
为了使问题更有趣，并涵盖一系列可能的 SVM 行为，让我们添加第三个浮点。因为(1，1)和(-1，-1)位于 y-x=0 的直线上，所以让第三点也位于这条直线上。让我们把它放在 x=y=u，同样，让我们给这个点一个正的标签(就像绿色的(1，1)点一样)。
只要 u>1，平面通过原点。但是当 u <1, the (u,u) point becomes the support vector and drives the separating purple plane downward.
Now, the intuition about support vectors tells us:
- If u> 1 时，最佳 SVM 线不会改变，因为支持向量仍然是(1，1)和(-1，-1)。
- 如果 u∈ (-1，1)，SVM 线随着 u 移动，因为支持向量现在从点(1，1)切换到(u，u)。
- 如果 u
Let’s see how the Lagrange multipliers can help us reach this same conclusion. From equation (14),
Eq (15)
Also, taking derivative of equation (13) with respect to b and setting to zero we get:
And for our problem, this translates to: α_0-α_1+α_2=0 (because the first and third points — (1,1) and (u,u) belong to the positive class and the second point — (-1,-1) belongs to the negative class).
Eq (16)
Next, equations 10-b imply simply that the inequalities should be satisfied. For our problem, we get three inequalities (one per data point). All specified by (per equation (7)):
for our three data points we get:
But from equation (15) we know that w_0=w_1. Further, the second point is the only one in the negative class. So, the inequality corresponding to it must be an equality. Using this and introducing new slack variables, k_0 and k_2 to convert the above inequalities into equalities (the squares ensure the three inequalities above are still ≥0):
Eq (17)
And finally, we have the complementarity conditions:
Eq (18)
From equation (17) we get: b=2w-1. Further, since we require α_0> 0，α_2>0，我们用α_0，α_2 来代替。从等式(15)和(16)我们得到:
等式(19)
将 b=2w-1 代入等式(17)的第一个，
等式(20)
并进入等式(17)的第三个，
等式(21)
现在，方程(18)到(21)很难手工求解。幸运的是，有一个用于求解多项式方程组的通用框架，称为“Buchberger 算法”,上面描述的方程基本上是一个多项式方程组。我在这里写了一篇关于 Buchberger 求解多项式方程组的算法的详细博客。它可以用来根据我们感兴趣的变量简化方程组(简化的形式称为“Groebner 基”)。该算法在 python 库 sympy 中实现。让我们看看它是如何工作的。
等式(18)至(21)及其 Groebner 基，其中 k2 是第一个要消除的变量，k0 是第二个。由于只有五个方程，但有六个变量，所以 Groebner 基的最后一个方程将只用 k0 和 k2 来表示。
上面代码中变量的顺序很重要，因为它告诉 sympy 它们的“重要性”。它试图让 Groebner 基末尾的方程，用末尾的变量来表示。在这种情况下，我们有六个变量，但只有五个方程。并且由于 k_0 和 k_2 是最后两个变量，所以基的最后一个方程将仅用它们来表示(如果有六个方程，则最后一个方程将仅用 k2 来表示)。我们得到:
这意味着 k0k 2 = 0，因此，它们中至少有一个是零。然而，我们知道它们不可能都为零(一般来说)，因为这意味着对应于(1，1)和(u，u)的约束都是紧的；这意味着它们都在离直线的最小距离处，这只有在 u=1 时才有可能。一个或另一个必须为 0 的事实是有意义的，因为(1，1)或(u，u)中恰好有一个必须最接近分隔线，使得相应的 k =0。
做一个类似的练习，但是用 u 和 k_0 表示的最后一个等式，我们得到:
这意味着要么 0 = 0，要么
类似地，根据 k_2 和 u 提取方程，我们得到:
这又意味着要么 k_2=0，要么
这意味着，如果 u>1，那么我们必须有 0 = 0，因为另一种可能性会使它成为虚数。因此，k_2 不可能是 0，而会变成(u-1)^.5.换句话说，对应于(1，1)的方程将变成等式，对应于(u，u)的方程将是“lose”(严格不等式)。这是有意义的，因为如果 u>1，(1，1)将是更接近超平面的点。
类似地，如果 u<1，k_2 将被迫变为 0，因此，k_0 将被迫取正值。我们看到了两点:(u，u)和(1，1)在 u 从小于 1 转变到大于 1 时切换作为支持向量的角色。另一方面，如果 u<0，则不可能找到均为非零实数的 k_0 和 k_2，因此方程没有实数解。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
