TowardsDataScience 2023 博客中文翻译（三百四十二）-CSDN博客

原文：TowardsDataScience

协议：CC BY-NC-SA 4.0

了解 LoRA — 低秩适配用于微调大型模型

原文：towardsdatascience.com/understanding-lora-low-rank-adaptation-for-finetuning-large-models-936bce1a07c6?source=collection_archive---------1-----------------------#2023-12-22

这一参数高效微调方法背后的数学原理

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传 Bhavin Jawade

关注发布于 Towards Data Science · 4 分钟阅读 · 2023 年 12 月 22 日

–

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

来源 — 图片由 DALLE-3 生成。提示：一个较小的机器人与一个较大的机器人握手。较小的机器人是紫色的，充满电，活跃而有活力，而较大的机器人则被冻结在冰中，呈灰色。

微调大型预训练模型在计算上具有挑战性，通常涉及调整数百万个参数。虽然这种传统的微调方法有效，但需要大量的计算资源和时间，这为将这些模型适应特定任务带来了瓶颈。LoRA 通过在微调过程中分解更新矩阵，提出了一个有效的解决方案。为了研究 LoRA，我们首先回顾一下传统的微调。

（ Δ W ）的分解

在传统的微调中，我们会修改预训练神经网络的权重以适应新任务。这种调整涉及改变网络的原始权重矩阵（ W ）。微调过程中对（ W ）所做的更改统称为（ Δ W ），以便更新后的权重可以表示为（ W + Δ W ）。

现在，LoRA 方法并不是直接修改（ W ），而是寻求分解（ Δ W ）。这种分解是减少与微调大模型相关的计算开销的关键步骤。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

传统的微调可以根据上述方法重新构想。在这里，W 是固定的，而 ΔW 是可训练的（图像来源于博客作者）

内在秩假设

内在秩假设表明，通过较低维度的表示可以捕捉到神经网络的显著变化。基本上，它认为（ Δ W ）的所有元素并不是同样重要的；相反，这些变化的一个较小子集可以有效地封装所需的调整。

引入矩阵（ A ）和（ B ）

在此假设基础上，LoRA 建议将（ Δ W ）表示为两个较小矩阵（ A ）和（ B ）的乘积，且秩较低。更新后的权重矩阵（ W’ ）因此变为：

[ W’ = W + BA ]

在这个方程中，（ W ）保持不变（即，在训练过程中不更新）。矩阵（ B ）和（ A ）的维度较低，其乘积（ BA ）表示（ Δ W ）的低秩近似。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

ΔW 被分解成两个矩阵 A 和 B，其中两个矩阵的维度都低于 d x d。（图像来源于博客作者）

低秩对可训练参数的影响

通过选择具有较低秩（ r ）的矩阵（ A ）和（ B ），可训练参数的数量显著减少。例如，如果（ W ）是一个（ d x d ）矩阵，传统上，更新（ W ）将涉及（ d² ）个参数。然而，使用大小为（ d x r ）和（ r x d ）的（ B ）和（ A ），总参数数量减少到（ 2dr ），当（ r << d ）时，这个数量要小得多。

通过低秩自适应（LoRA）方法实现的可训练参数数量的减少提供了几个重要的好处，特别是在微调大规模神经网络时：

减少内存占用：LoRA 通过降低需要更新的参数数量来减少内存需求，有助于管理大规模模型。
更快的训练和适应：通过简化计算需求，LoRA 加快了大型模型对新任务的训练和微调。
对较小硬件的可行性：LoRA 的较低参数数量使得可以在较不强大的硬件上，如普通 GPU 或 CPU，进行大型模型的微调。
扩展到更大的模型：LoRA 使得 AI 模型的扩展无需相应增加计算资源，使得管理不断增长的模型规模变得更为实际。

在 LoRA 的背景下，秩的概念在确定适应过程的效率和有效性方面起着关键作用。值得注意的是，论文指出矩阵A和B的秩可以非常低，有时低至 1。

尽管 LoRA 论文主要展示了在自然语言处理（NLP）领域的实验，但低秩适应的基本方法具有广泛的适用性，可能在训练不同领域的各种神经网络时有效。

结论

LoRA 通过将（Δ W）分解为低秩矩阵的乘积，有效地平衡了将大型预训练模型适应新任务与保持计算效率的需求。内在的秩概念对这一平衡至关重要，确保了模型学习能力的本质在参数大大减少的情况下得以保留。

参考文献：

[1] Hu, Edward J., et al. “Lora: Low-rank adaptation of large language models.” arXiv 预印本 arXiv:2106.09685 (2021).

理解马赛克图

原文：towardsdatascience.com/understanding-mosaic-plots-fcf148315f4b

PYTHON | DATA | VISUALISATION

一本全面的指南，讲述如何使用 statsmodels 和 Matplotlib 有效绘制多变量数据集

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传 David Farrugia

·发表于Towards Data Science ·7 分钟阅读·2023 年 6 月 13 日

–

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

照片来源：Dimitry B 于Unsplash

我们生活在一个数据泛滥的世界中——一个不断扩展的数字海洋。但在那片海洋中，有待发现的珍贵洞察就在其中。

寻找这些珍珠的关键是什么？数据可视化——将原始数据呈现成更易于理解和解读的视觉形式的过程。

借助数据可视化，你可以赋予原始数字生命，将它们转化为揭示隐藏模式、潜在趋势和关键联系的形式，这些联系可能在数据中隐匿不显。

在我们拥有的数据可视化工具中，有著名的 Matplotlib。

这个强大的 Python 库功能多样且强大。

在 Matplotlib 的工具包中隐藏着一个你可能未曾遇到过的宝石——马赛克图。

这些图表提供了一种强大的方法来可视化跨多个维度的分类数据。

想象一下能够俯瞰你的数据，理解不同变量如何相互作用和交集。马赛克图正是可以做到这一点，以直接且视觉上吸引人的方式呈现复杂关系。

无论是发现客户细分、理解用户行为，还是揭示人口统计趋势，马赛克图都是研究人员和数据科学家必备的绝佳工具。

在本文中，我们将探讨马赛克图的世界。我们将详细介绍如何在 Matplotlib 中创建马赛克图，并讨论如何解读它们，为你的数据可视化工具包增添一个额外的优势。

理解马赛克图

让我们首先将马赛克图比作某种织锦。每个复杂的部分代表数据中的一个独特类别，部分的大小反映了该类别的频率。

因此，织锦函数作为马赛克图，提供了一个视觉表现，帮助理解各种类别变量之间的互动。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

图片由 Tom Rumble 提供，来源于 Unsplash

马赛克图的独特之处在于它们能够同时处理多个维度。想象它就像是审视一个色彩斑斓的数据魔方。

从外部视角，可以观察到不同类别（或颜色，根据我们魔方的比喻）如何相互连接和融合。

想想，你可能在哪里使用这些图表？

想象一下你在 Netflix 等流媒体平台工作，任务是理解不同年龄组、性别和各种类型之间的相互作用。在这里，马赛克图大放异彩。它提供了一个视觉上的描绘，例如，18-25 岁的女性与喜剧的互动频率，与同龄男性进行对比，或者 35-45 岁年龄组对惊悚片的倾向。

马赛克图具有许多优点。

它们擅长处理多维度的类别数据，提供了数据的全景视角。
他们直观的颜色编码增强了视觉吸引力，加快了解读速度并促进快速决策。
它们揭示了原始数据表中可能隐藏的模式和关系。

尽管如此，了解它们的限制仍然至关重要。

如果类别数量过多或类别分布较为均匀，马赛克图可能会迅速变得复杂。
它们在处理定量数据时表现出限制，当精确的数值比较至关重要时，可能不是理想选择。

尽管有这些限制，马赛克图仍然是数据可视化工具中不可或缺的工具。

它们能够将枯燥的数字表转换为动态的趋势和关系表现。

因此，当未来面对多变量类别数据集时，可以考虑使用马赛克图来挖掘隐藏的洞察。

使用 Matplotlib 创建你的第一个马赛克图

初次接触马赛克图可能会感觉像迷宫一样，对吧？实际上，使用 Matplotlib，整个过程变得相对简单。

首先，我们需要一些数据。我们将使用开放的 Titanic 数据集作为示例（参见 www.openml.org/search?type=data&sort=runs&id=40945&status=active）。

我们的目标？弄清楚男性和女性之间的生存率如何分布。

让我们开始吧，好吗？

我们首先导入所需的库：

import seaborn as sns
import matplotlib.pyplot as plt
import pandas as pd
from scipy.io.arff import loadarff
from statsmodels.graphics.mosaicplot import mosaic

接下来我们加载数据集。数据集采用 ARFF 格式。我们可以通过以下方式将其加载到 pandas 中：

# Load the data
raw_data = loadarff('titanic.arff')
titanic = pd.DataFrame(raw_data[0])

我们调用 statsmodels 中的马赛克函数，并传入我们想要可视化的特征。

# Create the mosaic plot
mosaic(titanic, ['alive', 'sex'], gap=0.02)

# Customize the plot
plt.title('Survival by Sex on Titanic')
plt.xlabel('Survived')
plt.ylabel('Sex')

plt.show()

然后，我们得到如下结果：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

当然，这只是一个非常简单的示例。

让我们提升一下难度，好吗？

假设我们想创建一个马赛克图，表示泰坦尼克号上每个乘客舱的比例、他们的生存状态和性别。

我们还将为‘幸存’变量使用颜色编码，并添加一些标签到图表中。

这一次，我们还需要导入 numpy 库。

import numpy as np

我们将如下创建survived列：

# Convert 'survived' into string for clarity in the plot
titanic['survived'] = titanic['survived'].map({0:'Died', 1:'Survived'})

我们还定义了图表使用的属性。这些属性主要用于图表的颜色编码。

# Define properties function for colors
props = lambda key: {'color': '#1b9e77' if 'Survived' in key else '#d95f02'}

然后我们创建如下图表：

# Create a new figure with a defined size
fig, ax = plt.subplots(figsize=(10, 8))

# Create the mosaic plot
mosaic(titanic, ['class', 'sex', 'survived'], gap=0.02, properties=props, ax=ax)

# Customize the plot
plt.title('Survival by Class and Sex on Titanic', fontsize=15)
plt.ylabel('Class - Sex')
plt.xlabel('Proportion')

plt.show()

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

解析马赛克图

马赛克图分为三个主要的垂直部分，代表泰坦尼克号上的三个舱位：第一舱、第二舱和第三舱。

每个部分的宽度与每个舱位的乘客数量成比例。在这种情况下，似乎第三舱的乘客最多，其次是第一舱，最后是第二舱。

在每个舱位部分内，图表进一步水平划分为两个部分，分别代表男性和女性。这些部分的高度与每个舱位中男性和女性的数量成比例。例如，在第三舱中，男性比女性多。

每个性别部分再进一步划分为幸存和遇难两个部分。这些划分的宽度与每个舱位中每个性别的幸存者和遇难者数量成比例。

例如，在‘第一’舱的‘女性’部分，代表‘幸存’的绿色区域远大于代表‘遇难’的橙色区域。这表明第一舱的女性中有很大比例幸存。

颜色编码提供了额外的视觉提示：绿色表示幸存者，而橙色表示遇难者。

使用马赛克图的小贴士和技巧

保持简单 —— 马赛克图很容易变得令人困惑和杂乱。一次只处理少数类别。目标是传达洞察。
实践有效的颜色编码 —— 与任何数据可视化一样，颜色调色板有着显著的影响。马赛克图也不例外。
使用标签和注释 —— 始终标记图表的轴并提供描述性标题。背景信息是关键。

总结 remarks

就这样！

我们已经完成了马赛克图的探索。现在是时候思考一下我们在过程中学到了什么。

我们通过了解数据可视化为何如此重要来开启了这段旅程。我们了解了马赛克图作为观察多层次分类数据的有效方式，它帮助我们发现类别之间的详细相互作用。

然后我们利用 Titanic 数据集制作了我们的第一个马赛克图。我们还为图表添加了颜色、标签和标题，帮助我们用数据编织故事。

马赛克图的真正优势在于其将复杂数据转化为简单、直观故事的能力。但别忘了，保持简单很重要。过多的类别或混乱的颜色方案会将一个出色的图表变成视觉上的头痛。

现在，掌握了这项新知识，你已经准备好开始自己的数据可视化冒险了。也许你会使用马赛克图来了解你在网上商店的客户行为，或分析本地选举中的投票趋势。不论你处理的数据是什么，马赛克图都是你的空白画布，而你的发现就是你用来绘制它的颜色。

不要害羞，多多尝试，尝试不同的调整和各种数据集。实践出真知！

记住，每组数据都有一个故事要讲。所以，出去吧，让你的马赛克图编织出隐藏在数据中的故事。

祝你绘图愉快！

你喜欢这篇文章吗？每月$5，你可以成为会员，解锁对 Medium 的无限访问权限。你将直接支持我和你在 Medium 上的其他喜爱作者。对此，深表感谢！

[## 使用我的推荐链接加入 Medium - David Farrugia

获取对我所有⚡高级⚡内容以及 Medium 上的所有内容的独家访问权限。通过给我买杯咖啡来支持我的工作…

david-farrugia.medium.com](https://david-farrugia.medium.com/membership?source=post_page-----fcf148315f4b--------------------------------)

想要联系我吗？

我很想听听你对这个话题的看法，或者任何关于人工智能和数据的想法。

如果你希望联系我，可以发邮件到davidfarrugia53@gmail.com。

使用 Python 理解多项分布

原文：towardsdatascience.com/understanding-multinomial-distribution-using-python-f48c89e1e29f

多项分布背后的数学和直觉

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传 Reza Bagheri

·发布于Towards Data Science ·20 分钟阅读·2023 年 1 月 6 日

–

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

来源：pixabay.com/vectors/dice-game-die-luck-random-numbers-151867/

多项分布是二项分布的推广，用于在具有多个结果的实验中计算概率。本文对多项分布进行了直观的介绍，并讨论了其数学性质。此外，还将教你如何使用 Python 中的 SciPy 库来建模和可视化多项分布。

二项分布

由于多项分布是二项分布的推广，我们将在此简要回顾一下。有关单变量概率分布的详细讨论，请参阅另一篇文章，如果你对二项分布或随机变量和概率质量函数（PMF）等概念不熟悉，建议你先阅读那篇文章。

随机变量 是一个其值由随机实验的结果决定的变量。随机变量通常用大写字母表示，但我们使用小写字母表示它可以取的特定值。例如，我们可以定义一个随机变量 X 来表示掷硬币的结果。为此，我们需要给每个结果分配一个数值。我们可以用 1 表示得到正面，用 0 表示得到反面。现在，X=1 表示掷硬币的结果是正面，X=0 表示是反面。这样只能取特定值的随机变量称为离散随机变量。

我们将离散随机变量X的概率质量函数（PMF）定义为一个给出X等于某个特定值的概率的函数。从数学上讲，X的 PMF 定义为函数pₓ，使得

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

特定值X可以取的值用x表示，因此P(X=x)表示X=x的概率。我们知道，X所有可能值的概率之和应该等于 1，因此：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

具有伯努利分布的离散随机变量X，其参数为p，表示一个只有两个结果的随机实验，分别用 0 和 1 表示（我们用X~Bern(p)表示）。这里X=1（也称为‘成功’）发生的概率是p，而X=0（也称为‘失败’）发生的概率是 1-p。例如，我们可以用伯努利分布的随机变量X表示掷硬币的结果。现在我们可以假设X=1 和X=0 分别表示得到正面和反面，p是得到正面的概率。

假设我们有一个n次独立随机实验的序列，每次实验可以用一个伯努利分布参数p的随机变量表示。因此，每次实验只有两个结果，分别为 1（成功）和 0（失败），且X=1 发生的概率为p（这样的实验称为伯努利试验）。具有参数n和p的二项分布的离散随机变量X表示这一序列中的成功总数，我们可以写作X ~ Bin(n, p)。假设我们有一枚硬币，正面朝上的概率为p。我们可以用具有参数n和p的二项分布的随机变量X表示这一硬币在n次投掷中的正面总数（图 1）。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

图 1（作者提供的图片）

随机向量

多项分布是一种多变量分布。在统计学中，单变量分布是仅一个随机变量的概率分布。多变量分布是单变量分布的推广，适用于两个或更多随机变量。要理解这些分布，我们首先应讨论随机向量。随机向量是一个随机变量的向量。如果我们有n个随机变量X₁、X₂ …、Xₙ，我们可以将它们放在随机向量X中：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

我们用粗体大写字母表示随机向量。要表示随机向量的一个可能值（这也是一个向量），我们使用粗体小写字母，如x。例如：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

是一个包含两个元素的随机向量，且

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

是它的一个可能值。

联合 PMF

假设我们有一个随机向量

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

其中 X₁，X₂，…，Xₙ 是离散随机变量。X的联合概率质量函数（联合 PMF）定义为：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

其中 P(X₁=x₁，X₂=x₂，…，Xₙ=xₙ) 是 X₁=x₁，X₂=x₂，…，和 Xₙ=xₙ 同时发生的概率。如果我们知道随机向量 X 的联合 PMF，我们可以通过边缘化来推导其一个分量 Xᵢ 的分布。边缘 概率质量函数 可以从 X 的联合 PMF 导出，如下所示

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

这里 Rₓ 是 X 的支持集，表示随机向量 X 可以取的所有值的集合。因此，为了推导 Xᵢ 在点 x 的边缘概率 PMF，我们需要对 Rₓ 中 Xᵢ 等于 x 的所有向量的概率进行求和。

多项分布

多项分布是二项分布的一种推广。假设我们有 n 次独立试验。每次试验有 k 种不同的结果（k ≥ 2），第 i 种结果的概率是 pᵢ 和

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

向量 p 表示这些概率：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

设离散随机变量 Xᵢ 表示在n次试验中结果编号i出现的次数。随机向量 X 定义为

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

并且我们有 x₁+x₂ +…+xₖ=n。那么 X 就称为具有参数 n 和 p 的多项分布。X 的联合 PMF 定义如下：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

让我们看看是否可以证明这一点。首先，我们应该注意到如果 x₁+x₂ +…+xₖ≠n，那么事件 X₁=x₁，X₂=x₂，…，和 Xₖ=xₖ 是不可能的，因此其概率应为零。如果 x₁+x₂ +…+xₖ=n，事件 X₁=x₁，X₂=x₂，…，和 Xₖ=xₖ 的概率是

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

此外，获得 X₁=x₁，X₂=x₂，…，和 Xₖ=xₖ 的 n 次试验的不同方式的总数是

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

所以在n次试验中，X₁=x₁，X₂=x₂，…，和Xₖ=xₖ的总概率是

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

基于这个故事，多项分布可以用来描述一个 k 面的骰子。假设我们有一个 k 面的骰子，且获得面 i 的概率是 pᵢ。此外，设 Xᵢ 表示观察到面 i 的总次数。现在，如果我们掷骰子 n 次，随机向量 X 就具有参数 n 和 p 的多项分布（见图 2）。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

图 2（作者提供的图片）

我们也可以用另一种方式来表达多项式分布。假设我们有一个由 k 种不同类别的物品组成的总体（k ≥ 2）。总体中属于类别 i 的物品比例为 pᵢ，并且

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

现在我们从总体中随机选择 n 个物品（允许重复），我们假设随机变量 Xᵢ 代表选择的物品中属于类别 i 的数量（图 3）。由于这些物品是随机选择的且允许重复，因此选择是相互独立的，每个物品属于类别 i 的概率是 pᵢ。因此，我们可以将每次选择视为一个具有 k 种不同结果的独立试验，每个结果 i 的相应概率是 pᵢ。现在，如果我们定义向量 p 为

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

那么随机向量

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

具有参数 n 和 p 的多项式分布。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

图 3（图像来源：作者）

具有参数 n 和 p 的多项式分布的随机向量 X 可以表示为 n 个具有参数 1 和 p 的多项式分布的随机向量之和：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

证明（可选）：向量 Yᵢ 可以写成如下形式：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

由于 X₁+X₂ +…+Xₖ*=n，我们可以得出在右侧有 n 个向量。每个 Yᵢ 中有一个元素为 1，其余为 0。因此，每个 Yᵢ 可以表示一个具有 k 种不同结果的试验，其中等于 1 的元素表示观察到的结果。因此，每个 Yᵢ 具有参数为 1 和 p 的多项式分布。

我们可以合并多项式随机向量 X 中的多个元素来获得新的多项式随机向量。例如，设

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

现在我们可以合并前三个元素得到新的随机向量

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

而这个随机向量具有参数为 n 和的多项式分布

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

实际上，通过合并 X₁、X₂ 和 X₅，我们创建了一个新类别，该类别的比例是 p₁+p₂+p₅。因此，它可以在剩余类别 3 和 4 上创建一个多项式分布，同时这些类别的比例保持不变。

多项式分布是二项式分布的推广。如果 k=2，多项式分布简化为二项式分布。因此，如果我们有

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

那么 X₁ ~ Bin(n, p₁) 和 X₂ ~ Bin(n, p₂)。

证明（可选）：我们应该注意到 X₂= n-X₁ 和 p₂=1-p₁。因此，如果我们将结果 1 视为成功，将结果 2 视为失败，则每次试验只有两种结果（成功和失败），成功的概率是 p₁。因此，每次试验可以用具有参数 p₁ 的伯努利分布表示，X₁ 代表了 n 次伯努利试验中成功的次数，参数为 p₁。因此，我们得出 X₁ 具有参数 n 和 p₁ 的二项分布。类似地，如果我们将结果 2 视为成功，将结果 1 视为失败，则 X₂ 代表了 n 次伯努利试验中成功的次数，参数为 p₂，因此具有参数 n 和 p₂ 的二项分布。

我们可以将之前的结果扩展到 k>2 的多项分布。如果

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

那么每个 Xᵢ 的边际分布是具有参数 n 和 pᵢ 的二项分布：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

证明（可选）：我们假设结果 i 的成功概率为 pᵢ，所有其他结果的失败概率为 1-pᵢ。因此，Xᵢ 代表了 n 次伯努利试验中成功的次数，参数为 pᵢ，并且它具有参数 n 和 pᵢ 的二项分布。

此外，X 的某些元素的总和服从二项分布。从数学上讲，如果 *Xᵢ,*₁, *Xᵢ,*₂, …, Xᵢ,ₘ 是随机向量 X 的 m 个元素（m<k），它们在向量 p 中的对应概率是 *pᵢ,*₁, pᵢ,₂, …, pᵢ,ₘ，那么总和 *Xᵢ,*₁ + *Xᵢ,*₂ + …+ Xᵢ,ₘ 服从具有参数 n 和 *pᵢ,*₁ + *pᵢ,₂ + …+ pᵢ,ₘ 的二项分布。例如，对于 n=8 和 k=5，X₁+X₃+X₅ 具有参数 8 和 p₁ + p₃ + p₅ 的二项分布。

证明（可选）：如前所述，我们可以将多项式随机向量 X 中的多个元素合并，以获得一个新的多项式随机向量。因此，通过合并 *Xᵢ,*₁ + *Xᵢ,*₂ + …+ Xᵢ,ₘ，我们得到一个新的随机向量 Y，它具有多项分布：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

我们现在还展示了多项式随机向量每个元素的边际分布具有参数 n 和其对应概率的二项分布，因此我们得出结论 *Xᵢ,*₁ + *Xᵢ,*₂ + …+ Xᵢ,ₘ 具有参数 n 和 *pᵢ,*₁ + *pᵢ,*₂ + …+ pᵢ,ₘ 的二项分布。

我们可以使用 Python 中的 scipy 库生成多项分布。我们可以使用 scipy.stat 中的 multinomial 对象创建多项分布，并使用该对象的 pmf() 方法计算其联合 PMF。该对象接受 n 和 p 参数，其中 n 和 p 对应于 n 和 p。p 是一个类似数组的对象。p 的每个元素应在 [0, 1] 区间内，并且这些元素的和应为 1。如果它们的和不为 1，则 p 的最后一个元素不会被使用，而是用之前元素和的 1 减去替代（以使元素和为 1）。

首先，我们需要导入所有必需的库：

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import binom, multinomial
%matplotlib inline

列表 1 创建了以下多项分布

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

并计算其联合 PMF 在

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

# Listing 1

n = 5
p=[0.5, 0.3, 0.2]
mult = multinomial(n=n, p=p)
mult.pmf([3, 1, 1])

# Output
0.15

我们也可以直接将分布参数传递给 pmf() 函数：

multinomial.pmf([3, 1, 1], n=n, p=p)

# Output
0.15

方法 rvs() 可以用于生成随机变异。随机变异或简称变异是随机变量或随机向量的特定结果。使用此方法，我们可以从多项分布中抽取大小为 m 的随机样本，这意味着我们从具有多项分布的随机向量中生成 m 个随机变异。例如，要从方程 1 中的多项分布中抽取大小为 4 的随机样本，我们可以使用以下代码片段：

np.random.seed(2)
multinomial.rvs(n=n, p=p, size = 4)

# Output
array([[2, 3, 0],
       [3, 1, 1],
       [2, 2, 1],
       [2, 2, 1]])

它返回一个包含 4 个元素的 2d 数组。每个元素是一个数组，表示方程 1 中的随机向量可以取的可能值 X。

让我们看看如何可视化方程 1 定义的多项分布。在这里，我们希望收集X可能取的所有值。为此，我们只需计算元组 (X₁, X₂) 所有可能的值。一旦我们知道了 X₁ 和 X₂ 的值，就可以利用 X₃=n-X₁-X₂ 计算 X₃ 的值。由于 n=5，Xᵢ 可以取 0 到 5 之间的整数值。列表 2 使用 numpy 中的 meshgrid() 方法获取 (X₁, X₂) 的所有可能值及其对应的 X₃ 值。最后，(X₁, X₂, X₃) 的不同值被存储在数组 X_mat 中。我们得到了一些负值的 X₃。虽然这些负值是不可能取的，但我们无需丢弃它们，因为多项分布中任何 Xᵢ 为负值的联合 PMF 仅为零。

# Listing 2

n = 5
x1_array = np.arange(0, n+1)
X1_mat, X2_mat = np.meshgrid(x1_array, x1_array)
x3_array = n - X1_mat.flatten() - X2_mat.flatten()

X_mat = np.array([X1_mat.flatten(), X2_mat.flatten(), x3_array]).T
X_mat

# Output
array([[ 0,  0,  5],
       [ 1,  0,  4],
       [ 2,  0,  3],
       [ 3,  0,  2],
       [ 4,  0,  1],
       [ 5,  0,  0],
       [ 0,  1,  4],
       [ 1,  1,  3],
       ...

现在我们可以使用 X_mat 计算 X₁ 在 X₁=3 时的边际 PMF。为此，我们需要对所有值为 3 的 X 的联合 PMF 进行求和。

multinomial(n, p).pmf(X_mat[X_mat[:, 0]==3]).sum()

# Output
0.3125

我们还知道X₁的边际分布是参数n和p₁的二项分布。因此，我们可以使用参数n*=5 和p₁=0.5 的二项分布的 PMF。为了计算二项分布的 PMF，我们可以使用scipy.stat中的binom对象。我们计算X₁=3 处的 PMF 值，结果应该与之前的代码片段相同。

binom.pmf(k=3,n=n, p=p[0])

# Output
0.31249999999999983

结果几乎相同，但由于数值误差存在小的差异。我们还可以通过从多项分布中抽样来近似X₁的边际 PMF。列表 3 从方程 1 中的多项分布中抽取了大小为 100,000 的样本，并绘制了该样本中X₁值的条形图。结果如图 4 所示。正如图中所示，样本的条形图与参数n=5 和p₁=0.5 的二项分布的 PDF 匹配。

# Listing 3

np.random.seed(2)
n = 5
p=[0.5, 0.3, 0.2]
x = np.arange(n+1)
sample = multinomial.rvs(n=n, p=p, size=100000)
pmf_binomial = binom.pmf(k = x,n = n, p = p[0])
sample_marginal = sample[:, 0] 
values, counts = np.unique(sample_marginal, return_counts=True)
probs = counts / counts.sum()

fig = plt.figure(figsize=(8, 6))
plt.bar(values, probs, label='Marginal sample')
plt.plot(x, pmf_binomial, marker='o', color='red',
         label='Binomial distribution, n={}, p={}'.format(n, p[0]))
plt.xlabel('$x_1$', fontsize = 18)
plt.ylabel('Probability', fontsize=18)
plt.ylim([0, 0.45])
plt.legend(loc='best', fontsize = 15)
plt.show()

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

图 4

如果一个具有多项分布的随机变量只有 3 个元素（k=3），我们可以轻松地在 2d 或 3d 空间中绘制其联合 PMF。列表 4 绘制了图 5（顶部）中先前多项分布的联合 PMF。图是X₁与X₂的热图。我们不需要在此图中包含X₃，因为其值依赖于X₁与X₂（X₃=n-X₁-X）。X₁的二项分布 PMF 也绘制在图 5（底部）中。

# Listing 4

n = 5
p=[0.5, 0.3, 0.2]
x = np.arange(n+1)

x1_array = np.arange(0, n+1)
X1_mat, X2_mat = np.meshgrid(x1_array, x1_array)
x3_array = n - X1_mat.flatten() - X2_mat.flatten()
X_mat = np.array([X1_mat.flatten(), X2_mat.flatten(), x3_array]).T

pmf_mult = multinomial(n, p).pmf(X_mat)
pmf_binomial = binom.pmf(k = x, n = n, p = p[0])
pmf_grid = pmf_mult.reshape(6, 6)

f, (ax1, ax2) = plt.subplots(2, 1, figsize=(9, 11))
plt.subplots_adjust(hspace=0.4)

heatmap = ax1.pcolor(pmf_grid, cmap='coolwarm')
ax1.set_xlabel('$x_1$', fontsize = 18)
ax1.set_ylabel('$x_2$', fontsize = 18)

for x1_pos in range(pmf_grid.shape[0]):
    for x2_pos in range(pmf_grid.shape[1]):
        ax1.text(x1_pos + 0.5, x2_pos + 0.5,
                '%.5f' % pmf_grid[x2_pos, x1_pos],
                 horizontalalignment='center',
                 verticalalignment='center',
                 )

ax1.set_yticks(np.arange(pmf_grid.shape[0])+0.5)
ax1.set_yticklabels(np.arange(pmf_grid.shape[0]))

ax1.set_xticks(np.arange(pmf_grid.shape[1])+0.5)
ax1.set_xticklabels(np.arange(pmf_grid.shape[1]))
plt.colorbar(heatmap, ax=ax1)

ax2.bar(x, pmf_binomial,
        label='Binomial\ndisitrbution\nn={}, p={}'.format(n, p[0]))
ax2.bar_label(ax2.containers[0], fontsize = 13)
ax2.set_xlabel('$x_1$', fontsize = 18)
ax2.set_ylim([0, 0.4])
ax2.set_ylabel('Probability', fontsize=18)
ax2.legend(loc='best', fontsize = 14)

plt.show()

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

图 5（作者提供的图片）

热图中每列在X₁=x处包含所有值的联合 PMF，其中X₁等于x。因此，X₁在x处的边际 PMF 等于X₁=x列的总和（图 5）。类似地，X₂在x处的边际 PMF 等于X₂=x行的总和。

我们还可以使用 3d 条形图来可视化方程 1 中多项分布的联合 PMF。如前所述，我们只需要在图中包含X₁和X₂，因为X₃的值依赖于它们。列表 5 使用pmf()方法计算X_mat中数组的联合 PMF，并在图 6 中创建这些联合 PMF 值的 3d 条形图。此图中还显示了X₁和X₂的边际分布。

# Listing 5

n = 5
p=[0.5, 0.3, 0.2]

x1_array = np.arange(0, n+1)
X1_mat, X2_mat = np.meshgrid(x1_array, x1_array)
x3_array = n - X1_mat.flatten() - X2_mat.flatten()
X_mat = np.array([X1_mat.flatten(), X2_mat.flatten(), x3_array]).T

pmf_mult = multinomial(n, p).pmf(X_mat)
pmf_x1 = binom.pmf(k=x1_array,n=n, p=p[0])
pmf_x2 = binom.pmf(k=x1_array,n=n, p=p[1])

fig = plt.figure(figsize=(10, 10))
ax1 = fig.add_subplot(111, projection='3d')

x1 = X_mat[:, 0]
x2 = X_mat[:, 1]
z = np.zeros(len(x1))

width = 0.8
dx1 = np.repeat(width, len(x1))
dx2 = np.repeat(width, len(x1))

ax1.bar3d(x1-width/2, x2-width/2,
          z, dx1, dx2, pmf_mult, color='aqua')
ax1.bar3d(x1_array-width/2, np.repeat(6, n+1)-width/2,
          np.zeros(n+1), np.repeat(0.8, n+1),
          np.zeros(n+1), pmf_x1, color='blue', shade=False)
ax1.bar3d(np.repeat(0, n+1)-width/2, x1_array-width/2,
          np.zeros(n+1), np.zeros(n+1),
          np.repeat(0.8, n+1), pmf_x2,
          color='blue', shade=False)

ax1.set_xlabel('$x_1$', fontsize=20)
ax1.set_ylabel('$x_2$', fontsize=20)
ax1.set_zlabel("$p_\mathregular{X}(\mathregular{x})$", weight="bold",
               style="italic", fontsize=18, labelpad = 8)

ax1.text(4.2, 5, 0.25,"$p_{X_1}(x_1)$", fontsize= 18, color='b')
ax1.text(0.7, 1.7, 0.3, "$p_{X_2}(x_2)$", fontsize= 18, color='b')
ax1.text(1.3, 7, 0.25,"$X_1 \sim Bin(5, 0.5)$", fontsize= 18)
ax1.text(-3, 1.7, 0.3, "$X_2 \sim Bin(5, 0.3)$", fontsize= 18)

ax1.view_init(35, -45)
plt.show()

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

图 6

我们可以将方程 1 中的多项分布的前两个元素合并，得到新的随机向量。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

我们还知道这个随机向量具有参数n的多项分布。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

列表 6 通过从多项分布中采样来近似 X₁+X₂ 的边际 PMF。它从方程 86 中的多项分布中抽取一个 100,000 大小的样本，并绘制该样本中 X₁+X₂ 值的条形图。结果如图 7 所示，你可以看到条形图与参数 n=5 和 p₁+p₂ 的二项分布的 PDF 相匹配。

# Listing 6

np.random.seed(2)
n = 5
p=[0.5, 0.3, 0.2]
x = np.arange(n+1)
sample = multinomial.rvs(n=n, p=p, size=100000)
pmf_binomial = binom.pmf(k = x,n = n, p = p[0]+p[1])
sample_marginal = sample[:, 0] + sample[:, 1] 
values, counts = np.unique(sample_marginal, return_counts=True)
probs = counts / counts.sum()

fig = plt.figure(figsize=(8, 6))
plt.bar(values, probs, label='Marginal sample')
plt.plot(x, pmf_binomial, marker='o', color='red',
         label='Binomial distribution, n={}, p={}'.format(n, p[0]+p[1]))
plt.xlabel('$x_1+x_2$', fontsize = 18)
plt.ylabel('Probability', fontsize=18)
plt.ylim([0, 0.45])
plt.legend(loc='best', fontsize = 13)

plt.show()

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

图 7

列表 7 展示了如何从具有 k 个不同类别的项的总体中进行有放回采样可以产生多项分布（图 3）。这里总体由一个具有 3 个独特元素的列表表示：1、2 和 3。它们的比例分别是 0.5、0.3 和 0.2。我们可以从这个列表中随机选择 5 个元素，并计算 1、2 和 3 的出现次数（我们用 x₁、x₂ 和 x₃ 来表示）。现在这个向量

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

是方程 1 中多项分布的随机变量。我们也可以将其称为从该分布中抽取的大小为 1 的随机样本。请注意，我们也可以使用 rvs() 方法生成这个随机变量。现在我们可以重复相同的过程 n 次，从多项分布中抽取一个大小为 n 的随机样本。列表 7 在图 8（蓝色条形）中绘制了两个随机样本的条形图。这些样本的大小分别为 30 和 50000。图中还绘制了方程 1 中定义的多项分布的联合 PMF 的条形图（红色条形）。随着样本大小的增加，样本条形图的形状越来越接近多项分布联合 PMF 的条形图的形状。

# Listing 7

np.random.seed(1)
population = [1]*5 + [2]*3 + [3]*2
n = 5
xedges = np.arange(n+2)
yedges =  np.arange(n+2)
num_samples_list = [30, 50000]

x1_array = np.arange(0, n+1)
X1_mat, X2_mat = np.meshgrid(x1_array, x1_array)
x3_array = n - X1_mat.flatten() - X2_mat.flatten()
X_mat = np.array([X1_mat.flatten(), X2_mat.flatten(), x3_array]).T
pmf_mult = multinomial(n, p).pmf(X_mat)

x1 = X_mat[:, 0]
x2 = X_mat[:, 1]
z = np.zeros(len(x1))
width = 0.15
dx1 = np.repeat(width, len(x1))
dx2 = np.repeat(width, len(x1))

fig = plt.figure(figsize=(18, 10))
plt.subplots_adjust(wspace=0.1)
ax1 = fig.add_subplot(121, projection='3d')
ax2 = fig.add_subplot(122, projection='3d')
axs = [ax1, ax2]

for i, num_samples in enumerate(num_samples_list):
    samples = np.random.choice(population,
                 size=n*num_samples).reshape(num_samples, n)
    samples_count = np.stack(((samples==1).sum(axis=1),
                              (samples==2).sum(axis=1)), axis=-1)
    H, _, _ = np.histogram2d(samples_count[:, 0],
                             samples_count[:, 1],
                             bins=(xedges, yedges))
    H = H.T / H.sum()

    axs[i].bar3d(x1-2*width, x2-width, z, dx1, dx2, H.flatten(),
              color='aqua')
    axs[i].bar3d(x1+width/2, x2-width, z, dx1, dx2, pmf_mult,
              color='red')
    axs[i].set_xlabel('$x_1$', fontsize=20)
    axs[i].set_ylabel('$x_2$', fontsize=20)
    axs[i].set_zlabel("$p_\mathregular{X}(\mathregular{x})$",
                      weight="bold", style="italic",
                      fontsize=18, labelpad = 8)
    axs[i].set_title('Sample size={}'.format(num_samples), fontsize=20)
    axs[i].set_zlim([0, 0.2])
    axs[i].view_init(35, -135)

plt.show()

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

图 8

在列表 8 中，我们比较了在方程 1 中定义的多项分布的随机变量的联合 PMF 与具有 n 个随机向量的和的联合 PMF，这些随机向量具有参数 n=1 和 p。我们首先从方程 1 中定义的多项分布中抽取一个 5000000 大小的样本（sample1），该多项分布具有以下参数：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

然后我们从 5 个具有相同p和n=1 的多项分布中抽取一个 5000000 大小的样本。我们将这 5 个样本相加，得到一个 5000000 大小的样本并存储在sample2中。实际上，这个样本包含了随机向量

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

可以取。

# Listing 8

np.random.seed(50)
num_samples = 5000000
n = 5
p=[0.5, 0.3, 0.2]
sample1 = multinomial.rvs(n=n, p=p, size=num_samples)
sample2 = multinomial.rvs(n=1, p=p,
    size=n*num_samples).reshape(num_samples, n, len(p)).sum(axis = 1)

我们现在可以比较在sample1和sample2中获得特定值的概率。

(sample1 == [3, 2, 0]).all(axis=1).mean()

# output
0.1124632

(sample2 == [3, 2, 0]).all(axis=1).mean()

# Output
0.1125352

正如你所见，概率非常接近，微小的差异是由于样本大小的不同。随着样本大小趋于无穷大，这些概率趋于相同。

如前所述，如果我们有：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

然后每个Xᵢ具有参数n和p的二项分布。因此，利用二项分布的性质我们得出：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

我们定义随机向量X的均值为

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

我们还展示了X的一些元素的和具有二项分布。因此，对于每对i、j (i、j=1…k，i≠j)，Xᵢ + Xⱼ具有参数n和pᵢ + pⱼ的二项分布。因此，我们有：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

可以证明：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

由此可以得出：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

由于n、pᵢ和pⱼ是正数，我们得出Xᵢ和Xⱼ的协方差总是负数。让我们看看为什么相关性是负的。两个随机变量之间的负相关意味着当其中一个变量较高时，另一个变量倾向于较低，并且当一个变量增加时，另一个变量倾向于减少。我们知道X₁+X₂+ …+ Xₖ=n，所以假设X的所有其他分量（所有m ≠ i、j的Xm）都是常数，那么任何Xᵢ的增加都应该导致Xⱼ的减少，反之亦然。

我们定义X的协方差矩阵为：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

所以，协方差矩阵中的(i, j)元素是变量Xᵢ和Xⱼ的协方差，而第i个对角元素给出了Xᵢ的方差。根据协方差的定义，我们知道Cov(Xᵢ, Xⱼ)=Cov(Xⱼ, Xᵢ)。因此，协方差矩阵是对称矩阵，其中的(i, j)元素等于(j, i)元素。

我们可以使用 Python 轻松计算在方程 1 中定义的多项分布的均值和协方差矩阵。方法mean()返回方程 2 中定义的均值向量：

n = 5
p=[0.5, 0.3, 0.2]
multinomial(n, p).mean()

# Output
array([2.5, 1.5, 1\. ])

方法cov()返回协方差矩阵（方程 3）：

multinomial(n, p).cov()

# Output
array([[ 1.25, -0.75, -0.5 ],
       [-0.75,  1.05, -0.3 ],
       [-0.5 , -0.3 ,  0.8 ]])

在这篇文章中，我们讨论了多项分布背后的数学原理，并展示了如何在 Python 中实现它。多项分布在科学、工程和金融领域被广泛使用。它可以用于那些有超过两个可能结果的应用场景，并且系统不能通过成功-失败的描述进行建模。

希望你喜欢阅读这篇文章。如果你有任何问题或建议，请告诉我。本文中的所有代码列表可以从 GitHub 上作为 Jupyter 笔记本下载，链接为：github.com/reza-bagheri/probability_distributions/blob/main/multinomial_distribution.ipynb

理解 Naive Bayes 算法

原文：towardsdatascience.com/understanding-naive-bayes-algorithm-d753d3b76727

它是什么以及如何将其应用于实际场景

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传 Meagan Burkhart

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

–

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

照片由Google DeepMind提供，刊登在Unsplash

今年，我的决心是回到数据科学的基础。我每天都在处理数据，但如果你在完成重复的任务，很容易忘记一些核心算法的运作。我计划每周在 Towards Data Science 上深入探讨一个数据算法。这一周，我将介绍 Naive Bayes。

如何发音 Naive Bayes

只是为了澄清，你可以在这里学习如何发音 Naive Bayes。

既然我们知道怎么说它，我们来看看它的含义……

什么是 Naive Bayes 分类器？

这个概率分类器基于贝叶斯定理，可以总结如下：

当第二个事件已经发生时，第一个事件的条件概率是“事件 B 发生的情况下 A 的概率和 A 的概率除以事件 B 的概率”的乘积。

P(A|B) = P(B|A)P(A) / P(B)

一个常见的误解是贝叶斯定理和条件概率是同义词。

然而，有一个区别——贝叶斯定理使用条件概率的定义来找出所谓的“反向概率”或“逆概率”。

换句话说，条件概率是 A 发生的情况下 B 的概率。贝叶斯定理则找出 B 发生的情况下 A 的概率。

朴素贝叶斯算法的一个显著特征是其使用序列事件。简单来说，通过获取后续的附加信息，初始概率会进行调整。我们将这些称为先验概率/边际概率和后验概率。主要的结论是，通过了解另一个条件的结果，初始概率会发生变化。

朴素贝叶斯的应用

一个好的例子是医疗检测。例如，如果患者正在处理胃肠问题，医生可能会怀疑炎症性肠病（IBD）。这种情况的初始概率约为 1.3%。

然而，根据患者的症状，医生会要求进行C-反应蛋白（CRP）血液检查，这是一种炎症标志物。如果该标志物高于某个阈值（>50mg/l），则患有炎症性肠病的概率增加。

随后，医生进行结肠镜检查，结果显示有明显的炎症。现在，IBD 诊断的概率再次增加。

我们将通过本文的例子来理解不同的朴素贝叶斯分类器如何在医疗保健中应用。

朴素贝叶斯分类器的假设

正如其名字所示，朴素贝叶斯分类器与算法的特征是“朴素”或条件独立的相关。

另一个假设是模型中的特征对最终结果的贡献是相等的。然而，这些假设在现实世界中经常被违反。其目的是简化计算单一概率的数学模型。

朴素贝叶斯分类器的优缺点

朴素贝叶斯分类器的一个优点是它在小样本上表现良好，即使在违反一些上述严格假设的情况下。

根据具体问题的不同，朴素贝叶斯分类器有一些缺点。

例如，与回归模型不同，朴素贝叶斯并不打算衡量特征的重要性——主要是因为假设所有特征对结果的贡献是相等的。
此外，还存在一个被称为“零频率问题”的情况。当某个特征组合在训练数据中出现的概率为零时，所有概率的乘积也将等于零。为了克服这一问题，分析师通常会在值中加上 1，以防其等于零的可能性。

朴素贝叶斯分类器通常用于文本分类任务，包括识别“垃圾邮件”和标记情感。

朴素贝叶斯模型的不同类型

以下三种类型的朴素贝叶斯模型主要在关于条件概率分布的假设上有所不同。

1. 伯努利朴素贝叶斯

当变量有二进制分布时，使用伯努利朴素贝叶斯。这意味着所有特征都是分类的，可以取值 1（存在）或 0（不存在）。

在使用 Python 进行朴素贝叶斯时，我们可以使用 [sklearn.naive_bayes](https://scikit-learn.org/stable/modules/classes.html#module-sklearn.naive_bayes) 的 BernoulliNB 包将尚未是分类变量的变量二值化。如果所有特征的截断点相同，这会很有用。例如，可以使用 0。然而，也有其他方法可以对连续数据进行二值化，包括 Pandas 的 cut 函数，如下所示：

df['CRGlevel_high'] **=** pd.cut(x**=**df['CRG'], bins**=**[0, 50, np.inf], labels**=**['No', 'Yes'])

2. 高斯朴素贝叶斯分类器

关于高斯朴素贝叶斯，首先要注意的是模型中的特征是连续且正态分布的。基于这一假设，在对模型中的数据点进行分类时，算法将最终值分配给具有最大后验概率的类别。

高斯朴素贝叶斯可以应用于常见的 Iris 数据集。运行此算法的基本 Python 代码如下：

from sklearn.naive_bayes import GaussianNB
gnb = GaussianNB()

# Train the classifier on the training data
gnb.fit(X_train, y_train)

高斯 NB 模块的文档可以在这里找到。

关于高斯朴素贝叶斯，值得了解的一点是，如果算法的假设得到满足，则模型的功能与逻辑回归相同。

假设从 Bhowmik（2015）的一篇文章中总结如下：

Y 是布尔值，受伯努利分布控制
Xi ∼ N(µij , σi)，对所有的 i 适用
对于所有 i 和 k 6= i，Xi 和 Xk 在给定 Y 的条件下是条件独立的。

“如果朴素贝叶斯假设成立，NB 和 LR 产生的模型在渐近意义上是相同的。” 来源

在我们使用的示例中，可能需要考虑来自血液检查的多个值来确定一个人是否被诊断为 IBD。因此，我们可以输入原始的 CRG 值，以及其他 IBD 标志物，如 Calprotectin 和 SedRate，而不是二值化连续变量。

3. 多项式朴素贝叶斯

当特征值表示相对频率而非二进制字段时，我们可以使用多项式朴素贝叶斯。

这对于文本分类问题很有用，这些问题通常依赖于 token 的频率。在我们开始讨论的 IBD 示例中，我们可以考虑使用电子健康记录（EHR）笔记来应用多项式朴素贝叶斯。

根据与患者症状相关的笔记，模型可能会使用“炎症”、“腹泻”、“便秘”等关键字的频率来判断患者是否可能被诊断为溃疡性结肠炎（UC）或克罗恩病。

当然，仅凭文本并不能确保做出诊断，但一个潜在的使用案例是使用 NB 分类器根据患者的症状识别出患有炎症性肠病（IBD）风险的患者，并鼓励医疗从业者订购额外的检查，以排除差异诊断或确认诊断。

我在医疗保健中的建议应用

总之，朴素贝叶斯和贝叶斯定理做出了一些假设，这些假设在现实世界中可能并不适用。然而，当你处理一个小数据集并寻找一个好的基线模型时，NB 可以很好地工作，因为它计算简单且易于理解。

基于我的研究，医疗保健中的潜在应用似乎很有前景。例如，在排除差异诊断时，可以使用朴素贝叶斯分类器来确定需要进行多少测试才能得出最终诊断。由于每个特征被期望是独立的并且贡献相等，这可能有助于指导从业者确定哪种测试组合是“足够的”以达成结论性诊断。

我预期这可以减少不必要的测试，并缩短最终诊断的时间，这样可以为医生节省时间，为健康保险计划减少费用，并为患者减少接受治疗和进行侵入性检查所需的时间。

了解 NeRFs

原文：towardsdatascience.com/understanding-nerfs-2a082e13c6eb

场景表示的重大突破

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传 Cameron R. Wolfe, Ph.D.

·发布于 Towards Data Science ·11 min read·2023 年 4 月 28 日

–

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

(照片由 nuddle 提供，来源于 Unsplash)

正如我们通过 DeepSDF [2] 和 SRNs [4] 等方法所见，将 3D 对象和场景编码在前馈神经网络的权重中是一种内存高效的隐式 3D 数据表示，既准确又高分辨率。然而，到目前为止，我们看到的方法还无法以足够的保真度捕捉现实和复杂的场景。相反，离散表示（例如三角网格或体素网格）在内存分配充足的情况下，能产生更准确的表示。

这在 Neural Radiance Fields (NeRFs) [1] 的提议下发生了改变，NeRFs 使用前馈神经网络来建模场景和物体的连续表示。NeRFs 使用的表示称为辐射场，与先前提案有些不同。特别是，NeRFs 将五维坐标（即空间位置和视角方向）映射到体积密度和视角相关的 RGB 颜色。通过在不同视角和位置上累积这些密度和外观信息，我们可以渲染出逼真的新视角场景。

像SRNs [4]一样，NeRFs 可以仅使用一组图像（以及它们相关的相机姿态）对底层场景进行训练。与之前的方法相比，NeRF 渲染在质量和数量上都更好。值得注意的是，NeRFs 甚至可以捕捉复杂效果，如物体表面的视角依赖反射。通过在前馈神经网络的权重中隐式建模场景，我们在不需要过多内存成本的情况下实现了离散场景表示的准确性。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

（摘自[1]）

为什么这篇论文很重要？ 本帖子是我关于 3D 形状和场景深度学习系列的一部分。NeRFs 在这一领域是一个革命性的提案，因为它们能够从任意视角实现极其准确的 3D 重建。NeRFs 生成的场景表示质量非常高，我们将在本帖的其余部分看到这一点。

背景

理解 NeRFs 所需的大部分背景概念已在之前的帖子中涵盖，包括：

前馈神经网络 [link]
表示 3D 对象 [link]
离散表示的问题 [link]

在介绍 NeRFs 的工作原理之前，我们只需要覆盖几个更多的背景概念。

位置编码

NeRFs 不是直接将[x, y, z]坐标作为输入传递给神经网络，而是将这些坐标转换为高维位置嵌入。我们在之前关于transformer 架构的帖子中讨论过位置嵌入，因为位置嵌入用于向自注意力模块提供令牌排序和位置的概念。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

（摘自[1]）

简而言之，位置嵌入将一个标量数作为输入（例如，坐标值或表示序列中位置的索引），并产生一个高维向量作为输出。我们可以在训练期间学习这些嵌入，或者使用固定函数生成它们。对于 NeRFs，我们使用上面显示的函数，该函数将一个标量p作为输入，并产生一个2L维的位置编码作为输出。

其他内容

在本概述中，我们可能会遇到一些其他（可能）不熟悉的术语。让我们现在快速澄清这些术语。

端到端训练。 如果我们说一个神经架构可以“端到端”学习，这仅仅意味着系统的所有组件都是可微分的。因此，当我们为某些数据计算输出并应用我们的损失函数时，我们可以通过整个系统（即，端到端）进行微分，并用梯度下降进行训练！

并非所有系统都可以进行端到端训练。例如，如果我们正在建模表格数据，我们可能会执行特征提取过程（例如，一热编码），然后在这些特征上训练机器学习模型。由于特征提取过程是手工设计的且不可微分，我们不能对系统进行端到端训练！

朗伯特反射。 在阅读 NeRFs 之前，这个术语对我完全陌生。朗伯特反射指的是物体表面的反射特性。如果一个物体的表面是哑光的，并且从不同角度看不会改变，我们称这个物体是朗伯特型的。另一方面，一个“光亮”的物体，其反射光线的方式根据观察角度不同而变化，则被称为非朗伯特型的。

建模 NeRFs

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

（来自 [1]）

使用 NeRF 渲染场景视角的高层次过程如下：

使用射线行进方法生成场景的 3D 点和视角方向样本。
将点和视角方向作为输入提供给前馈神经网络，以生成颜色和密度输出。
执行体积渲染，将场景中的颜色和密度累积到 2D 图像中。

我们现在将详细解释这个过程的每个组件。

辐射场。 如前所述，NeRFs 模型是一个 5D 向量值（即，意味着函数输出多个值）函数，称为辐射场。该函数的输入是一个 [x, y, z] 空间位置和一个 2D 视角方向。视角方向有两个维度，对应于表示 3D 空间中方向的两个角度；见下文。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

3D 空间中的方向可以用两个角度表示。

实际上，视角方向仅表示为一个 3D 笛卡尔单位向量。

这个函数的输出有两个组成部分：体积密度和颜色。颜色只是一个 RGB 值。然而，这个值是视角依赖的，这意味着给定不同的视角方向作为输入，颜色输出可能会变化！这种特性使 NeRFs 能够捕捉反射和其他视角依赖的外观效果。相比之下，体积密度仅依赖于空间位置，捕捉不透明度（即，光线通过该位置时积累的程度）。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

NeRF 使用前馈神经网络建模辐射场（来自[1]）

神经网络。 在[1]中，我们用前馈神经网络建模辐射场，该网络接收 5D 输入并训练以生成相应的颜色和体积密度作为输出；见上文。然而，请记住，颜色是视角相关的，而体积密度则不是。为了解决这个问题，我们首先将输入的 3D 坐标通过几个前馈层，这些层产生体积密度和特征向量作为输出。然后，将这个特征向量与视角方向连接，并通过一个额外的前馈层来预测视角相关的 RGB 颜色；见下文。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

NeRF 的前馈架构。

体积渲染（简要说明）。 体积渲染是一个太复杂的话题，无法在这里深入讨论，但我们应该了解以下内容：

它可以从离散数据样本（例如颜色和密度值）生成底层场景的图像。
它是可微的。

对于那些对体积渲染的更多细节感兴趣的人，请查看这里的解释和[1]的第四部分。

大局观。 NeRF 使用前馈网络生成关于场景几何和外观的相关信息，这些信息沿着众多不同的相机光线（即从特定相机视点向场景中的某个方向移动的 3D 空间中的一条线）进行传递，然后使用渲染将这些信息汇聚成一张 2D 图像。

值得注意的是，这两个组件都是可微的，这意味着我们可以端到端地训练整个系统！给定一组带有对应相机姿态的图像，我们可以通过仅生成/渲染已知视角，并使用(随机)梯度下降来最小化 NeRF 输出与实际图像之间的误差，来训练 NeRF 以生成新的场景视点；见下文。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

（来自[1]）

一些额外细节。 我们现在了解了 NeRF 的大部分组件。然而，迄今为止我们描述的方法在[1]中实际上被证明是低效的，并且通常不能很好地表示场景。为了改进模型，我们可以：

用位置嵌入替换空间坐标（包括空间位置和视角方向）。
采用层次采样方法进行体积渲染。

通过使用位置嵌入，我们将前馈网络的输入（即空间位置和视角方向坐标）映射到更高维度。先前的工作表明，与直接使用空间或方向坐标作为输入相比，这种方法使神经网络更好地建模场景的高频（即变化很大/很快）特征[5]。这使得 NeRF 的输出质量大大提高；见下文。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

（来自 [1]）

NeRF 使用的分层采样方法使渲染过程更高效，只对可能影响最终渲染结果的地点和视角方向进行采样（并通过前馈神经网络传递）。这样，我们仅在需要的地方评估神经网络，避免在空白或遮挡区域浪费计算。

给我们展示一些结果吧！

NeRFs 被训练为一次仅表示一个场景，并在多个数据集上进行评估，这些数据集包括合成和真实对象。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

（来自 [1]）

如上表所示，NeRFs 显著超越了 SRNs [4] 和 LLFF [6] 等替代方法。除了定量结果之外，与其他方法的输出进行视觉对比也非常有用。首先，我们可以立即看出，使用位置编码和以视角依赖的方式建模颜色是非常重要的；见下文。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

（来自 [1]）

我们会立即注意到的一个改进是，NeRFs — 因为它们以视角依赖的方式建模颜色 — 能够捕捉场景中的复杂反射（即非朗伯特特征）和视角依赖的模式。此外，NeRFs 还能够以惊人的精度建模底层几何体的复杂方面；见下文。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

（来自 [1]）

NeRF 场景表示的质量在以视频形式查看时最为明显。正如下面的视频所示，NeRFs 以令人印象深刻的准确性和一致性建模了底层场景，并在不同视角之间保持一致。

要获取更多可以使用 NeRF 生成的逼真场景视角的示例，我强烈建议查看这里链接的项目网站！

主要结论

正如我们在评估中所见，NeRFs 在场景表示质量上取得了巨大的突破。因此，这项技术在人工智能和计算机视觉研究社区中获得了极大的关注。由于其场景表示的质量，NeRF 的潜在应用（例如虚拟现实、机器人等）几乎是无穷无尽的。主要结论列在下方。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

（来自 [1]）

NeRFs 捕捉复杂细节。 使用 NeRFs，我们能够捕捉场景中的细粒度细节，例如船只中的装配材料；见上文。除了几何细节，NeRFs 还能处理非朗伯特效应（即反射和基于视角的颜色变化），因为它们以视角依赖的方式建模颜色。

我们需要智能采样。 到目前为止，我们见过的所有 3D 场景建模方法都使用神经网络在 3D 空间上建模一个函数。这些神经网络通常会在所考虑的空间体积内的每个空间位置和方向上进行评估，如果处理不当，可能会非常昂贵。对于 NeRF，我们使用一种分层采样方法，仅评估可能影响最终渲染图像的区域，这大大提高了采样效率。类似的方法也被先前的工作采用；例如， ONets [3] 使用基于 octree 的分层采样方法更高效地提取对象表示。

位置嵌入效果很好。 到目前为止，我们见过的大多数场景表示方法都直接将坐标值作为输入传递给前馈神经网络。通过 NeRF，我们发现将这些坐标进行位置嵌入要好得多。特别是，将坐标映射到更高维度似乎允许神经网络捕捉场景几何和外观的高频变化，这使得结果场景渲染更加准确且视角一致。

仍然节省内存。 NeRF 隐式地建模了基础场景的连续表示。这种表示可以在任意精度下进行评估，并且具有固定的内存成本——我们只需存储神经网络的参数！因此，NeRF 可以在不占用大量内存的情况下提供准确的高分辨率场景表示。

“至关重要的是，我们的方法克服了在高分辨率下建模复杂场景时离散体素网格的高昂存储成本。” — 来源于 [1]

局限性。 尽管 NeRF 在先进技术方面取得了显著进展，但它们并不完美——表示质量仍有改进空间。然而，NeRF 的主要局限性在于它们每次只能建模一个场景，并且训练成本高（即，每个场景在单个 GPU 上需要 2 天）。未来在这一领域的进展如何找到更高效的生成 NeRF 质量场景表示的方法将会非常有趣。

结束语

非常感谢你阅读这篇文章。我是 Cameron R. Wolfe，Rebuy 的人工智能总监，也是莱斯大学的博士生。我研究深度学习的经验和理论基础。你也可以查看我在 medium 上的其他文章！如果你喜欢这篇文章，请关注我的 twitter 或订阅我的 Deep (Learning) Focus 新闻通讯，在这里我通过对该领域流行论文的易懂概述，帮助读者深入理解深度学习研究中的主题。

参考文献

[1] Mildenhall, Ben 等人。“Nerf：将场景表示为用于视图合成的神经辐射场。” ACM 通讯 65.1 (2021): 99–106。

[2] Park, Jeong Joon 等人。“Deepsdf：学习用于形状表示的连续符号距离函数。” IEEE/CVF 计算机视觉与模式识别会议论文集。2019 年。

[3] Mescheder, Lars 等人。“占用网络：在函数空间中学习 3d 重建。” IEEE/CVF 计算机视觉与模式识别会议论文集。2019 年。

[4] Sitzmann, Vincent、Michael Zollhöfer 和 Gordon Wetzstein。“场景表示网络：连续的 3d 结构感知神经场景表示。” 神经信息处理系统进展 32 (2019)。

[5] Rahaman, Nasim 等人。“神经网络的谱偏差。” 国际机器学习会议。PMLR，2019 年。

[6] Mildenhall, Ben 等人。“局部光场融合：带有规定性采样指南的实用视图合成。” ACM 图形学交易 (TOG) 38.4 (2019): 1–14。

理解机器学习中的噪声数据和不确定性

原文：towardsdatascience.com/understanding-noisy-data-and-uncertainty-in-machine-learning-4a2995a84198

你的机器学习模型无法正常工作的真正原因

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传 Harrison Hoffman

·发表于 Towards Data Science ·阅读时间：9 分钟·2023 年 1 月 23 日

–

人工智能和机器学习领域比以往任何时候都更火热。随着像 Chat GPT 和 Stable Diffusion 这样的模型席卷全球，人工智能和机器学习的炒作重新回归并吸引了大众。面对这些炒作，我们必须提醒自己，机器学习成功的关键在于高质量数据。

在缺乏优质训练数据的情况下，监督机器学习模型没有任何用处。不幸的是，大多数现实世界的数据科学项目失败是因为在数据源质量完全了解之前对模型性能有不切实际的期望。本文将尝试提供对噪声数据的直观理解以及为什么机器学习模型无法有效工作的原因。我们将探讨监督学习和确定性函数的本质、不同类型的模型不确定性，并讨论减少这种不确定性和管理期望的方法。

监督学习表述

从根本上讲，监督机器学习就是函数近似。我们向模型提供一些输入（X）和目标（y），并期望模型通过优化目标函数来近似（或学习）将 X 映射到 y 的函数。更正式地说，监督学习的表述大致如下：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

监督学习的简要表述。图像来源：作者。

给定一个数据集（X，y），我们假设 X 和 y 之间的关系至少部分是确定性的。随机误差项 epsilon 的方差在函数 f(X)的确定性水平中发挥了至关重要的作用。epsilon 的方差越高，X 和 y 之间的关系就越随机，越不容易预测。与监督学习的常见教学方式相反，值得注意的是，epsilon 可以（并且通常是）X 的函数。这意味着噪声量和相应的模型不确定性在 X 的不同区域中有所不同。

为了更好地理解确定性和随机误差的作用，我们将深入研究确定性函数的性质，检查它们的特征以及在机器学习中的影响。

确定性函数

确定性函数的输出完全由输入指定。这意味着对于定义域中的每个输入，函数都输出一个唯一的量。这些是我们从第一门代数课程开始就见过的大多数函数：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

确定性函数示例。图像由作者提供。

在上述每个示例中，传递到函数中的每个输入都产生一个单一的输出。注意，这并不意味着函数对每个可能的数字都定义。例如，函数（4）仅在 x1 > 0 且 x1 和 x2 的和不为零时定义。

从视觉上看，在一维中，确定性函数通过“垂直线测试”——这意味着如果我们在定义域中的任何点画一条垂直线，它最多会与函数相交一次：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

确定性函数的垂直线测试。图像由作者提供。

与此相比，纯非确定性函数则不同。这些函数对于单一输入可能有多个输出。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

非确定性函数的垂直线测试。图像由作者提供。

在上图中，垂直线 x = 4 穿过曲线的-2 和 2 处。如果有人问我们 f(4)是什么，我们必须说是-2 或 2。在最极端的情况下，非确定性函数可能会对单一输入输出无限多个值。

尽管上述非确定性函数在可视化和量化上很简单，但机器学习模型永远无法学习它。为了解释这一点，假设我们按如下方式对函数进行采样：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

非确定性函数的样本。图像由作者提供。

我们将使用上述样本作为训练集，并拟合一个决策树回归器：

import numpy as np
from sklearn.tree import DecisionTreeRegressor

# Create the training set for the function y = +- sqrt(x)
y_train = np.random.normal(0, 20, 100)
x_train = y_train**2

# Fit a decision tree regressor to the training data
model = DecisionTreeRegressor()
model.fit(x_train.reshape(-1, 1), y_train)

# Create the test set for the function y = +- sqrt(x)
y_test = np.random.normal(0, 20, 100)
x_test = y_test**2

# Make predictions on the test set
preds = model.predict(x_test.reshape(-1, 1))

我们可以可视化模型在测试集上的表现：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

决策树回归器在测试集上的预测叠加。图像由作者提供。

决策树在测试集上本质上是在猜测。由于领域中的每个非零输入映射到两个不同的值，机器学习模型无法知道在测试集上预测哪个值。由于基础函数的性质，训练数据中的噪声会导致任何机器学习模型做出具有高不确定性的预测。

不确定性来源

现在我们已经了解了确定性和非确定性/随机函数之间的区别，让我们探讨从学习这些函数中可能产生的不确定性类型。

数据（随机）不确定性

数据的不确定性，也称为随机不确定性，源于观察数据集的固有复杂性。在分类设置中，这通常表现为重叠的类别。例如，以下散点图描绘了二维空间中重叠的类别：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

带有重叠类别的分类数据。图片来源：作者。

对于 (x1, x2) 邻域中的大多数观测值，没有明显的方法来区分两个类别。虽然确实存在一个足够复杂的分类器可以完美分类上述数据集中的例子，但同样的分类器在测试集上的表现可能不会比随机猜测更好。

在回归设置中，数据不确定性通常源于基础数据中的附加噪声。数据不确定性会出现在类似于以下数据集的情况中：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

噪声回归数据。图片来源：作者。

在这种情况下，x 是唯一考虑的特征，对于 x 的任何邻域，y 都可能取很多不同的值。例如，如果我们基于数据来预测 x = 0 时 y 的值，没有明显的答案显现出来。一般来说，附加噪声分布的方差越大，数据不确定性就会越高。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

带有附加噪声的二次函数样本。图片来源：作者。

如前所述，噪声水平和数据不确定性在领域的不同区域中常常有所不同。数据不确定性通常在领域的边缘区域较高，如下所示：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

噪声和数据不确定性随 x 变化。附加噪声项的方差是 x 的函数。图片来源：作者。

数据的不确定性通常表现为模型中遗漏了关键特征。为了说明这一点，考虑以下数据集，我们希望使用 x1 来预测 y：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

像愤怒的青蛙一样的噪声数据。图片来源：作者。

如果仅使用 x1 作为特征来预测 y，没有任何机器学习模型能够表现良好。例如，当 x1 = -0.49 时，在这个数据集中 y 取 14 个不同的值。在没有其他特征的情况下，预测 y 看起来无望。然而，上述数据集实际上是以下函数的一个样本：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

一个方程。图像由作者提供。

如果我们从 x2 取样，并在三维中可视化 x1、x2 和 y，则会出现一个清晰的表面：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

一个三维表面。图像由作者提供。

互动表面。图像由作者提供。

x1 和 y 之间的关系噪声较大，任何试图学习这种关系的模型都会导致高数据不确定性。然而，(x1, x2) 和 y 之间的关系没有噪声，因为它完全由一个封闭形式的方程决定。只要有足够的训练数据，足够非线性的模型可以学习这种关系，从而没有数据不确定性。

许多金融应用面临数据不确定性，因为许多未观察到的人类行为会影响市场。例如，无数因素影响一个人是否会违约。即使两个贷款接收者有相同的信用档案并在特征空间中看起来相同（大多数金融机构有大量的特征），仍然有可能因为生活紧急情况而导致接收者 A 违约，而接收者 B 没有违约。除非违约模型接触到预测生活紧急情况的数据，否则这些情况下总会存在数据不确定性。

知识（认知）不确定性

知识（认知）不确定性来自于在领域的某些区域稀疏采样的数据。考虑以下数据集，我们希望建模 x 和 y 之间的关系：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

高知识不确定性。图像由作者提供。

当 x 在 1 到 3 之间或在 4 到 7 之间取值时，这个数据集没有关于 y 的信息。机器学习模型在预测这些区域的 y 值时不得不做出盲目猜测。幸运的是，与数据不确定性不同，知识不确定性可以通过采样缺失数据的区域来减少。如果我们对领域进行完全采样，x 和 y 之间的关系将变得清晰：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

通过采样消除知识不确定性。图像由作者提供。

知识不确定性在实际应用中很常见，因为大多数数据集不是从均匀分布中采样的，导致领域中的稀疏区域。一般来说，领域采样越均匀，建模过程中的知识不确定性就会越小。

处理噪声数据和不确定性的方法

现在我们对噪声数据和不确定性的性质有了一些直觉，让我们探讨一些可以采取的实际措施来应对这个问题。

1. 停止尝试寻找更好的模型

快速失败并知道何时转移是无价的数据科学技能。我们生活在一个几乎对任何机器学习问题都有最先进、高性能模型的时代。任何在模型选择或超参数调整过程中遇到次优结果的数据科学家应该寻求其他地方的性能提升。换句话说，模型可能并不需要改进，数据才需要。

2. 获取更多数据

正如我们之前所见，表现不佳的模型通常是由于数据不完整。即使是一个看似不重要的特征，在与更多特征一起使用时也可以提供显著的预测能力。因此，数据科学家必须确保他们对问题及成功建模所需的数据源有全面的理解。获取正确的数据和工程化良好的特征可以说是数据科学家工作中最困难但最关键的部分。数据科学家应该能够与利益相关者一起权衡获取更多数据的成本和收益，因为这最终是带来最大商业价值的地方。

3. 量化模型预测的不确定性

对于那些噪声数据不可避免的高风险机器学习问题，模型量化预测的不确定性通常是有价值的。经过校准的分类器提供的预测概率的不确定性估计比预测类别可能更有用。在回归中，像分位数回归、符合预测、贝叶斯神经网络和自然梯度提升等方法可以用于生成不确定性估计。简而言之，拥有一个在做出预测时能够说“我不知道”的模型可以非常有价值。

4. 管理期望

作为数据科学家，重要的是要记住机器学习是实现特定目标的工具，应该在清楚理解基本业务问题的情况下使用。全面理解问题以及机器学习在解决问题中的角色，可以更准确地设定模型性能的期望，从而最终更有效地使用这一工具。这将有助于开发出一个针对业务特定需求的机器学习模型，提供最佳结果。机器学习是一种手段。通常情况下，不需要一个完美的模型来解决问题，但正确利用不完美的模型则是必要的。

成为会员： https://harrisonfhoffman.medium.com/membership

喜欢我的文章？请给我买杯咖啡： https://www.buymeacoffee.com/HarrisonfhU

参考文献

使用 Catboost 估计不确定性 — https://catboost.ai/en/docs/references/uncertainty

理解预测性维护 — 数据采集与信号去噪

原文：towardsdatascience.com/understanding-predictive-maintenance-data-acquisition-and-signal-denoising-0309b1b103a7

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传 Marcin Stasko

·发表于 Towards Data Science ·阅读时长 10 分钟·2023 年 11 月 8 日

–

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

图片由 Michael Dziedzic 提供，来源于 Unsplash

文章目的

我想开始一系列文章，给你提供预测性维护的实际操作经验，并让你更容易入门信号处理。在这篇文章中，我们将重点关注数据采集和信号清理。如果你对某些部分感兴趣，我会考虑进一步详细讲解。在这篇文章的下一部分，我为你准备了一些实际练习。你可以使用我准备的代码进行自己的实验，通过实践来学习。

预测性维护的数据科学

数据科学中的预测性维护就像是为机器提供一种超级智能的护理方式。我们不是在机器坏了之后才进行修理，而是利用高级计算机程序和过去的数据来预测机器可能出现的问题。这有点像为机器提供了一个水晶球！通过这种方式，公司可以节省资金并延长重要机器的使用寿命。这种方法包括实时监控机器，收集数据，并使用智能计算机程序来告诉我们何时需要进行维护。因此，我们可以在问题出现之前进行修复，就像是在机器生病之前为其进行健康检查！

数据采集与处理

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

图片由 Mika Baumeister 提供，来源于 Unsplash

一切都始于数据。我们需要深入了解一些通信理论原则，如 Shannon–Hartley 定理和 Nyquist 率，以确保传感器数据的准确和高效传输。

Shannon–Hartley 定理

Shannon–Hartley 定理就像是关于通过通信频道传输多少信息的规则手册。它告诉我们频道的宽度，即频道能处理的数据量，至关重要。因此，在选择用于监控机器或传感器的设备或工具之前，我们需要确保频道足够宽，以处理我们想要的所有数据而不会丧失质量。

为了更好地利用频道，我们可以聪明地使用频道。这就像是找到最佳的方式来安排数据，以便我们能最有效地利用频道空间。这被称为优化频谱效率。因此，在选择传感器时，我们应该选择那些在可用频道空间内表现优秀的传感器。

定理应用

评估频道带宽 在选择传感器之前，仔细检查通信频道可以处理多少数据。如果频道宽度不足，考虑升级频道或寻找需要较少带宽的传感器。

优化传感器选择

选择那些在使用可用频道空间方面效率高的传感器。一些传感器可能在传输数据时更有效，不占用频道过多的空间。

Nyquist 率

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

图片由 Jair Lázaro 拍摄，来源于 Unsplash

Nyquist 率就像是我们应该用传感器收集数据的速度限制。如果我们收集数据的速度太慢，我们可能会错过关于事物运作的重要细节。因此，当我们设置系统来监控事物时，需要确保我们以能够准确展示发生情况的速度收集数据，而不会造成混淆。

我们必须确保传感器以遵循 Nyquist 率的速度收集数据。这意味着数据收集速度要足够快，以捕捉所有重要细节，而不会遗漏任何重要信息。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

Nyquist 定律

确保传感器能够以至少是最快变化（fmax) 两倍的速度（fs) 测量。这可以防止遗漏重要细节，并确保测量的准确性。选择速度与所需测量匹配的传感器。

Nyquist 的超级英雄——抗混叠滤波器

为了避免混淆或错误，我们使用像抗锯齿滤波器这样的特殊滤波器。这些滤波器帮助我们的传感器专注于信息的重要部分，并去除任何额外的噪音或混淆的细节。这就像使用放大镜来清晰地查看事物一样。因此，当我们使用传感器监控机器时，我们应该确保这些滤波器到位，以便获得最佳和最准确的信息。那么，什么是锯齿效应和抗锯齿呢？

锯齿效应

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

锯齿效应和抗锯齿的比较来源：www.zilliondesigns.com/

想象一下你正试图拍摄一辆快速移动的汽车。然而，你的相机每秒只拍摄一张照片。如果汽车移动得非常快，当你拍摄下一张照片时，它可能已经在完全不同的位置了。因此，当你查看这些照片时，汽车看起来像是处于奇怪的位置，因为你错过了它在每张照片之间的移动。这种事物看起来与实际情况不同的奇怪效果被称为锯齿效应。

抗锯齿

现在，想象你有一台神奇的相机，可以非常快速地拍摄大量照片。它不是每秒拍一张，而是在同一秒内拍摄许多张照片，捕捉汽车移动的每一个细节。当你查看这些照片时，你会看到汽车移动的平滑而准确的表现。确保捕捉到所有细节并避免锯齿效应的过程称为抗锯齿。

简单来说，抗锯齿就像拥有一台超级快速的相机，以确保事物看起来如实，无任何奇怪的失真或模糊。它有助于创建清晰而准确的图像，尤其是在事物快速移动或变化迅速时。

好的，但为什么在信号处理中这很重要？这些动画将向你展示我们如何使用它来去噪实际信号，并捕捉这些信号的“核心”。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

信号处理中的工作原理来源：siemens.com

实践经验

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

shraga kopstein 的照片，来源于 Unsplash

让我们深入代码，动手实践吧！首先，我们将生成一个示例信号，并设置一个绘图函数用于可视化。在这个例子中，我们将使用基本的正弦波，并尝试调整其参数。

import numpy as np
import matplotlib.pyplot as plt

def generate_signal(frequency, duration, sampling_rate):
    t = np.linspace(0, duration, int(sampling_rate * duration), endpoint=False)
    signal = np.sin(2 * np.pi * frequency * t)
    return t, signal

def plot_signals(t, original_signal, filtered_signal, title):
    plt.figure(figsize=(10, 6))
    plt.plot(t, original_signal, label='Original Signal', linewidth=2)
    plt.plot(t, filtered_signal, label='Filtered Signal', linestyle='dashed', linewidth=2)
    plt.title(title)
    plt.xlabel('Time (s)')
    plt.ylabel('Amplitude')
    plt.legend()
    plt.grid(True)
    plt.show()

接下来，我们将创建一个使用简单滤波器（如巴特沃斯滤波器）进行信号抗锯齿的示例。

将巴特沃斯滤波器视为信号中不同频率的门控器。它允许低频通过，同时减少高频。调整如滤波器阶数和截止频率等参数可以微调其选择性。这有助于我们防止诸如混叠等问题，确保信号的更清晰表示。让我们在代码中查看实际效果！

from scipy.signal import butter, lfilter

def apply_antialiasing(signal, cutoff_frequency, sampling_rate, order=4):
    # Design a low-pass Butterworth filter (maximal flat magnitute)
    nyquist = 0.5 * sampling_rate
    # Nyquist law in practice
    normal_cutoff = cutoff_frequency / nyquist
    b, a = butter(order, normal_cutoff, btype='low', analog=False)

    # Apply the filter to the signal
    filtered_signal = lfilter(b, a, signal)
    return filtered_signal

让我们进行一次实验

 frequency = 30.0          # Frequency of the signal
duration = 1.0            # Duration of the signal in seconds
sampling_rate = 100.0     # Sampling rate in Hz
cutoff_frequency = 20.0   # Cutoff frequency of the anti-aliasing filter

# Generate a signal
t, original_signal = generate_signal(frequency, duration, sampling_rate)

# Apply anti-aliasing filter
filtered_signal = apply_antialiasing(original_signal, cutoff_frequency, sampling_rate)

# Plot the original and filtered signals
plot_signals(t, original_signal, filtered_signal, 'Original and Filtered Signals')

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

这个实验的输出

当涉及到采样时，我们的纯正弦波看起来并不像我们期望的那样“美丽”。某些部分被截断了，这是由于采样过程造成的。

我鼓励你复制这段代码并尝试调整参数。这非常有趣！

我们也可以尝试不同类型的滤波器。

from scipy.signal import butter, cheby1, cheby2, ellip, lfilter
def apply_filter(signal, cutoff_frequency, sampling_rate, filter_type='butter', order=4):
    nyquist = 0.5 * sampling_rate
    normal_cutoff = cutoff_frequency / nyquist

    if filter_type == 'butter':
        b, a = butter(order, normal_cutoff, btype='low', analog=False)
    elif filter_type == 'cheby1':
        b, a = cheby1(order, 5, normal_cutoff, btype='low', analog=False)
    elif filter_type == 'cheby2':
        b, a = cheby2(order, 40, normal_cutoff, btype='low', analog=False)
    elif filter_type == 'ellip':
        b, a = ellip(order, 5, 40, normal_cutoff, btype='low', analog=False)

    filtered_signal = lfilter(b, a, signal)
    return filtered_signal

我们将稍微修改我们的绘图函数，以便将所有结果汇集在一个图中。

 def plot_signals_subplots(t, original_signal, filtered_signals, titles):
    num_filters = len(filtered_signals)

    fig, axes = plt.subplots(num_filters + 1, 1, figsize=(12, 2 * (num_filters + 1)))

    # Plot original signal
    axes[0].plot(t, original_signal, label='Original Signal', linewidth=2, alpha=0.7)
    axes[0].set_title('Original Signal')
    axes[0].set_xlabel('Time (s)')
    axes[0].set_ylabel('Amplitude')
    axes[0].legend()
    axes[0].grid(True)

    # Plot filtered signals
    for i, (filtered_signal, filter_type) in enumerate(zip(filtered_signals, filter_types)):
        label = f'Filtered Signal ({filter_type})'
        axes[i + 1].plot(t, filtered_signal, label=label, linestyle='dashed', linewidth=2)
        axes[i + 1].set_title(f'Filtered Signal ({filter_type})')
        axes[i + 1].set_xlabel('Time (s)')
        axes[i + 1].set_ylabel('Amplitude')
        axes[i + 1].legend()
        axes[i + 1].grid(True)

    plt.tight_layout()
    plt.show()

# Generate a signal (exactly the same as previous)
t, original_signal = generate_signal(frequency, duration, sampling_rate)

# Apply different filters
filter_types = ['butter', 'cheby1', 'cheby2', 'ellip']
filtered_signals = [apply_filter(original_signal, cutoff_frequency, sampling_rate, f) for f in filter_types]

# Plot the original and filtered signals
plot_signals(t, original_signal, filtered_signals, 'Original and Filtered Signals')

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

在美丽的纯正弦波中，某些效果可能不容易看出。让我们加入一些噪声。

def generate_signal_with_noise(frequency, duration, sampling_rate, noise_amplitude=0.1):
    t = np.linspace(0, duration, int(sampling_rate * duration), endpoint=False)
    signal = np.sin(2 * np.pi * frequency * t)

    # Add noise to the signal
    noise = noise_amplitude * np.random.normal(size=len(signal))
    signal_with_noise = signal + noise

    return t, signal_with_noise

再次运行实验。

# Generate a noisy signal
t, original_signal = generate_signal(frequency, duration, sampling_rate)
t, original_signal_with_noise = generate_signal_with_noise(frequency, duration, sampling_rate, noise_amplitude=0.8)
# Apply different filters
filter_types = ['butter', 'cheby1', 'cheby2', 'ellip']
filtered_signals = [apply_filter(original_signal_with_noise, cutoff_frequency, sampling_rate, f) for f in filter_types]

# Plot the original and filtered signals
plot_signals_subplots(t, original_signal_with_noise, filtered_signals, 'Original and Filtered Signals')

现在，我们可以观察去噪的效果，并了解各种滤波器的工作原理。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

在这个实验中，我们精确识别了我们生成的函数（一个纯正弦波）。在研究实际数据时，我们旨在确定描述现象的函数，例如故障信号或持续增加的设备磨损。噪声在我们的数据中总是存在，理解如何处理它至关重要。

季节性分解

季节性分解是信号处理中的一种强大技术，用于理解信号的不同组成部分并研究其季节性模式。在这个背景下，信号可以表示各种现象，如经济数据、环境变量甚至电信号。季节性分解的目标是将复杂信号拆解为其基本部分，使分析和解释变得更容易。

一个信号通常由三个主要组件组成：趋势、季节性成分和残差（或噪声）。趋势代表信号的长期行为或整体模式，季节性成分捕捉在规律间隔内发生的重复模式，而残差是无法用趋势或季节性解释的剩余变异或噪声。

通过将信号分解为这些组件，分析师可以洞察潜在的模式和趋势，帮助他们做出更明智的决策。这个过程在金融、气候科学和制造业等领域尤为重要，因为理解季节性变化对于准确预测和有效决策至关重要。

在这次季节性分解的探索中，我们将深入探讨信号分解的方法，逐一检查每个组件，并学习如何调查和解释从这种分析方法中获得的见解。

from statsmodels.tsa.seasonal import seasonal_decompose    

# Generate a signal with noise
t, original_signal_with_noise = generate_signal_with_noise(frequency, duration, sampling_rate, noise_amplitude=0.8)

# Perform seasonal decomposition
decomposition = seasonal_decompose(original_signal_with_noise, period=25)  # Adjust the period as needed

# Get the trend, seasonal, and residual components
trend = decomposition.trend
seasonal = decomposition.seasonal
residual = decomposition.resid

# Plot signals in separate subplots
plot_signals_subplots(t, original_signal_with_noise, [trend, seasonal, residual], 'Original Signal with Noise and Decomposed Components')

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

通过添加有趣的信号来使我们的正弦波更有趣。我将调整信号生成以模拟设备的磨损，想象它像机器的零部件随着时间的推移逐渐磨损。这可能会导致由于摩擦增加而产生更多的振动。

def generate_signal_with_wear_and_noise(frequency, duration, sampling_rate, wear_slope=0.02, noise_amplitude=0.1):
    t = np.linspace(0, duration, int(sampling_rate * duration), endpoint=False)
    signal = np.sin(2 * np.pi * frequency * t)

    # Simulate equipment wearing with a linear trend
    wear = wear_slope * t
    signal_with_wear = signal + wear

    # Add noise to the signal
    noise = noise_amplitude * np.random.normal(size=len(signal))
    signal_with_wear_and_noise = signal_with_wear + noise

    return t, signal_with_wear_and_noise

开始实验。

# Generate a signal with equipment wearing and noise
    t, original_signal_with_wear_and_noise = generate_signal_with_wear_and_noise(
        frequency, duration, sampling_rate, wear_slope=0.5, noise_amplitude=0.1
    )

    # Perform seasonal decomposition
    decomposition = seasonal_decompose(original_signal_with_wear_and_noise, period=25)  # Adjust the period as needed

    # Get the trend, seasonal, and residual components
    trend = decomposition.trend
    seasonal = decomposition.seasonal
    residual = decomposition.resid

    # Plot signals in separate subplots with component names
    plot_signals_subplots(t, original_signal_with_wear_and_noise, [trend, seasonal, residual], 'Signal with Wearing and Noise, and Decomposed Components')

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

哇！在我们原始的噪声信号中，我们找不到任何东西。然而，在将其分解成部分后，我们发现了趋势组件，现在已经清晰可见。这是一个信号。

接下来做什么？

在下一篇文章中，我们将探讨为什么检查平稳性很重要，什么是单位根，以及进行特征工程。跟随我，不要错过哦 😃

理解预测性维护 — 单位根和稳态

原文：towardsdatascience.com/understanding-predictive-maintenance-unit-roots-and-stationarity-f05322f7b6df

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传 Marcin Stasko

·发表于Towards Data Science ·阅读时间 13 分钟·2023 年 11 月 13 日

–

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

照片由艾玛·戈塞特拍摄，来源于Unsplash

文章目的

在这篇文章中，我们将深入探讨单位根和稳态的关键概念。请做好准备，我们将探讨为何检查稳态至关重要，什么是单位根，以及这些元素如何在我们的预测维护工具箱中发挥重要作用。我们还将掌握混沌理论！

这篇文章是“理解预测性维护”系列的一部分。我计划以类似风格创建整个系列。

点击此链接查看整个系列。请关注我，以确保不会错过新的文章。

数据稳态 — 捉迷藏的分析游戏

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

照片由米切尔·洛拍摄，来源于Unsplash

是否曾经想过你的数据是否在玩捉迷藏？直接切入主题——我们在谈论稳态。它不仅仅是一个花哨的术语，它是理解你的时间依赖数据究竟有多稳定和可预测的秘密武器。请做好准备，我们将探讨为何数据稳态在建模和预测中是游戏规则的改变者。

稳态的关键规则

常数均值：一个平稳的时间序列应表现出一致的平均值。如果均值发生变化，可能表明过程的基本行为发生了变化。
常数方差：时间序列的方差，代表数据点的分布，应该保持不变。方差的波动会使得准确预测变得困难。
常数自相关：自相关测量时间序列与其滞后值之间的相关性。在平稳序列中，自相关的强度和模式应该保持一致。

仅仅是统计性质的“稳定性”。

为什么平稳性如此重要

想象一下你的预测模型就像是专家导航员在数据的海洋中航行。为了顺利导航，它们更喜欢平静的水域——这就是平稳性的作用。平稳数据就像是一片宁静的海洋，模式保持一致。但是，如果你的数据是一片风雨交加、波涛汹涌的海洋（非平稳），准确预测就会变得非常困难。这就是为什么我们需要发现这些风暴，并将数据转变成一个平静的池塘，以便进行有效的时间序列分析。

现实世界的影响

数据的平稳性不仅仅是技术问题；它无处不在，影响着从金融到天气预测的决策。在金融领域，精确性对于风险和回报的估计至关重要，假设平稳性就像拥有一个可靠的指南针。气候科学家依赖平稳模型来预测长期天气模式——这就像拥有一个可靠的天气应用程序来预测地球的未来。

迈向深刻分析的旅程

使我们的数据平稳不仅仅是一个技术任务；它是迈向清晰的冒险。这就像将混乱的宝藏地图转变成清晰的指南，帮助分析师和决策者理解一切。在动态的时间依赖数据世界中，平稳性成为我们可靠的地图，引导我们理解表面下的模式，使我们在数据的旅程中更加顺利。

好了，现在我们了解了为什么平静的数据很酷，接下来让我们学习如何让它平静下来。但是等等，在我们动手写代码之前，让我给你介绍一下“单位根”。把它们当作影响我们数据行为的特殊成分。了解单位根就像拥有一个秘密配方，将我们波动的、混乱的数据变成一个平滑的池塘，准备好让我们潜入其中探索。所以，准备好迎接我们旅程的下一部分吧！

单位根——数据历史书中的顽皮时间旅行者

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

Andy Beales 拍摄于 Unsplash

单位根是时间序列分析中的基本概念，在理解现实世界数据的行为和特征中发挥着关键作用。在这次探讨中，我们将深入了解单位根是什么，它们在真实数据分析中为何重要，以及它们如何影响预测性维护领域。当然，我们会在动手实验部分做一些实验。

什么是单位根？

时间序列变量中的单位根意味着一个随机过程，其中变量在任何给定时间的值受其过去值的影响。形式上，单位根表明非平稳性，表明该变量不会随着时间的推移回到一个恒定的均值。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

单位根的数学解释

单位根的存在将持久性引入时间序列，导致建模和预测中的挑战。增强型迪基-富勒（ADF）测试及其他统计方法用于检测单位根的存在，提供了非平稳性的定量度量。

单位根就像数据的讲述者，编织出超越个体时刻的叙事，创造出连续的故事情节。它们标志着历史影响的持续性，将记忆的元素引入我们数据集的数字结构中。

想象你的数据集是一部历史小说，每个数据点代表着正在展开的故事中的一章。在这个背景下，单位根就是那些反复出现的主题和角色，它们在叙事中留下了不可磨灭的印记，以微妙而一致的影响引导情节发展。

为什么这对我们很重要？

理解单位根对时间序列分析师和建模师来说是基础性的。非平稳数据带来了挑战，因为传统模型通常假设数据是平稳的，以便进行准确的预测。分析师必须通过采用如差分等变换来处理单位根，以诱导平稳性并促进模型开发。

在预测维护场景中，单位根在确保预测模型准确性方面发挥着至关重要的作用。嵌入单位根的长期影响可能显著影响预测的可靠性，因此其识别和缓解对于有效的维护策略至关重要。

在我们进行这次技术探索时，我们将深入探讨单位根测试方法，解读结果，并探索处理非平稳时间序列数据的策略。单位根的理论基础为我们在分析旅程中的实际应用提供了坚实的基础。

增强型迪基-富勒（ADF）帮助我们

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

由 MD_JERRY 提供的照片，来源于 Unsplash

想象你有一行蚂蚁朝某个方向移动。ADF 测试检查这些蚂蚁是否有目的地行进（平稳）或是随机散布在各处（非平稳）。

ADF 测试涉及一些数学运算，但我们可以简化它：

零假设 (*H0*)：这就像默认假设。ADF 的零假设是数据有单位根，这意味着它是非平稳的。就像说蚂蚁在随机游走。

H0: 数据具有单位根（非平稳）
替代假设 (*H1*): 这是我们试图证明的假设。替代假设是数据是平稳的，就像蚂蚁在清晰的直线上行走一样。

H1: 数据是平稳的
测试统计量: ADF 测试计算一个称为测试统计量的数字。如果这个数字很小，表明数据可能是平稳的。

p 值: 这是一个概率分数。如果 p 值 很小（低于某个阈值，例如 0.05），我们就拒绝零假设并接受替代假设，认为我们的数据可能是平稳的。

这并不复杂，只需运行测试并检查 p 值

from statsmodels.tsa.stattools import adfuller
# Perform the Augmented Dickey-Fuller (ADF) test for stationarity
adf_statistic, adf_p_value, adf_lags,
adf_nobs, adf_critical_values, adf_reg_results = adfuller(stationary_series)

# Check if the series is stationary based on the p-value
is_stationary = adf_p_value < 0.05  # Using a significance level of 0.05

你通常会像这样使用 adf:

# What youy will probably will use most of the time
_, adf_p_value, _, _, _, _= adfuller(stationary_series)

但我会解释这些变量背后的含义

adf_statistic: ADF 测试中的统计量，指示了针对非平稳性零假设的证据强度。
adf_p_value: 与零假设相关的 p 值。较低的 p 值表示更强的反对非平稳性的证据。
adf_lags: 测试中使用的滞后数。
adf_nobs: ADF 测试中使用的观察数量。
adf_critical_values: 在不同显著性水平下的测试统计量的临界值。
adf_reg_results: 回归结果，提供了关于测试中执行的线性回归的额外信息。

虽然混沌可能看起来令人生畏，但我们可以通过理解和利用其模式将其转变为我们的盟友。在数据和分析的领域，混沌可以是一种强大的力量，当正确引导时，提供洞察、预测和更清晰的前进道路。这一切都是关于将不可预测性转变为优势，让混沌成为我们探索和理解旅程中的战略伙伴。

你的随机性有多“随机”？

让我们先生成一个简单的平稳序列，但这里有个提醒：“随机”并非全都相同。随机性主要有两种——真正的随机和伪随机。你可能更常接触伪随机，因为这是计算机的首选。

在计算中，生成真正的随机数是一个挑战，因为计算机是确定性的机器。伪随机数，顾名思义，并非真正随机，而是由模拟随机性的算法生成。这些算法以一个称为种子的初始值开始，并使用它生成看起来随机的数字序列。

种子

种子是伪随机数生成中的一个关键元素。它作为算法的起点。如果你使用相同的种子，每次都会得到相同的伪随机数序列。这种确定性在你希望重现的场景中是有利的。例如，如果你运行一个涉及随机性的模拟或实验，设置种子可以让你重现确切的随机数序列。

另一方面，更改种子会导致不同的伪随机数序列。这一特性常用于引入模拟中的变异性或为使用随机性的算法提供不同的初始条件。

总之，伪随机数是由算法生成的，而种子是这些算法的起点。控制种子可以让你控制伪随机数的序列，在计算机生成的随机性中提供确定性与变异性之间的平衡。

生成我们的伪随机分布的时间到了。

import numpy as np
import matplotlib.pyplot as plt
np.random.seed(1992) # WOW this is our deterministic seed.

def generate_stationary_series_pseudorandom(size=100):
    stationary_series = np.random.randn(size)
    return stationary_series

我们可以使用真正的随机性吗？

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

图片由 aj_aaaab 提供，来源于 Unsplash

现在我们可能会惊讶于即使是我们大多数时间所接触的随机性也是确定性随机。但是我们能否创造真正的随机性，确保没有确定性在其背后？

好消息！我们可以利用真正的物理现象——大气噪声。还记得你电视屏幕上的那些闪烁的黑白点吗？那就是我们的大气噪声，我们将利用它来产生真正的随机性。所以，你的电视不仅仅是用来看节目的；它是你摆脱确定性世界的门票。

import requests

def generate_stationary_series_random(size=100):
    # Fetch truly random values from random.org atmospheric noise API
    response = requests.get(f'https://www.random.org/integers/?num={size}&min=-10000&max=10000&col=1&base=10&format=plain&rnd=new')
    if response.status_code == 200:
        stationary_series = [int(value) for value in response.text.strip().split('\n')]
        return stationary_series
    else:
        raise Exception(f"Failed to fetch random values. Status code: {response.status_code}")

使用这个函数，我们可以生成真正的随机性，万岁！

平稳性检验

首先，让我们生成序列。

# Generate series
stationary_series_pseudorandom = generate_stationary_series_pseudorandom()
stationary_series_random = generate_stationary_series_pseudorandom()

创建结果图。

titles = [
    'stationary_series_pseudorandom',
    'stationary_series_random'
]

plot_multiple_series(stationary_series_random, stationary_series_pseudorandom, 
                     titles=titles)

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

嗯，“壮观”

_, adf_p_value, _, _, _, _= adfuller(stationary_series_pseudorandom)
print(f'PseudoRandom adf p-value: {adf_p_value}')
_, adf_p_value, _, _, _, _= adfuller(stationary_series_random)
print(f'TrueRandom adf p-value: {adf_p_value}')

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

结果

当 p 值非常小（<0.05）时，它提供了反对零假设的证据，表明你的数据很可能是平稳的。

所以，在这种情况下，p 值远小于 0.05，你可以有信心地说，“是的，我们的数据是平稳的。”

现在，让我们花一点时间来处理这些数据。我们的伪随机数的 p 值大约比真正的随机数小 200 万倍。

为什么会发生这种情况？伪随机数是由算法生成的，这引入了一定程度的确定性。这些算法可能会无意中在数据中引入模式或结构。另一方面，真正的随机数据，如大气噪声，更有可能表现出纯随机性的特征。ADF 检验专注于检测指示非平稳性的模式，可能会在真正的随机数据中发现较少的此类模式，从而导致相对较高的 p 值。

实践经验

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

图片由 Eddie Kopp 提供，来源于 Unsplash

现在是时候通过代码来动手了。我们将进行一些实验，帮助你熟悉文章中的概念。我建议你复现这些实验。在我们深入探讨平稳性之前，我想问你一个问题。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

照片由 Ian Barsby 拍摄，来自 Unsplash

现在我们将添加几个示例，展示如何使这些数据变为非平稳数据，我们将打破平稳性的关键规则。解释完毕后，我们将绘制所有图表。

线性趋势（非恒定均值）

def generate_non_stationary_linear_trend(size=100):
    time = np.arange(size)
    linear_trend = 0.5 * time
    non_stationary_series = np.random.randn(size) + linear_trend
    return non_stationary_series

引入线性趋势以违反恒定均值规则意味着在时间上添加系统性的增加或减少。在非平稳的线性趋势序列中，值随着时间线性增加。这违反了恒定均值规则，因为序列的平均值在变化，表明过程的基本行为发生了变化。在这种情况下，单位根有助于线性趋势的持续性，导致变量在任何给定时间的值受到其过去值的影响。

正弦幅度（非恒定方差）

def generate_non_stationary_sin_amplitude(size=100):
    time = np.arange(size)
    amplitude = 0.5 + 0.02 * time
    sin_amplitude_component = amplitude * np.sin(2 * np.pi * time / 10)
    non_stationary_series = np.random.randn(size) + sin_amplitude_component
    return non_stationary_series

添加一个幅度逐渐增加的正弦组件违反了恒定方差规则。在非平稳的季节性组件序列中，正弦组件的幅度随时间线性增长。这导致数据点的分布波动，使得方差变得不恒定。单位根有助于季节性组件的持续性，影响方差随着幅度的变化而变化。

指数增长（非恒定自相关）

def generate_non_stationary_exponential_growth(size=100, growth_rate=0.05):
    time = np.arange(size)
    exponential_growth_component = np.exp(growth_rate * time)
    non_stationary_series = np.random.randn(size) + exponential_growth_component
    return non_stationary_series

引入指数增长模式违反了恒定自相关规则。非平稳的扩展幅度序列表现出指数增长，导致自相关模式随着值的增加而变化。单位根在时间序列中引入了持续性，导致建模和预测中的挑战。单位根的存在意味着非平稳性，表明变量随时间不会恢复到恒定均值。

开始实验

执行代码生成时间序列并绘制结果。

# Example usage
stationary_series_pseudorandom = generate_stationary_series_pseudorandom()
non_stationary_linear_trend_series = generate_non_stationary_linear_trend()
non_stationary_sin_amplitude_series = generate_non_stationary_sin_amplitude()
non_stationary_exponential_growth_series = generate_non_stationary_exponential_growth()

# Visualize the examples
plot_multiple_series(stationary_series_pseudorandom, 
                     non_stationary_linear_trend_series, 
                     non_stationary_sin_amplitude_series, 
                     non_stationary_exponential_growth_series,
                     titles=[
                         'Stationary series',
                         'Linear Trend (Non-Constant Mean)',
                         'Sinusoidal Amplitude (Non-Constant Variance)',
                         'Exponential Growth (Non-Constant Autocorrelation)'
                     ])

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

在探索性数据分析中，发现线性趋势或指数增长相对简单，因为这些模式展示了明显的视觉线索。然而，在处理正弦幅度时，区分平稳状态和非平稳状态变得具有挑战性。仅通过查看数据，很难判断幅度是平稳还是非平稳。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

这个案例将展示统计测试的力量。我们手中有强大的工具。

_, adf_p_value_stationary, _, _, _, _ = adfuller(stationary_series_pseudorandom)
_, adf_p_value_linear_trend, _, _, _, _ = adfuller(generate_non_stationary_linear_trend())
_, adf_p_value_sin_amplitude, _, _, _, _ = adfuller(generate_non_stationary_sin_amplitude())
_, adf_p_value_exponential_growth, _, _, _, _ = adfuller(generate_non_stationary_exponential_growth())

# Print the results
print(f'PseudoRandom ADF P-value (Stationary Series): {adf_p_value_stationary}')
print(f'PseudoRandom ADF P-value (Linear Trend): {adf_p_value_linear_trend}')
print(f'PseudoRandom ADF P-value (Sinusoidal Amplitude): {adf_p_value_sin_amplitude}')
print(f'PseudoRandom ADF P-value (Exponential Growth): {adf_p_value_exponential_growth}')

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

测试结果表明，只有平稳序列在 ADF 测试中是平稳的。

ADF 测试在平稳和非平稳时间序列之间提供了明确的区分。在第一种情况下，我们可以自信地拒绝原假设，表明时间序列是平稳的。然而，对于其他情况，我们必须接受原假设，得出数据是非平稳的结论。具体来说，在正弦幅度的情况下，尽管非平稳性在视觉上显而易见，但 ADF 测试通过不允许我们拒绝原假设来确认我们的观察结果。

实践变换

现在，让我们玩一玩变换，尝试将我们的非平稳时间序列转变为平稳序列——就像做一点反向工程。在现实场景中，确定所需的确切变换通常是一个试错过程。我建议进行探索性数据分析，绘制时间序列，并进行经验尝试。如果某个变换使序列平稳，你不仅达到了平稳性，还可以获得有关数据特征的宝贵见解。

def make_linear_trend_stationary(series):
    # Subtract the linear trend to make the mean constant.
    time = np.arange(len(series))
    linear_trend = 0.5 * time # Somehow we have found this trend :)
    stationary_series = series - linear_trend
    return stationary_series

def make_sin_amplitude_stationary(series):
    # Apply differencing to stabilize and make the variance constant.
    diff_series = np.diff(series)
    return diff_series

def make_exponential_growth_stationary(series, epsilon=1e-8):
    # Add a small constant to avoid zero or negative values
    series = np.where(series <= 0, epsilon, series)

    # Add a small constant to avoid non-finite values
    series += epsilon

    # Apply the log for stabilization
    series = np.log(series)

    # Take the first difference to remove the exponential growth
    stationary_series = np.diff(series)

    return stationary_series

在定义了我们的变换函数后，是时候将它们付诸实践了。让我们将这些变换应用于我们的非平稳时间序列，看看是否能够成功引入平稳性。

# Apply transformations to make non-stationary examples stationary
stationary_linear_trend = make_linear_trend_stationary(generate_non_stationary_linear_trend())
stationary_sin_amplitude = make_sin_amplitude_stationary(generate_non_stationary_sin_amplitude())
stationary_exponential_growth = make_exponential_growth_stationary(generate_non_stationary_exponential_growth())

# Perform ADF test for the transformed series
adf_p_value_stationary_linear_trend = adfuller(stationary_linear_trend)[1]
adf_p_value_stationary_sin_amplitude = adfuller(stationary_sin_amplitude)[1]
adf_p_value_stationary_exponential_growth = adfuller(stationary_exponential_growth)[1]

# Print the results
print(f'ADF P-value (Stationary Linear Trend): {adf_p_value_stationary_linear_trend}')
print(f'ADF P-value (Stationary Sinusoidal Amplitude): {adf_p_value_stationary_sin_amplitude}')
print(f'ADF P-value (Stationary Exponential Growth): {adf_p_value_stationary_exponential_growth}')

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

现在我的数据平稳了，太棒了！

数据的样子如下：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

好消息！由于我们的数据现在是平稳的，我们可以在每种情况下自信地拒绝原假设。现在，为了增加一点趣味，我将接受用给定种子反向工程你的随机生成迭代的挑战。让我们看看我能否揭开这个谜团！😄

在这个链接中查看整个系列。通过关注我，确保你不会错过新文章。

理解预测性维护——波数据：特征工程（第一部分）

原文：towardsdatascience.com/understanding-predictive-maintenance-wave-data-feature-engineering-part-1-fd22b6a19323

开始学习波数据信号处理所需的所有信息

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传 Marcin Stasko

·发布于 Towards Data Science ·阅读时间 16 分钟·2023 年 11 月 21 日

–

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

照片由 Lukas Tennie 提供，发布在 Unsplash 上。

文章目的

我们即将深入探讨一些有趣的内容——波数据信号处理。这在预测性维护中非常重要，但在其他领域也是如此。我将在这一系列中逐步讲解，使其易于理解。如果你有任何想法，请随时分享！

本文是《理解预测性维护》系列的一部分。

查看整个系列的链接。通过关注我，确保你不会错过新的文章。所有没有说明文字的图片均由我创作。

为什么我们的分析在时域中进行很重要？

信号处理中的时域分析是一种方法，重点关注信号在时间上的行为和特征。与频域分析不同，后者探讨信号成分的频率内容，时域分析提供了信号在不同时间间隔内变化的洞察。这种方法使我们能够观察信号表现出的变化、模式和趋势，提供了关于系统或过程的动态和时间方面的宝贵信息。

为什么在预测性维护中它如此重要？

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

照片由 James Lewis 提供，发布在 Unsplash 上。

通过将这一分析技术应用于设备数据，维护专业人员可以识别和分析机械性能中的时间模式。监测变化有助于及早发现异常或偏离预期行为，从而及时干预以解决潜在问题，防止其升级。这种前瞻性的维护方法提高了设备的可靠性，减少了停机时间，并最终有助于更加经济高效的运营过程。

理解信号的时间域特征使得行业能够超越被动维护实践。通过时间域分析获得的见解支持预测性维护，使组织能够根据设备的实际状态安排维护活动，而不是基于任意的时间间隔。这不仅优化了资源利用，还延长了机械的使用寿命，从而带来显著的成本节省和改进的整体运营性能。

振动数据——预测性维护的核心

理解振动数据在预测性维护中至关重要，原因有几个。首先，异常振动通常是机械故障的早期指示。通过持续监测和分析振动数据，维护团队可以在故障发生之前检测到异常。其次，振动分析提供了潜在问题的具体性质的见解，从而允许进行有针对性和及时的干预。最后，通过利用振动数据，预测性维护策略可以摆脱基于时间的例行维护，转向更高效的基于状态的维护方法，从而优化设备性能并减少停机时间。

特征工程理论

在实践部分，我将为每个特征提供代码示例，并附上解释以便于实际应用。让我们探索特征工程背后的理论及其在振动数据分析中的应用。

时间域特征

在这一类别中，我们为每个振动信号计算统计度量，如均值、标准差、偏度和峰度。此外，我们还深入探讨诸如均方根（RMS）和峰值因子等指标，以提供信号能量和峰值特征的整体度量。

分布统计度量 计算每个振动信号的统计度量，如均值、标准差、偏度和峰度。
RMS（均方根） 提供信号整体能量的度量。
峰值因子 峰值与RMS值的比率。

频域特征

转换到频域时，我们采用如快速傅里叶变换 (FFT)的技术来转换时域信号。提取的特征包括主频率、谱熵和谱峭度。功率谱密度 (PSD)提供了关于功率分布和谐波关系的见解。

FFT (快速傅里叶变换) 将时域信号转换为频域。从结果频谱中提取特征，如主频率、谱熵和谱峭度。
功率谱密度 (PSD) 描述信号的功率如何在频率上分布。

时频特征

探索时频域涉及技术，如小波变换和短时傅里叶变换 (STFT)，提供信号的动态表示，并捕捉频率内容随时间的变化。

小波变换 提供信号的时频表示。从小波系数中提取特征。
短时傅里叶变换 (STFT) 表示信号频率内容如何随时间变化。

包络分析

解调技术，如Hilbert 变换或小波变换，用于提取信号包络。分析包络内的特征增加了另一层理解。

解调 使用Hilbert 变换或小波变换提取信号的包络。分析包络的特征。

窗口上的统计测量

滚动统计，通过固定大小的窗口计算，允许捕捉趋势和模式。此外，高阶统计矩在窗口上，称为波形矩，提供了宝贵的见解。

滚动统计 在固定大小的窗口上计算统计量，捕捉趋势和模式。
波形矩 高阶统计矩在窗口上。

重现图

深入研究重现图并利用重现定量分析 (RQA)可以辨别数据结构中的模式，为振动信号提供独特的视角。

重现定量分析 (RQA) 分析重现图的结构以捕捉数据中的模式。

特定领域特征

特定领域特征，如峰值特征和形状特征，旨在识别和分析振动信号中的峰值和整体波形形状。

峰值特征 识别和分析振动信号中的峰值。
形状特征 提取与信号波形形状相关的特征。

尽管这些示例并未涵盖所有可能性，但其中一些可能对您的需求有用。 😃

实践经验

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

图片由Amauri Mejía提供，来自Unsplash

现在是时候通过代码亲自操作了。我们将进行一些实验，帮助你熟悉文章中的概念。我建议你复现这些实验。

为实验创建信号

我们需要模拟振动信号并增加更多现实感，以复现设备磨损

def generate_vibration_signal(duration, sampling_rate, frequency, amplitude, noise_level, max_wear, wear_threshold):
    t = np.linspace(0, duration, int(sampling_rate * duration), endpoint=False)

    # Generate a sinusoidal signal
    signal = amplitude * np.sin(2 * np.pi * frequency * t)

    # Add random noise to simulate real-world conditions
    noise = np.random.normal(0, noise_level, signal.shape)
    signal_with_noise = signal + noise

    # Simulate equipment wear
    wear = np.linspace(0, max_wear, len(t))
    wear[wear > wear_threshold] = 0  # Reset wear if it exceeds the threshold
    signal_with_wear = signal_with_noise + wear

    return t, signal_with_wear

在这段代码中，磨损在达到特定值后会重置——模拟设备更换

让我们生成信号并绘图

# Parameters
duration = 20         # seconds
sampling_rate = 20    # Hz
frequency = 5         # Hz (vibration frequency)
amplitude = 1.0       # Min Max range
noise_level = 0.3     # Noise factor to increase reality
max_wear = 1          # Maximum wear before reset
wear_threshold = 0.5  # Wear threshold for reset

# Generate synthetic vibration signal with wear and threshold
time, vibration_signal = generate_vibration_signal(duration, sampling_rate, frequency, amplitude, noise_level, max_wear, wear_threshold)

# Plot the signal
plt.plot(time, vibration_signal)
plt.title('Synthetic Vibration Signal with Equipment Wear')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
plt.show()

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

带磨损效应的噪声信号用于下一步实验

本文的任务是向你介绍这些酷炫的功能。我们在这里不会构建整个过程——今天没有管道给我们。那是另一个文章的故事！现在，让我们深入探索从信号中创建特征的有趣世界。准备好迎接特征盛宴了吗？出发吧！🚀

窗口化还是不窗口化，这是个问题。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

由我使用 DiscoDiffusion 模型生成

时间序列窗口化就像是在连续时间线中查看快照，这非常有用，特别是在预测维护中。想象一下你在看一部电影，但不是看完整部电影，而是每隔几分钟暂停一次并拍一张照片。这些照片就是你的“窗口”。这些窗口帮助我们理解事物如何随时间变化。在机器和设备的世界中，了解它们过去的行为有助于我们预测未来可能发生的情况。

使用这些窗口的一个大优点是它们使理解发生了什么变得更容易。这就像将一个大故事拆分成更小的章节。每个窗口就是一个章节，通过查看它们，我们可以发现该时间段内发生的任何奇怪或有趣的事情。这种详细的视角帮助我们找出机器可能出现磨损或故障的原因。此外，这些窗口帮助我们处理信息获取频率的变化，并处理数据中的任何异常，确保我们的预测是可靠的。

但当然，这并非全是阳光和彩虹。选择这些窗口的正确大小有点棘手。如果它们太大或太小，我们可能会错过重要细节或添加不必要的噪声。这就像为相机选择正确的镜头——你想捕捉到恰到好处的量。此外，决定这些窗口是否应该重叠也是一个难题。重叠的窗口提供更多的上下文，但重叠过多可能会使数据变得重复。这就像试图平衡在书籍每一章中包含多少背景故事。找到这个甜蜜点对确保我们关于机器维护的预测准确无误至关重要。

窗口化示例

df_windowed = pd.DataFrame({'time': time, 'vibration_signal': vibration_signal})

# Make some experiments
window_size = int(2)  

# Apply mean windowing using the 'rolling' function
df_windowed['mean_amplitude'] = df_windowed['vibration_signal'].rolling(window=window_size, min_periods=1).mean()

# Plot the original signal and the mean windowed signal
plt.plot(df_windowed['time'], df_windowed['vibration_signal'], label='Original Signal')
plt.plot(df_windowed['time'], df_windowed['mean_amplitude'], label=f'Mean Window ({window_size} samples)')
plt.title('Synthetic Vibration Signal with Mean Windowing')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
plt.show()

我打算举办一个不同大小的窗口派对，向你展示它们如何改变事物。这就像 Mean Windowing 是我们的酷 DJ，在窗口范围内旋转平均值。让我们看看数据舞池如何跟随不同窗口大小的节拍！

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

平均窗口大小 = 2

窗口大小为 2 时，很难看到任何清晰的模式；这引入了太多噪音。我们需要增加窗口大小，以更好地了解数据中发生的情况。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

平均窗口大小 = 200

现在，窗口太大了，这不好，因为我们丢失了很多数据细节。我们需要一个合适的窗口，以捕捉所有重要信息。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

平均窗口大小 = 20

窗口大小为 20 时，数据模式变得非常明显，这让我们能够识别在信号生成过程中引入的合成“磨损效应”。在训练模型时，进行试错以找出最佳窗口大小是至关重要的。本文中，我将使用窗口大小 20 来生成特征。

时间域特征

分布统计度量

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

分布统计（Latex 编译）

让我们玩一下分布测量。当然，我可能会跳过mean和standard deviation，因为它们很明显，但我尝试让这个文章有点趣味，以便深入解释每个内容。

均值

想象一下你和你的朋友们在吃披萨派对上。每个人都喜欢各种配料的披萨。这个mean就像是计算每个人披萨上意大利辣肠片的平均数量。如果一个朋友有很多，而另一个朋友只有几个，平均值帮助你知道每个人的意大利辣肠片数量。就像找到披萨的和谐！

标准差

现在，让我们谈谈standard deviation。想象一群猫。一些猫非常放松和懒惰，而另一些猫则非常活跃和好动。这个standard deviation就像是测量每只猫的能量水平如何偏离或不同于所有猫的平均能量水平。如果标准差很高，那你就有一群既懒又活跃的猫。如果标准差很低，大多数猫的能量水平相似——也许是一个悠闲的猫派对！

偏度

让我们用水果篮的场景来更清楚地理解positive和negative skewness之间的区别。

Positive Skewness (向右偏斜) 想象你的朋友们正在填充一个水果篮。大多数朋友决定添加各种水果，但有几个朋友特别热衷，添加了额外的香蕉、苹果和橙子。由于这种额外的水果热情，水果篮的秋千向右倾斜。这就是positive skewness，表示右侧有更多的兴奋。
Negative Skewness (Light to the Left) 现在，假设几个朋友决定保持轻松，只在篮子里添加几颗葡萄和浆果。这种轻盈的水果方式使跷跷板向左倾斜。这就是 negative skewness，表示向左方向的轻微偏斜。

峰度

现在，想象你在坐过山车。一些过山车很疯狂，充满了曲折，而其他的则比较温和。Kurtosis 是我们过山车的评论员，评估车程的刺激程度。Positive kurtosis 意味着过山车有急转弯和意外的回旋，而 negative kurtosis 表示平稳、更温和的骑行。Kurtosis 就是我们统计主题公园的刺激因子！

偏度比较

import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import skew, kurtosis

# Set a random seed for reproducibility
np.random.seed(1992)

# Generate synthetic datasets with varying skewness and kurtosis
# Normal distribution
normal_data = np.random.normal(loc=170, scale=5, size=1000)

# Positively skewed distribution                
skewed_data = np.random.gamma(shape=2, scale=5, size=1000)

# Negatively skewed distribution                 
negative_skewed_data = -np.random.gamma(shape=2, scale=5, size=1000)      

# Calculate mean and median for each dataset
normal_mean, normal_median = np.mean(normal_data), np.median(normal_data)
skewed_mean, skewed_median = np.mean(skewed_data), np.median(skewed_data)
negative_skewed_mean, negative_skewed_median = np.mean(negative_skewed_data), np.median(negative_skewed_data)

# Plot the distributions
plt.figure(figsize=(12, 6))

plt.subplot(1, 3, 1)
plt.hist(normal_data, bins=30, color='blue', alpha=0.7)
plt.axvline(x=normal_mean, color='red', linestyle='--', label=f'Mean: {normal_mean:.2f}')
plt.axvline(x=normal_median, color='green', linestyle='--', label=f'Median: {normal_median:.2f}')
plt.legend()
plt.title('Normal Distribution')

plt.subplot(1, 3, 2)
plt.hist(skewed_data, bins=30, color='orange', alpha=0.7)
plt.axvline(x=skewed_mean, color='red', linestyle='--', label=f'Mean: {skewed_mean:.2f}')
plt.axvline(x=skewed_median, color='green', linestyle='--', label=f'Median: {skewed_median:.2f}')
plt.legend()
plt.title('Positively Skewed Distribution')

plt.subplot(1, 3, 3)
plt.hist(negative_skewed_data, bins=30, color='green', alpha=0.7)
plt.axvline(x=negative_skewed_mean, color='red', linestyle='--', label=f'Mean: {negative_skewed_mean:.2f}')
plt.axvline(x=negative_skewed_median, color='green', linestyle='--', label=f'Median: {negative_skewed_median:.2f}')
plt.legend()
plt.title('Negatively Skewed Distribution')

plt.tight_layout()
plt.show()

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

偏度比较图与均值和中位数

Normal 分布的数据在两侧均匀分布，均值与中位数紧密对齐。

Positive / Right-skewed 分布在右侧有较长或较胖的尾部，表示左侧的数据点更多。均值通常大于中位数。

Negative / Left-skewed 分布在左侧有较长或较胖的尾部，表明右侧的数据点更多。均值通常小于中位数。

峰度比较

# Leptokurtic distribution (heavier tails)
heavy_tails_data = np.random.exponential(scale=10, size=1000)

# Platykurtic distribution (lighter tails)
light_tails_data = np.random.uniform(low=160, high=180, size=1000)        

# Calculate mean and median for each dataset
normal_mean, normal_median = np.mean(normal_data), np.median(normal_data)
heavy_tails_mean, heavy_tails_median = np.mean(heavy_tails_data), np.median(heavy_tails_data)
light_tails_mean, light_tails_median = np.mean(light_tails_data), np.median(light_tails_data)

# Plot the distributions with mean and median
plt.figure(figsize=(12, 6))

plt.subplot(1, 3, 1)
plt.hist(normal_data, bins=30, color='blue', alpha=0.7)
plt.axvline(x=normal_mean, color='red', linestyle='--', label=f'Mean: {normal_mean:.2f}')
plt.axvline(x=normal_median, color='green', linestyle='--', label=f'Median: {normal_median:.2f}')
plt.legend()
plt.title('Normal (esokurtic) Distribution')

plt.subplot(1, 3, 2)
plt.hist(heavy_tails_data, bins=30, color='red', alpha=0.7)
plt.axvline(x=heavy_tails_mean, color='red', linestyle='--', label=f'Mean: {heavy_tails_mean:.2f}')
plt.axvline(x=heavy_tails_median, color='green', linestyle='--', label=f'Median: {heavy_tails_median:.2f}')
plt.legend()
plt.title('Leptokurtic Distribution (Heavier Tails)')

plt.subplot(1, 3, 3)
plt.hist(light_tails_data, bins=30, color='green', alpha=0.7)
plt.axvline(x=light_tails_mean, color='red', linestyle='--', label=f'Mean: {light_tails_mean:.2f}')
plt.axvline(x=light_tails_median, color='green', linestyle='--', label=f'Median: {light_tails_median:.2f}')
plt.legend()
plt.title('Platykurtic Distribution (Lighter Tails)')

plt.tight_layout()
plt.show()

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

峰度图比较

现在，让我们计算统计数据进行比较。

# Calculate skewness and kurtosis for each dataset
normal_skewness = skew(normal_data)
normal_kurtosis = kurtosis(normal_data)

skewed_skewness = skew(skewed_data)
skewed_kurtosis = kurtosis(skewed_data)

negative_skewness = skew(negative_skewed_data)
negative_kurtosis = kurtosis(negative_skewed_data)

heavy_tails_skewness = skew(heavy_tails_data)
heavy_tails_kurtosis = kurtosis(heavy_tails_data)

light_tails_skewness = skew(light_tails_data)
light_tails_kurtosis = kurtosis(light_tails_data)

# Print the calculated values
print("Normal Distribution:")
print(f"Skewness: {normal_skewness}, Kurtosis: {normal_kurtosis}\n")

print("Positively Skewed Distribution:")
print(f"Skewness: {skewed_skewness}, Kurtosis: {skewed_kurtosis}\n")

print("Negatively Skewed Distribution:")
print(f"Skewness: {negative_skewness}, Kurtosis: {negative_kurtosis}\n")

print("Leptokurtic Distribution (Heavier Tails):")
print(f"Skewness: {heavy_tails_skewness}, Kurtosis: {heavy_tails_kurtosis}\n")

print("Platykurtic Distribution (Lighter Tails):")
print(f"Skewness: {light_tails_skewness}, Kurtosis: {light_tails_kurtosis}\n")

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

计算的统计结果输出

正态分布

Skewness -0.0237（略微负偏）
Kurtosis 0.1356（平峰度，比正常分布更平坦）

正偏分布

Skewness 1.3753（强烈正偏）
Kurtosis 2.7358（峰态分布，尾部比正常分布更重）

负偏分布

Skewness -1.3357（强烈负偏）
Kurtosis 2.4060（峰态分布，尾部比正常分布更重）

峰态分布（更重的尾部）

Skewness 1.8463（正偏）
Kurtosis 4.4461（高峰度，非常重的尾部）

平峰度分布（更轻的尾部）

Skewness -0.0243（略微负偏）
Kurtosis -1.1587（平峰度，尾部比正常分布轻得多）

使用滚动窗口应用统计

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

滚动窗口解释（Latex 编译）

def calculate_rolling_statistics(signal, window_size):

    df = pd.DataFrame({'signal': signal})
    rolling_stats = df['signal'].rolling(window=window_size, min_periods=1)

    mean_values = rolling_stats.mean()
    std_dev_values = rolling_stats.std()
    skewness_values = rolling_stats.apply(skew, raw=True)
    kurtosis_values = rolling_stats.apply(kurtosis, raw=True)

    return mean_values, std_dev_values, skewness_values, kurtosis_values

window_size = 20  

# Calculate rolling statistics
rolling_means, rolling_std_devs, rolling_skewness, rolling_kurtosis = calculate_rolling_statistics(vibration_signal, window_size)

我会将结果绘制出来，为了更好的可视化，我会分开 skewness 和 kurtosis

# Plot the signal and rolling statistics
plt.figure(figsize=(12, 6))

# Plot Rolling Mean, Rolling Mean + Std Dev, Rolling Mean - Std Dev
plt.subplot(2, 1, 1)
plt.plot(time[:len(rolling_means)], vibration_signal[:len(rolling_means)], label='Vibration Signal')
plt.plot(time[:len(rolling_means)], rolling_means, label='Rolling Mean')
plt.plot(time[:len(rolling_means)], rolling_means + rolling_std_devs, label='Rolling Mean + Std Dev', linestyle='--')
plt.plot(time[:len(rolling_means)], rolling_means - rolling_std_devs, label='Rolling Mean - Std Dev', linestyle='--')
plt.title(f'Synthetic Vibration Signal with Rolling Mean and Standard Deviation (Window Size = {window_size})')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()

# Plot Rolling Skewness and Rolling Kurtosis
plt.subplot(2, 1, 2)
plt.plot(time[:len(rolling_means)], vibration_signal[:len(rolling_means)], label='Vibration Signal')
plt.plot(time[:len(rolling_means)], rolling_skewness, label='Rolling Skewness', linestyle='--')
plt.plot(time[:len(rolling_means)], rolling_kurtosis, label='Rolling Kurtosis', linestyle='--')
plt.title(f'Synthetic Vibration Signal with Rolling Skewness and Kurtosis (Window Size = {window_size})')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()

plt.tight_layout()
plt.show()

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

滚动应用与原始信号

观察结果可以明显看出，滚动窗口作为一种有效的去噪技术。我建议尝试不同的窗口大小，因为数据科学通常涉及通过试验和错误进行实证探索。

幕后 **.apply()** 的故事

这里有一点额外的内容，展示函数和窗口如何在后台一起工作。我会仅仅做一次，以给你一个关于当应用函数在处理窗口时如何工作的直观了解，使用我们的简单示例。

def calculate_rolling_statistics_behind_scenes(signal, window_size):
    mean_values = np.convolve(signal, np.ones(window_size)/window_size, mode='valid')
    std_dev_values = np.array([np.std(signal[i-window_size+1:i+1]) for i in range(window_size-1, len(signal))])
    skewness_values = np.array([skew(signal[i-window_size+1:i+1]) for i in range(window_size-1, len(signal))])
    kurtosis_values = np.array([kurtosis(signal[i-window_size+1:i+1]) for i in range(window_size-1, len(signal))])

    return mean_values, std_dev_values, skewness_values, kurtosis_values

mean_values 这是通过使用 np.convolve 函数计算的，该函数执行卷积操作。在这种情况下，它通过用一个窗口与信号进行卷积来计算滚动平均值。mode='valid' 参数确保卷积仅在完整窗口可以适合而无需零填充的地方进行。
std_dev_values 这是通过使用列表推导式遍历信号来计算的。对于信号中的每个位置 i，它计算子数组 signal[i-window_size+1:i+1] 的标准差。这代表了滚动标准差。
skewness_values 类似于标准差，它是通过使用列表推导式遍历信号来计算的。对于每个位置 i，它计算子数组 signal[i-window_size+1:i+1] 的偏度。
kurtosis_values 再次，类似于标准差和偏度，它是通过使用列表推导式遍历信号来计算的。对于每个位置 i，它计算子数组 signal[i-window_size+1:i+1] 的峰度。

RMS（均方根）

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

RMS 方程（Latex 编译）

Root Mean Square (RMS) 就像是数学中的超级英雄。它接收一组数字，对每一个数字进行平方，计算平均值，然后取平方根。这个过程会给出一个单一的数字，代表原始数字的典型大小或强度。它在各个领域中都是一个方便的工具，从测量机械振动到评估信号强度。

def calculate_rolling_rms(signal, window_size):
    df = pd.DataFrame({'signal': signal})
    rolling_stats = df['signal'].rolling(window=window_size, min_periods=1)

    rms_values = np.sqrt(rolling_stats.apply(lambda x: np.mean(x**2), raw=True))

    return rms_values

window_size = 20
rolling_rms = calculate_rolling_rms(vibration_signal, window_size)

plt.plot(time[:len(rolling_rms)], vibration_signal[:len(rolling_rms)], label='Vibration Signal')
plt.plot(time[:len(rolling_rms)], rolling_rms, label='Rolling RMS', linestyle='--')
plt.title(f'Synthetic Vibration Signal with Rolling RMS (Window Size = {window_size})')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()

plt.tight_layout()
plt.show()

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

RMS 揭示了一个单一而强大的指标，捕捉了潜在的信号功率，同时有效地减少了噪声。这个图表作为一个视觉证明，展示了 RMS 在各种应用中提高信号清晰度和精度的实际效果。

峰值因子

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

峰值因子方程（Latex 编译）

Crest Factor (CF) 就像信号分析世界中的副手 Root Mean Square (RMS)。虽然 RMS 给出了整体强度，Crest Factor 则突出显示了峰值。它是最高点与 RMS 值的比率，告诉你信号的尖锐程度或“峰值”程度。可以将 CF 视为帮助你理解数据中尖锐峰值的超级英雄伙伴，无论是在声音波、电子信号还是其他波动测量中。它们一起组成了揭示数据中隐藏秘密的动态组合。

def calculate_crest_factor_and_peak(signal, window_size):
    df = pd.DataFrame({'signal': signal})
    rolling_stats = df['signal'].rolling(window=window_size, min_periods=1)

    peak_values = rolling_stats.apply(lambda x: np.max(np.abs(x)), raw=True)
    rms_values = np.sqrt(rolling_stats.apply(lambda x: np.mean(x**2), raw=True))

    crest_factor_values = peak_values / rms_values

    return crest_factor_values, peak_values

# Calculate rolling Crest Factor and Peak values
rolling_crest_factor, rolling_peak_values = calculate_crest_factor_and_peak(vibration_signal, window_size)

创建图表

# Plot the vibration signal, rolling Crest Factor, and Peak values

plt.plot(time[:len(rolling_crest_factor)], vibration_signal[:len(rolling_crest_factor)], label='Vibration Signal')
plt.plot(time[:len(rolling_crest_factor)], rolling_crest_factor, label='Rolling Crest Factor', linestyle='--')
plt.plot(time[:len(rolling_peak_values)], rolling_peak_values, label='Rolling Peak Values', linestyle='-.')
plt.title(f'Synthetic Vibration Signal with Rolling Crest Factor and Peak Values (Window Size = {window_size})')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()

plt.tight_layout()
plt.show()