使用 Sigma 规则进行异常检测：构建您自己的 Spark 流处理检测
轻松在 Spark 流处理管道中部署 Sigma 规则：一种支持即将发布的 Sigma 2 规范的未来-proof 解决方案
·
Follow 发表在 Towards Data Science · 13 分钟阅读 · 2023 年 6 月 12 日
–
由 Dana Walker 拍摄于 Unsplash
在我们之前的文章中，我们详细阐述并设计了一个名为 flux-capacitor 的有状态函数。
flux-capacitor 有状态函数可以记住日志事件之间的父子（和祖先）关系。它还可以记住在特定时间窗口内在同一主机上发生的事件，Sigma 规范称之为 时间临近相关性。
要深入了解 flux-capacitor 的设计，请参阅 第一部分，第二部分，第三部分，第四部分 和 第五部分。不过，你不需要理解功能的实现就可以使用它。
在本文中，我们首先展示一个执行离散检测的 Spark 流处理作业。离散检测是一个 Sigma 规则，它使用单个日志行（单个事件）的特征和值。
然后，我们利用 flux-capacitor 功能来处理日志事件之间的有状态父子关系。flux-capacitor 还能检测在特定时间窗口内在同一主机上发生的多个事件，这些在即将发布的 Sigma 规范中被称为 时间临近相关性。这些 Spark 流处理作业的完整演示可以在我们的 git repo 中找到。
离散检测
执行离散测试相当简单，这要归功于 Spark 中开箱即用的所有内置功能。Spark 支持读取流数据源、写入接收器、检查点、流流连接、窗口聚合等功能。有关所有可能功能的完整列表，请参阅全面的 Spark 结构化流处理编程指南。
这里是一个高层次的图示，展示了一个 Spark 流处理作业，它从“start-process”窗口事件的 Iceberg 表中消费事件（1）。一个经典的例子可以在 Windows Security Logs (事件 ID 4688) 中找到。
离散检测的拓扑结构
源表（1）名为process_telemetry_table。Spark 作业读取所有事件，检测异常事件，标记这些事件，并将其写入名为tagged_telemetry_table的表（3）。被认为异常的事件也会写入一个包含警报的表（4）。
定期我们轮询一个包含我们想要应用的 Sigma 规则自动生成的 SQL 的 git 仓库（5）。如果 SQL 语句发生变化，我们重新启动流处理作业以将这些新的检测添加到管道中。
以这个 Sigma 规则为例：
来自proc_creation_win_rundll32_sys.yml 在 Sigma HQ 的截图
detection部分是 Sigma 规则的核心，包括一个condition和一个或多个命名测试。selection1和selection2是命名的布尔测试。Sigma 规则的作者可以为这些测试命名有意义的名称。condition是用户可以在最终评估中组合测试的地方。有关编写 Sigma 规则的更多细节，请参见Sigma 规范。
从现在开始，我们将这些命名的布尔测试称为标签。
Spark 流处理作业的内部工作分为 4 个逻辑步骤：
-
读取源表
process_telemetry_table -
执行模式匹配
-
评估最终条件
-
写入结果
模式匹配步骤包括评估在 Sigma 规则中找到的标签，最终条件评估评估condition。
在该图的右侧，我们展示了在此处理阶段该行的样子。蓝色的列表示从源表中读取的值。模式匹配步骤添加了一列名为Sigma tags，这是所有执行的测试及其通过或失败情况的映射。灰色列包含最终的 Sigma 规则评估。最后，棕色列是在 foreachBatch 函数中添加的。生成了一个 GUID，从 Sigma 标签映射中提取出为真的规则名称，并从规则名到规则类型的查找映射中检索检测action。这为生成的警报提供了上下文。
此图描绘了事件的属性如何组合成标签、最终评估和最终上下文信息。
现在让我们看看实际的 pyspark 代码。首先，我们使用readStream函数将 spark 连接到源表，并指定从中读取冰山表的名称。load函数返回一个数据帧，我们用它来创建一个名为process_telemetry_view的视图。
spark
.readStream
.format("iceberg")
.option("stream-from-timestamp", ts)
.option("streaming-skip-delete-snapshots", True)
.option("streaming-skip-overwrite-snapshots", True)
.load(constants.process_telemetry_table)
.createOrReplaceTempView("process_telemetry_view")
process_telemetry_view中的数据如下所示：
+-------------------+---+---------+---------------------+
|timestamp |id |parent_id|Commandline |
+-------------------+---+---------+---------------------+
|2022-12-25 00:00:01|11 |0 | |
|2022-12-25 00:00:02|2 |0 |c:\win\notepad.exe |
|2022-12-25 00:00:03|12 |11 | |
|2022-12-25 00:00:08|201|200 |cmdline and args |
|2022-12-25 00:00:09|202|201 | |
|2022-12-25 00:00:10|203|202 |c:\test.exe |
+-------------------+---+---------+---------------------+
在此视图上，我们应用了一个 模式匹配 步骤，该步骤由 Sigma 编译器生成的 SQL 语句组成。patern_match.sql 文件如下所示：
select
*,
-- regroup each rule's tags in a map (ruleName -> Tags)
map(
'rule0',
map(
'selection1', (CommandLine LIKE '%rundll32.exe%'),
'selection2', (CommandLine LIKE '%.sys,%' OR CommandLine LIKE '%.sys %'),
)
) as sigma
from
process_telemetry_view
我们使用 spark.sql() 将此语句应用于 process_telemetry_view 视图。
df = spark.sql(render_file("pattern_match.sql"))
df.createOrReplaceTempView("pattern_match_view")
请注意，每个 Sigma 规则中找到的标签的结果存储在一个布尔值映射中。sigma 列保存了每个 Sigma 规则中找到的每个标签的结果。通过使用 MapType，我们可以轻松引入新的 Sigma 规则，而不会影响表的模式。添加新规则只是简单地在 sigma 列（一个 MapType）中增加一项新条目。
+---+---------+---------------------+----------------------------------+
|id |parent_id|Commandline |sigma
+---+---------+---------------------+----------------------------------+
|11 |0 | |{rule0 -> {
selection1 -> false,
selection2 -> false
},
}
同样，评估最终条件 步骤应用了 Sigma 规则中的条件。condition 被编译成一个 SQL 语句，该语句使用 map, map_filter, map_keys 来构建一个名为 sigma_final 的列。此列包含所有条件评估为真规则的名称。
select
*,
map_keys( -- only keep the rule names of rules that evaluted to true
map_filter( -- filter map entries keeping only rules that evaluated to true
map( -- store the result of the condition of each rule in a map
'rule0',
-- rule 0 -> condition: all of selection*
sigma.rule0.selection1 AND sigma.rule0.selection2)
)
, (k,v) -> v = TRUE)) as sigma_final
from
pattern_match_view
自动生成的语句使用 spark.sql() 应用。
df = spark.sql(render_file("eval_final_condition.sql"))
这是添加了新列 sigma_final 后的结果，列出了触发的规则数组。
+---+---------+-------------------------------------+-------------+
|id |parent_id|sigma | sigma_final |
+---+---------+-------------------------------------+-------------+
|11 |0 |{rule0 -> { | [] |
selection1 -> false,
selection2 -> false
}
}
我们现在准备开始对我们的数据帧进行流处理作业。请注意，我们将回调函数 for_each_batch_function 传递给 foreachBatch。
streaming_query = (
df
.writeStream
.queryName("detections")
.trigger(processingTime=f"{trigger} seconds")
.option("checkpointLocation", get_checkpoint_location(constants.tagged_telemetry_table) )
.foreachBatch(foreach_batch_function)
.start()
)
streaming_query.awaitTermination()
for_each_batch_function 在每个微批次中被调用，并接收评估后的 batchdf 数据帧。for_each_batch_function 将 batchdf 的全部内容写入 tagged_telementry_table，并对任何评估为真的 Sigma 规则写入警报。
def foreach_batch_function(batchdf, epoch_id):
# Transform and write batchDF
batchdf.persist()
batchdf.createOrReplaceGlobalTempView("eval_condition_view")
run("insert_into_tagged_telemetry")
run("publish_suspected_anomalies")
spark.catalog.clearCache()
详细信息请参见我们的 Git 仓库中的 insert_into_tagged_telemetry.sql 和 publish_suspected_anomalies.sql。
如上所述，使用 Spark 内置功能编写流式异常检测处理离散测试相对简单。
基于过去事件的检测
迄今为止，我们展示了如何使用离散 Sigma 规则检测事件。在本节中，我们利用 flux-capacitor 函数来启用缓存标签和测试过去事件的标签。正如我们在之前的文章中讨论的，flux-capacitor 使我们能够检测父子关系以及过去事件的仲裁特征序列。
这些类型的 Sigma 规则需要同时考虑当前事件和过去事件的标签。为了执行最终规则评估，我们引入了 Time travel tags 步骤，以检索事件的所有过去标签，并将它们与当前事件合并。这正是 flux-capacitor 函数设计的目的，它缓存和检索过去的标签。现在，过去标签和当前标签在同一行上，Eval final condition 可以像我们在上述离散示例中所做的那样进行评估。
现在检测结果如下：
flux-capacitor 接收 Pattern Match 步骤生成的 Sigma tags。flux-capacitor 存储这些标签以备后用。红色的列具有与我们之前使用的 Sigma tags 列相同的模式。然而，它结合了当前标签和从内部状态检索到的过去标签。
借助 flux-capacitor 函数，添加过去标签的缓存和检索变得容易。在我们的 Spark 异常检测中，我们是这样应用 flux-capacitor 函数的。首先，将 Pattern Match 步骤生成的数据框传递给 flux_stateful_function，该函数返回另一个数据框，其中包含过去的标签。
flux_update_spec = read_flux_update_spec()
bloom_capacity = 200000
# reference the scala code
flux_stateful_function = spark._sc._jvm.cccs.fluxcapacitor.FluxCapacitor.invoke
# group logs by host_id
jdf = flux_stateful_function(
pattern_match_df._jdf,
"host_id",
bloom_capacity,
flux_update_spec)
output_df = DataFrame(jdf, spark)
为了控制 flux_stateful_function 的行为，我们传入一个 flux_update_spec。flux-capacitor 规范是由 Sigma 编译器生成的 yaml 文件。规范详细说明了哪些标签应该被缓存和检索，以及它们应该如何处理。action 属性可以设置为 parent、ancestor 或 temporal。
让我们使用 Sigma HQ 的一个具体例子 proc_creation_win_rundll32_executable_invalid_extension.yml
来自 Sigma HQ github 的截图
再次强调，检测的核心由标签和一个最终的 condition 组成，它将所有这些标签组合在一起。然而，请注意，这条规则（我们称之为规则 1）涉及对 CommandLine 的测试，还测试了父进程 ParentImage。ParentImage 不是在启动进程日志中找到的字段。相反，它指的是父进程的 Image 字段。
如之前所见，这个 Sigma 规则将被编译成 SQL 用于评估标签，并将它们组合成最终条件。
为了传播父标签，Sigma 编译器还生成了一个 flux-capacitor 规范。规则 1 是一个 parent 规则，因此规范必须指定什么是父字段和子字段。在我们的日志中，这些对应于 id 和 parent_id。
规范还指定了哪些 tags 应该被 flux-capacitor 函数缓存和检索。以下是自动生成的规范：
rules:
- rulename: rule1
description: proc_creation_win_run_executable_invalid_extension
action: parent
tags:
- name: filter_iexplorer
- name: filter_edge_update
- name: filter_msiexec_system32
parent: parent_id
child: id
请注意规则 0 不包含在 flux-capacitor 函数中，因为它没有时间标签。
标记传播示例
为了更好地理解 flux-capacitor 的作用，你可以在流分析之外使用该函数。这里我们展示了一个简单的 祖先 示例。我们想要传播 pf 标签。例如，pf 可能表示包含 rundll32.exe 的 CommandLine。
spec = """
rules:
- rulename: rule2
action: ancestor
child: pid
parent: parent_pid
tags:
- name: pf
"""
df_input = spark.sql("""
select
*
from
values
(TIMESTAMP '2022-12-30 00:00:05', 'host1', 'pid500', '', map('rule1', map('pf', true, 'cf', false))),
(TIMESTAMP '2022-12-30 00:00:06', 'host1', 'pid600', 'pid500', map('rule1', map('pf', false, 'cf', false))),
(TIMESTAMP '2022-12-30 00:00:07', 'host1', 'pid700', 'pid600', map('rule1', map('pf', false, 'cf', true)))
t(timestamp, host_id, pid, parent_pid, sigma)
""")
打印数据框 df_input 我们看到 pid500 启动并具有带有 pf 特性的 CommandLine。然后 pid500 启动了 pid600。之后 pid600 启动了 pid700。Pid700 具有子特性 cf。
+-------------------+------+----------+--------------+-------------------------------------+
|timestamp |pid |parent_pid|human_readable|sigma |
+-------------------+------+----------+--------------+-------------------------------------+
|2022-12-30 00:00:05|pid500| |[pf] |{rule2 -> {pf -> true, cf -> false}} |
|2022-12-30 00:00:06|pid600|pid500 |[] |{rule2 -> {pf -> false, cf -> false}}|
|2022-12-30 00:00:07|pid700|pid600 |[cf] |{rule2 -> {pf -> false, cf -> true}} |
+-------------------+------+----------+--------------+-------------------------------------+
Sigma 规则是 pf 和 cf 的组合。为了将 pf 标签带回当前行，我们需要对 pf 标签应用时间旅行。将 flux-capacitor 函数应用于 df_input 数据框。
jdf = flux_stateful_function(df_input._jdf, "host_id", bloom_capacity, spec, True)
df_output = DataFrame(jdf, spark)
我们获得了 df_output 数据框。注意 pf 标签如何在时间中传播。
+-------------------+------+----------+--------------+------------------------------------+
|timestamp |pid |parent_pid|human_readable|sigma |
+-------------------+------+----------+--------------+------------------------------------+
|2022-12-30 00:00:05|pid500| |[pf] |{rule2 -> {pf -> true, cf -> false}}|
|2022-12-30 00:00:06|pid600|pid500 |[pf] |{rule2 -> {pf -> true, cf -> false}}|
|2022-12-30 00:00:07|pid700|pid600 |[pf, cf] |{rule2 -> {pf -> true, cf -> true}} |
+-------------------+------+----------+--------------+------------------------------------+
本笔记本 TagPropagationIllustration.ipynb 包含了更多类似的父子和时间邻近的示例。
构建带有上下文的警报
flux-capacitor 函数缓存所有过去的标签。为了节省内存，它使用布隆过滤器段缓存这些标签。布隆过滤器具有极小的内存占用，查询和更新都很快。然而，它们确实引入了可能的假阳性。因此，我们的检测可能实际上是一个假阳性。为了解决这个问题，我们将怀疑的异常放入队列 (4) 进行重新评估。
为了消除假阳性，第二个名为 警报生成器 的 Spark 流作业读取怀疑的异常 (5) 并检索需要重新评估规则的事件 (6)。
例如，在父子 Sigma 规则的情况下，警报生成器将读取怀疑的异常 (5)，检索子进程事件。接下来，在 (6) 中，它将检索该子事件的父进程。然后，利用这两个事件重新评估 Sigma 规则。然而，这次 flux-capacitor 配置为将标签存储在哈希映射中，而不是布隆过滤器中。这消除了假阳性，并且作为额外的好处，我们获得了所有参与检测的事件。我们将此警报及证据行（父子事件）存储到警报表中 (7)。
带有状态检测的拓扑（时间）
警报生成器 处理的量仅为 (2) 流检测 的一小部分。由于在 (5) 中读取的低量，历史搜索到标记的遥测 (6) 是可能的。
想要更深入地了解，请查看 流检测 的 Spark 作业 streaming_detections.py 和 警报生成器 streaming_alert_builder.py。
性能
为了评估这个概念验证的性能，我们在拥有 16 个 CPU 和 64G 内存的机器上进行了测试。我们编写了一个简单的数据生成器，每秒创建 5,000 个合成事件，并在 30 天内进行了实验。
Spark Streaming Detections 作业在一台机器上运行。该作业配置为每分钟触发一次。每个微批次（触发）读取 300,000 个事件，平均需要 20 秒完成。该作业可以轻松跟上输入事件的速率。
Spark Streaming Detections
Spark Alert Builder 在单台机器上运行，并配置为每分钟触发一次。这个作业完成需要 30 到 50 秒的时间。该作业对 tagged_telemetry_table 的组织方式非常敏感。这里我们看到维护作业的效果，每小时组织并排序最新数据。因此，每小时，Spark Alert Builder 的微批处理执行时间恢复到 30 秒。
Spark Streaming Alert Builder
表维护
我们的 Spark 流作业每分钟触发一次，因此每分钟产生小型数据文件。为了在这个表中实现快速搜索和检索，定期对数据进行压缩和排序非常重要。幸运的是，Iceberg 提供了内置程序来组织和维护您的表。
例如，这个脚本 maintenance.py 每小时运行一次，用于对 Iceberg 的 tagged_telemetry_table 中新增加的文件进行排序和压缩。
CALL catalog.system.rewrite_data_files(
table => 'catalog.jc_sched.tagged_telemetry_table',
strategy => 'sort',
sort_order => 'host_id, has_temporal_proximity_tags',
options => map('min-input-files', '100',
'max-concurrent-file-group-rewrites', '30',
'partial-progress.enabled', 'true'),
where => 'timestamp >= TIMESTAMP \'2023-05-06 00:00:00\' '
)
每天结束时，我们还会重新对这个表进行排序，以在长时间搜索期间（几个月的数据）实现最大的搜索性能。
CALL catalog.system.rewrite_data_files(
table => 'catalog.jc_sched.tagged_telemetry_table',
strategy => 'sort',
sort_order => 'host_id, has_temporal_proximity_tags',
options => map('min-input-files', '100',
'max-concurrent-file-group-rewrites', '30',
'partial-progress.enabled', 'true',
'rewrite-all', 'true'),
where => 'timestamp >= TIMESTAMP \'2023-05-05 00:00:00\' AND timestamp < TIMESTAMP \'2023-05-06 00:00:00\' '
)
我们还执行另一个维护任务，即从流表中删除旧数据。这些表仅用作生产者和消费者之间的缓冲区。因此，每天我们会对流表进行老化处理，保留最近 7 天的数据。
delete from catalog.jc_sched.process_telemetry_table
where
timestamp < current_timestamp() - interval 7 days
最后，每天我们执行标准的 Iceberg 表维护任务，如过期快照和删除孤立文件。我们在所有表上运行这些维护作业，并在Airflow上安排这些作业的时间。
结论
在本文中，我们展示了如何构建一个 Spark 流异常检测框架，通用地应用 Sigma 规则。新的 Sigma 规则可以轻松添加到系统中。
这个概念验证在合成数据上进行了广泛测试，以评估其稳定性和可扩展性。它显示出很大的潜力，将在生产系统上进一步评估。
除非另有说明，所有图片均为作者所有
使用 Sigma 规则进行异常检测（第一部分）：利用 Spark SQL 流处理
Sigma 规则用于检测网络安全日志中的异常。我们使用 Spark 结构化流处理来大规模评估 Sigma 规则。
·
关注 发表在 Towards Data Science ·8 min 阅读·2023 年 1 月 24 日
–
摄影：Tom Carnegie，来源于 Unsplash，加拿大最高法院
数据草图的兴起
数据草图是一个总括性术语，涵盖了使用理论数学、统计学和计算机科学的各种数据结构和算法，以解决集合基数、分位数、频率估计等问题，并具有数学上证明的误差范围。
数据草图比传统方法快几个数量级，它们需要更少的计算资源，有时是解决大数据问题的唯一可行解决方案。要了解更多关于数据草图的信息，请查看 Apache Data Sketch 项目。
草图实现了可以从中提取信息的算法。
在单次传输的数据流中，这也被称为“单次触摸”处理。
Spark 大量利用数据草图，例如：维度缩减、局部敏感哈希、计数最小草图。
在这一系列文章中，我们将带你深入了解高性能欺诈检测系统的设计。通过实际例子，我们评估并对比了传统算法与基于数据草图的算法的性能。
什么是 Sigma 规则
Sigma 是一种通用签名格式，允许你在日志事件中进行检测。规则易于编写，适用于任何类型的日志。最重要的是，Sigma 规则是抽象的，不绑定于任何特定的 SIEM，使得 Sigma 规则可共享。
一旦网络安全研究人员或分析师开发了检测方法，他们可以使用 Sigma 描述并与他人分享他们的技术。以下是来自 Sigma HQ 的一段话：
规则的核心是检测部分。当条件评估为 true 时，意味着我们进行了检测。条件由命名表达式组成。例如，这里声明了选择和过滤表达式。这些表达式对日志的属性进行测试。在这种情况下是 web 日志。
Sigmac 生成 SQL
Sigmac 编译器用于将抽象的 Sigma 规则转换为具体的形式，以便实际的 SIEM 或处理平台进行评估。Sigmac 具有许多后端，能够将规则转换为 QRadar、ElasticSearch、ArcSight、FortiSIEM 和通用 SQL。
使用 SQL sigmac 后端，我们可以将上述规则转换为：
SELECT
*
FROM
(
SELECT
(
(cs-uri-query LIKE '%cmd=read%'
OR cs-uri-query LIKE '%connect&target%'
OR cs-uri-query LIKE '%cmd=connect%'
OR cs-uri-query LIKE '%cmd=disconnect%'
OR cs-uri-query LIKE '%cmd=forward%')
AND (cs-referer IS NULL
AND cs-USER-agent IS NULL
AND cs-METHOD LIKE 'POST'))
AS web_webshell_regeorg,
*
FROM
test_webserver_logs
)
WHERE
web_webshell_regeorg = TRUE
这些 SQL 语句通常由调度器在特定触发间隔（例如 1 小时）下调用。每小时，检测系统会搜索最新的事件。
然而，一些 Sigma 规则应用了时间聚合。例如，通过全球目录进行枚举计算一段时间窗口内事件的发生次数。
detection:
selection:
EventID: 5156
DestPort:
- 3268
- 3269
timeframe: 1h
condition: selection | count() by SourceAddress > 2000
使用上述批处理模型，这些类型的查询会一遍又一遍地重新处理相同的事件。特别是当相关窗口很大时。此外，如果我们试图通过将触发频率提高到每5 分钟来减少检测延迟，我们将引入更多对相同事件的重新处理。
理想情况下，为了减少对相同事件的反复处理，我们希望异常检测能够记住上一个处理的事件是什么，以及迄今为止计数器的值。这正是 Spark Structured Streaming 框架所提供的功能。流式查询每分钟触发一次微批处理（可配置）。它读取新事件，更新所有计数器并将其持久化（用于灾难恢复）。
在这种模型中，每个事件只评估一次。提高触发频率不会像无状态批处理模型那样产生相同的成本。而且由于事件只评估一次，复杂的检测（如正则表达式匹配）不会产生膨胀的成本。
使用 Spark Streaming 运行检测
Spark Structured Streaming 可以轻松评估 sigmac 编译器生成的 SQL。首先，我们通过连接到我们最喜欢的队列机制（EventHubs，Kafka）来创建一个流数据框。在本例中，我们将从一个 Iceberg 表中readStream，该表中事件会被增量插入。了解有关 Iceberg 流能力的更多信息，请点击这里。
# current time in milliseconds
ts = int(time.time() * 1000)
# create a streaming dataframe for an iceberg table
streamingDf = (
spark.readStream
.format("iceberg")
.option("stream-from-timestamp", ts)
.option("streaming-skip-delete-snapshots", True)
.load("icebergcatalog.dev.events_table")
)
# alias the dataframe to a table named "events"
streamingDf.createOrReplaceTempView("events")
注意，我们将流数据框别名为视图名称events。这样做是为了在 SQL 语句中引用此流数据框，即：select * from events。我们现在要做的就是配置 sigmac 编译器，以便对events表生成 SQL 语句。例如，生成的 sql 文件可能如下所示：
SELECT
(cs-uri-query LIKE '%cmd=read%'
OR cs-uri-query LIKE '%connect&target%'
OR cs-uri-query LIKE '%cmd=connect%'
OR cs-uri-query LIKE '%cmd=disconnect%'
OR cs-uri-query LIKE '%cmd=forward%')
AND (cs-referer IS NULL
AND cs-USER-agent IS NULL
AND cs-METHOD LIKE 'POST'))
AS web_webshell_regeorg,
-- another detection here
cs-uri-query LIKE '%something%'
AS detection2
*
FROM
events
在我们的分析中，我们加载生成的 SQL 并要求 Spark 从中创建一个hitsDf数据框。
# load auto-generated SQL statement
with open('./generated_sql_statement.sql', 'r') as f:
detections_sql = f.read()
hitsDf = spark.sql(detections_sql)
我们通过调用writeStream启动流查询，并配置查询以每分钟触发一次微批处理。此流查询将无限期运行，将检测结果写入我们选择的接收端。这里我们只是将结果写入控制台接收端，但我们也可以写入另一个 Iceberg 表。或者，我们可以使用[forEachBatch](https://spark.apache.org/docs/latest/structured-streaming-programming-guide.html#using-foreach-and-foreachbatch)执行一些任意的 python 代码，例如，将通知推送到 REST 端点。或者我们甚至可以同时做这两件事。
# start a streaming query printing results to the console
query = (
hitsDf.writeStream
.outputMode("append")
.format("console")
.trigger(processingTime="1 minute")
.start()
)
父进程挑战
到目前为止，我们已经看到了如何检测离散事件中的异常。然而，Sigma 规则可以将事件与之前的事件相关联。一个经典的例子是在 Windows 安全日志 (事件 ID 4688) 中找到的。在这个日志源中，我们可以找到有关进程创建的信息。该日志中的一个关键部分是启动此进程的进程。你可以使用这些 Process ID 来确定程序在运行时做了什么等。
以这个 Sigma 规则为例：Rundll32 执行没有 DLL 文件。
detection:
selection:
Image|endswith: '\rundll32.exe'
filter_empty:
CommandLine: null
filter:
- CommandLine|contains: '.dll'
- CommandLine: ''
filter_iexplorer:
ParentImage|endswith: ':\Program Files\Internet Explorer\iexplore.exe'
CommandLine|contains: '.cpl'
filter_msiexec_syswow64:
ParentImage|endswith: ':\Windows\SysWOW64\msiexec.exe'
ParentCommandLine|startswith: 'C:\Windows\syswow64\MsiExec.exe -Embedding'
filter_msiexec_system32:
ParentImage|endswith: ':\Windows\System32\msiexec.exe'
ParentCommandLine|startswith: 'C:\Windows\system32\MsiExec.exe -Embedding'
filter_splunk_ufw:
ParentImage|endswith: ':\Windows\System32\cmd.exe'
ParentCommandLine|contains: ' C:\Program Files\SplunkUniversalForwarder\'
filter_localserver_fp:
CommandLine|contains: ' -localserver '
condition: selection and not 1 of filter*
在原始遥测中，一个事件只知道父级 Process ID。然而，规则中提到 ParentImage 和 ParentCommandLine。规则基本上假设已经进行了连接。
幸运的是，Spark Structured Streaming 支持 流-流连接。为了检索 ParentImage 和 ParentCommandLine，我们对进程日志进行自连接。我们将 current 侧与 parent_of_interest 侧连接。连接条件如下：
current.ParentProcessID = parent_of_interest.ProcessID
左侧：为每个检测规则设置标志
我们使用 c 作为当前进程的约定，使用 r1 作为规则 1。
因此，在 Rundll32 执行没有 DLL 文件（规则 1）中，filter_empty 被
%%sparksql --view current --output skip
select
*,
ID,
CommandLine,
ImagePath,
-- rule 1
ImagePath ilike '%\\\\rundll32.exe' as cr1_selection,
Commandline is null as cr1_filter_empty,
Commandline ilike '%.dll%' OR Commandline = '' as cr1_filter,
Commandline ilike '% -localserver %' as cr1_filter_localserver_fp
from
events
右侧：在 parents_of_interest 上过滤消息
对于应用于父进程的条件，我们也进行相同的操作。然而，在这种情况下，我们还会过滤表格。这意味着被过滤的父进程的所有标志必然都设置为 FALSE。通过过滤，我们大大减少了执行流式连接时需要缓存的 parents_of_interest 键的数量。
%%sparksql --output skip --view parents_of_interest
select
*
from (
select
host_id as parent_host_id,
ID as parent_id,
ImagePath as parent_imagepath,
CommandLine as parent_commandline,
-- rule 1
(ImagePath ilike '%:\Program Files\Internet Explorer\iexplore.exe'
AND CommandLine ilike '%.cpl%')
as pr1_filter_iexplorer,
(ImagePath ilike '%:\Windows\SysWOW64\msiexec.exe'
AND CommandLine ilike 'C:\Windows\syswow64\MsiExec.exe -Embedding%')
as pr1_filter_msiexec_syswow64,
(ImagePath ilike '%:\Windows\System32\msiexec.exe' AND
CommandLine ilike 'C:\Windows\system32\MsiExec.exe -Embedding%')
as pr1_filter_msiexec_system32
from
events
)
where
pr1_filter_iexplorer
OR pr1_filter_msiexec_syswow64
OR pr1_filter_msiexec_system32
将当前与其父级连接
我们对父级侧进行左连接。由于父级侧被过滤，因此可能找不到对应的父进程 ID。当未找到父进程时，列将具有标志设置为 NULL。我们使用 coalesce 将这些父级标志的值设置为 FALSE。pr3_selection_atexec 是一个父标志，因此我们像这样应用 coalesce：
coalesce(pr3_selection_atexec, FALSE)
我们还组合来自当前和父级的条件。例如，selection_atexec 条件由父级和子级条件组成。
selection_atexec:
ParentCommandLine|contains:
- 'svchost.exe -k netsvcs'
- 'taskeng.exe'
CommandLine|contains|all:
- 'cmd.exe'
- '/C'
- 'Windows\Temp\'
- '&1'
因此我们像这样组合它们：
cr3_selection_atexec AND coalesce(pr3_selection_atexec, FALSE)
as r3_selection_atexec,
r3_selection_atexec 是规则 3 中 selection_atexec 的最终标志。
%%sparksql --view joined --output skip
select
--rule1
cr1_selection as r1_selection,
cr1_filter_empty as r1_filter_empty,
cr1_filter as r1_filter,
cr1_filter_localserver_fp as r1_filter_localserver_fp,
coalesce(pr1_filter_iexplorer, FALSE) as r1_filter_iexplorer,
coalesce(pr1_filter_msiexec_syswow64, FALSE) as r1_filter_msiexec_syswow64,
coalesce(pr1_filter_msiexec_system32, FALSE) as r1_filter_msiexec_system32,
parent_host_id,
parent_id,
parent_imagepath,
parent_commandline
from
current as c
left join parents_of_interest as p
on c.ParentProcessID = p.parent_id
最后我们应用 Sigma 规则条件
例如规则 1 的条件是：
condition: selection and not 1 of filter*
我们只需应用这个条件，并将结果命名为 rule1。
r1_selection AND NOT (r1_filter_empty
OR r1_filter
OR r1_filter_localserver_fp
OR r1_filter_iexplorer
OR r1_filter_msiexec_syswow64
OR r1_filter_msiexec_system32)
as rule1,
以下是完整语句中的条件。
%%sparksql --output json
select
*
from (
select
*,
-- rule 1 -> condition: selection and not 1 of filter*
r1_selection AND NOT (r1_filter_empty
OR r1_filter
OR r1_filter_localserver_fp
OR r1_filter_iexplorer
OR r1_filter_msiexec_syswow64
OR r1_filter_msiexec_system32)
as rule1,
from
joined
)
where
rule1 = TRUE
请注意，sigmac 编译器不会生成这种类型的 SQL。然而，我们计划编写一个自定义的 sigma 编译器以生成上述 SQL 语句。
执行这些 SQL 语句与我们最初的示例没有区别。
Spark Structured Streaming 可以在微批处理之间保持和持久化状态。在文档中，Spark 称之为窗口分组聚合。相同的原则适用于流-流连接。你可以配置 Spark 来缓存聚合，或者在这种情况下，缓存parents_of_interests的行在一个窗口中。
然而，这种方式的扩展性如何？我们可以在 Spark 的状态存储窗口中保留多少行parents_of_interests？
在我们下一篇文章中，我们将回答这些问题。为了不遗漏，请关注我们并订阅邮件获取这些故事。敬请期待！
除非另有说明，所有图像均由作者提供。
使用 Sigma 规则的异常检测（第二部分） Spark 流-流连接
一类 Sigma 规则检测时间相关性。我们评估了 Spark 的有状态对称流-流连接在执行时间相关性时的可扩展性。
·
关注 发表在Towards Data Science ·7 分钟阅读·2023 年 2 月 2 日
–
由 Naveen Kumar 拍摄，来自 Unsplash
跟进我们之前的文章，我们评估了 Spark 将一个开始进程事件与其父开始进程事件进行连接的能力。
在这篇文章中，我们评估了 Spark 流-流连接的扩展性。具体来说，它能在连接窗口中容纳多少事件。
在我们的研究中，我们评估了几种方法：
完全连接
完整的流-流连接需要缓存连接右侧的所有先前事件（父项启动过程）。由于只有这些父项启动过程事件的子集是感兴趣的，因此不需要所有过去的父项启动过程详细信息。例如，一个 Sigma 规则可能会指定一个包含字符串 .cpl 的父项 CommandLine——所有其他事件可以被忽略。
加入感兴趣的父项
感兴趣的父项是将过滤条件应用于连接的右侧得到的结果。这可以大大减少需要记住的父项数量。连接完成后，我们对当前处理和父项处理应用条件。
与感兴趣的父项特征连接
更好的解决方案是存储我们在右侧评估的条件，并丢弃所有其他属性——CommandLine、Image 等。这样我们只保留有限数量的布尔标志，而不是可能很长的字符串。在下图中，Features 是 Sigma 过滤器表达式名称的映射，值是测试结果。例如：
features = {
'rule3_selection_atexec' -> false,
'rule2_selection' -> true
}
在我们的研究中，我们很快意识到减少 Spark 需要存储的状态量至关重要。因此，我们选择只保留感兴趣的父项。这些是我们正在寻找的特征的父项。我们丢弃所有其他父项，只保留这些父项的最小信息集：连接键、时间戳和特征标志。
模拟测试框架
为了评估 Spark 流-流连接的性能，我们创建了一个 cause 和 effect 事件的模拟流。在我们的实验中，我们通过设置水印为零来禁用迟到支持。
cause = (
spark.readStream
.format("rate")
.option("rowsPerSecond", rate)
.load()
.withWatermark("timestamp", "0 seconds")
.withColumn("name", F.lit(name))
)
这些cause和effect流事件是通过 Spark 流-流连接操作进行合并的：
cause = cause.select('cause_timestamp', 'cause_key', 'cause_load')
effect = effect.select(
'effect_timestamp',
'effect_key',
'effect_load',
'host_id',
'id',
'name',
'value')
joindf = (
effect.join(
cause,
F.expr(f"""
effect_key = cause_key
and effect_timestamp >= cause_timestamp
and effect_timestamp < (cause_timestamp + interval {window} seconds)
"""),
"left"
)
joindf
.writeStream
.format("iceberg")
.outputMode("append")
.trigger(processingTime="1 minutes")
.queryName(name)
.option("checkpointLocation", checkpoint)
.toTable(output_table_name)
注意连接表达式 effect_key = cause_key 和窗口子句，表示效果时间必须在 cause 之后，但不能比 window 秒更久。
Linxiao Ma 在他的文章 Spark Structured Streaming Deep Dive (7) — Stream-Stream Join 中详细解释了，在这些条件下，原因事件会被缓存到 Spark 的有状态存储中，最多 window 秒。然而，effect 事件不需要被存储。对于每个在流-流连接中经过的 effect 事件，都会进行查找以找到相应的 cause。对于每个进入连接的 effect，都会写入一行 cause+effect。
选择合适的状态存储
Spark 有两个状态存储实现。最初的是名为 HDFSBackedStateStore 的，它是一个由 HDFS 文件支持的简单内存哈希表。最新的状态存储基于 RocksDB。RocksDB 是一个可嵌入的键值持久存储用 C++编写。RocksDB 的状态部分保存在内存中，部分保存在本地磁盘上。在每个检查点，Spark 将更改的文件副本保存到中央位置（数据湖）。
当你需要存储大量键时，Spark 推荐使用RocksDB。根据DataBricks，一个大型 Spark 工作节点可以缓存多达 1 亿个键。
由于我们的流流连接会缓存大量的感兴趣的父项行，我们决定在评估中使用 RocksDB 状态存储。
.config("spark.sql.streaming.stateStore.providerClass",
"org.apache.spark.sql.execution.streaming.state.RocksDBStateStoreProvider")
我们所有的实验都在一个具有 48G RAM 和 16 CPU 的单个 Spark 工作节点上进行。我们模拟了来自 50,000 个主机的日志。
我们的测试框架非常灵活。它允许我们更改许多参数，例如每秒事件数、每个事件的大小、连接的时间窗口、键大小、时间窗口中事件的分布等。
Spark 分区的影响
在我们的第一个实验中，我们在 10,000 秒（约 2.77 小时）的窗口中连接效果和原因。我们模拟了每个感兴趣的父项会有 12 个布尔标志。我们设置了每秒 10,000 个事件的速率。在这里，我们展示了 Spark 分区（单个任务）数量变化的影响。
更改分区数量对性能没有影响。执行一个微批次的时间大约是 225 秒。记住我们每 60 秒触发一次.trigger(processingTime="1 minutes")。Spark 将立即开始下一个微批次。因此，事件处理的延迟最多为 225 秒。
窗口大小的影响
在第二个实验中，我们调整了流流连接窗口的大小（时间）。在每秒 5,000 个事件的速率下，作业不稳定。每个微批次的执行时间越来越长。我们正在落后。
如果我们将窗口减少到 18 小时，并将速率降低到每秒 2,500 个事件，作业会稳定下来，并在每个微批次约 300 秒时达到稳定。
然而，实际上，我们不会保留每一个父事件。我们只会保留“感兴趣的父事件”。这些事件中有一个或多个为真的 Sigma 规则表达式。重要的是要衡量 Spark 存储父事件的能力。我们可以轻松计算：2,500 事件/秒 x 64,000 秒。Spark 可以缓存 1.6 亿个“感兴趣的父事件”。我们的实验结果确认了 Databricks 关于 RocksDB StateStore 的声明，即每台机器可以处理 1 亿个键。如果我们假设这些事件来自 50,000 个主机，那么 Spark 可以每台主机保存 3,200 个“感兴趣的父事件”。
有趣的是，Spark 存储了特征标志和键。它需要存储 causes 的键，以便将其与 effect 键连接。
我们可以做的另一个有趣的观察是，通过键查找我们检索到什么？我们检索到一个包含布尔标志的事件（即一行）。实际上，这些标志通常是互斥的。也就是说，在某一行中，只有一个标志可能为真，而所有其他标志都为假。Spark 存储 cause 键和所有相关的标志，无论它们是真还是假。
Bloom 过滤器连接
是否有更好的方法来跟踪特征标志？是的，答案是 Bloom 过滤器。
Bloom 过滤器是一种概率数据结构，可以存储一个键并测试键的存在性。Bloom 过滤器对键进行哈希，并利用哈希结果在位数组中设置一些位。
Bloom 过滤器非常紧凑。为了节省空间，你可能会付出假阳性的代价。然而，一旦检测完成，触发的 Sigma 规则可以重新评估以确认正确性。
我们可以使用 Bloom 过滤器来执行上述连接。假设我们使用一个复合键（parent_key + feature_id），其中 feature_id 是给定 Sigma 过滤器表达式的名称。应用于父进程的过滤器表达式存储在 Bloom 过滤器中，但仅在它们为真的情况下。测试复合键的存在性，如果键在 Bloom 中，则返回真；如果不在，则返回假。
Bloom 可以存储一定数量的键。超出这个数量后，假阳性会急剧增加。通过只存储真实的特征，我们可以延长 Bloom 过滤器的有效性。
因此，连接被建模为在 Bloom 过滤器中的查找。
在我们的下一篇文章中，我们将构建一个自定义的 Spark 有状态连接函数，利用 Bloom 过滤器。
除非另有说明，否则所有图像均由作者提供
使用 Sigma 规则进行异常检测（第三部分）基于布隆过滤器的时间相关性
一个基于布隆过滤器的定制状态映射函数能否超越通用的 Spark 流-流连接？
·
关注 发表在 Towards Data Science ·6 分钟阅读·2023 年 2 月 14 日
–
由 Kalpaj 拍摄于 Unsplash，Peggys Cove，加拿大新斯科舍省
这是我们系列文章的第 3 篇。请参考第一部分和第二部分以获得一些背景信息。
Spark 的 flatMapGroupsWithState 函数允许用户在分组数据上应用自定义代码，并提供支持来持久化用户定义的状态。
在本文中，我们将实现一个有状态的函数，该函数检索父进程的标签（特征）。解决方案的关键是创建一个由进程 ID（下图中的 e_key）和我们想要记住的标签（特征）名称组成的复合键。例如，在第 1 行中，我们创建一个布隆键 300.pr1_filter_msiexec_syswow64。我们将此键存储在布隆过滤器中。请注意，由于布隆过滤器的特性，键实际上并不存储在布隆中。相反，键被哈希，并且该哈希用于在布隆的位数组中打开一些位。因此，如果标签为 false，我们不将其存储在布隆中。只有标记为 true 的才会存储。
在第二行中，我们展示了如何使用父进程的进程 ID (pe_key 而非 e_key) 创建一个检索布隆键。因此，查找键为 300.pr1_filter_msiexec_syswow64。如果查询布隆键返回 false，我们可以确定它从未存储在布隆中。然而，如果查询返回 true，我们知道该键确实存储在布隆中（但存在假阳性的可能性）。
我们使用布隆返回的结果来更新当前行的 pr1_filter_msiexec_syswow64。这有效地执行了我们上一篇文章中的连接操作。使用当前行的标签和检索到的当前行父进程的标志，我们最终可以评估完整的 Sigma 规则条件。
测试工具
为了评估这种方法的性能，我们在 Scala 中构建了一个初步的 flatMapGroupsWithState 原型，并使用了 Guava 布隆过滤器库。
为了估算布隆过滤器的大小，我们使用了这个 实用的在线计算器。容量越大（可以存储的标签数量），布隆过滤器的大小就越大。例如，容量为 200,000 的布隆过滤器的大小约为 234KiB，假阳性率为 1%。
我们为每个主机使用一个布隆过滤器。继续模拟 50,000 个主机，这大约给我们提供了 12GiB 的状态。
我们最初使用 RocksDB 状态存储进行所有实验，因为在进行流流连接时表现非常好。然而，正如经常发生的那样，我们进行了一个意外的发现。旧版 HDFS 状态存储在这种用例中实际上比 RocksDB 表现更好。以下是结果的简要总结：
使用 RocksDB 的执行时间不断增加，而 HDFS 在每个微批次约 200 秒时保持稳定。
RocksDB 性能不佳的原因是它重新组织本地状态文件。在每次触发时，我们加载和修改每个 bloom，每个 bloom 的大小为 200KiB。然后，RocksDB 需要对其 SST 文件进行压缩和重组。RocksDB 似乎更适合处理大量较小的值。而 HDFS 则更适合处理较大的值，但在处理非常大量的数据时表现不佳。
与流-流连接的结果比较
现在让我们将这些结果与流-流连接的结果进行比较。在流-流连接中，我们使用了每秒 2,500 个事件的速率，而使用 bloom 时的速率为 5,000。
此外，在流-流连接中，我们可以保持平均 3,200 个“兴趣的父标签”。即与 3,200 个事件相关联的标签。如果我们假设每个事件只有一个真实标签，这将产生每个主机 3,200 个标签。
使用 bloom，我们可以每台主机保持约 100,000 个标签。这大约是原来的 30 倍！
使用自定义的 flatMapGroupsWithState，我们仅读取一次输入流（例如来自 Kafka 或 Iceberg 表）。而使用连接时，我们必须读取两次流，一次用于连接的每一侧。
自定义的 flatMapGroupsWithState 更加灵活
此外，bloom 方法适用于更有趣的用例。例如，假设你不仅想记住你父级的标签，还想记住你父级的父级的标签。这使用流-流连接可能会相当困难。你要连接多少次？
通过简单的调整，我们可以支持传播祖先标签。在上面的示例中，在第 2 行，我们提取了 bloom 键 300.pr1_filter_msiexec_syswow64 和列中的值 pr1_filter_msiexec_syswow64。为了支持祖先标签，我们将该标签重新放入 bloom 中，但使用第 2 行的 e_key 值。因此，我们存储了 bloom 键 600.pr1_filter_msiexec_syswow64。
现在，假设我们得到一个子进程 900（进程 300 的孙子）。如前所述，我们使用我们的父键 pe_key（600）进行提取，形成 bloom 键 600.pr1_filter_msiexec_syswow64，从而有效地提取我们祖父的标签。
超越父子关系
到目前为止，我们集中于父子关系。然而，通过存储和检索标签，我们还可以支持 时间邻近关联。时间邻近关联是即将发布的 Sigma 规范中的一项新功能。规范中给出的示例如下：
action: correlation
type: temporal
rules:
- recon_cmd_a
- recon_cmd_b
- recon_cmd_c
group-by:
- ComputerName
- User
timespan: 5m
ordered: false
该规则规定，3 个事件（recon_cmd_a/b/c）必须发生在一个时间窗口内。这些事件不必按任何特定顺序出现。
正如我们之前所做的，我们可以将这些事件表示为标签，并对每个主机（每个 ComputerName）使用一个布隆过滤器。我们使用列User，而不是使用列e_key/pe_key（父/子）。列User为我们存储和从布隆过滤器中检索的复合键提供了上下文。
在第一行，我们看到recon_cmd_b为真，因此我们使用User作为上下文将其存储在布隆过滤器中，从而将其存储在Alice.recon_cmd_b下。
在第二行，我们看到recon_cmd_a为真，因此我们将其存储在Alice.recon_cmd_a下。在每一行中，我们总是尝试提取键：Alice.recon_cmd_a、Alice.recon_cmd_b和Alice.recond_cmd_c。因此，recom_cmd_b被更新为真。
对于第三行，我们推送Alice.recond_cmd_c并提取其他两个。
Sigma 规范支持这种关联的ordered版本。支持有序时间接近性也是可能的。为了支持有序事件，我们有条件地提取标签，而不是始终提取。
如果我们假设上述规则是有序的，我们将把程序从始终存储和提取所有键改为有条件地存储recon_cmd_b，仅当recon_cmd_a已经被看到时，并且永远不存储recon_cmd_c。
在我们的下一篇文章中，我们将改进我们的概念验证。我们将使其更加灵活，并能够支持上述描述的附加用例。
除非另有说明，否则所有图片均由作者提供
使用 Sigma 规则进行异常检测（第四部分）：Flux 电容器设计
我们实现了一个 Spark 结构化流式状态映射函数，以处理网络安全日志中的时间相关性。
·
关注 发布于 数据科学前沿 · 6 分钟阅读 · 2023 年 3 月 2 日
–
图片来自 Pixabay 的 Robert Wilson
这是我们系列文章的第 4 篇。有关更多背景信息，请参考 第一部分 、第二部分 和 第三部分。
在本文中，我们将详细介绍自定义 Spark flatMapWithGroupState 函数的设计。我们选择用 Scala 编写这个函数，因为 Spark 本身是用 Scala 编写的。我们将这个函数命名为 Flux Capacitor。电容器积累电荷并随后释放它们。类似地，我们的 flatMapWithGroupState 将积累标签（评估为真/假的 Sigma 表达式），然后释放它们。
我们的 Flux Capacitor 函数易于配置，允许用户指定每个标签的存储和检索方式。在我们的实现中，我们将标签的更新和检索行为封装在一个标签评估器类层次结构中。
标签评估器
除非另有说明，所有标签将由过渡评估器进行评估。这个评估器是一个无操作的评估器，它仅仅将当前标签值传递出去。可缓存评估器是一个能够在布隆过滤器中存储和检索标签的基类。这个基类有两个模板方法，makeBloomPutKey 和 makeBloomGetKey，让子类决定标签的存储和检索方式。默认标签使用相同的键在布隆过滤器中存取值。这个键只是规则和标签名称的拼接。当用户将以下规范传递给 Flux Capacitor 时，将创建默认标签评估器对象来处理标签的存储和检索。
rules:
- rulename: rule5
description: un-ordered set of events by host
action: temporal
ordered: false
tags:
- name: recon_cmd_a
- name: recon_cmd_b
- name: recon_cmd_c
我们设计的 Flux Capacitor 规范类似于 Sigma 规则。然而，这个规范仅描述了 Flux Capacitor 中标签的处理方式。
有序评估器通过引入“依赖标签”来专业化默认评估器。依赖标签必须首先被评估，只有当依赖标签为真时，有序标签才会被评估。当用户指定 ordered: true 时，将使用此评估器。
层级结构的另一侧处理父子关系。与默认标签相反，父评估器的 put 键由 规则名称 + 标签名称 + 父 ID 组成，而 get 键由 规则名称 + 标签名称 + 当前 ID 组成。父子关系的规范允许用户指定哪个列名包含父 ID，哪个列名包含当前消息 ID。在这个例子中，parent_id 和 id 是表示父子关系的列。
rules:
- rulename: rule4
description: propagate to child only
action: parent
parent: parent_id
child: id
tags:
- name: parent_process_feature1
正如我们在上一篇文章中所解释的，祖先评估器通过将自身存储在布隆过滤器中，从而沿父子层级传播自身，专业化了父评估器。
规则评估器
到目前为止，我们只应用了一个 Sigma 规则。我们这样做是为了保持简单。然而，实际上我们希望同时应用多个 Sigma 规则。当我们处理日志源并解析其内容时，我们希望应用所有适用的 Sigma 规则。
为了保持组织性，我们在MapType类型的列中表示每个规则及其相关的标签。这个映射中的键是规则名称，值是标签名称到标签值（true/false）的映射。
select
-- regroup each rule's tags in a map (ruleName -> Tags)
map(
'rule1', map('cr1_selection', cr1_selection,
'cr1_filter_empty', cr1_filter_empty,
'cr1_filter', cr1_filter,
'cr1_filter_localserver_fp', cr1_filter_localserver_fp,
'pr1_filter_iexplorer', pr1_filter_iexplorer,
'pr1_filter_msiexec_syswow64', pr1_filter_msiexec_syswow64,
'pr1_filter_msiexec_system32', pr1_filter_msiexec_system32),
'rule2', map('cr2_filter_provider', cr2_filter_provider,
'cr2_filter_git', cr2_filter_git,
'pr2_selection', pr2_selection,
'pr2_filter_git', pr2_filter_git),
'rule3', map('cr3_selection_other', cr3_selection_other,
'cr3_selection_atexec', cr3_selection_atexec,
'pr3_selection_other', pr3_selection_other,
'pr3_selection_atexec', pr3_selection_atexec)
) as sigma
通过使用映射，我们可以保持输入和输出行的模式不变。任何稍后引入的新 Sigma 规则将简单地添加到这个映射中，而不会影响整体输入和输出行模式。
我们的 Flux Capacitor 规范反映了这一点。规范适用于一系列规则名称。每个规则规范都说明了如何更新该规则的标签。
rules:
- rulename: rule1
description: ordered list of events by host
action: temporal
ordered: true
tags:
- name: recon_cmd_a
- name: recon_cmd_b
- name: recon_cmd_c
- rulename: rule3
description: propagate to all decendents
action: ancestor
tags:
- name: ancestor_process_feature1
parent: parent_id
child: id
在我们的实现中，每个规则规范由一个Rule类处理。Rules类依次评估每个Rule。
映射函数
Spark Dataframe flatMapGroupsWithState在每个组上调用用户提供的函数，同时在调用之间维护用户定义的每组状态。该函数的签名接受一个键、一个输入行迭代器（每组）和一个状态，并返回一个输出行迭代器。我们的 Flux Capacitor 函数如下：
case class FluxCapacitorMapFunction(
val tagCapacity: Int,
val specification: String
) {
def processBatch(
key: String,
rows: Iterator[Row],
state: GroupState[FluxState]
): Iterator[Row] = {
val bloom =
if (state.exists()) {
state.get()
} else {
BloomFilter.create(Funnels.stringFunnel(), tagCapacity, 0.01)
}
val rulesConf = RulesConf.load(specification)
val rules = new RulesAdapter(new Rules(rulesConf, bloom))
val outputRows = rows.map(row => rules.evaluateRow(row))
state.update(bloom)
outputRows.iterator
}
第一次调用函数时，我们将没有状态。在这种情况下，我们创建一个假阳性概率为 1%的 Guava 布隆过滤器。
然后我们从 YAML 字符串加载规范。稍后我们将看到用户如何将这个规范传递给函数。规范被解析并提供给Rules类。Rules类创建相应的标签评估器。
接下来，我们将evaluateRow应用于这一组输入行，生成一个输出行的列表。
调用evaluateRow会修改布隆过滤器（一些标签将被存储在布隆中）。因此，我们使用state.update(bloom)来持久化布隆。
最后，我们返回一个输出行的迭代器。
应用 Flux Capacitor
从用户的角度来看，使用 Flux Capacitor 进行异常检测是非常简单的。
假设用户已经创建了一个 dataframe taggedEventsDF。在后续的文章中，我们将展示如何利用 sigma 编译器创建这个 dataframe。
val taggedEventsDF = ...
val specification = """
rules:
- rulename: rule1
description: ordered list of events by host
action: temporal
ordered: true
tags:
- name: recon_cmd_a
- name: recon_cmd_b
- name: recon_cmd_c
- rulename: rule3
description: propagate to all decendents
action: ancestor
tags:
- name: ancestor_process_feature1
parent: parent_id
child: id
"""
val tagCapacity = 200000
val flux = new FluxCapacitorMapFunction(tagCapacity, specification)
一旦创建了 FluxCapacitorMapFunction，它需要传递给 Spark 的flatMapGroupsWithState。用户还需要指定哪个列包含主机 ID。
// create encoders to serialize/deserialize the output rows and the bloom
val outputEncoder = RowEncoder(taggedEventsDF.schema).resolveAndBind()
val stateEncoder = Encoders.javaSerialization(BloomFilter.class)
// we associate a bloom with each host
var groupKeyIndex = df.schema.fieldIndex("host_id")
taggedEventsDF
.groupByKey(row => row.getString(groupKeyIndex))
.flatMapGroupsWithState(
Append,
GroupStateTimeout.NoTimeout()
)(flux.processBatch)(stateEncoder, outputEncoder)
我们为每个主机关联一个布隆过滤器，因此我们的分组键groupByKey是host_id列。
Spark 需要知道如何序列化和反序列化状态。我们为输出行和布隆过滤器（状态）创建 Spark 编码器。
我们指定了一种附加模式，而不是更新或完成。在附加模式下，我们输出的每一行都被添加到结果表中，这在异常检测中是所期望的。
最后，我们将我们的flux.processBatch传递给flatMapGroupsWithState函数。
结论
在本文中，我们展示了如何实现一个可配置且可重用的有状态映射函数，能够处理多种用例：时间邻近相关性（有序和无序）以及父子关系（包括祖先）。
目前，我们的状态永远不会超时，我们始终使用相同的布隆过滤器，因此布隆过滤器最终会被填满。在我们的下一篇文章中，我们将用一个会忘记的布隆过滤器来解决这个问题，并优化 Flux Capacitor 的性能。
除非另有说明，所有图片均由作者提供
使用 Sigma 规则进行异常检测（第五部分）：Flux Capacitor 优化
为了提升性能，我们实现了一种遗忘布隆过滤器和一个定制的 Spark 状态存储提供者
·
关注 发表于 Towards Data Science ·8 min read·Mar 17, 2023
–
照片来自 Shippagan, NB, Canada 的 Leora Winter，Unsplash
这是我们系列文章的第 5 篇。请参阅第一部分，第二部分，第三部分和第四部分获取一些背景信息。
在我们之前的文章中，我们展示了使用布隆过滤器所获得的性能提升。我们还展示了如何利用布隆过滤器实现时间接近相关性、父子和祖先关系。
到目前为止，我们一直在每个主机上使用一个布隆过滤器。最终，布隆过滤器将被标签填满，并会产生大量假阳性。使用这个在线布隆过滤器计算器，我们可以看到获得假阳性的概率。注意到假阳性率在超过 200,000 个标签后迅速增加。（这个图表的 n=200,000 和 p=1%）
图片由作者提供
健忘布隆过滤器
我们需要一种方法来处理非常旧的标签。我们需要一个健忘的布隆过滤器。正如 Redis Labs 在这篇优秀论文中解释的Age-Partitioned Bloom Filter，有许多方法可以实现健忘的布隆过滤器。我们将使用最基本的方法：
基于分段的方法使用几个不相交的分段，这些分段可以单独添加和退役。最简单且多次提到的方法是使用一系列普通的布隆过滤器，每代一个，当使用中的布隆过滤器满了时，添加一个新的并退役最旧的一个。
我们选择使用 10 代。因此，每台主机使用 10 个布隆过滤器。每个布隆过滤器最多可以容纳 20,000 个标签。
我们使用“活动”布隆过滤器来插入新标签。当“活动”布隆过滤器满了时，我们创建一个新的。当我们达到 10 个布隆过滤器时，我们丢弃最旧的布隆过滤器。
我们通过测试“活动”布隆过滤器来查询标签。如果没有找到标签，我们测试下一个（更旧的）布隆过滤器，直到到达末尾。
请注意，对于每个我们想要测试的标签，我们可能会在 10 个不同的布隆过滤器中执行 10 次测试。每次测试都有一定的概率报告假阳性。因此，通过使用 10 个布隆过滤器，我们将机会提高了 10 倍。为了降低假阳性的几率，我们使用了假阳性率为 1/1000 的布隆过滤器，而不是 1/100 的。实际上，我们将展示我们甚至可以使用 1/10000 的假阳性率。
为了适应多个布隆过滤器，我们不再在状态存储中保存布隆对象：
val stateEncoder = Encoders.javaSerialization(BloomFilter.class)
相反，我们将持久化一个 FluxState 对象，其中包含一个布隆过滤器列表：
val stateEncoder = Encoders.product[FluxState]
FluxState 包含以下字段：
case class FluxState(
var version: Long = 0,
var active: Int = 0,
var serializedBlooms: List[Array[Byte]] = List()
) extends TagCache {
出于性能原因，我们自己序列化 bloom 过滤器。由于我们知道这些对象的大小，我们可以通过预分配序列化缓冲区来优化序列化。serializedBlooms 字段保存序列化的 blooms。active 字段跟踪此列表中活动 bloom 的索引。我们稍后会解释版本号的使用。这就是我们序列化 blooms 的方式：
val padding = 4000
val n = tagCapacity / NUM_BLOOMS
// Formula taken from https://hur.st/bloomfilter
// m = ceil( (n * log(p)) / log(1 / pow(2, log(2))))
val mBits = Math.ceil(
(n * Math.log(desiredFpp)) / Math.log(1 / Math.pow(2, Math.log(2))))
val numBytes = (mBits / 8).toInt + padding
val byteArrayOut = new ByteArrayOutputStream(numBytes)
val store = new ObjectOutputStream(byteArrayOut)
store.writeObject(bloom)
store.close
byteArrayOut.toByteArray()
高效的检查点
我们将大型 bloom 分割成 10 个较小的 bloom。由于 bloom 过滤器的特性，10 个 20,000 标签的 bloom 所占用的空间大致与一个 200,000 标签的大 bloom 相同，大约为 200KiB。
Spark HDFS 状态存储提供者将所有 FluxState 对象保存在内存中。如果我们假设有一个 50,000 主机的集群，这将导致大约 10GiB 的 RAM 实际上，HDFS 状态存储的内存使用量被测量为 25GiB。
图片来源：作者
它的使用量更高的原因是 HDFS 状态存储默认会保留状态的 2 份副本。我们可以通过 spark.sql.streaming.maxBatchesToRetainInMemory 将其更改为只存储一份副本。这将内存使用量降低到大约 12GiB 的 RAM，这与我们的估计相符。
作为检查点的一部分，Spark 会将所有状态写入数据湖，并在每个微批处理完成后执行此操作。Spark 花费大量时间来持久化 12 GiB 的状态，并且反复进行此操作。
然而，在每个微批处理中，我们只修改 10 个 blooms 中的 1 个（活动 bloom）。其他 9 个 blooms 可能会被查询但保持不变。默认的 HDFS 状态存储提供者并不知晓哪个 bloom 被更改，它只是持久化 FluxState 对象。如果状态存储提供者知道哪个 bloom 是活动 bloom，它可以更高效，只检查点修改过的活动 bloom。这可能会将序列化减少到 12GiB 的 1/10。
自定义状态存储提供者
HDFSBackedStateStoreProvider.scala 类处理 put 和 get 状态请求。它将这些键/值对保存在内存映射中，并将这些键/值持久化到数据湖中。
通过对 put 和 get 方法进行一些轻微的修改，我们可以专门化 HDFSBackedStateStoreProvider 类的行为，使其能够意识到我们遗忘 bloom 过滤器策略。
基本的想法是将每个 bloom 段存储在一个单独的键下。我们不将整个状态存储在键“windows_host_abc”下，而是将每个段存储在“windows_host_abc_segment1”、“windows_host_abc_segment2”、“windows_host_abc_segment3”等下。
put函数接收要存储的键和值。键将是主机 ID，值将是一个 FluxState 对象。Spark 在调用此方法之前，将键和值都编码为 UnsafeRow：
override def put(key: UnsafeRow, value: UnsafeRow): Unit = {
require(value != null, "Cannot put a null value")
verify(state == UPDATING, "Cannot put after already committed or aborted")
//val keyCopy = key.copy()
val valueCopy = value.copy()
// Create the key to store this bloom using the bloom id (active)
val newRowKey = makeRowKey(key.getString(0), getActiveValue(valueCopy))
// Store key/value as normal
mapToUpdate.put(newRowKey, valueCopy)
writeUpdateToDeltaFile(compressedStream, newRowKey, valueCopy)
}
我们的put函数与原始的完全相同。唯一的变化是键。我们将“活跃的”布隆过滤器索引附加到原始键上。
需要注意的是，我们还修改了我们的 FluxState 类，只序列化“活跃的”布隆过滤器，而不是所有 10 个布隆过滤器。
def toState(): FluxState = {
version += 1
// we only store the active blooms, all other blooms were un-changed
val approxsize = tagCapacity / NUM_BLOOMS * 3
val byteArrayOut = new ByteArrayOutputStream(approxsize)
val store = new ObjectOutputStream(byteArrayOut)
store.writeObject(getActiveBloom())
store.close
serializedBlooms = List(byteArrayOut.toByteArray())
}
原始的get方法是这样的：
override def get(key: UnsafeRow): UnsafeRow = {
mapToUpdate.get(key)
}
我们修改了get方法以收集主机 ID 的 10 个布隆过滤器段。首先，我们通过迭代 10 个布隆索引来构建一个 FluxState 列表。然后，我们创建一个新的 FluxState 来保存所有的布隆过滤器。我们通过版本号来确定哪个是活跃的布隆过滤器。
override def get(key: UnsafeRow): UnsafeRow = {
// CCCS: The list of blooms for this key are put
// in separate entries in the store
// we will find all of these entries and create a FluxState
val groupKey = key.getString(0)
val fluxStates = Range(0, NUM_BLOOMS)
.map(bloomIndex => mapToUpdate.get(makeRowKey(groupKey, bloomIndex)))
.filter(_ != null)
.map(value => deserializeFluxState(value))
if(fluxStates.length > 0) {
makeRowValue(coalesceFluxStates(fluxStates))
}
else {
// else we found none of the blooms
null
}
}
调整假阳性概率
现在我们将布隆过滤器分成了 10 部分，我们可以查询 10 个布隆过滤器，因此可能会有更多的假阳性。为了解决这个问题，我们将假阳性概率降低到 1/10000，从而使总体假阳性概率为 1/1000。这比我们之前的实验少了十倍的假阳性机会。然而，由于我们只序列化“活跃的”布隆过滤器，总体性能要好得多。
结果
以前，当我们序列化所有段时，我们可以在每秒 5,000 个事件的速度下实现每主机 100,000 个标签的容量。
采用分段的方法，其中我们只序列化活跃的布隆过滤器，我们可以在每秒 5,000 个事件的速度下实现每主机 300,000 个标签的容量。或者，我们可以减小布隆过滤器的大小以容纳更多的事件：200,000 个标签 @ 8,000 个事件每秒。
作者提供的图片
在所有以前的实验中，我们一直在用随机标签对新创建的布隆过滤器进行“填充”。我们这样做是为了防止 HDFSBackStore 在将其状态保存到数据湖时压缩布隆过滤器。一个空的布隆过滤器几乎会压缩到零，而一个容量满的布隆过滤器（具有最大熵）几乎是不可压缩的。当我们首次启动实验时，由于压缩效果惊人，性能表现非常好。要看到标签在布隆过滤器中的效果需要很长时间。为了解决这个问题，我们将所有布隆过滤器填充到 95%的容量。换句话说，我们一直在测量最坏情况。
然而，在实际操作中，布隆过滤器会慢慢填满。一些布隆过滤器的填充速度比其他的快。从统计上讲，我们不可能在同一时间点上所有 50,000 个布隆过滤器都达到 95%的容量。使用随机填充可以进行更现实的模拟。
def createBloom() = {
val bloomCapacity: Int = tagCapacity / NUM_BLOOMS
val bloom = BloomFilter.create(
Funnels.stringFunnel(),
bloomCapacity,
desiredFpp)
val prep = (bloomCapacity * Math.random()).toInt
(1 to prep).map("padding" + _).foreach(s => bloom.put(s))
bloom
}
由于可压缩的布隆过滤器，我们可以以每秒 10,000 个事件的速度运行该模拟，并且每个主机的总体容量为 400,000 个标签，总计 20 亿个标签在单个 Spark 工作节点上。这远远超过了我们在流流连接中能够实现的 1 亿个标签。
从 bloom 中存储和检索标签非常迅速。平均而言，一台机器可以每秒进行约 200,000 次测试。将标签存储在 bloom 中成本稍高，但一台机器仍能每秒存储 20,000 个标签。这意味着我们可以同时支持大量 Sigma 规则。
这里是我们实验中使用的不同策略的总结。
作者提供的图像
结论
结果清晰地显示了与自定义状态存储结合使用的bloom 策略的性能改进，该存储仅保存“活跃”的 bloom 段。bloom 策略还比流-流连接方法更通用，因为它可以处理如“祖先”和时间接近性（有序和无序）等用例。概念验证可以在这里找到github.com/cccs-jc/flux-capacitor。如果你对这个 flux-capacitor 功能有新的用例，我们很乐意听到。
异常根本原因分析 101
如何找到每个指标异常的解释
·
关注 发表在 Towards Data Science ·14 min 阅读·2023 年 6 月 28 日
–
照片由 Markus Winkler 提供，发布于 Unsplash
我们使用指标和关键绩效指标（KPIs）来监控产品的健康状况：确保一切稳定或产品按预期增长。但有时，指标会突然变化。转化率可能在一天内上升 10%，或收入可能在几个季度内略微下降。在这种情况下，企业不仅需要了解发生了什么，还需要理解原因以及应采取的措施。这时分析师就发挥了重要作用。
我第一次从事的数据分析角色是 KPI 分析师。异常检测和根本原因分析是我近三年来的主要关注点。我发现了数十个 KPI 变化的关键驱动因素，并制定了一种处理这些任务的方法。
在这篇文章中，我想与大家分享我的经验。这样，下次你遇到意外的指标行为时，你将有一个可以遵循的指南。
应该关注什么？
在进行分析之前，让我们定义我们的主要目标：我们希望实现什么。那么，我们的异常根本原因分析的目的是什么？
最直接的答案是理解指标变化的关键驱动因素。不言而喻，从分析师的角度来看，这是一个正确的答案。
但我们从商业角度来看。花费资源进行此研究的主要原因是为了最小化对客户的潜在负面影响。例如，如果由于昨天发布的新版本应用中的一个漏洞导致转化率下降，那么今天找出这个问题会比一个月后当数百名客户已经流失时要好。
我们的主要目标是最小化对客户的潜在负面影响。
作为分析师，我喜欢为我的工作任务设置优化指标。最小化潜在不利影响听起来像是一种合适的心态，帮助我们专注于正确的事情。
因此，牢记主要目标，我会尝试找出以下问题的答案：
-
这是一个真正影响客户行为的问题，还是仅仅是一个数据问题？
-
如果我们的客户行为确实发生了变化，我们能做些什么？不同选项的潜在效果是什么？
-
如果这是一个数据问题，我们能否使用其他工具监控相同的过程？我们如何修复损坏的过程？
步骤 1：自己动手做
根据我的经验，最好的第一步是重现受影响的客户旅程。例如，假设 iOS 上的电商应用订单数量减少了 10%。在这种情况下，值得尝试购买一些东西，并仔细检查是否存在产品问题：按钮不可见、横幅无法关闭等。
另外，记得查看日志以确保信息正确捕捉。客户体验可能一切正常，但我们可能会丢失关于购买的数据。
我相信这是开始异常调查的一个重要步骤。首先，经过 DIY，你将更好地理解受影响的客户旅程部分：步骤是什么，数据如何记录。其次，你可能会找到根本原因，从而节省几个小时的分析时间。
提示： 如果异常幅度较大，更可能重现问题，这意味着问题影响了许多客户。
步骤 2：检查数据
正如我们之前讨论的，首先，了解客户是否受到影响，或者只是数据异常，这一点非常重要。
我强烈建议你检查数据是否是最新的。你可能会看到昨天的收入减少了 50%，因为报告只捕捉了当天的前半部分。你可以查看原始数据或与数据工程团队沟通。
如果没有已知的数据相关问题，你可以使用不同的数据源来双重检查指标。在许多情况下，产品有客户端数据（例如，Google Analytics 或 Amplitude）和后端数据（例如，应用日志、访问日志或 API 网关日志）。因此，我们可以使用不同的数据源来验证 KPI 动态。如果你在一个数据源中看到异常，那么你的问题可能与数据有关，并不会影响客户。
另一点需要记住的是时间窗口和数据延迟。有一次，一位产品经理找我说激活出现了问题，因为从注册到首次成功操作（例如在电子商务中是购买）的转化率在过去三周里一直在下降。然而，这实际上是一种日常情况。
作者基于合成数据的示例
下降的根本原因是时间窗口。我们跟踪注册后前 30 天内的激活情况。因此，注册超过 4 周的用户有整整一个月的时间来进行首次操作。但最近一批用户只有一周的时间来转化，因此他们的转化率预期会低得多。如果你想比较这些用户群体的转化情况，可以将时间窗口更改为一周或等待。
如果数据出现延迟，你可能会看到最近几天有类似的下降趋势。例如，我们的移动分析系统在设备使用 Wi-Fi 网络时会批量发送事件。因此，平均而言，从所有设备收集所有事件需要 3 到 4 天。所以在过去的 3 到 4 天看到较少的活跃设备是很正常的。
对于这种情况，好的做法是从图表中去掉最后一个句点。这将防止你的团队基于数据做出错误的决策。然而，人们仍然可能会不小心遇到这些不准确的指标，因此在深入进行根本原因分析之前，你应该花些时间理解这些方法论上准确的指标。
第三步：全景视图
下一步是更全面地查看趋势。首先，我倾向于缩小视野，查看较长时间的趋势，以获得整体图景。
例如，看看购买数量。订单数量已经稳定增长了好几个星期，预计在 12 月底（圣诞节和新年期间）会有所下降。但接着，在 5 月初，KPI 显著下降并持续减少。我们应该开始感到恐慌了吗？
作者基于合成数据的示例
实际上，大多数情况下，没有必要惊慌。我们可以查看过去三年的指标趋势，并注意到每年夏季购买数量都会减少。这是季节性因素的影响。许多产品在夏季的参与度较低，因为客户会去度假。然而，这种季节性模式并不是普遍存在的：例如，旅游或夏季节庆网站可能会有相反的季节性趋势。
作者基于合成数据的示例
再看一个例子——另一个产品的活跃客户数量。我们可以看到自六月以来有所下降：月活跃用户曾经为 380K — 400K，现在只有 340–360K（下降约 -10%）。我们已经检查过以往几年的夏季没有发生这样的变化。我们是否应该得出结论认为我们的产品出现了问题？
作者基于合成数据的示例
等等，还不行。在这种情况下，放大查看也可能有帮助。考虑到长期趋势，我们可以看到最近三周的数值接近于二月和三月的数值。真正的异常是从四月初到五月中旬的 1.5 个月高客户数量。我们可能错误地得出 KPI 下降的结论，但它只是回到了正常水平。考虑到那是 2020 年春季，网站上的高流量很可能是由于 COVID 隔离：客户在家里，在线时间增加。
作者基于合成数据的示例
你初步分析的最后一点是确定 KPI 变化的确切时间。在某些情况下，变化可能会在 5 分钟内突然发生。而在其他情况下，可能只是趋势的微小变化。例如，活跃用户曾经每周增长 +5%，但现在仅为 +3%。
尝试尽可能准确地确定变化点（甚至到分钟级别）是值得的，因为这将帮助你后续选择最可信的假设。
指标变化的速度可以给你一些线索。例如，如果转化率在 5 分钟内发生变化，那不可能是由于新应用版本的发布（通常需要几天时间才能更新应用），更可能是由于后端的变化（例如 API）。
步骤 4：获取背景信息
理解整个背景（发生了什么）对我们的调查可能至关重要。
我通常检查以查看整体情况的内容：
-
内部变化。不用说，内部变化会影响 KPI，因此我通常会查看所有的发布、实验、基础设施事件、产品变化（例如，新设计或价格变动）以及供应商更新（例如，升级到我们用于报告的最新版本 BI 工具）。
-
外部因素可能会根据你的产品有所不同。在金融科技领域，货币汇率可能会影响客户行为，而重大新闻或天气变化可能会影响搜索引擎市场份额。你可以为你的产品头脑风暴类似的因素。试着在思考外部因素时富有创意。例如，我们曾发现网站流量的下降是由于我们最重要地区的网络问题。
-
竞争对手活动。尝试了解你的主要竞争对手是否正在做某些事情——大规模的营销活动、产品不可用的事件或市场关闭。最简单的方法是查看 Twitter、Reddit 或新闻上的提及。此外，还有许多监控服务问题和故障的网站（例如，DownDetector或DownForEveryoneOrJustMe），你可以在这些网站上检查竞争对手的健康状况。
-
客户反馈。你可以通过客户支持团队了解产品的问题。因此，不要犹豫，问问他们是否有新的投诉或特定类型的客户联系量增加。然而，请记住，可能只有少数人会联系客户支持（特别是如果你的产品对日常生活并非必需）。例如，多年前，我们的搜索引擎在~100K 老版本 Opera 浏览器用户中完全崩溃。问题持续了几天，但不到十个客户联系了支持。
由于我们已经定义了异常时间，因此很容易获取所有发生在附近的事件。这些事件就是你的假设。
提示： 如果你怀疑内部变更（发布或实验）是 KPI 下降的根本原因。最佳实践是恢复这些变更（如果可能的话），然后尝试理解确切的问题。这将帮助你减少对客户的潜在负面影响。
步骤 5：切片与切块
此时，你希望已经对异常发生时的周围情况有所了解，并对根本原因有了一些假设。
我们先从更高的层次看待异常。例如，如果在美国客户的 Android 设备上出现了异常，值得检查 iOS、网页以及其他地区的客户。这样你将能够适当地理解问题的规模。
之后，是时候深入挖掘并尝试定位异常（尽可能狭窄地定义受 KPI 变化影响的段或多个段）。最直接的方法是查看你产品在不同维度上的 KPI 趋势。
这样有意义的维度列表可能会根据你的产品有所不同，所以值得与团队进行头脑风暴。我建议查看以下因素组：
-
技术特性：例如，平台、操作系统、应用版本；
-
客户特征：例如，新客户或现有客户（群体），年龄，地区；
-
客户行为：例如，产品功能的采用，实验标志，营销渠道。
在按不同维度拆分 KPI 趋势时，最好只查看足够显著的细分。例如，如果收入下降了 10%，没有必要查看对总收入贡献不到 1%的国家。小组中的指标往往更具波动性，因此不显著的细分可能会增加过多噪音。我更倾向于将所有小片段分组到其他组中，以避免完全丧失该信号。
例如，我们可以查看按平台划分的收入。不同平台的绝对数值可能差异很大，因此我将所有系列标准化到第一个点，以比较随时间的动态。有时，最好对前 N 个点进行平均化。例如，将前七天的平均值来捕捉每周季节性。
这就是你可以在 Python 中做到的。
import plotly.express as px
norm_value = df[:7].mean()
norm_df = df.apply(lambda x: x/norm_value, axis = 1)
px.line(norm_df, title = 'Revenue by platform normed on 1st point')
图表告诉我们整个故事：在 5 月之前，不同平台的收入趋势非常接近，但随后 iOS 出现了变化，iOS 收入下降了 10-20%。所以 iOS 平台主要受到这一变化的影响，而其他平台则相对稳定。
基于合成数据的作者示例
步骤 6：理解你的指标
在确定受到异常影响的主要细分后，让我们尝试分解我们的 KPI。这可能会给我们更好的了解情况。
我们通常在分析中使用两种类型的 KPI：绝对数值和比率。所以让我们讨论每种情况的分解方法。
我们可以通过标准化来分解绝对数值。例如，让我们看看在服务中总共花费的时间（内容产品的标准 KPI）。我们可以将其分解为两个单独的指标。
然后我们可以查看两个指标的动态。在下面的例子中，我们可以看到活跃客户数保持稳定，而每位客户花费的时间下降了，这意味着我们并没有完全失去客户，但由于某些原因，他们开始在我们的服务上花费更少的时间。
基于合成数据的作者示例
对于比率指标，我们可以分别查看分子和分母的动态。例如，让我们使用 30 天内从注册到首次购买的转化率。我们可以将其分解为两个指标：
-
在注册后 30 天内完成购买的客户数量（分子），
-
注册人数（分母）。
在下面的例子中，转化率在四月份从 43.5%下降到 40%。注册人数和转化客户数都增加了。这意味着有更多的客户转化率较低。这可能是由于不同的原因：
-
新的营销渠道或质量较低的用户的营销活动；
-
数据中的技术变化（例如，我们改变了地区定义，现在我们考虑了更多的客户）；
-
网站上的欺诈或机器人流量。
作者基于合成数据的示例
提示： 如果我们看到转换用户的下降，而总用户数量保持稳定，这将表明产品或数据中关于转换的事实存在问题。
对于转换，将其转化为漏斗也可能会有所帮助。例如，在我们的案例中，我们可以查看以下步骤的转换：
-
完成注册
-
产品目录
-
将商品添加到购物车
-
下订单
-
成功支付。
每一步的转换动态可以显示客户旅程中发生变化的阶段。
第 7 步：得出结论
根据上述所有分析阶段，你应该对当前情况有一个比较全面的了解：
-
具体改变了什么；
-
哪些细分市场受到影响；
-
周围发生了什么。
现在是总结的时候了。我倾向于以结构化的方式记录所有信息，描述测试过的假设和我们得出的结论，以及当前对主要根本原因的理解和下一步（如果需要）。
提示： 记录所有测试过的假设（不仅仅是已证明的假设）是值得的，因为这将避免重复不必要的工作。
现在最重要的是验证我们的主要根本原因是否能够完全解释 KPI 变化。我通常在没有已知影响的情况下对情况进行建模。
例如，在从注册到首次购买的转换案例中，我们可能发现了欺诈攻击，并知道如何通过 IP 地址和用户代理来识别机器人流量。因此，我们可以查看在已知主要根本原因——欺诈流量——影响之外的转换率。
作者基于合成数据的示例
如你所见，欺诈流量仅解释了大约 70%的下降，还有其他因素可能影响 KPI。这就是为什么最好双重检查你是否发现了所有重要因素。
有时，证明你的假设可能是具有挑战性的，例如价格或设计的变化，你无法进行适当的 A/B 测试。我们都知道，相关性并不意味着因果关系。
在这种情况下检查假设的可能方法：
-
查看过去类似情况，例如价格变化是否与 KPI 有类似的相关性。
-
尝试识别行为发生变化的客户，比如那些开始在我们的应用中花费更少时间的客户，并进行调查。
在此分析之后，你仍然可能对效果产生疑问，但这可能会增加你对找到正确答案的信心。
提示： 如果你陷入困境，调查也可能有帮助：你已经检查了所有假设但仍未找到解释。
如何为下一次根本原因分析做好准备？
在详细调查结束时，是时候考虑如何让下次变得更轻松、更好。
经过多年的异常调查经验，我总结了以下最佳实践：
-
拥有特定于您的产品的清单是非常有帮助的——这可以为您和您的同事节省数小时的工作。值得整理出假设列表和检查这些假设的工具（例如仪表盘链接、关于竞争对手的外部信息来源等）。请记住，编写清单不是一次性活动：一旦遇到新的异常类型，您应该将新知识添加到清单中，以保持其最新。
-
另一个有价值的文档是包含所有重要事件的变更日志，例如价格变动、竞争产品的推出或新功能发布。变更日志将帮助您在一个地方找到所有重要事件，而无需翻阅多个聊天记录和维基页面。记住更新变更日志可能会很有挑战性。您可以将其纳入分析值班职责，以确立明确的责任。
-
在大多数情况下，您需要来自不同人员的输入，以了解情况的整体背景。预先准备好的工作组和 KPI 异常调查渠道可以节省宝贵的时间，并保持所有利益相关者的更新。
-
最后但同样重要的是，为了最小化对客户的潜在负面影响，我们应该建立监控系统，以便尽快了解异常情况并开始寻找根本原因。因此，请腾出时间来建立和改进您的警报和监控系统。
TL;DR
我希望您牢记的关键消息：
-
处理根本原因分析时，您应该专注于最小化对客户的潜在负面影响。
-
尽量发挥创造力并广泛考虑：了解您产品内部发生的情况、基础设施的状况以及潜在的外部因素。
-
深入挖掘：从不同角度查看您的指标，尝试检查不同的细分市场并分解您的指标。
-
做好准备：如果您已经有产品清单、变更日志和工作组进行头脑风暴，那么处理这样的研究会容易得多。
非常感谢您阅读这篇文章。我希望现在您不会在面对根本原因分析任务时感到困惑，因为您已经有了手头的指南。如果您有任何后续问题或评论，请随时在评论区留言。
另一种（符合性）预测概率分布的方法
原文：
towardsdatascience.com/another-conformal-way-to-predict-probability-distributions-fcc63e78680d
使用 Catboost 进行符合性多分位回归
·发表于 Towards Data Science ·阅读时间 11 分钟·2023 年 3 月 8 日
–
德克萨斯州。图像来源：作者。
在 上一篇文章 中，我们探索了 Catboost 的多分位损失函数的能力，该函数允许使用单一模型预测多个分位数。这种方法优雅地克服了传统分位回归的一项限制：后者需要为每个分位数开发一个单独的模型，或将整个训练集存储在模型中。然而，分位回归还有另一个缺点，我们将在本文中讨论：预测的分位数可能存在偏差，无法保证校准和覆盖。本文将演示如何通过符合性多分位回归来克服这一问题。我鼓励尚未跟进本系列的人在阅读之前回顾以下文章：
探索 Catboost 的多分位回归
towardsdatascience.com [## 理解机器学习中的噪声数据和不确定性
实际原因是你的机器学习模型未能正常工作
## 如何使用“符合分位数回归”预测风险比例区间 如何预测风险比例区间与符合分位数回归
这个算法——由斯坦福学者在 2019 年发布——将分位数回归与符合预测结合在一起。这里…
回顾：为何选择多分位数回归？
多分位数回归使我们能够使用一个模型预测多个目标分位数。因为没有计算限制要求每个分位数一个模型，也没有像 KNN 或分位数回归森林那样需要在模型中存储整个训练集的限制，我们可以更有效地预测更多分位数，并更好地了解条件目标分布的样貌。
使用传统的分位数回归，生成 95%的预测区间需要一个用于 2.5 分位数的模型，一个用于 97.5 分位数的模型，可能还需要一个用于期望值或 50 分位数的模型。每个模型的单次预测结果大致如下：
预测的 CDF 样本用于单个测试样本（三个独立的分位数模型）。图像由作者提供。
假设这些分位数经过校准，它们揭示了一些见解。首先，给定特征，目标小于或等于 3.6 的概率大约是 0.50 或 50%。类似地，给定特征，目标值在 3.25 和 4.38 之间的概率大约是 0.95 或 95%。
虽然模型的输出很好且符合我们的要求，但我们可能希望动态调整风险容忍度。例如，如果我们需要更保守的 99%预测区间怎么办？类似地，如果我们更愿意承担风险，接受 90%或 80%预测区间怎么办？如果我们想知道“给定特征，目标大于 y1 的概率是多少？”我们可能还想问“给定特征，目标在 y1 和 y2 之间的概率是多少？”多分位数回归通过预测尽可能多的分位数来帮助回答这些问题：
预测的 CDF 样本用于单个测试样本（一个多分位数模型）。图像由作者提供。
能够准确预测的分位数越多，风险容忍度可以随时调整的余地就越大，我们可以更好地回答关于条件目标分布的一般概率问题。
请注意，单一决策树模型已被用于生成多个分位数预测。然而，这依赖于树将所有目标值存储在叶子节点中。在预测时，指定一个分位数，并从叶子节点中的数据中经验性地计算出来，这要求模型存储整个训练集。这也意味着深度树可能在叶子节点中只有很少的样本可供使用。
Catboost 本质上有所不同，因为它仅在终端节点中存储指定的分位数的数量。此外，损失函数被优化以预测每个指定的分位数。我们还享受了 Catboost 提供的性能提升，这是其底层架构带来的优势。
分位回归的问题
在传统的和多分位回归中，没有总是有统计保证 分位数是无偏的。这意味着，对于训练来预测目标分布的第 95 分位数的模型，并不能保证 95%的观察值实际上会小于或等于预测值。这在需要准确概率表示来做出关键决策的高风险应用中是一个问题。
分位回归也可能产生过于保守的预测区间，进而导致信息不足。一般来说，预测区间应尽可能窄，同时保持所需的覆盖水平。
符合性多分位回归
符合性分位回归的想法是调整预测的分位数，以准确反映所需的风险容忍度和区间长度。这是通过一个“校准”步骤来完成的，该步骤计算“符合性得分”以纠正预测的分位数。有关符合性分位回归的更多细节可以在这篇论文和这篇文章中找到。对于符合性多分位回归，我们将利用以下定理：
左尾和右尾符合性分位回归。来源。
如果这看起来过于抽象，请不要担心，步骤实际上很简单：
-
创建训练集、校准集和测试集。在训练集上拟合多分位模型以预测所有感兴趣的分位数。
-
在校准集上进行预测。对于每个校准实例和预测的分位数，计算预测分位数与对应目标值之间的差异。这些就是一致性分数。
-
对于每个测试示例和预测的分位数（假设为 q），从模型预测的分位数中减去对应于分位数 q 的一致性分数的 1-q 分位数。这些就是新的预测分位数。
我们可以在 Python 类中实现这个逻辑：
import numpy as np
import pandas as pd
from catboost import CatBoostRegressor, CatBoostError
from typing import Iterable
class ConformalMultiQuantile(CatBoostRegressor):
def __init__(self, quantiles:Iterable[float], *args, **kwargs):
"""
Initialize a ConformalMultiQuantile object.
Parameters
----------
quantiles : Iterable[float]
The list of quantiles to use in multi-quantile regression.
*args
Variable length argument list.
**kwargs
Arbitrary keyword arguments.
"""
kwargs['loss_function'] = self.create_loss_function_str(quantiles)
super().__init__(*args, **kwargs)
self.quantiles = quantiles
self.calibration_adjustments = None
@staticmethod
def create_loss_function_str(quantiles:Iterable[float]):
"""
Format the quantiles as a string for Catboost
Paramters
---------
quantiles : Union[float, List[float]]
A float or list of float quantiles
Returns
-------
The loss function definition for multi-quantile regression
"""
quantile_str = str(quantiles).replace('[','').replace(']','')
return f'MultiQuantile:alpha={quantile_str}'
def calibrate(self, x_cal, y_cal):
"""
Calibrate the multi-quantile model
Paramters
---------
x_cal : ndarray
Calibration inputs
y_cal : ndarray
Calibration target
"""
# Ensure the model is fitted
if not self.is_fitted():
raise CatBoostError('There is no trained model to use calibrate(). Use fit() to train model. Then use this method.')
# Make predictions on the calibration set
uncalibrated_preds = self.predict(x_cal)
# Compute the difference between the uncalibrated predicted quantiles and the target
conformity_scores = uncalibrated_preds - np.array(y_cal).reshape(-1, 1)
# Store the 1-q quantile of the conformity scores
self.calibration_adjustments = \
np.array([np.quantile(conformity_scores[:,i], 1-q) for i,q in enumerate(self.quantiles)])
def predict(self, data, prediction_type=None, ntree_start=0, ntree_end=0, thread_count=-1, verbose=None, task_type="CPU"):
"""
Predict using the trained model.
Parameters
----------
data : pandas.DataFrame or numpy.ndarray
Data to make predictions on
prediction_type : str, optional
Type of prediction result, by default None
ntree_start : int, optional
Number of trees to start prediction from, by default 0
ntree_end : int, optional
Number of trees to end prediction at, by default 0
thread_count : int, optional
Number of parallel threads to use, by default -1
verbose : bool or int, optional
Verbosity, by default None
task_type : str, optional
Type of task, by default "CPU"
Returns
-------
numpy.ndarray
The predicted values for the input data.
"""
preds = super().predict(data, prediction_type, ntree_start, ntree_end, thread_count, verbose, task_type)
# Adjust the predicted quantiles according to the quantiles of the
# conformity scores
if self.calibration_adjustments is not None:
preds = preds - self.calibration_adjustments
return preds
示例：超导性数据集
我们将对 超导性数据集 上的 超导体 数据集进行一致性多分位回归。该数据集提供了 21,263 个包含 81 个 超导体 特征及其 临界温度 （目标）。数据被划分为 ~64% 用于训练，~16% 用于校准，20% 用于测试。
# Dependencies
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from catboost import CatBoostRegressor, CatBoostError
from sklearn.model_selection import train_test_split
from typing import Iterable
pd.set_option('display.max.columns', None)
sns.set()
# Read in superconductivity dataset
data = pd.read_csv('train.csv')
# Predicting critical temperature
target = 'critical_temp'
# 80/20 train/test split
x_train, x_test, y_train, y_test = train_test_split(data.drop(target, axis=1), data[target], test_size=0.20)
# Hold out 20% of the training data for calibration
x_train, x_cal, y_train, y_cal = train_test_split(x_train, y_train, test_size=0.20)
print("Training shape:", x_train.shape) # Training shape: (13608, 81)
print("Calibration shape:", x_cal.shape) # Calibration shape: (3402, 81)
print("Testing shape:", x_test.shape) # Testing shape: (4253, 81)
我们将指定一组要预测的分位数。为了展示多分位回归的强大功能，该模型将预测从 0.005 到 0.99 的 200 个分位数——这在实际中可能有些过度。接下来，我们将 拟合 一致性多分位模型，进行 未经校准的预测，在 校准集 上 校准 模型，并进行 校准后的预测。
# Store quantiles 0.005 through 0.99 in a list
quantiles = [q/200 for q in range(1, 200)]
# Instantiate the conformal multi-quantile model
conformal_model = ConformalMultiQuantile(iterations=100,
quantiles=quantiles,
verbose=10)
# Fit the conformal multi-quantile model
conformal_model.fit(x_train, y_train)
# Get predictions before calibration
preds_uncalibrated = conformal_model.predict(x_test)
preds_uncalibrated = pd.DataFrame(preds_uncalibrated, columns=[f'pred_{q}' for q in quantiles])
# Calibrate the model
conformal_model.calibrate(x_cal, y_cal)
# Get calibrated predictions
preds_calibrated = conformal_model.predict(x_test)
preds_calibrated = pd.DataFrame(preds_calibrated, columns=[f'pred_{q}' for q in quantiles])
preds_calibrated.head()
结果预测应如下所示：
前五个观察值的前几个预测分位数。图片来源于作者。
在测试集上，我们可以测量未经校准和校准的预测如何与它们所表示的左尾概率对齐。例如，如果分位数已校准，则 40% 的目标值应小于或等于预测分位数 0.40，90% 的目标值应小于或等于预测分位数 0.90 等。下面的代码计算期望左尾概率与实际左尾概率之间的平均绝对误差（MAE）：
# Initialize an empty DataFrame
comparison_df = pd.DataFrame()
# For each predicted quantile
for i, quantile in enumerate(quantiles):
# Compute the proportion of testing observations that were less than or equal
# to the uncalibrated predicted quantile
actual_prob_uncal = np.mean(y_test.values <= preds_uncalibrated[f'pred_{quantile}'])
# Compute the proportion of testing observations that were less than or equal
# to the calibrated predicted quantile
actual_prob_cal = np.mean(y_test.values <= preds_calibrated[f'pred_{quantile}'])
comparison_df_curr = pd.DataFrame({
'desired_probability':quantile,
'actual_uncalibrated_probability':actual_prob_uncal,
'actual_calibrated_probability':actual_prob_cal}, index=[i])
comparison_df = pd.concat([comparison_df, comparison_df_curr])
comparison_df['abs_diff_uncal'] = (comparison_df['desired_probability'] - comparison_df['actual_uncalibrated_probability']).abs()
comparison_df['abs_diff_cal'] = (comparison_df['desired_probability'] - comparison_df['actual_calibrated_probability']).abs()
print("Uncalibrated quantile MAE:", comparison_df['abs_diff_uncal'].mean())
print("Calibrated quantile MAE:", comparison_df['abs_diff_cal'].mean())
# Uncalibrated quantile MAE: 0.02572999018133225
# Calibrated quantile MAE: 0.007850550660662823
未校准的分位数平均偏差约为 0.026，而校准后的分位数偏差为 0.008。因此，校准后的分位数与期望的左尾概率更加一致。
实际 vs 预测分位数的左尾概率。图片来源于作者。
这可能看起来不像校准中有戏剧性的变化，但通过分析实际与期望覆盖率，可以更清楚地看出未经校准模型中的误差：
coverage_df = pd.DataFrame()
for i, alpha in enumerate(np.arange(0.01, 0.41, 0.01)):
lower_quantile = round(alpha/2, 3)
upper_quantile = round(1 - alpha/2, 3)
# Compare actual to expected coverage for both models
lower_prob_uncal = comparison_df[comparison_df['desired_probability'] == lower_quantile]['actual_uncalibrated_probability'].values[0]
upper_prob_uncal = comparison_df[comparison_df['desired_probability'] == upper_quantile]['actual_uncalibrated_probability'].values[0]
lower_prob_cal = comparison_df[comparison_df['desired_probability'] == lower_quantile]['actual_calibrated_probability'].values[0]
upper_prob_cal = comparison_df[comparison_df['desired_probability'] == upper_quantile]['actual_calibrated_probability'].values[0]
coverage_df_curr = pd.DataFrame({'desired_coverage':1-alpha,
'actual_uncalibrated_coverage':upper_prob_uncal - lower_prob_uncal,
'actual_calibrated_coverage':upper_prob_cal - lower_prob_cal}, index=[i])
coverage_df = pd.concat([coverage_df, coverage_df_curr])
coverage_df['abs_diff_uncal'] = (coverage_df['desired_coverage'] - coverage_df['actual_uncalibrated_coverage']).abs()
coverage_df['abs_diff_cal'] = (coverage_df['desired_coverage'] - coverage_df['actual_calibrated_coverage']).abs()
print("Uncalibrated Coverage MAE:", coverage_df['abs_diff_uncal'].mean())
print("Calibrated Coverage MAE:", coverage_df['abs_diff_cal'].mean())
# Uncalibrated Coverage MAE: 0.03660674817775689
# Calibrated Coverage MAE: 0.003543616270867622
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(coverage_df['desired_coverage'],
coverage_df['desired_coverage'],
label='Perfect Calibration')
ax.scatter(coverage_df['desired_coverage'],
coverage_df['actual_uncalibrated_coverage'],
color='orange',
label='Uncalibrated Model')
ax.scatter(coverage_df['desired_coverage'],
coverage_df['actual_calibrated_coverage'],
color='green',
label='Calibrated Model')
ax.set_xlabel('Desired Coverage')
ax.set_ylabel('Actual Coverage')
ax.set_title('Desired vs Actual Coverage')
ax.legend()
plt.show()
实际覆盖率与期望覆盖率。图像来源：作者。
未校准模型往往过于保守，覆盖了比期望更多的示例。而校准模型则与每个期望的覆盖率几乎完美对齐。
此外，校准模型生成的预测区间的平均长度少于未校准模型。因此，校准模型的覆盖率更好，预测区间也更具信息性。
按期望覆盖率的平均预测区间长度。图像来源：作者。
可能会有人问，如果我们允许未校准模型将校准集作为训练数据会发生什么。这在实践中是合理的，因为我们不会无缘无故地丢弃好的训练数据。以下是结果：
# Fit a model using the training and calibration data
regular_model = ConformalMultiQuantile(iterations=100,
quantiles=quantiles,
verbose=10)
regular_model.fit(pd.concat([x_train, x_cal]), pd.concat([y_train, y_cal]))
# Fit a model on the training data only
conformal_model = ConformalMultiQuantile(iterations=100,
quantiles=quantiles,
verbose=10)
conformal_model.fit(x_train, y_train)
# Get predictions before calibration
preds_uncalibrated = regular_model.predict(x_test)
preds_uncalibrated = pd.DataFrame(preds_uncalibrated, columns=[f'pred_{q}' for q in quantiles])
# Calibrate the model
conformal_model.calibrate(x_cal, y_cal)
# Get calibrated predictions
preds_calibrated = conformal_model.predict(x_test)
preds_calibrated = pd.DataFrame(preds_calibrated, columns=[f'pred_{q}' for q in quantiles])
comparison_df = pd.DataFrame()
# Compare actual to predicted left-tailed probabilities
for i, quantile in enumerate(quantiles):
actual_prob_uncal = np.mean(y_test.values <= preds_uncalibrated[f'pred_{quantile}'])
actual_prob_cal = np.mean(y_test.values <= preds_calibrated[f'pred_{quantile}'])
comparison_df_curr = pd.DataFrame({
'desired_probability':quantile,
'actual_uncalibrated_probability':actual_prob_uncal,
'actual_calibrated_probability':actual_prob_cal}, index=[i])
comparison_df = pd.concat([comparison_df, comparison_df_curr])
comparison_df['abs_diff_uncal'] = (comparison_df['desired_probability'] - comparison_df['actual_uncalibrated_probability']).abs()
comparison_df['abs_diff_cal'] = (comparison_df['desired_probability'] - comparison_df['actual_calibrated_probability']).abs()
print("Uncalibrated quantile MAE:", comparison_df['abs_diff_uncal'].mean())
print("Calibrated quantile MAE:", comparison_df['abs_diff_cal'].mean())
# Uncalibrated quantile MAE: 0.023452756375340143
# Calibrated quantile MAE: 0.0061827359227361834
即使在训练数据少于未校准模型的情况下，校准模型输出的分位数也更好。更重要的是，当我们将预测的分位数的期望值与目标值进行比较时，这些模型的表现相似：
from sklearn.metrics import r2_score, mean_absolute_error
print(f"Uncalibrated R2 Score: {r2_score(y_test, preds_uncalibrated.mean(axis=1))}")
print(f"Calibrated R2 Score: {r2_score(y_test, preds_calibrated.mean(axis=1))} \n")
print(f"Uncalibrated MAE: {mean_absolute_error(y_test, preds_uncalibrated.mean(axis=1))}")
print(f"Calibrated MAE: {mean_absolute_error(y_test, preds_calibrated.mean(axis=1))} \n")
# Uncalibrated R2 Score: 0.8060126144892599
# Calibrated R2 Score: 0.8053382438575666
# Uncalibrated MAE: 10.622258046774979
# Calibrated MAE: 10.557269513856014
最终思考
机器学习中没有万无一失的解决方案，符合性分位回归也不例外。支撑符合性预测理论的粘合剂是数据的可交换性假设。例如，如果数据的分布随时间漂移（这在许多实际应用中通常发生），那么符合性预测无法再提供强有力的概率保证。虽然有绕过这一假设的方法，但这些方法最终取决于数据漂移的严重程度和学习问题的性质。将宝贵的训练数据用于校准也可能不是最佳选择。
一如既往，机器学习从业者负责理解数据的性质并应用适当的技术。感谢阅读！
成为会员： https://harrisonfhoffman.medium.com/membership
参考文献
-
Catboost 损失函数 —
catboost.ai/en/docs/concepts/loss-functions-regression#MultiQuantile -
符合性分位回归 —
arxiv.org/pdf/1905.03222.pdf -
符合性预测超越可交换性 —
arxiv.org/pdf/2202.13415.pdf -
超导数据集 —
archive.ics.uci.edu/ml/datasets/Superconductivty+Data -
如何使用保形分位回归预测风险比例区间 —
towardsdatascience.com/how-to-predict-risk-proportional-intervals-with-conformal-quantile-regression-175775840dc4 -
如何使用机器学习保形预测预测完整概率分布 —
valeman.medium.com/how-to-predict-full-probability-distribution-using-machine-learning-conformal-predictive-f8f4d805e420
蚁群优化算法的实际应用
原文：
towardsdatascience.com/ant-colony-optimization-in-action-6d9106de60af
一只滑雪的蚂蚁。图片由作者使用 Dall·E 创建。
使用 ACO 在 Python 中解决优化问题并提升结果
·发布于Towards Data Science ·10 分钟阅读·2023 年 9 月 20 日
–
欢迎回来！在我的 上一篇文章****中，我介绍了蚁群优化（ACO）的基本原理。在这一部分，我们将深入探讨如何从头开始实现 ACO 算法，以解决两种不同类型的问题。
我们将要解决的问题是旅行商问题（TSP）和二次分配问题（QAP）。为什么是这两个问题？因为 TSP 是一个经典挑战，而 ACO 恰好是一种有效的算法，用于找到图中最具成本效益的路径。另一方面，二次分配问题代表了与优化物品排列相关的不同问题类别，在这篇文章中，我旨在展示 ACO 也可以成为解决此类分配相关问题的有价值工具。这种多功能性使得 ACO 算法适用于广泛的问题。最后，我将分享一些更快速获得改进解决方案的技巧。
旅行商问题
TSP（旅行商问题）描述起来很简单，但在寻找解决方案时可能会面临重大挑战。基本定义是：你的任务是发现一条访问图中所有节点的最短路径。这个问题属于NP 难题的范畴，这意味着如果你尝试探索所有可能的路径，找到解决方案可能需要不切实际的时间。相反，更有效的方法是寻求在合理时间内的高质量解决方案，这正是我们通过 ACO（蚁群优化算法）要实现的目标。
问题定义
使用以下代码，我们可以创建一个具有给定节点数量的 TSP 实例：
import itertools
import math
import random
from typing import Tuple
import networkx as nx
import networkx.algorithms.shortest_paths.dense as nxalg
class TSP:
"""
Creates a TSP problem with a certain number of nodes
"""
def __init__(self, nodes: int = 30, dimensions: Tuple[int, int] = (1000, 1000), seed: int = 5):
if seed:
random.seed(seed)
graph = nx.Graph()
nodes_dict = dict()
for i in range(nodes):
nodes_dict[i] = (random.randint(0, dimensions[0]), random.randint(0, dimensions[1]))
graph.add_node(i)
for i, j in itertools.permutations(range(nodes), 2):
graph.add_edge(i, j, weight=self.calculate_distance(nodes_dict[i], nodes_dict[j]))
self.graph = graph
self.nodes = nodes_dict
self.distance_matrix = nxalg.floyd_warshall_numpy(graph)
@staticmethod
def calculate_distance(i, j):
"""
Calculate the Euclidian distance between two nodes
"""
return int(math.sqrt((i[0] - j[0])**2 + (i[1] - j[1])**2))
我们将用来演示 ACO 的 TSP 示例是以下的（默认设置）：
访问所有节点并返回起始节点。图片由作者提供。
对于这个问题的最优解（用混合整数规划计算）如下：
旅行推销员问题的最优解。图片由作者提供。
这条路径的距离是 4897。
用 ACO 解决 TSP
下一步是用蚁群优化解决这个问题，看看我们能接近最优解到什么程度。如果你不熟悉 ACO，想了解算法如何工作，可以阅读我之前的文章。然后，你可以回到这里查看 ACO 的实际应用。
ACO 的代码：
import datetime
import json
import logging
import random
import time
import matplotlib.pyplot as plt
import numpy as np
from problem import TSP
class AntColonyOptimization:
"""
Ant colony optimization algorithm for finding the shortest route in a graph.
Parameters:
m = number of ants
k_max = number of iterations
alpha = pheromone importance
beta = distance importance
rho = pheromone evaporation rate
Q = pheromone deposit
tau = pheromone
eta = distance
"""
def __init__(self, problem, **kwargs):
self.graph = problem.graph
self.nodes = list(problem.nodes)
self.coordinates = list(problem.nodes.values())
self.n = len(self.nodes)
self.distance_matrix = problem.distance_matrix
self.m = kwargs.get("m", 100)
self.k_max = kwargs.get("k_max", 50)
self.alpha = kwargs.get("alpha", 1)
self.beta = kwargs.get("beta", 5)
self.rho = kwargs.get("rho", 0.9)
self.Q = kwargs.get("Q", 1)
self.time_limit = kwargs.get("time_limit", 5)
# initialization of tau and eta
self.tau = np.full(self.distance_matrix.shape, 0.1)
self.eta = 1 / (self.distance_matrix + 1e-10)
self.history = []
def ant_colony_optimization(self):
"""
Ant colony optimization algorithm
"""
start_time = time.time()
x_best, y_best = [], float("inf")
for _ in range(self.k_max):
self.edge_attractiveness()
self.tau *= (1-self.rho)
for _ in range(self.m):
x_best, y_best = self.ant_walk(x_best, y_best)
if time.time() - start_time > self.time_limit:
logging.info("Time limit reached. Stopping ACO.")
return x_best, y_best
return x_best, y_best
def edge_attractiveness(self, plot: bool = False):
"""
Calculate edge attractiveness
tau = pheromone
eta = distance
alpha = pheromone importance
beta = distance importance
"""
self.A = (self.tau ** self.alpha) * (self.eta ** self.beta)
def ant_walk(self, x_best, y_best, plot: bool = True):
"""
Ant walk
"""
x = [0] # Start at first node
while len(x) < self.n:
i = x[-1]
neighbors = [j for j in range(self.n) if j not in x and self.distance_matrix[i][j] > 0]
if len(neighbors) == 0:
return x_best, y_best
p = [self.A[(i, j)] for j in neighbors]
sampled_neighbor = random.choices(neighbors, weights=p)[0]
x.append(sampled_neighbor)
x.append(0)
y = self.score(x)
self.history.append(y)
for i in range(1, self.n):
self.tau[(x[i-1], x[i])] += self.Q / y
if y < y_best:
logging.info("Better ACO solution found. Score: %.2f", y)
return x, y
return x_best, y_best
def score(self, x):
"""
Score a solution
"""
y = 0
for i in range(len(x) - 1):
y += self.distance_matrix[x[i]][x[i + 1]]
return y
让我们一步步分析代码中的关键部分：
-
第一步是初始化
__init__。我们定义问题，这里是一个 TSP 实例，并根据需要提供超参数。 -
ant_colony_optimization部分包含核心执行。在指定的迭代次数k_max内，算法致力于提升当前的最佳解。它涉及部署多个蚂蚁m，每只蚂蚁在图中移动。 -
ant_walk提供了单只蚂蚁旅程的模拟。在 while 循环中，蚂蚁的路径通过选择基于其吸引力A的下一条边来构建。边的吸引力使用edge_attractiveness方法计算，考虑了信息素矩阵tau、alpha、距离矩阵和 beta。每只蚂蚁走完后，信息素矩阵会更新。
运行算法以处理问题实例，只需执行以下操作：
problem = TSP()
aco = AntColonyOptimization(problem)
best_solution, best_score = aco.ant_colony_optimization()
现在，ACO 与最优解相比表现如何？你可以通过 GIF 可视化进展，GIF 中的每个图像显示当前的最佳路径，让你观察到随时间的改进。
路径的改进。GIF 由作者提供。
最终解决方案的得分是 4944，非常接近最优解（误差小于 1%）！观察解决过程也很有趣：
解决过程。每个点表示一次蚂蚁行走。图片由作者提供。
在这个图表中，x 轴表示蚂蚁的数量，而 y 轴表示蚂蚁旅程中覆盖的距离。水平红线表示最优解决方案的分数，红点象征蚂蚁发现了新的改进解决方案。值得注意的是，通常需要几只蚂蚁才能发现更好的解决方案。然而，最后一个红点非常接近最优解决方案。还有一些策略可以提高 ACO 的性能，我将在接下来的部分中详细说明。
分配问题
TSP 是一个路径规划问题，而 ACO 最初是为了应对路径规划问题而设计的。虽然这已经是很久以前的事了，但在此期间，人们发现了使用 ACO 解决不同类型挑战的方法。其中一个值得强调的例子是其在分配问题中的应用。
分配问题是将“某物”分配到“其他物”的问题。一个例子是二次分配问题（QAP）。假设你有一组位置，并且你想将一组设施分配到这些位置。目标是确定最佳分配，以最小化总成本。将设施f分配到位置l的成本由以下公式决定：
-
设施之间的流量：设施之间有一定的流量或互动。这表示在设施之间转移了多少“物品”或“活动”。
-
位置之间的距离：每对位置都有一个对应的距离，代表在这些位置之间运输或操作所需的成本或努力。
将特定设施分配到特定位置的成本由设施之间的流量和位置之间的距离决定。具体来说，计算一对设施的成本是将它们的流量与它们被分配到的地点之间的距离相乘。
要找到特定分配的总成本，你需要将所有可能的设施对的成本加总，这些成本基于它们的位置信息。
一个五个设施放置在五个位置上的最优解决方案的示例。蓝色表示流量。点击放大。图片由作者提供。
当你想将 ACO 应用于分配问题而不是路径规划问题时，问题的表述会发生怎样的变化？我会把编码留给你，但我将从蚂蚁的视角提供对路径规划问题和分配问题之间差异的直观理解。
TSP（旅行商问题）主要是找到访问各个地点的最佳顺序。另一方面，QAP（设施布局问题）则将焦点转向决定物品或设施的放置位置。在 TSP 的 ACO（蚁群优化）中，蚂蚁学会偏向于特定的访问序列。相比之下，在处理 QAP 时，蚂蚁则倾向于选择特定的设施放置于特定的位置。在这种情况下，信息素轨迹表示将设施分配到特定位置的期望程度。
你可以这样想象：每只蚂蚁都有一个可用位置（QAP 中的位置）列表，在这些位置上分配物品（QAP 中的设施）。蚂蚁不断重复这些步骤，直到确定所有物品的最佳安排。
简而言之，可以将其视为蚂蚁协作确定最有效的物品分配方式，依据他们所获得的关于物品到位置分配的最佳效果的知识。
提高解决方案质量
有几种策略可以在更短的时间内获得更好的解决方案。以下是三条宝贵的建议，它们可以显著影响你的结果：
超参数搜索
如果你处理的是多个相同类型的问题，强烈建议进行超参数搜索。诸如迭代次数、蚂蚁数量、alpha、beta、rho 和 Q 等参数可以对算法性能产生重大影响。
让我们看一个例子。在下面的图表中，我们测试了两个不同的 alpha 值，而其他参数保持不变。线形图显示了 100 次蚂蚁运行的移动平均值，而点表示某次运行找到的最佳解决方案。
使用橙色设置（alpha = 2）的算法不仅发现了一个优越的解决方案（由橙点表示），而且比使用蓝色设置（alpha = 1）的算法更快速地完成了这一任务。
为了进一步强调超参数调整的影响，考虑一个包含 100 个节点的 TSP 问题。如果我们在超参数上进行随机搜索（10 次迭代）并绘制移动平均图，我们会观察到性能的显著变化。
在顶部，橙色和黄色线条产生了次优结果，而其他线条则接近于最佳解决方案。作为背景，需要指出的是，OR 工具求解器花费了 5 分钟找到最优解，而 ACO 算法则在最多 5 秒内完成运行。
预热程序
在让蚂蚁开始之前，你可以选择对信息素矩阵进行预热。这个准备步骤通常减少了找到最佳解决方案所需的迭代次数。在这篇论文中，提出了一种预热程序，并比较了其效果。
探索与开发
类似于许多算法，如果它们的探索不够充分，可能会遇到停滞。为了解决这个问题，你可以使用最大-最小蚂蚁系统。MMAS 通过将高信息素值分配给未探索的路径，鼓励蚂蚁探索这些路径。当停滞发生时，这些路径的轨迹会被重置为这个较高的值。与原始算法相比，MMAS 的一个额外优势是只有全局最佳路径或迭代中的最佳路径才能增强其信息素轨迹。这些调整有助于在探索和利用之间实现更有效的平衡。
照片由Shardar Tarikul Islam提供，来自Unsplash
结论
蚁群优化算法是一个有趣的算法，它能够在问题规模增加的情况下仍然找到高质量的解决方案。这使它成为一个强大的算法。
为了进一步提升算法的性能并发现更优解，考虑融入一些关键策略，如超参数搜索、预热过程以及平衡探索和利用的技术。这些调整可以显著改善你的结果。
如果你对 ACO 产生了兴趣并想进一步阅读，我推荐这篇文章。其中一位作者是 1992 年提出 ACO 的人。它提供了对这一卓越优化技术的宝贵见解和全面理解。
相关内容
图论及其应用简介
解释、参数、优点、缺点及使用案例
基于蚂蚁行为的较少为人知的启发式算法简介
towardsdatascience.com
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
