TowardsDataScience 2023 博客中文翻译（二百五十二）

原文：TowardsDataScience

协议：CC BY-NC-SA 4.0

粒子群优化：搜索过程可视化

原文：towardsdatascience.com/particle-swarm-optimization-search-procedure-visualized-4b0364fb3e5a

直觉 + 数学 + 代码，为从业者提供

詹姆斯·柯博士

·发表于数据科学前沿 ·阅读时间 9 分钟·2023 年 12 月 1 日

–

图像由 DALL·E 3 根据提示“绘制一幅科幻主题的图像，描绘一群无人机在城市景观中寻找目标”创建

人类喜欢模仿自然界中的许多事物。

我们模仿青蛙游泳。我们通过在飞机上安装翅膀来模仿鸟类以提供升力。我们模仿起重机/蛇/螳螂来进行武术。我们模仿白蚁来建造具有高效温控的结构（参见东门大厦）。

这也扩展到了数学算法中，例如你可能听说过蜜蜂算法、蚂蚁算法、布谷鸟搜索算法和萤火虫算法。我之前还谈过进化算法，它遵循自然选择的原则。

今天，我将讲解 PSO——粒子群优化。在本文的最后，你将获得能够实现解决方案的代码，并生成一个 gif 来可视化搜索过程。

使用案例

在高维空间中寻找最优解是困难的。刚接触机器学习的学生可能在第一周就听说过“维度诅咒”这个术语。

高维空间不仅仅是一个抽象的数学概念。考虑一个供应链问题。一家公司必须决定在哪里建立生产工厂、仓库、分销中心和零售店。为了简单起见，我们假设每种设施只有一个。这已经使我们要寻找的解决方案是 8 维的——(x₁, y₁, x₂, y₂, x₃, y₃, x₄, y₄)。

由 DALL·E 3 根据提示“绘制一个包含四个关键设施——生产工厂、仓库、配送中心和零售店——的全景视图”创建的图像。

绝对位置(x₁, y₁, x₂, y₂, x₃, y₃, x₄, y₄)影响每个设施之间的相对位置。绝对位置和相对位置都会影响运营成本以及预期收益，从而影响利润。虽然不完全准确，但我们可以大致说搜索空间和对目标函数的影响是连续的。

问题的公式化

实际上，将输入（候选解）映射到输出（目标）的底层函数是一个黑箱，无法用数学方法表示。

如果可以的话，可以直接获得解析解。然而，在处理黑箱时，我们将进行采样。一个天真的方法是执行网格搜索。在本文末尾，你将掌握更好的工具，更重要的是，理解为什么它有效。

让我们将供应链问题和黑箱转化为以下数学方程式（以便于操作）。同时，我们将导入所有需要的库。

import math
import numpy as np
import random
import matplotlib.pyplot as plt
from tqdm import tqdm

def blackbox(x1, y1, x2, y2, x3, y3, x4, y4):
    return math.sqrt((x1+2)**2 + (y1-3)**2) * \
            math.sin(
                math.atan2(2*(y2-4), 3*(x2+1)) * (
                    math.cos(x3-x1) + math.sin(y3-y1)
                    + math.cos(3*(y4+3))
                    + math.sin(2*(x4-2))
                ) 
            )

注意一切是如何被故意纠缠在一起的。这是因为如果使用一个简单的组件加法函数（例如，下方所示），可以执行 O(n²)的网格搜索而不是 O(n⁸)。这也与现实不一致，因为我们不能天真地解决生产设施和零售店位置独立的最佳仓库位置问题。

def blackbox_addition(x1, y1, x2, y2, x3, y3, x4, y4):
    return np.sin(x1) * np.cos(y1) \
            + np.cos(x2) * np.sin(y2) \
            + np.sin(x3) * np.sin(y3) \
            + np.cos(x4) * np.cos(y4)

问题的目标是找到(x₁, y₁, x₂, y₂, x₃, y₃, x₄, y₄)，使得blackbox给出的输出尽可能高，而不知道底层方程。为了使问题更具现实性，我们在采样过程中加入噪声。当然，我们可以简单地在每个点重复测量并取平均值，但一个鲁棒的搜索算法即使没有这样做也应能表现良好。

class Task:
    def __init__(self):
        pass

    def score(self, x_arr, with_noise=True):
        if with_noise:
            return blackbox(*x_arr) + 0.1*np.random.randn()
        else:
            return blackbox(*x_arr)

解决方案的公式化

让我们从 PSO 的基本构建块——粒子开始。一个粒子是一个候选解向量，携带其个人最佳解和一组向量，这些向量决定了个人影响和社会影响的程度。在每次迭代中，每个粒子i根据其速度探索新的解：

所有图像均为作者提供，除非另有说明。

该速度受其惯性控制，同时也受到个人最佳解和群体全局最佳解的倾斜影响。[1]

每个粒子的速度由三个组件决定。

w是粒子继续沿其先前轨迹的倾向，而ϕ₁和ϕ₂可以看作是对个人和全球最佳的学习率。拥有惯性是个好主意，因为如果我们之前已经朝着一个假设的良好方向前进，那么继续保持是合理的。事实上，这种动量概念也存在于 Adam 优化器中，用于在梯度下降过程中更新神经网络的权重。

个人最佳，甚至全球最佳，是基于迄今为止发现的内容。可能有些方面（即维度）较好，至少相对而言，但不是全部。然而，我们不知道哪些维度应该保留，因此只能尝试各种组合。秉持进化算法的精神，重要的是最佳粒子而不是平均种群，因此探索时只有上升潜力。

将这个概念转化为数学方程，我们有向量U₁和U₂从均匀分布中随机抽取，乘以各自的学习率ϕ₁和ϕ₂。每个粒子对其个人最佳和全球最佳有不同维度上的亲和力。综合来看，我们可以为Particle类写出如下代码。

class Particle:
    def __init__(self, idx, n_dim, param_limits=(-1,1)):
        self.idx = idx
        self.low = param_limits[0]
        self.high = param_limits[1]
        self.pos = np.random.uniform(
            low=self.low, high=self.high, size=(n_dim,)
        )
        self.w = np.random.uniform(0.1, 0.5)
        self.r1 = np.random.uniform(0, 1, size=(n_dim,))
        self.r2 = np.random.uniform(0, 1, size=(n_dim,))
        self.vel = np.random.normal(size=n_dim)
        self.pbest = np.zeros(n_dim,)
        self.pbest_fitness = -1e8

    def update(self, gbest, task):
        self.vel = self.w * self.vel \
                + self.r1 * (self.pbest - self.pos) \
                + self.r2 * (gbest - self.pos)

        self.vel = np.clip(self.vel, 0.2*self.low, 0.2*self.high)
        self.pos = np.clip(self.pos + self.vel, self.low, self.high)

        fitness = task.score(self.pos)
        if fitness > self.pbest_fitness:
            self.pbest_fitness = fitness
            self.pbest = deepcopy(self.pos)
        return self.pos

注意，我将ϕ₁U₁和ϕ₂U₂分别合并为r₁和r₂，以便（非常轻微地）节省计算。实际上，我将学习率规定为 1 的数量级，以节省这个简单问题的时间。你也可以使用random.uniform(0,0.1)或你认为合适的其他范围。

接下来，一个群体包含大量粒子。每个粒子是独立创建的，并且对个人最佳和全球最佳的不同维度有其独特的亲和力。在每次迭代中，所有粒子根据其‘速度’搜索附近的解决方案，并跟踪其速度，如本节开头所示。

class Swarm:
    def __init__(self, n_dim, n_population, param_limits):
        self.n_dim = n_dim
        self.particles = [Particle(i+1, n_dim, param_limits) for i in range(n_population)]
        self.reset()

    def reset(self):
        self.gbest = np.zeros(self.n_dim,)
        self.gbest_fitness = -1e8
        self.history = {}

    def update_all_particles(self, task):
        for particle_Obj in self.particles:
            particle_Obj.update(self.gbest, task)

    def solve(self, task, num_generation):
        self.reset()
        for t in tqdm(range(num_generation)):
            self.update_all_particles(task)
            particles_pbest_fitness = [task.score(p.pos) for p in self.particles]
            self.history[t] = dict([(p, deepcopy(p.pos)) for p in self.particles])
            if np.max(particles_pbest_fitness) > self.gbest_fitness:
                self.gbest_fitness = np.max(particles_pbest_fitness)
                self.gbest_index = np.argmax(particles_pbest_fitness)
                self.gbest_particle = self.particles[self.gbest_index]
                self.gbest = self.gbest_particle.pos                
        return self.gbest

实现可以如下完成。这里，种群包括 2000 个粒子作为示例。每个粒子都有一个噪声测量值，但噪声会从最终解决方案中去除，以获得算法性能的准确度量。注意，gbest不需要在每次迭代中更新。这是因为没有保证每次解决方案都会改进。

position_min, position_max = -5, 5
swarm = Swarm(
    n_dim=8, n_population=2000, param_limits=(position_min, position_max)
)
solution = swarm.solve(task, num_iteration=200)

print(solution)
print(task.score(solution, with_noise=False))

结果

仅需几秒钟即可找到解决方案。为了公平起见，我重复了十次实验，每次循环都删除swarm并重新初始化所有内容。结果不可避免地有所不同，因为每个粒子都被随机初始化，侧重于不同的组件。

进行 10 次 PSO，每次有 200 次迭代。

如果没有 PSO 或其他优秀的搜索算法，你的替代方案是进行暴力网格搜索。让我们看看这种方法与 PSO 的比较如何。

best_score = -np.inf
solution = None
intervals = 11

for x1 in tqdm(np.linspace(position_min, position_max, intervals)):
    for y1 in np.linspace(position_min, position_max, intervals):
        for x2 in np.linspace(position_min, position_max, intervals):
            for y2 in np.linspace(position_min, position_max, intervals):
                for x3 in np.linspace(position_min, position_max, intervals):
                    for y3 in np.linspace(position_min, position_max, intervals):
                        for x4 in np.linspace(position_min, position_max, intervals):
                            for y4 in np.linspace(position_min, position_max, intervals):
                                sample = task.score(np.array([x1,y1,x2,y2,x3,y3,x4,y4]))
                                if sample > best_score:
                                    best_score = sample
                                    solution = np.array([x1,y1,x2,y2,x3,y3,x4,y4])

print("solution: ", solution)
print(task.score(solution, with_noise=False))

即使在每个维度上仅进行 11 个点（即分成十分之一）的极其粗略搜索，也需要进行 10⁸次黑箱函数计算。在我的计算机上，这需要超过 20 分钟。请注意，这是在“最佳性能”模式下，并且背景中没有其他程序运行；甚至没有音乐（我在等待期间去喝了咖啡）。在现实中，每次计算可能涉及一个需要更长时间的模拟。

完成八维网格搜索所需的时间，以及对应的最佳解决方案和目标函数值。

仅使用极少的计算预算，PSO 解决了问题（中位数 10.61，最大 10.63），效果与蛮力网格搜索（最终值 10.61）一样好。值得注意的是，我们这里只涉及 8 维。如果涉及 1000 维等问题，网格搜索显然不可行，而 PSO 的好处将得到放大。

可视化

让我们看看粒子在迭代中的分布情况。尽管在 8 维中可视化事物是不可能的，我们可以每次关注两个维度，同时忽略其他组件。

在下面的图中，热图展示的是在仅改变两个维度的同时，保持其他 6 个维度为该特定迭代的swarm.gbest_particle值的目标值。因此，每个快照实际上仅适用于一个粒子（且没有实际方法同时呈现 2000 个热图）。尽管如此，它仍然提供了一些想法，胜于空白背景。

x_range = np.linspace(-5, 5, 100)
y_range = np.linspace(-5, 5, 100)
X, Y = np.meshgrid(x_range, y_range)
vectorized_blackbox = np.vectorize(blackbox)

for t in [1,2,3,4,5,6,8,10,12,15,20,30,50,100,200]:
    data = np.array(list(swarm.history[t-1].values()))
    ref = data[swarm.gbest_particle.idx-1,:]
    heatmap = {
        1: vectorized_blackbox(X, Y, ref[2], ref[3], ref[4], ref[5], ref[6], ref[7]),
        2: vectorized_blackbox(ref[0], ref[1], X, Y, ref[4], ref[5], ref[6], ref[7]),
        3: vectorized_blackbox(ref[0], ref[1], ref[2], ref[3], X, Y, ref[6], ref[7]),
        4: vectorized_blackbox(ref[0], ref[1], ref[2], ref[3], ref[4], ref[5], X, Y)
    }

    fig, axs = plt.subplots(2, 2, figsize=(10, 10))
    for i in range(2):
        for j in range(2):
            c = axs[i,j].contourf(X, Y, heatmap[2*i+j+1], vmin=0, vmax=12, levels=50, alpha=0.5)
            axs[i,j].scatter(data[:, 2*i+j], data[:, 2*i+j+1], s=2, alpha=0.4)
            axs[i,j].scatter(data[swarm.gbest_particle.idx-1, 2*i+j], data[swarm.gbest_particle.idx-1, 2*i+j+1], c='w', s=5, alpha=1)
            axs[i,j].set_xlabel("x_%d"%(2*i+j+1))
            axs[i,j].set_ylabel("y_%d"%(2*i+j+1))
            axs[i,j].set_xlim([-5, 5])
            axs[i,j].set_ylim([-5, 5])
    fig.colorbar(c, ax=axs)
    plt.suptitle("Distribution of particles at iteration %d"%t)
    plt.savefig("plot_%03d.png"%t)

使用imageio库创建 gif，如之前共享的，我们可以将所有图形组合成以下 gif。

白点表示swarm.gbest_particle的位置，而其他 1999 个粒子由较小的蓝点表示。

可能会很诱人急于得出结论，声称 PSO 代理是“愚蠢的”，只需稍微朝地图上更亮的区域移动即可。这样的论点是有缺陷的。我们看到的热图只是多维解决空间的一部分，正如 gif 中所示，当其他维度取不同的值时，它可以有很大变化。这不仅仅是贪婪地逐个优化组件维度的问题。如果采用了不同的表示组合，比如说x₁与y₃，或x₄与y₂，情况会有所不同。没有一种万全的技术可以同时让所有人满意。

最终，结果说明了一切，PSO 在几秒钟内解决问题的事实证明了其有效性。

结论

在这篇文章中，我们探讨了如何高效地实现 PSO 来寻找未知黑箱的最优解，并了解了实现这一目标的幕后过程。

我设想继续进行一系列以自然为灵感的算法。毕竟，自然是美丽的。数学是美丽的。自然中的数学？这必定是最棒的。

在我的下一篇文章中，我将解释 ABC（人工蜂群）的实现方法，并与 PSO 进行比较，探讨 ABC 在何种问题上表现更佳。敬请关注。

参考文献

[1] A. E. Eiben 和 J. E. Smith, 进化计算导论 (2015), Springer-Verlag Berlin Heidelberg

PatchTST：时间序列预测中的突破

原文：towardsdatascience.com/patchtst-a-breakthrough-in-time-series-forecasting-e02d48869ccc

从理论到实践，了解 PatchTST 算法，并在 Python 中与 N-BEATS 和 N-HiTS 一起应用。

Marco Peixeiro

·发表于 Towards Data Science ·阅读时间 10 分钟·2023 年 6 月 20 日

–

图片由 Ray Hennessy 提供，来源于 Unsplash

基于变换器的模型已成功应用于许多领域，如自然语言处理（例如 BERT 或 GPT 模型）和计算机视觉等。

然而，在时间序列方面，最先进的结果主要由 MLP 模型（多层感知器）如 N-BEATS 和 N-HiTS 实现。最近的一篇论文甚至表明，简单的线性模型在许多基准数据集上优于复杂的基于变换器的预测模型（见 Zheng 等人, 2022）。

尽管如此，一种新的基于变换器的模型已经被提出，并在长期预测任务中取得了最先进的结果：PatchTST。

PatchTST 代表补丁时间序列变换器，它首次由 Nie、Nguyen 等人在 2023 年 3 月提出，详细介绍在他们的论文中：时间序列值得 64 个词：使用变换器进行长期预测。他们提出的方法在与其他基于变换器的模型比较时，取得了最先进的结果。

在这篇文章中，我们首先直观地探索 PatchTST 的内部工作原理，而不使用任何方程式。接着，我们将该模型应用于预测项目，并将其性能与 MLP 模型（如 N-BEATS 和 N-HiTS）进行比较和评估。

当然，关于 PatchTST 的更多细节，请参考 原始论文。

使用我的 免费时间序列备忘单 在 Python 中学习最新的时间序列分析技术！获取统计和深度学习技术的实现，全部使用 Python 和 TensorFlow！

让我们开始吧！

探索 PatchTST

如前所述，PatchTST 代表补丁时间序列 transformer。

正如名字所示，它利用了补丁和 transformer 架构。它还包括通道独立性来处理多变量时间序列。下面展示了总体架构。

PatchTST 模型架构。我们看到该模型利用通道独立性来处理多变量时间序列。在 transformer 主干中，我们还看到使用了补丁（由矩形表示）。此外，模型有两个版本：监督学习和自监督学习。图片由 Nie Y.、Nguyen N.、Sinthong P.、Kalagnanam J. 提供，来源于 “一个时间序列值 64 个词：使用 Transformer 进行长期预测”。

从上图中需要收集大量信息。这里的关键要素是，PatchTST 使用通道独立性来预测多变量时间序列。然后，在其 transformer 主干中，该模型使用补丁，这些补丁由小的垂直矩形表示。此外，该模型有两个版本：监督学习和自监督学习。

让我们更详细地探讨 PatchTST 的架构和内部工作原理。

通道独立性

在这里，多变量时间序列被视为多通道信号。每个时间序列基本上是一个包含信号的通道。

PatchTST 模型概述。在这里，我们特别强调通道独立性的实现，其中每个输入 token 到 Transformer 主干仅包含来自一个通道或一个时间序列的信息。图片由 Nie Y.、Nguyen N.、Sinthong P.、Kalagnanam J. 提供，来源于 “一个时间序列值 64 个词：使用 Transformer 进行长期预测”。

在上图中，我们看到如何将多变量时间序列分离成单独的序列，并将每个序列作为输入 token 传递给 Transformer 主干。然后，对每个序列进行预测，结果被连接起来形成最终预测。

补丁

大多数基于 Transformer 的预测模型的工作集中在构建新的机制来简化原始注意力机制。然而，它们仍然依赖于逐点注意力，这在时间序列的情况下并不理想。

在时间序列预测中，我们希望提取过去时间步和未来时间步之间的关系以进行预测。使用逐点注意力时，我们试图从单个时间步中检索信息，而不考虑该点周围的内容。换句话说，我们孤立一个时间步，而不查看之前或之后的点。

这就像试图理解一个词的意义而不看句子中的其他词一样。

因此，PatchTST 利用分块来提取时间序列中的局部语义信息。

分块如何工作

每个输入序列被分成块，这些块只是来自原始序列的较短序列。

PatchTST 的 Transformer 主干。这里，我们看到输入时间序列（图底部）经过分块，结果是多个块（垂直矩形），然后被送到 Transformer 编码器。图片由 Nie Y.、Nguyen N.、Sinthong P.、Kalagnanam J. 提供，来自 A Time Series is Worth 64 Words: Long-Term Forecasting with Transformers。

这里，块可以是重叠的也可以是不重叠的。块的数量取决于块的长度 P 和步幅 S。这里，步幅类似于卷积，它只是分隔连续块开始的时间步数。

可视化分块。这里，我们有一个 15 个时间步的序列，分块长度为 5，步幅也为 5，结果是三个块。图片由作者提供。

在上图中，我们可以可视化分块的结果。这里，我们有一个序列长度 (L) 为 15 个时间步，分块长度 (P) 为 5，步幅 (S) 为 5。结果是序列被分为 3 个块。

分块的优势

通过分块，模型可以通过查看时间步的组而不是单个时间步来提取局部语义信息。

它还有一个额外的好处，就是大大减少了送入 Transformer 编码器的标记数量。这里，每个块成为输入到 Transformer 的输入标记。这样，我们可以将标记数量从 L 减少到大约 L/S。

这样，我们大大减少了模型的空间和时间复杂度。这反过来意味着我们可以给模型输入更长的序列，以提取有意义的时间关系。

因此，通过分块，模型变得更快、更轻，并且可以处理更长的输入序列，这意味着它可能会对序列学习得更多并做出更好的预测。

Transformer 编码器

一旦序列被分块，它会被送入 Transformer 编码器。这是经典的 Transformer 架构，没有任何修改。

然后，将输出送入线性层，并进行预测。

使用表示学习改进 PatchTST

论文的作者建议通过使用表示学习来改进模型。

在 PatchTST 中可视化自监督表示学习。在这里，模型将随机掩盖补丁并学习重建它们。图片由 Nie Y.、Nguyen N.、Sinthong P.、Kalagnanam J.提供，来自A Time Series is Worth 64 Words: Long-Term Forecasting with Transformers。

从上图中，我们可以看到 PatchTST 可以利用自监督表示学习来捕捉数据的抽象表示。这可能会导致预测性能的潜在改善。

在这里，过程相当简单，因为随机补丁将被掩盖，意味着它们将被设置为 0。这在上图中由空白的垂直矩形表示。然后，模型被训练以重建原始补丁，这就是图顶端输出的灰色垂直矩形。

现在我们对 PatchTST 的工作原理有了很好的理解，让我们将其与其他模型进行测试，看看其表现如何。

使用 PatchTST 进行预测

在论文中，PatchTST 与其他基于 Transformer 的模型进行了比较。然而，最近也发布了基于 MLP 的模型，如 N-BEATS 和 N-HiTS，它们在长期预测任务中也展示了最先进的性能。

本节的完整源代码可在GitHub上找到。

在这里，让我们应用 PatchTST，以及 N-BEATS 和 N-HiTS，并评估它在这两个基于 MLP 的模型中的表现。

对于这个练习，我们使用 Exchange 数据集，这是研究中用于长期预测的常见基准数据集。该数据集包含 1990 年至 2016 年期间相对于美元的八个国家的每日汇率。该数据集通过 MIT 许可证提供。

初始设置

让我们首先导入所需的库。在这里，我们将使用neuralforecast，因为它们有一个现成的 PatchTST 实现。对于数据集，我们使用datasetsforecast库，其中包含所有流行的数据集，用于评估预测算法。

import torch
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from neuralforecast.core import NeuralForecast
from neuralforecast.models import NHITS, NBEATS, PatchTST

from neuralforecast.losses.pytorch import MAE
from neuralforecast.losses.numpy import mae, mse

from datasetsforecast.long_horizon import LongHorizon

如果你已经安装了 CUDA，那么neuralforecast将自动利用你的 GPU 来训练模型。在我这边，我没有安装 CUDA，这就是为什么我没有进行广泛的超参数调整或在非常大的数据集上训练的原因。

完成这一步后，让我们下载 Exchange 数据集。

Y_df, X_df, S_df = LongHorizon.load(directory="./data", group="Exchange")

在这里，我们看到得到三个 DataFrames。第一个包含每个国家的每日汇率。第二个包含外生时间序列。第三个包含静态外生变量（如日期、月份、年份、小时或我们已知的任何未来信息）。

对于这个练习，我们只使用Y_df。

然后，确保日期具有正确的类型。

Y_df['ds'] = pd.to_datetime(Y_df['ds'])

Y_df.head()

Exchange 数据集的前五行。作者提供的图片。

在上图中，我们看到有三列。第一列是唯一标识符，当使用neuralforecast时，需要有一个 id 列。然后，ds列有日期，y列有汇率。

Y_df['unique_id'].value_counts()

显示每个唯一 id 的观测数量。图片由作者提供。

从上图可以看出，每个唯一 id 对应一个国家，我们每个国家有 7588 个观测值。

现在，我们定义验证集和测试集的大小。在这里，我选择了 760 个时间步作为验证集，1517 个时间步作为测试集，如[datasets](https://github.com/Nixtla/datasetsforecast/blob/main/datasetsforecast/long_horizon.py) 库所规定。

val_size = 760
test_size = 1517

print(n_time, val_size, test_size)

然后，让我们绘制其中一个序列，看看我们在处理什么。在这里，我决定绘制第一个国家的序列（unique_id = 0），但可以自由绘制其他序列。

u_id = '0'

x_plot = pd.to_datetime(Y_df[Y_df.unique_id==u_id].ds)
y_plot = Y_df[Y_df.unique_id==u_id].y.values

x_plot

x_val = x_plot[n_time - val_size - test_size]
x_test = x_plot[n_time - test_size]

fig, ax = plt.subplots(figsize=(12,8))
ax.plot(x_plot, y_plot)
ax.set_xlabel('Date')
ax.set_ylabel('Exhange rate')
ax.axvline(x_val, color='black', linestyle='--')
ax.axvline(x_test, color='black', linestyle='--')

plt.text(x_val, -2, 'Validation', fontsize=12)
plt.text(x_test,-2, 'Test', fontsize=12)

plt.tight_layout()

第一个国家从 1990 年到 2016 年的每日汇率。图片由作者提供。

从上图中，我们看到数据相当嘈杂，没有明显的季节性。

建模

探索了数据后，让我们开始使用neuralforecast进行建模。

首先，我们需要设置预测期。在这种情况下，我使用 96 个时间步，因为这个预测期也在PatchTST 论文中使用。

然后，为了公平评估每个模型，我决定将输入大小设置为预测期的两倍（即 192 个时间步），并将最大训练轮数设置为 50。所有其他超参数保持默认值。

horizon = 96

models = [NHITS(h=horizon,
               input_size=2*horizon,
               max_steps=50),
         NBEATS(h=horizon,
               input_size=2*horizon,
               max_steps=50),
         PatchTST(h=horizon,
                 input_size=2*horizon,
                 max_steps=50)]

然后，我们初始化NeuralForecast对象，指定我们要使用的模型和预测频率，在这里是每日。

nf = NeuralForecast(models=models, freq='D')

我们现在准备进行预测。

预测

为了生成预测，我们使用cross_validation方法利用验证集和测试集。它将返回一个包含所有模型预测值及相关真实值的 DataFrame。

preds_df = nf.cross_validation(df=Y_df, val_size=val_size, test_size=test_size, n_windows=None)

预测 DataFrame 的前五行。图片由作者提供。

如您所见，对于每个 id，我们有来自每个模型的预测以及y列中的真实值。

现在，为了评估模型，我们需要将实际值和预测值的数组重塑为形状为(series 的数量, 窗口数量, 预测期)。

y_true = preds_df['y'].values
y_pred_nhits = preds_df['NHITS'].values
y_pred_nbeats = preds_df['NBEATS'].values
y_pred_patchtst = preds_df['PatchTST'].values

n_series = len(Y_df['unique_id'].unique())

y_true = y_true.reshape(n_series, -1, horizon)
y_pred_nhits = y_pred_nhits.reshape(n_series, -1, horizon)
y_pred_nbeats = y_pred_nbeats.reshape(n_series, -1, horizon)
y_pred_patchtst = y_pred_patchtst.reshape(n_series, -1, horizon)

完成后，我们可以选择性地绘制模型的预测。在这里，我们绘制了第一个序列的第一个窗口中的预测。

fig, ax = plt.subplots(figsize=(12,8))

ax.plot(y_true[0, 0, :], label='True')
ax.plot(y_pred_nhits[0, 0, :], label='N-HiTS', ls='--')
ax.plot(y_pred_nbeats[0, 0, :], label='N-BEATS', ls=':')
ax.plot(y_pred_patchtst[0, 0, :], label='PatchTST', ls='-.')
ax.set_ylabel('Exchange rate')
ax.set_xlabel('Forecast horizon')
ax.legend(loc='best')

plt.tight_layout()

第一个序列在第一个窗口中的每日汇率预测。图片由作者提供。

这个图有点让人失望，因为 N-BEATS 和 N-HiTS 的预测结果与实际值差异很大。然而，尽管 PatchTST 也有偏差，但似乎最接近实际值。

当然，我们必须对此持保留态度，因为我们只是可视化了一个序列中的一个预测窗口。

评估

所以，让我们评估每个模型的表现。为了复制论文中的方法，我们使用 MAE 和 MSE 作为性能指标。

data = {'N-HiTS': [mae(y_pred_nhits, y_true), mse(y_pred_nhits, y_true)],
       'N-BEATS': [mae(y_pred_nbeats, y_true), mse(y_pred_nbeats, y_true)],
       'PatchTST': [mae(y_pred_patchtst, y_true), mse(y_pred_patchtst, y_true)]}

metrics_df = pd.DataFrame(data=data)
metrics_df.index = ['mae', 'mse']

metrics_df.style.highlight_min(color='lightgreen', axis=1)

所有模型的性能。在这里，PatchTST 实现了最低的 MAE 和 MSE。图片由作者提供。

在上表中，我们看到 PatchTST 是冠军模型，因为它实现了最低的 MAE 和 MSE。

当然，这不是最彻底的实验，因为我们只使用了一个数据集和一个预测范围。尽管如此，看到一个基于 Transformer 的模型能够与最先进的 MLP 模型竞争仍然很有趣。

结论

PatchTST 是基于 Transformer 的模型，通过补丁来提取时间序列数据中的局部语义意义。这使得模型的训练速度更快，并且具有更长的输入窗口。

相较于其他基于 Transformer 的模型，它已经达到了最先进的性能。在我们的练习中，我们看到它的表现也优于 N-BEATS 和 N-HiTS。

虽然这并不意味着它比 N-HiTS 或 N-BEATS 更好，但它仍然是一个有趣的选项，尤其是在长期预测时。

感谢阅读！希望你喜欢这篇文章，并且学到了一些新东西！

想要掌握时间序列预测？那就看看Python 应用时间序列预测。这是唯一一个通过 15 个引导实践项目使用 Python 实现统计学、深度学习和最先进模型的课程。

干杯 🍻

支持我

喜欢我的工作吗？通过给我买杯咖啡来表示支持，这是鼓励我的简单方式，我可以享受一杯咖啡！如果你愿意，请点击下面的按钮 👇

参考文献

时间序列的价值为 64 个词：使用 Transformer 进行长期预测 作者：Nie Y., Nguyen N. 等。

Neuralforecast 作者：Olivares K., Challu C., Garza F., Canseco M., Dubrawski A.

多变量时间序列的主成分分析：动态高维数据的预测

原文：towardsdatascience.com/pca-for-multivariate-time-series-forecasting-dynamic-high-dimensional-data-ab050a19e8db

噪声和序列相关性存在下的系统预测

Marco Cerliani

·发表在数据科学前沿 ·5 分钟阅读·2023 年 1 月 31 日

–

图片由Viva Luna Studios提供，来源于Unsplash

多步骤的多变量时间序列预测被认为是一个复杂的预测任务。我们必须考虑到输入和输出的高维度；必须充分处理横截面和时间依赖关系；最后但同样重要的是，我们必须确保长时间准确性的可接受水平。

如今，处理时间和数量维度巨大数据的分析应用非常普遍。因此，所有基于这些系统构建的解决方案必须能够操作大数据集。在物联网（IoT）时代，处理大量时间序列是很常见的，这些时间序列在大多数情况下展示了强相关模式。这些动态在电信、工业制造、金融、电网等领域非常常见。

假设你是一个数据科学家。我们负责开发一个预测分析应用，提供多步骤的预测，针对由相关和噪声传感器组成的物联网系统。多变量预测在文献中广泛讨论。从 VAR（向量自回归模型）等统计技术到更复杂的基于深度学习的方法，存在许多可用于预测任务的解决方案。然而，现实世界比预期的更复杂和残酷。实时管理大量高频传感器需要开发兼具适当准确度和合理响应延迟的解决方案。

在这篇文章中，我们尝试开发一个用于多变量和多步骤传感器预测的预测应用程序，该应用程序可以在接近实时模式下使用。这是通过将降维与适合多变量上下文的预测技术相结合来实现的。所提出的方法在经济预测文献中很受欢迎，被称为动态因子建模[1]。换句话说，我们将最喜欢的预测算法叠加在降维技术（如 PCA）的结果之上，以预测未来的系统动态。

实验设置

在本文范围内，我们生成了多个合成时间序列。这些序列可以根据其正弦动态分成两组。一切都很完美，只是信号被噪声掩盖了。

数据中的合成动态 [image by the author]

我们在时间序列中加入了大量噪声，以复制现实世界系统的混沌行为，使预测任务变得更加困难。

带有噪声的模拟序列预测 [image by the author]

拥有多个时间序列可供使用，我们的目标是预测它们未来的多个步骤。为了以接近实时的模式进行多变量预测任务，我们应该在预测准确性和推理过程的持续时间之间找到一个折衷。

让我们看看如何解决这个问题。

多变量动态预测

动态因子建模（DFM）是一种源自经济文献的多变量预测技术[1]。DFM 的基本思想是少量系列可以解释大量变量的时间行为。如果我们可以获得这些因素的准确估计，整个预测任务可以通过使用估计的动态因素来简化，而不是使用所有系列。

动态因子建模估计流程 [image by the author]

使用 DFM 获得的预测质量取决于两个主要方面：因素估计的优良性和因素预测的准确性。估计动态因素的方式有很多种。在机器学习生态系统中最常见的方法是通过正交旋转（PCA）获得一组数据的主成分[2]。

DFM 也是一种与模型无关的技术。换句话说，任何降维和任何预测策略都可以用于进行预测。对于我们的实验，我们使用标准 PCA 和朴素直接预测。下面是进行 DFM 估计和预测的代码快照。

scaler_pca = make_pipeline(StandardScaler(), PCA(n_components=2))

X_train_factors = scaler_pca.fit_transform(X_train)

forecaster = ForecastingChain(
    Ridge(),
    n_estimators=test_size,
    lags=range(1,25),
    use_exog=False,
    n_jobs=-1,
).fit(None, X_train_factors)

y_pred_factors = forecaster.predict(np.arange(test_size))

y_pred = scaler_pca.steps[0][-1].inverse_transform(
    scaler_pca.steps[-1][-1].inverse_transform(y_pred_factors)
)

最后的步骤是通过简单地对我们所有可用的系列采用多变量直接预测来解决相同的任务。无论结果如何，这种方法都不可持续，因为它需要对所有可用系列进行滞后特征计算。这可能需要处理大量的滞后变量，这使得大多数系统由于物理（内存）和时间限制而无法进行预测。

我们采用了时间交叉验证策略来验证这两种方法的结果并存储性能数据。

两种方法预测结果的视觉比较 [作者提供的图像]

DFM 超越了幼稚的多变量直接预测。在较短时间内（推断/估计时间依赖于系统变量的数量）实现了更高的准确性。检查生成的预测，我们可以观察到 DFM 能够正确区分并再现原始系列中存在的双正弦波动态。

多变量预测性能比较 [作者提供的图像]

总结

在这篇文章中，我们提出了动态因子建模的实际应用。它被证明是一种有效的多变量时间序列预测建模方法。它特别适用于预测高维数据，同时也显示出可能的高噪声水平。正如往常一样，不存在适用于所有情况的完美预测技术。作为数据科学家，我们有责任尝试以前未曾了解的技术。只有通过持续的自我学习，我们才能选择和区分适用于日常任务的最佳解决方案。

查看我的 GitHub 仓库

保持联系: Linkedin

参考文献

[1] M. Forni, M. Hallin, M. Lippi 和 L. Reichlin, “广义动态因子模型”，《美国统计学会杂志》，第 100 卷，第 471 期，第 830–840 页。

[2] G. Bontempi, Y. -A. Le Borgne 和 J. de Stefani, “动态因子机器学习方法用于多变量和多步预测”，2017 年 IEEE 国际数据科学与高级分析会议 (DSAA)，东京，日本，2017 年，第 222–231 页。

PCA/LDA/ICA：组件分析算法比较

原文：towardsdatascience.com/pca-lda-ica-a-components-analysis-algorithms-comparison-c5762c4148ff

回顾这些著名算法的概念和差异

Yoann Mocquin

·发表于Towards Data Science ·8 分钟阅读·2023 年 2 月 19 日

–

在深入比较算法之前，让我们先独立回顾一下它们。

注意，本文章并不旨在深入解释每个算法，而是比较它们的目标和结果。

如果你想了解更多关于 PCA 和 ZCA 的区别，请查看我之前基于 numpy 的文章：

PCA-whitening vs ZCA-whitening：一个 numpy 2D 可视化

数据白化的过程包括一种变换，使得变换后的数据具有单位矩阵作为…

[towardsdatascience.com

PCA：主成分分析

• PCA 是一种无监督的线性降维技术，旨在找到一组新的正交变量，这些变量捕捉数据中最重要的变异来源。

• 它广泛用于特征提取和数据压缩，并且可以用于探索性数据分析或作为机器学习算法的预处理步骤。

• 生成的组件根据它们解释的方差量进行排名，可以用来可视化和解释数据，以及用于聚类或分类任务。

LDA：线性判别分析

• LDA 是一种监督式线性降维技术，旨在找到一组新的变量，这些变量最大限度地提高类之间的分离，同时最小化每个类内的变化。

• PCA 广泛用于特征提取和分类，并且可以在保留类间判别信息的同时减少数据的维度。

• 结果组件按其区分能力排序，可以用来可视化和解释数据，以及进行分类或回归任务。

ICA：独立成分分析

• ICA 是一种无监督的线性降维技术，旨在找到一组统计上独立且非高斯的变量。

• 它广泛用于信号处理和源分离，可以用来提取数据中其他技术无法访问的潜在变异源。

• 结果组件按其独立性排序，可以用来可视化和解释数据，以及进行聚类或分类任务。

鸢尾花数据集的结果

让我们使用 sklearn 比较它们在著名的鸢尾花数据集上的结果。首先，使用每个 4 个数值特征的 pairplot 绘制鸢尾花数据集，并将颜色作为类别特征：

import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris

# Load the iris dataset
iris = load_iris()
data = iris.data
target = iris.target
target_names = iris.target_names

# Convert the iris dataset into a pandas DataFrame
iris_df = sns.load_dataset('iris')
iris_df['target'] = target

# Generate the pairplot∑
sns.pairplot(data=iris_df, hue='target', palette=['navy', 'turquoise', 'darkorange'], markers=['o', 's', 'D'],
             plot_kws=dict(s=25, alpha=0.8, edgecolor='none'), diag_kws=dict(alpha=0.8, edgecolor='none'))

# Set the title and adjust plot spacing
plt.suptitle('Iris Pairplot')
plt.subplots_adjust(top=0.92)

plt.show()

作者提供的图像：鸢尾花数据集的 pairplot

现在我们可以计算每个变换并绘制结果。注意我们只使用了 2 个组件，因为 LDA 最多需要 (N-1) 个组件，其中 N 是类别的数量（这里等于 3，因为有 3 种鸢尾花类型）。

from sklearn.datasets import load_iris
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA, FastICA
import matplotlib.pyplot as plt

# Load the Iris dataset
iris = load_iris()
X = iris.data
y = iris.target
target_names = iris.target_names

# Standardize the data
scaler = StandardScaler()
X_std = scaler.fit_transform(X)

# Apply LDA with 2 components
lda = LinearDiscriminantAnalysis(n_components=2)
X_lda = lda.fit_transform(X_std, y)

# Apply PCA with 2 components
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X_std)

# Apply ICA with 2 components
ica = FastICA(n_components=2)
X_ica = ica.fit_transform(X_std)

# Plot the results
plt.figure(figsize=(15, 5))

plt.subplot(1, 3, 1)
for target, color in zip(range(len(target_names)), ['navy', 'turquoise', 'darkorange']):
    plt.scatter(X_lda[y == target, 0], X_lda[y == target, 1], color=color, alpha=.8, lw=2,
                label=target_names[target])
plt.legend(loc='best', shadow=False, scatterpoints=1)
plt.title('LDA')
plt.xlabel('LD1')
plt.ylabel('LD2')

plt.subplot(1, 3, 2)
for target, color in zip(range(len(target_names)), ['navy', 'turquoise', 'darkorange']):
    plt.scatter(X_pca[y == target, 0], X_pca[y == target, 1], color=color, alpha=.8, lw=2,
                label=target_names[target])
plt.legend(loc='best', shadow=False, scatterpoints=1)
plt.title('PCA')
plt.xlabel('PC1')
plt.ylabel('PC2')

plt.subplot(1, 3, 3)
for target, color in zip(range(len(target_names)), ['navy', 'turquoise', 'darkorange']):
    plt.scatter(X_ica[y == target, 0], X_ica[y == target, 1], color=color, alpha=.8, lw=2,
                label=target_names[target])
plt.legend(loc='best', shadow=False, scatterpoints=1)
plt.title('ICA')
plt.xlabel('IC1')
plt.ylabel('IC2')

plt.show()

这段代码加载了鸢尾花数据集，应用了 LDA、PCA 和 ICA，每种方法使用 2 个组件，然后使用不同的颜色绘制每个类别的结果。

注意，通常在应用 PCA、ICA 或 LDA 之前对数据进行标准化是一个好习惯。标准化很重要，因为这些技术对输入特征的尺度很敏感。标准化数据可以确保每个特征具有零均值和单位标准差，这使所有特征处于相同的尺度上，避免了某一特征对其他特征的支配。

由于 LDA 是一种有监督的降维技术，它以类别标签作为输入。相比之下，PCA 和 ICA 是无监督的技术，这意味着它们仅使用输入数据而不考虑类别标签。

LDA 的结果可以解释为将数据投影到一个最大化类别分离的空间，而 PCA 和 ICA 的结果可以解释为将数据投影到一个分别捕捉最重要的变异源或独立源的空间。

作者提供的图像：LDA、PCA 和 ICA 在鸢尾花数据集上的比较

注意，尽管 ICA 的目的不是分离类别，但它仍然显示了类别之间的分离：这是因为类别在输入数据集中已经相当有序。

让我们暂时放下 LDA，关注 PCA 和 ICA 之间的区别——因为 LDA 是一种有监督的技术，专注于分离类别并强制最大组件数量，而 PCA 和 ICA 侧重于创建一个与输入矩阵形状相同的新矩阵。

让我们查看 PCA 和 ICA 的 4 个组件的输出：

左：PCA 的配对图 / 右：ICA 的配对图（图片由作者提供）

我们还可以比较每个转换数据的相关性矩阵：注意到两种方法都产生了不相关的向量（换句话说，转换后的数据特征是正交的）。这是因为这是 PCA 算法中的一个约束——每个新向量必须正交于之前的向量——而 ICA 算法的一个结果——这意味着原始数据集是混合在一起的独立信号，必须被重构。

左：ICA 的相关性热图 / 右：PCA 的相关性热图（图片由作者提供）

所以 PCA 和 ICA 似乎给出了具有类似属性的结果：这是由于以下两个原因：

• 独立性在这两种算法中都“编码”了

• 鸢尾花数据集展示了良好的类分离

这就是为什么我们需要另一个例子，更适合 ICA 的原因。

另一个例子：

让我们看另一个例子：我们首先生成一个合成数据集，其中包含两个独立的源，一个正弦波和一个方波，它们被线性组合在一起以创建一个混合信号。

实际的、真正的、独立的信号如下：

它们被线性组合在一起，作为 2 个线性组合：

让我们看看 PCA 和 ICA 在这个新数据集上的表现：

注意 PCA 如何创建了一个新的组件，它作为输入的线性组合展示了大量的方差，但这绝对不匹配原始数据：这确实不是 PCA 的目的。

相反，ICA 在恢复原始数据集方面表现非常好，无论方差组成如何。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import FastICA

# Generate a synthetic dataset with two independent sources
np.random.seed(0)
n_samples = 2000
time = np.linspace(0, 8, n_samples)

s1 = np.sin(2 * time) # Source 1: sine wave
s2 = np.sign(np.sin(3 * time)) # Source 2: square wave

S = np.c_[s1, s2]
S += 0.2 * np.random.normal(size=S.shape) # Add noise to the sources
S /= S.std(axis=0) # Standardize the sources

# Mix the sources together to create a mixed signal
A = np.array([[0.5, 0.5], [0.2, 0.8]]) # Mixing matrix
X = np.dot(S, A.T) # Mixed signal

# Standardize the data
X = (X - np.mean(X, axis=0)) / np.std(X, axis=0)

# Use PCA to reduce the dimensionality of the data
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X)

# Use ICA to separate the sources from the mixed signal
ica = FastICA(n_components=2)
X_ica = ica.fit_transform(X) # Estimated sources

# Plot the results
plt.figure()

models = [X, S, X_pca, X_ica]
names = ['Observations (mixed signal)',
         'True Sources',
         'PCA features', 'ICA estimated sources']
colors = ['red', 'steelblue']

for ii, (model, name) in enumerate(zip(models, names), 1):
    plt.subplot(4, 1, ii)
    plt.title(name)
    for sig, color in zip(model.T, colors):
        plt.plot(sig, color=color)

plt.tight_layout()
plt.show()

结论

PCA、LDA 和 ICA 算法可能看起来像是彼此的定制版本，但它们的目的实际上并不相同。总结如下：

• PCA 旨在创建新的组件，这些组件包含输入数据的最大方差

• LDA 旨在创建基于类别特征分离簇的新组件

• ICA 旨在恢复在输入数据集中线性组合在一起的原始特征

希望你能更好地理解这些算法之间的差异，并能够在未来快速识别你需要的算法。

如果你喜欢这个故事，确保关注我并帮助我实现 100 个订阅者的目标 😃

查看下面的一些其他故事：

[## 使用 numpy 的有限差分法实现 300 倍更快的分辨率

有限差分法是一种强大的技术来解决复杂问题，numpy 使它变得快速！

使用 numpy 提升有限差分方法的分辨率 300 倍 [## 使用 numpy 提升有限差分方法的分辨率 300 倍

容器方法。

深入了解 seaborn 调色板 [## 深入了解 seaborn 调色板

被 seaborn 调色板淹没了吗？

PCA-whitening 与 ZCA-whitening：一个 numpy 2d 可视化 [## PCA-whitening 与 ZCA-whitening：一个 numpy 2d 可视化

数据白化的过程包括一种变换，使得变换后的数据具有单位矩阵作为…

包裹 numpy 的数组 [## 包裹 numpy 的数组

PCA 与自编码器在小数据集上的降维比较

原文：towardsdatascience.com/pca-vs-autoencoders-for-a-small-dataset-in-dimensionality-reduction-67b15318dea0

神经网络和深度学习课程：第四十五部分

Rukshan Pramoditha

·发表在 Towards Data Science ·阅读时间 8 分钟·2023 年 2 月 16 日

–

图片由 Robert Katzki 提供，来源于 Unsplash

小数据集上，通用机器学习算法是否能超越神经网络？

一般而言，深度学习算法如神经网络需要大量数据才能实现合理的性能。因此，像自编码器这样的神经网络可以从我们用来训练模型的非常大的数据集中获益。

有时，当通用机器学习算法在非常小的数据集上训练时，它们的表现可能会超越神经网络算法。

自编码器也可以用于降维应用，尽管它们在图像去噪、图像生成、图像着色、图像压缩、图像超分辨率等其他流行应用中被广泛使用。

早些时候，我们通过在非常大的 MNIST 数据集上训练模型，比较了自编码器在降维中的表现与 PCA 的表现。在那里，自编码器模型轻松超越了 PCA 模型 [参考¹]，因为 MNIST 数据集大且具有非线性特征。

参考¹：自编码器如何在降维中超越 PCA

自编码器在处理大规模和非线性数据时表现良好。

尽管自编码器是一种神经网络，但仍然可以使用较小的数据集，只要选择正确的模型架构和超参数值，就能够取得一定的性能。

自编码器功能强大且灵活，足以捕捉数据中的复杂和非线性模式。

今天，我们将通过在非常小的葡萄酒数据集上训练模型，比较自编码器（神经网络模型）在降维方面与 PCA（通用机器学习算法）的表现。

葡萄酒数据集有 178 个样本和 13 个特征。与 MNIST 数据集的 60,000 个样本和 784 个特征相比，这个数据集非常小！

加载葡萄酒数据集

葡萄酒数据集已经预加载到 Scikit-learn 中，可以按如下方式加载。

# Load the Wine dataset
from sklearn.datasets import load_wine

wine = load_wine()
X = wine.data
y = wine.target

print("Wine dataset size:", X.shape)

葡萄酒数据集的形状（作者提供的图片）

使用葡萄酒数据运行 PCA

首先，我们需要通过使用所有组件运行 PCA 来选择最佳组件数量 [ref²]。

ref²: 进行主成分分析 (PCA) 以实现降维的 3 个简单步骤

# Feature scaling
from sklearn.preprocessing import StandardScaler

X_scaled = StandardScaler().fit_transform(X)

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA

# Select the best number of components
# by running PCA with all components
pca = PCA(n_components=None)
pca.fit(X_scaled)

# Plot cumulative explained variances
exp_var = pca.explained_variance_ratio_ * 100
cum_exp_var = np.cumsum(exp_var)

plt.bar(range(1, 14), exp_var, align='center',
        label='Individual explained variance')

plt.step(range(1, 14), cum_exp_var, where='mid',
         label='Cumulative explained variance', color='red')

plt.ylabel('Explained variance percentage')
plt.xlabel('Principal component index')
plt.xticks(ticks=list(range(1, 14)))
plt.legend(loc='best')
plt.tight_layout()

plt.savefig("cumulative_explained_variance_plot.png")

累计解释方差图（作者提供的图片）

通过分析上述累计解释方差图，我决定保留前七个组件，这些组件捕获了数据中约 90% 的方差。因此，这七个组件将准确地代表原始葡萄酒数据集。

使用选定的组件重新运行 PCA

要将降维应用于葡萄酒数据集，我们需要使用选定的组件重新运行 PCA 并应用转换。

# Run PCA again with selected (7) components
pca = PCA(n_components=7)
X_pca = pca.fit_transform(X_scaled)
print("PCA reduced wine dataset size:", X_pca.shape)

应用 PCA 后的葡萄酒数据集的形状（作者提供的图片）

现在，数据集中只有七个组件（特征）。因此，数据的维度已经减少！

使用自编码器进行降维

以下代码定义了一个自编码器，其中编码器部分可以用于获得数据的低维（编码）表示。

import numpy as np
from tensorflow.keras import Model, Input
from tensorflow.keras.layers import Dense

# Build the autoencoder
input_dim = X.shape[1]
latent_vec_dim = 7

input_layer = Input(shape=(input_dim,))

# Define encoder
x = Dense(8, activation='relu')(input_layer)
x = Dense(4, activation='relu')(x)
encoder = Dense(latent_vec_dim, activation="tanh")(x)

# Define decoder
x = Dense(4, activation='relu')(encoder)
x = Dense(8, activation='relu')(x)
decoder = Dense(input_dim, activation="sigmoid")(x)

autoencoder = Model(inputs=input_layer, outputs=decoder)

# Compile the model with optimizer and loss function
autoencoder.compile(optimizer="adam", loss="mse")

# Train the model with 100 epochs
autoencoder.fit(X_scaled, X_scaled, epochs=100, verbose=0,
                batch_size=16, shuffle=True)

# Use the encoder part to obtain the lower dimensional,
# encoded representation of the data
encoder_model = Model(inputs=input_layer, outputs=encoder)
X_encoded = encoder_model.predict(X_scaled)

# Print the shape of the encoded data
print("Autoencoder reduced wine dataset size:", X_encoded.shape)

自编码数据后的葡萄酒数据集的形状（作者提供的图片）

自编码器模型具有输入层、编码层和解码层。输入维度是葡萄酒数据集中的输入特征数，即 13。潜在向量维度为 7，等于我们在 PCA 中先前选择的组件数量。

所有层都使用 Keras 功能 API 方法 [ref³] 连接。然后，通过连接输入层和解码器部分创建整个自编码器模型。

ref³: 构建 Keras 模型的两种不同方式：顺序 API 和功能 API

然后，我们使用 Adam 优化器和均方误差（mse）损失函数编译整个自编码器。

模型在标准化（缩放）的葡萄酒数据上训练 100 个周期，批量大小为 16。

潜在向量表示了输入数据中最重要的特征的低维形式 [ref⁴]。因此，在训练完整个自动编码器之后，我们可以使用其编码器部分来获得数据的低维（编码）表示。

ref⁴: 深度学习中的自动编码器简介

encoder_model 是通过连接输入层和编码器部分创建的。然后，我们可以在缩放后的 Wine 数据上调用它的 predict() 方法，以获得 Wine 数据集的低维（编码）表示，这由变量 X_encoded 表示。由于**“潜在向量维度”**被设置为 7，因此编码表示具有 7 个特征，数据的维度已被降低！

通过可视化数据来比较 PCA 和自动编码器模型

高维数据的可视化可以通过降维实现 [ref⁵]。因此，降维对于数据可视化极为有用。

ref⁵: 降维的 11 种不同用途

仅使用两个组件（维度），我们绘制了 PCA 和自动编码器模型返回的 Wine 数据集输出。

import matplotlib.pyplot as plt

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))

# Plot PCA output
ax1.scatter(X_pca[:, 0], X_pca[:, 1], c=y, s=25, cmap='plasma')
ax1.set_title('PCA')
ax1.set_xlabel('Component 1')
ax1.set_ylabel('Component 2')

# Plot autoencoder output
ax2.scatter(X_encoded[:, 0], X_encoded[:, 1], c=y, s=25, cmap='plasma')
ax2.set_title('Autoencoder')
ax2.set_xlabel('Dimension 1')
ax2.set_ylabel('Dimension 2')

plt.savefig("Output.png")

Wine 数据的二维表示：PCA 与自动编码器（图像作者）

如你所见，使用 PCA 获得的 Wine 数据的二维表示显示了三种葡萄酒类别之间的明显分离，但在自动编码器输出中分离效果不够好。

PCA 在小数据集如 Wine 数据集上效果很好。

注意： 由于算法的随机特性以及超参数的值，如隐藏层的数量和隐藏单元、每层使用的激活函数类型、损失函数类型、优化器类型、训练轮数和批量大小，自动编码器的输出（右图）可能会有显著变化！但是，三个类别的分离可能不如 PCA 输出的效果好。

PCA 与自动编码器：哪种在降维方面更好？

选择 PCA 和自动编码器进行降维取决于以下因素。

• 数据集的大小

• 数据集的复杂性（线性还是非线性，图像还是数值数据）

• 分析的目标

• 计算资源的可用性

• 可解释性

PCA 在小数据集上效果很好。它也可以用于较大的数据集。然而，自动编码器在非常大的数据集上效果非常好。

PCA 对线性数据效果很好，因为它是一种线性降维技术。对于非线性数据效果不佳。相反，自动编码器可以轻松捕捉数据中的复杂和非线性模式。因此，它们在非线性数据上表现良好。

自动编码器在图像数据上表现良好。PCA 在数值数据上表现良好。

无法确定自编码器模型潜在向量中每个组件（特征）的重要性。但是在 PCA 中，我们可以为此创建累积解释方差图。

由于自编码器是神经网络，它的架构可能变得复杂。此外，它还需要大量的数据。因此，自编码器比 PCA 需要更多的计算资源。

自编码器是黑箱模型。因此，它们难以解释。我们不知道它们如何从我们提供的数据中选择重要特征。因此，这些模型的解释非常困难。

结论

PCA 和自编码器都可以用于执行降维。PCA 是一个通用的机器学习算法，而自编码器是一种神经网络架构，需大数据集和大量计算资源。

自编码器在处理大规模和非线性数据时表现良好。它们足够强大和灵活，可以捕捉数据中的复杂和非线性模式，但在较小的数据集上可能无法超越通用机器学习算法如 PCA！

如果你真的想在执行降维时考虑类可分性，LDA（线性判别分析）是最佳选择。请阅读下面的完整指南。

LDA 在分类数据集中比 PCA 更有效

今天的文章到此结束。

如果你有任何问题或反馈，请告诉我。

那么 AI 课程怎么样？

加入我的神经网络和深度学习课程，这是 Medium 上的首个课程（作者截图）

支持我作为作者

我希望你喜欢阅读这篇文章。如果你愿意支持我作为作者，请 注册会员 以获得对 Medium 的无限制访问权限。每月只需$5，我将从你的会员费中获得一部分。

[## 通过我的推荐链接加入 Medium - Rukshan Pramoditha

阅读 Rukshan Pramoditha 的每一篇故事（以及 Medium 上成千上万的其他作者的故事）。您的会员费用直接…

rukshanpramoditha.medium.com

加入我的私人邮件列表

不要再错过我的精彩故事。通过 订阅我的邮件列表，你将直接在我发布故事时收到。

非常感谢你的持续支持！下篇文章见。祝大家学习愉快！

葡萄酒数据集信息

• 数据集来源： 你可以从这里下载原始数据集。

• 数据集许可证： 该数据集在CC BY 4.0（知识共享署名 4.0）许可证下提供。

• 引用： Lichman, M. (2013). UCI 机器学习库 [archive.ics.uci.edu/ml]。加州尔湾：加利福尼亚大学信息与计算机科学学院。

Rukshan Pramoditha

2023 年 2 月 16 日

Pearson、Spearman 和 Kendall 相关系数，手动计算

原文：towardsdatascience.com/pearson-spearman-and-kendall-correlation-coefficients-by-hand-d2e1676ca73f

学习如何手动计算 Pearson、Spearman 和 Kendall 相关系数，以评估两个变量之间的关系

Antoine Soetewey

·发表于 Towards Data Science ·15 分钟阅读·2023 年 9 月 5 日

–

图片由 S O C I A L . C U T 提供

介绍

在统计学中，相关性用于评估两个变量之间的关系。

在之前的帖子中，我们展示了如何 在 R 中计算相关性并执行相关性测试。在本帖中，我们展示了如何手动计算 Pearson、Spearman 和 Kendall 相关系数，并在两种不同的场景下（即有平局和无平局）进行计算。

数据

为了展示有平局和无平局的方法，我们考虑了两个不同的数据集，一个带有平局，另一个没有平局。

带有平局

对于带有平局的场景插图，假设我们有以下大小为 5 的样本：

表格由作者提供

图表由作者提供

如我们所见，由于变量 x 中有两个相同的观测值，因此存在一些平局。

无平局

对于不需要处理平局的情况，我们将考虑以下大小为 3 的样本：

表格由作者提供

图表由作者提供

手动计算的相关系数

三种最常见的相关性方法是：1

1. Pearson，用于两个 定量连续 变量之间具有线性关系

2. Spearman，用于两个定量变量，如果关系部分线性，或用于一个定性序数变量和一个定量变量

3. Kendall，通常用于两个定性序数变量

每种方法将在接下来的章节中介绍。

请注意，本帖的目的是展示如何手动计算三种相关系数和在两种不同情况下的操作；我们并不关心验证基本假设。

Pearson

幸运的是，无论是否存在重复值，计算 Pearson 相关系数的过程是相同的，所以我们不区分这两种情况。

有无重复值的情况

Pearson 相关系数，样本情况下用 r 表示，可以如下计算2

其中：

• n 是样本大小

• xi 和 yi 是观测值

• ¯x¯ 和 ¯y¯ 是 x 和 y 的样本均值

• sx 和 sy 是 x 和 y 的样本标准差

我们逐步展示如何计算。

步骤 1.

如你所见，在计算相关系数之前，我们首先需要计算两个变量的均值和标准差。

x 的均值计算如下：

x 的标准差计算如下：

这些公式也可以类似地用于变量 y。

我们从计算两个变量的均值开始：

步骤 2.

我们现在需要计算每个变量的标准差。为了简化计算，最好使用表格，从观测值开始：

作者提供的表格

计算标准差时，我们需要每个观测值与其均值之间平方差的总和，即：

我们首先在表格中创建两列新数据，标记为 x-xbar 和 y-ybar，对应于：

作者提供的表格

我们取这两列新数据的平方，标记为 (x-xbar)² 和 (y-ybar)²：

作者提供的表格

然后我们对这两列求和，结果是：

标准差因此为：

步骤 3.

我们现在需要

因此，我们在表格中添加一列新数据，对应于 xiyi，标记为 x*y：

作者提供的表格

步骤 4。

我们取最后一列的总和，结果是：

步骤 5。

最后，可以通过将上述找到的值代入初始公式来计算皮尔逊相关系数：

请注意，还有其他公式可以计算皮尔逊相关系数。例如，

所有公式当然会给出完全相同的结果。

供您参考，将皮尔逊相关系数 r 平方得到决定系数 R2，适用于简单的线性回归。

斯皮尔曼

现在我们介绍斯皮尔曼相关系数，也称为斯皮尔曼等级相关系数。这个系数实际上与皮尔逊系数相同，只是计算基于等级值而不是原始观察值。

再次介绍如何一步步手动计算，但这次我们区分了两种情况：

1. 如果有并列

2. 如果没有并列

含并列

斯皮尔曼相关系数（含并列），记作 rs，定义如下：

其中：

这里是如何手动计算的。

步骤 1。

如前所述，斯皮尔曼系数基于等级。因此，我们首先需要分别为两个变量添加观察值的等级（从最低到最高）。

对于 y，我们看到：

• -3 是最小值，所以我们将其分配为等级 1。

• -1 是第二小的值，所以我们将其分配为等级 2。

• 然后是 0，所以我们将其分配为等级 3。

• 然后是 1，所以我们将其分配为等级 4。

• 最终，2 是最大值，所以我们将其分配为等级 5。

对于 x 也是如此，只是这里有两个观察值为 5（因此等级中会有并列）。在这种情况下，我们取平均等级：

• -1 是最小值，所以我们将其分配为等级 1。

• 2 是第二小的值，所以我们将其分配为等级 2。

• 然后是 3，所以我们将其分配为等级 3。

• 最终，两个最大值属于等级 4 和 5，所以我们将它们都分配为等级 4.5。

我们在表中包括等级，表示为Rx和Ry：

作者提供的表格

步骤 2。

从那里开始，计算方法类似于皮尔逊系数，只是我们处理的是等级而不是初始观察值。为了避免剩余步骤中的任何混淆，我们从表格中删除初始观察值，仅保留等级：

作者提供的表格

我们从等级的平均值开始：

对于标准差，我们使用与皮尔逊系数相同的表格：

作者提供的表格

我们求出最后两列的总和，这给出：

标准差为：

步骤 3.

我们现在需要

因此，我们在表格中添加了一列，对应于 RxiRyi，用Rx*Ry表示：

按作者分类的表格

步骤 4.

我们求出最后一列的总和，这给出

步骤 5.

最后，通过将所有值代入初始公式来计算斯皮尔曼相关系数：

无平局

当每个变量中的所有初始值都不同，即所有排名都是不同的整数时，就没有平局。在这种特定情况下，斯皮尔曼系数可以用以下简化公式计算：

例如，假设样本中没有平局的数据介绍如下：

按作者分类的表格

步骤 1.

每个变量中的所有观察值都不同，因此所有排名都是不同的整数，并且没有平局：

按作者分类的表格

步骤 2.

我们只需要计算每行两个排名之间的差异，用Rx-Ry表示：

按作者分类的表格

求这些差异的平方，记作(Rx-Ry)²：

按作者分类的表格

然后求出最后一列的总和，这给出

步骤 3.

最后，我们可以填入初始公式以找到斯皮尔曼系数：

Kendall

Kendall 系数相关性，也称为 Kendall 的τ系数，与斯皮尔曼系数类似，只是它通常在样本较小且有许多排名平局时更受欢迎。

在这里我们也区分两种情况：

• 如果没有平局

• 如果有平局

与斯皮尔曼系数不同，我们首先说明没有平局的情况，然后是有平局的情况。

无平局

当没有平局时，Kendall 系数，记作τa，定义如下：

其中：

• C 是一致对的数量

• D 是不一致对的数量

这个系数，也称为 Kendall tau-a，不对平局进行任何调整。

让我们看看使用无平局的数据（与斯皮尔曼相关性无平局的相同数据）的一致对和不一致对：

按作者分类的表格

步骤 1.

我们首先计算每个变量（Rx 和 Ry）的排名，就像我们为斯皮尔曼系数做的那样：

按作者分类的表格

步骤 2.

然后我们在两个变量 Rx 或 Ry 中任意选择一个作为参考变量。假设我们在这里选择 Rx 作为参考变量。

我们按照这个参考变量对数据集进行排序：

表格由作者提供

从现在开始，我们只看第二个变量（非 参考水平的那个，这里是 Ry）的排名，因此为了避免在剩下的步骤中产生混淆，我们只保留 Ry 列：

表格由作者提供

步骤 3.

我们现在逐行检查 Ry，并检查表中它下面的行是否更小或更大。

在我们的表格中，Ry 的第一行是 1。我们看到它下面的值是 3，大于 1。由于 3 大于 1，这被称为和谐配对。我们将其写入表中：

表格由作者提供

下一行，2，也大于 1，因此它也是一个和谐配对。我们也将其写入表中：

表格由作者提供

我们重新从 Ry 的第二行开始，即 3。再次，我们查看下面的行，检查它是否更大或更小。在这里，下面的行是 2，小于 3，因此我们有一个 discordant 配对。我们将此信息添加到表中：

表格由作者提供

步骤 4.

现在我们计算和谐配对和 discordant 配对的总数。总共有：

• 2 对和谐配对，所以 C = 2，和

• 1 对 discordant 配对，所以 D = 1。

步骤 5.

最后，我们将刚刚找到的值代入初始公式中：

或者，我们也可以使用以下公式

其中 n 仍然是样本大小，两者都给出完全相同的结果。

有重复值的情况

我必须承认，有重复值的过程比没有重复值的情况稍微复杂一些。

Kendall tau-b 系数，即对重复值进行调整的系数，定义如下

其中：

请注意，这里 C2n 中的字母 C 表示“组合”，而不是“和谐”。

让我们用皮尔逊相关系数和有重复的斯皮尔曼相关系数的数据集来说明这个场景和公式，即：

表格由作者提供

步骤 1.

与没有重复的情况类似，有重复的情况下，这个过程也是基于排名的，因此我们首先为每个变量添加排名，通常表示为Rx和Ry：

表格由作者提供

步骤 2.

必须选择两个变量中的一个（x 或 y）作为参考变量。这次，我们不是任意选择，而是选择没有重复值的那个。在我们的例子中，变量 x 有重复值，而变量 y 没有任何重复值。因此 y 将是我们的参考变量。

然后我们：

• 按参考变量排序数据集，这里是 Ry，并

• 我们只保留必要的列，以避免在剩余步骤中出现任何混淆：

作者提供的表格

步骤 3.

我们检查它是否是一个一致对或不一致对，与我们处理无联系时的方式相同，但这次我们不计算联系。

列 Rx 中的第一行是 1。我们将表中所有低于该值的行与之比较。表中 1 以下的所有行都更大，因此我们写下它们是一致对：

作者提供的表格

我们对每一行重复相同的过程。

例如，对于列 Rx 的第二行，我们有值 2。再次，表中所有低于该值的行都更大，因此我们写下它们是一致对：

作者提供的表格

对于列 Rx 的第三行，我们有值 4.5。表中紧接其下的行（= 3）小于 4，因此我们写下 D 表示不一致对：

作者提供的表格

现在，如你所见，最后一行也是 4.5。由于它等于我们正在比较的值，它既不是一致对也不是不一致对，因此我们在表格中写下“T”表示“联系”：

作者提供的表格

最后，列 Rx 中的第四行是 3，我们将其与 4.5 进行比较，发现 4.5 更大，因此它是一个一致对：

作者提供的表格

步骤 4.

然后我们汇总一致对和不一致对的数量：

• C = 8

• D = 1

步骤 5.

我们现在拥有计算 τb 分子的所有信息，但仍需找到

以计算分母。

如上所述，C2n 是可能对的总数，因此对应于两个值的组合数量。这个对的数量可以通过组合公式找到，即：

其中 n 是样本大小。

在我们的例子中，我们的样本大小为 5，所以我们有：3

此外，nx 和 ny 分别是 x 和 y 变量中具有联系的可能对的数量。

当查看上述表格中的 Ry 和 Rx 时：

作者提供的表格

我们看到有：

• 变量 x 中的 2 个相同值

• 变量 y 中的 0 个相同值

这意味着具有联系的对的数量等于：4

• 对于变量 x（唯一的可能联系对是 {4.5, 4.5}）

• 0 对于变量 y（由于所有等级都是独特的，因此不存在成对的并列情况）

因此，

现在我们拥有计算 τb 所需的所有信息！

步骤 6。

将上述找到的值代入初始公式中，我们得到

在 R 中验证

为了完整性，我们使用 R 验证每个系数和场景的结果。

Pearson：

x <- c(-1, 3, 5, 5, 2)
y <- c(-3, 1, 0, 2, -1)

cor(x, y, method = "pearson")

## [1] 0.8769051

带并列情况的 Spearman：

cor(x, y, method = "spearman")

## [1] 0.8207827

不带并列情况的 Spearman：

x2 <- c(3, 5, 2)
y2 <- c(5, 1, -1)

cor(x2, y2, method = "spearman")

## [1] 0.5

不带并列情况的 Kendall：

cor(x2, y2, method = "kendall")

## [1] 0.3333333

带并列情况的 Kendall：

cor(x, y, method = "kendall")

## [1] 0.7378648

我们确实手动计算的结果与 R 中得到的结果相同（任何差异都是由于四舍五入造成的）！

结论

记住，相关系数（无论是 Pearson、Spearman 还是 Kendall）的范围是 -1 到 1，其中 0 表示没有相关性，绝对值越接近 1，相关性越强。

广义上讲，正相关意味着一个变量的高值与另一个变量的高值相关（反之亦然）。负相关意味着一个变量的高值与另一个变量的低值相关。

最后但同样重要的是，记住从样本中计算出的相关系数得出的结论不能在没有适当统计检验的情况下推广到总体，即相关性检验。

感谢阅读。

希望这篇文章帮助你手动计算 Pearson、Spearman 和 Kendall 相关系数（包括和不包括并列情况）。

和往常一样，如果你对本文涉及的主题有任何问题或建议，请在评论中添加，这样其他读者也能从讨论中受益。

1. 供参考，Pearson 相关系数被视为参数检验，而 Spearman 和 Kendall 相关系数被视为非参数检验。↩︎

2. 在这里，我们假设我们有一个样本而不是总体。如果你拥有的观察值代表整个总体，则相关系数记作 ρ，公式略有不同。如果需要，请参见样本和总体之间的差异的回顾。↩︎

3. 记住，n 的阶乘，记作 n!，是 n! = n×(n−1)×⋯×1。例如，5! = 5×4×3×2×1=120。↩︎

4. 如果你有许多相同的值，并希望通过公式而非手动计数来计算可能的配对数量，你可以再次使用组合的公式，即 C2n = n! / x!(n−x)!，其中 n 是给定组中相同值的数量。请记住，组合的数量必须对所有组进行求和。在我们的例子中，我们只有一个包含 2 个值的组，因此 nx = C22 = 2! / 2!(2−2)! = 1。↩︎

相关文章

• R 中的相关系数和相关性检验

• R 中的描述性统计

• 什么是生存分析？手动示例和 R 中的示例

• R 中的单样本 Wilcoxon 检验

• 手动假设检验

最初发布于 https://statsandr.com 于 2023 年 9 月 5 日。

人员分析是新的热门话题，了解它的原因在于此

原文：towardsdatascience.com/people-analytics-is-the-new-big-thing-heres-why-you-should-know-about-it-61cda54da5e7

关注工作中的数据与分析

拉希·德赛

·发布于 Towards Data Science ·5 分钟阅读·2023 年 7 月 7 日

–

人员分析可以定义为处理人力资源数据并将其转化为数据驱动的可操作见解，以改善业务运作方式的过程。

José Martín Ramírez Carrasco 在 Unsplash上的照片

人员分析（也称为人才分析或劳动力分析）是一种深度数据驱动、基于证据的数据分析流派，其中专业人员负责使用数据科学概念、实施数据模型和可视化数据，以推荐可操作的见解。简单来说，人员分析是从提取到传达给业务的端到端数据生命周期，但针对的是人力资源数据。

在当今这个时代，人员分析对于一个组织至关重要，因为现在，比以往任何时候都更能影响一个组织的成败。人员分析与业务能够通过内部人员对业务目标、使命和愿景的影响而取得的成就相关，并评估人力资源举措的有效性。

人员分析并不是一个小的、利基的团队；相反，它涵盖了整个组织，并在推动劳动力战略方面发挥着关键作用。

为什么选择人员分析？

超过 70%的公司现在表示他们认为人员分析是一项高优先级任务，但仅有 8%的组织报告说他们拥有可用的数据，而只有 9%的公司相信他们对推动其绩效的人才维度有很好的理解——这为更多公司考虑人员分析提供了广阔的机会！

当今组织面临着劳动力、工作场所和工作世界的根本变化。在这个变化的时代，未来的组织即将到来，人才能够获取和员工体验是业务关键绩效指标的重要贡献因素。

数据与人员之间的相互作用因多种原因显得很有前景，但主要因为——

• 人力资源部门已经拥有用于人员分析的原材料——人力资源跟踪、存储和维护一切数据。

• 很多数据随着时间的推移并不会变化——有关个人的事实如种族、性别和工作历史并不会经常变化。

• 人力资本是企业的一项增值资产，它既可以为你节省成本，也可能使你付出代价。

对于任何擅长分析人力资源数据的人来说，优化人才管理和组织绩效、研究人类心理，并将其与员工体验相关联，人员分析是适合你的！

一天的生活：人员分析师

人员分析师负责理解、分析、执行和解决业务和人力资源问题，并为现有和新员工发现机会。

人员分析师的工作日程与其他数据分析师的生活日常大致相同，不同之处在于所处理的数据——人员分析师处理的人力资源数据，包括但不限于——

• 人才获取

• 离职率（员工离开公司的比率）

• 员工名册

• 员工情绪

• 员工留存和参与度

人员分析师负责领导数据模型、工具和程序的设计、开发和实施，进行统计分析、定量和定性研究，以及预测建模，以推动对人力资源影响和业务洞察的发现和表达。

人员分析作为一种先进的分析支持/资源，帮助领导者头脑风暴并制定有效且高效的人才策略。

成为人员分析师需要什么？

• 基本分析工具包——熟练掌握从数据提取到数据可视化的整个数据生命周期，使用 SQL、Python、R、SAS、Qlik、Tableau、Power BI、Microsoft Excel、Alteryx 等工具，并对数据库架构有良好的理解

• 商业敏锐性

• 良好的沟通和人际交往技能（加上一点幽默感）

• 数据讲故事——创建和进行信息性演示，以推荐和影响决策

• 能够通过外部视角解决复杂问题，跳出框框，创新日常流程

• 能够在最复杂的情况下独立工作和团队协作，且只需最少的指导

• 能够将数据翻译为业务——运行预测模型很简单，但让业务理解其含义才是真正的挑战

• 咨询技能——问题定义、假设构建、项目管理、解决方案开发、变更管理、利益相关者管理等

• 网络建设——从现实世界获取脉搏检查

• 有心理学背景？赚取那笔钱吧！

人员分析的作用是什么？

• 人员概况： 为相关受众提供关键人力资源指标的高层次概览，如员工数量、人才组合、人口统计数据、员工流失、人才获取渠道等

• 流失与离职： 提供关于员工流失、内部流动、角色转移（横向和晋升）、新招聘/高绩效员工流失、主动或被动离职的见解

• 人员统计细分： 提供有关员工的清晰洞察，包括领导者、部门、经理级别、职位名称等，涵盖内部和承包商员工

• 人员利用： 展示现有员工的才能组合，包括内部和外部招聘人员的人口统计信息、晋升和离职情况，并确保公平的薪酬待遇

• 人才获取渠道： 展示候选人在申请过程中的各个阶段的才能组合

• 公平工作环境： 工作于员工保留和参与，解决需要了解驱动员工体验的因素，以创造一个愉快的工作环境

• 改进的招聘机会： 分析候选人在招聘过程中从申请到录用信的时间分布，识别机会和改进领域

• 成本效益： 通过预测谁将辞职来节省资金

• 员工情绪： 通过脉搏调查测量员工情绪，并将假设与调查结果进行对比，识别机会和改进领域

以上就是我对这篇博客的所有内容。感谢阅读！希望你觉得这篇文章很精彩。在评论中告诉我你在讲故事方面的经验，你在数据领域的旅程，以及你对 2023 年的期望！

如果你喜欢阅读这样的故事，可以通过这个链接考虑注册成为 Medium 会员**。**

祝数据探测愉快！

拉希是来自芝加哥的数据高手，她热衷于分析数据并创建数据故事来传达洞察。她是一名全职的医疗数据分析顾问，周末时会喝杯咖啡，写博客谈论数据……

感知机：第一个神经网络模型

原文：towardsdatascience.com/perceptrons-the-first-neural-network-model-8b3ee4513757

概述及在 Python 中的实现

Dr. Roi Yehoshua

·发表于 Towards Data Science ·14 分钟阅读·2023 年 3 月 28 日

–

Hal Gatewood 在 Unsplash 上的照片

感知机是神经网络（NN）的最早计算模型之一，它们构成了今天更复杂、更深层网络的基础。理解感知机模型及其理论将为你理解神经网络中的许多关键概念提供良好的基础。

背景：生物神经网络

生物神经网络（例如我们大脑中的网络）由大量称为神经元的神经细胞组成。

每个神经元通过称为树突的纤维接收来自其邻近神经元的电信号（冲动）。当其接收到的信号总和超过某个阈值时，神经元会通过长纤维称为轴突发出自己的信号，这些轴突连接到其他神经元的树突上。

两个神经元之间的连接称为突触。平均而言，每个神经元连接到大约 7000 个突触，这表明我们大脑网络的高度连通性。当我们学习两个概念之间的新联想时，表示这些概念的神经元之间的突触强度会增强。这种现象被称为赫布规则（1949），其表述为“一起发火的细胞会一起连接”。

生物神经元 (公共图像 在维基共享资源下自由授权)

感知机模型

感知机模型由弗兰克·罗森布拉特于 1957 年提出，是生物神经元的简化模型。

感知机有 m 个二进制输入，表示为 x₁, …, xₘ，这些输入代表来自其邻近神经元的信号，并输出一个二进制值 o，表示感知机是否“激发”。

感知机模型

每个输入神经元 xᵢ 通过一个强度由权重 wᵢ 表示的连接与感知机相连。权重较高的输入对感知机的输出有更大的影响。

感知机首先计算其输入信号的加权和，通过将每个输入乘以其对应的权重。这个加权和通常称为净输入，表示为 z：

感知机的净输入

如果净输入超过某个预定义的阈值 θ，则感知机激发（输出为 1），否则不激发（输出为 0）。换句话说，感知机当且仅当：

我们的目标是找到权重 w₁, …, wₘ 和阈值 θ，使得感知机能正确地将其输入 x₁, …, xₘ（代表数据中的特征）映射到期望的输出 y（代表标签）。

为了简化学习过程，我们添加了一个特殊的输入神经元，称为偏置神经元，它总是输出值 1。这个神经元通常表示为x₀，其连接权重表示为b或w₀。

因此，感知机的净输入变为：

包括偏置的净输入

这种公式允许我们像学习输入信号的权重一样学习正确的阈值（偏置）。

在向量形式中，我们可以将 z 写作输入向量 x = (x₁, …, xₘ)ᵗ 与权重向量 w = (w₁, …, wₘ)ᵗ 的点积加上偏置：

净输入的向量形式

只有当净输入非负时，感知机才会激发，即：

更一般地，感知机对净输入应用激活函数 f(z)，以生成其输出。感知机中最常用的两个激活函数是：

1. 阶跃函数（也称为海维赛德函数）是一个对于负输入值为 0，对于非负输入值为 1 的函数：

阶跃函数

2. 符号函数是一个对于负输入值为-1，对于非负输入值为 1 的函数：

符号函数

在更复杂的网络中使用其他类型的激活函数，例如多层感知机（MLPs）。在本文的其余部分，我将假设感知机使用步进函数。

总结一下，感知机的计算包括两个步骤：

1. 输入值 x₁, …, xₘ 乘以其对应的权重 w₁, …, wₘ，然后加上偏置 b，得到感知机的净输入 z = wᵗx + b**。**

2. 在生成二进制输出（0/1 或 -1/+1）的净输入上应用激活函数 f(z)。

我们可以将整个计算写成一个方程：

其中 f 是选择的激活函数，o 是感知机的输出。

使用感知机实现逻辑门

为了演示感知机的工作原理，让我们尝试构建计算逻辑函数“与”和“或”的感知机。

作为提醒，逻辑“与”函数有两个二进制输入，当两个输入都为真时返回真（1），否则返回假（0）。

“与”函数的真值表

实现“与”函数的感知机有两个二进制输入和一个偏置。我们希望这个感知机只有在两个输入都“激活”时才会“触发”。例如，可以通过为两个输入选择相同的权重来实现这一点，例如 w₁ = w₂ = 1，然后选择偏置在-2, -1)范围内。这样可以确保当两个神经元都激活时，净输入 2 + b 将是非负的，但当只有一个神经元激活时，净输入 1 + b 将是负的（当没有神经元激活时，净输入 b 也为负）。

计算逻辑“或”函数的感知机

验证你是否理解这个感知机是如何工作的！

作为练习，尝试构建一个用于 NAND 函数的感知机，其真值表如下所示：

NAND 函数的真值表

感知机作为线性分类器

感知机是一种线性分类器，因为它将输入空间划分为由下述超平面分隔的两个区域

分隔超平面的方程

权重向量w与这个超平面正交，从而决定其方向，而偏置 b 定义了它离原点的距离。

每个超平面上方的例子（wᵗx + b> 0**））被感知机分类为正例，而每个超平面下方的例子（wᵗx + b< 0**））被分类为负例。

作为线性分类器的感知器

其他线性分类器包括逻辑回归和线性支持向量机（SVM）。

线性分类器只能学习线性可分的问题，即正例和负例之间的决策边界是一个线性表面（一个超平面）。

例如，以下数据集是非线性可分的，因此感知器无法正确分类此数据集中的所有示例：

非线性可分数据集

感知器学习规则

感知器具有一个简单的学习规则，如果数据是线性可分的，它可以找到分离超平面。

对于每一个被感知器误分类的训练样本（即oᵢ ≠ yᵢ），我们将以下更新规则应用于权重向量：

感知器学习规则

其中α是学习率（0 < α ≤ 1），它控制每次更新中权重调整的大小。

换句话说，我们将每个连接权重wⱼ的感知器在这个例子上的误差（真实标签yᵢ和输出oᵢ之间的差异）乘以相应输入xᵢⱼ的值和学习率。

这个学习规则尝试做的是减少感知器输出oᵢ和真实标签yᵢ之间的差距。为了理解它为什么有效，让我们检查感知器误分类的两种可能情况：

1. 真实标签是yᵢ = 1，但感知器的预测是oᵢ =0，即wᵗxᵢ + b < 0。在这种情况下，我们希望增加感知器的净输入，使其最终变为正值。

   为此，我们将量(yᵢ — oᵢ)xᵢ = xᵢ添加到权重向量中（乘以学习率）。这会增加正值输入的权重（其中xᵢⱼ > 0），同时减少负值输入的权重（其中xᵢⱼ < 0）。因此，感知器的整体净输入会增加。

2. 真实标签是yᵢ = 0，但感知器的预测是oᵢ = 1，即wᵗxᵢ + b > 0。类似于前面的情况，这里我们希望减少感知器的净输入，使其最终变为负值。这是通过将量(yᵢ — oᵢ)xᵢ = -xᵢ添加到权重向量中实现的，因为这会减少正值输入的权重（其中xᵢⱼ > 0），而增加负值输入的权重（其中xᵢⱼ < 0）。因此，感知器的整体净输入会减少。

该学习规则会依次应用于所有训练样本（顺序可以是任意的）。通常需要对整个训练集进行多次迭代（称为一个周期），才能找到正确的权重向量（即分隔正负样本的超平面的向量）。

根据感知器收敛定理，如果数据是线性可分的，则重复应用感知器学习规则最终将收敛于分隔超平面的权重（在有限的步骤中）。有兴趣的读者可以在这篇论文中找到该定理的正式证明。

感知器学习算法

实际上，感知器学习过程可能需要较长时间才能收敛（即，在训练集上达到零错误）。此外，数据本身可能不是线性可分的，在这种情况下，算法可能永远不会终止。因此，我们需要通过一些预定义的参数来限制训练周期的数量。如果感知器在达到该数量之前在训练集上达到了零错误，我们可以提前停止训练。

感知器学习算法的总结如下伪代码：

请注意，权重通常初始化为小的随机值以打破对称性（如果所有权重相等，则感知器的输出对于每个输入都是恒定的），而偏置则初始化为零。

示例：学习多数函数

例如，我们来看看感知器学习算法如何用于学习三个二进制输入的多数函数。多数函数是一个在其输入中有一半或更多为真时返回真（1），否则返回假（0）的函数。

感知器的训练集包括所有 8 种可能的二进制输入：

多数函数的训练集

在这个例子中，我们将假设初始权重和偏置为 0，学习率为 α = 0.5。

让我们跟踪第一次训练周期中的权重更新。呈现给感知器的第一个样本是x = (0, 0, 0)ᵗ。在这种情况下，感知器的净输入为：z = wᵗx + b = 0 × 0 + 0 × 0 + 0 × 0 + 0 = 0。因此，它的输出是 o = 1（记住，当输入≥0 时，阶跃函数的输出为 1）。然而，在这种情况下，目标标签是 y = 0，所以感知器的错误为 y - o = -1。

根据感知器学习规则，我们通过将 α(y - o)xᵢ = -0.5xᵢ 加到每个权重 wᵢ 上来更新权重。由于在这种情况下所有输入都是 0，除了偏置神经元（x₀ = 1），我们只更新偏置为-0.5，而不是 0。

我们对其他 7 个训练样本重复相同的过程。下表显示了每个样本之后的权重更新：

第一轮训练

在第一轮中，感知器发生了 4 次错误。第一轮后的权重向量为w = (0, 0.5, 1)ᵗ，偏置为 0。

在第二轮训练中，我们得到以下权重更新：

第二轮训练

这次感知器只发生了三次错误。第二轮后的权重向量为w = (0.5, 0.5, 1)ᵗ，偏置为-0.5。

第三轮的权重更新如下：

第三轮训练

在这一轮中，第二个示例更新后，感知器已收敛到解决此分类问题的权重向量：

w = (0.5, 0.5, 0.5)ᵗ 和 b = -1。由于所有权重相等，感知器仅在至少两个输入为 1 时激发，此时它们的加权和至少为 1，即大于或等于偏置的绝对值，因此感知器的净输入是非负的。

Python 中的感知器实现

现在我们来实现 Python 中的感知器学习算法。

我们将通过扩展sklearn.base.BaseEstimator类来实现它作为自定义 Scikit-Learn 估计器。这将允许我们像使用其他 Scikit-Learn 估计器一样使用它（例如，将其添加到一个pipeline中）。

自定义估计器需要实现**fit()和predict()方法，并在init()**方法中设置所有超参数。

我会首先展示这个类的完整代码，然后逐步讲解。

from sklearn.base import BaseEstimator

class Perceptron(BaseEstimator):
    def __init__(self, alpha, n_epochs):
        self.alpha = alpha        # the learning rate
        self.n_epochs = n_epochs  # number of training iterations

    def fit(self, X, y):
        (n, m) = X.shape  # n is the number of samples, m is the number of features

        # Initialize the weights to small random values
        self.w = np.random.randn(m)
        self.b = 0

        # The training loop
        for epoch in range(self.n_epochs):
            n_errors = 0  # number of misclassification errors

            for i in range(n):
                o = self.predict(X[i])
                if o != y[i]:
                    # Apply the perceptron learning rule
                    self.w += self.alpha * (y[i] - o) * X[i]
                    self.b += self.alpha * (y[i] - o)
                    n_errors += 1

            # Compute the accuracy on the training set
            accuracy = 1 - (n_errors / n)
            print(f'Epoch {epoch + 1}: accuracy = {accuracy:.3f}')

            # Stop the training when there are no more errors
            if n_errors == 0:
                break

    def predict(self, X):
        z = X @ self.w + self.b
        return np.heaviside(z, 1)  # the step function

类的构造函数初始化了模型的两个超参数：学习率（alpha）和训练轮数（n_epochs）。

def __init__(self, alpha, n_epochs):
    self.alpha = alpha
    self.n_epochs = n_epochs

**fit()**方法在给定数据集X（带标签y）上运行学习算法。我们首先通过查询X的形状来找出数据集中有多少样本和特征：

(n, m) = X.shape

n 是训练样本的数量，m 是特征的数量。

接下来，我们使用标准正态分布（均值为 0，标准差为 1）初始化权重向量，将偏置设为 0：

self.w = np.random.randn(m)
self.b = 0

我们现在运行训练循环进行n_epochs次迭代。在每次迭代中，我们遍历所有训练样本，对于每个样本，我们通过调用**predict()**方法并将其输出与真实标签进行比较，检查感知器是否正确分类：

for i in range(n):
    o = self.predict(X[i])
    if o != y[i]:

如果感知器错误地分类了样本，我们将权重更新规则应用于权重向量和偏置，然后将误分类错误的数量增加 1：

self.w += self.alpha * (y[i] - o) * X[i]
self.b += self.alpha * (y[i] - o)
n_errors += 1

当轮次终止时，我们报告感知器在训练集上的当前准确性，如果错误数量为 0，我们终止训练循环：

accuracy = 1 - (n_errors / n)
print(f'Epoch {epoch + 1}: accuracy = {accuracy:.3f}')

if n_errors == 0:
    break

predict() 方法非常直接。我们首先计算感知机的净输入，即输入向量与权重的点积加上偏置：

z = X @ self.w + self.b

最后，我们使用 NumPy 的 heaviside() 函数将步进函数应用于净输入并返回输出：

return np.heaviside(z, 1)

np.heaviside() 的第二个参数指定了函数在 z = 0 时的值。

现在让我们在由 Scikit-Learn 的 make_blobs() 函数生成的数据集上测试我们的实现。

我们首先生成一个包含 100 个随机点的数据集，并将其分为两组：

from sklearn.datasets import make_blobs

X, y = make_blobs(n_samples=100, n_features=2, centers=2, cluster_std=0.5)

我们将 cluster_std 设置为 0.5（而不是默认的 1），以确保数据是线性可分的。

让我们绘制数据集：

import seaborn as sns

sns.scatterplot(x=X[:, 0], y=X[:, 1], hue=y, style=y, markers=('s', 'o'), 
                palette=('r', 'b'), edgecolor='black')
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')

blobs 数据集

现在我们实例化 Perceptron 类并将其拟合到数据集上：

perceptron = Perceptron(alpha=0.01, n_epochs=10)
perceptron.fit(X, y)

训练过程中的输出是：

Epoch 1: accuracy = 0.250
Epoch 2: accuracy = 0.950
Epoch 3: accuracy = 1.000

感知机在经过三次训练周期后已经收敛。

我们可以使用以下函数绘制感知机找到的决策边界和两个类别区域：

def plot_decision_boundary(model, X, y):
    # Retrieve the model parameters
    w1, w2, b = model.w[0], model.w[1], model.b

    # Calculate the intercept and slope of the separating line
    slope = -w1 / w2
    intercept = -b / w2

    # Plot the line
    x1 = X[:, 0]
    x2 = X[:, 1]
    x1_min, x1_max = x1.min() - 0.2, x1.max() + 0.2
    x2_min, x2_max = x2.min() - 0.5, x2.max() + 0.5
    x1_d = np.array([x1_min, x1_max])
    x2_d = slope * x1_d + intercept

    # Fill the two classification areas with two different colors
    plt.plot(x1_d, x2_d, 'k', ls='--')
    plt.fill_between(x1_d, x2_d, x2_min, color='blue', alpha=0.25)
    plt.fill_between(x1_d, x2_d, x2_max, color='red', alpha=0.25)
    plt.xlim(x1_min, x1_max)
    plt.ylim(x2_min, x2_max)

    # Draw the data points
    sns.scatterplot(x=x1, y=x2, hue=y, style=y, markers=('s', 'o'), 
                    palette=('r', 'b'), edgecolor='black')
    plt.xlabel('$x_1$')
    plt.ylabel('$x_2$')

plot_decision_boundary(perceptron, X, y)

感知机找到的分离超平面

Scikit-Learn 提供了自己的 Perceptron 类，实现了类似的算法，但提供了更多选项，如正则化和提前停止。

感知机模型的局限性

尽管感知机模型显示出一些初步的成功，但很快就意识到感知机无法学习一些简单的函数，如 XOR 函数：

XOR 问题无法通过感知机解决

XOR 问题不是线性可分的，因此线性模型如感知机无法解决它。

这一发现使得神经网络领域停滞了很多年（这一时期被称为“AI 冬天”），直到意识到通过将多个感知机堆叠在层中可以解决更复杂和非线性的问题，例如 XOR 问题。

多层感知机（MLPs）在 这篇文章 中有详细介绍。

最后的备注

除非另有说明，所有图片均为作者提供。

你可以在我的 GitHub 上找到这篇文章的代码示例：github.com/roiyeho/medium/tree/main/perceptrons

感谢阅读！