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

原文：TowardsDataScience

协议：CC BY-NC-SA 4.0

我的 #30DayMapChallange 第三周

原文：towardsdatascience.com/my-third-week-of-the-30daymapchallange-b42efbe38b51

我对 #30DayMapChallange 第三周的个人看法，这是一项旨在每年 11 月每天设计主题地图的社交挑战。

Milan Janosov

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

–

自 2019 年以来，地理信息系统（GIS）和空间分析社区每年 11 月都非常繁忙——这要归功于一个有趣的挑战，称为 #30DayMapChallange。每年，这个挑战都有一个主题日程，提出一个主题作为当天地图可视化的主要方向。虽然预设的每日主题无疑对创造力构成了限制，但它们也帮助参与者找到共同兴趣，分享数据源，并在视觉和技术上表达个人风格。

在这里，我想简要回顾一下我在这个挑战的第三周，详细介绍并展示我创建的不同地图——通常使用 Python 和各种空间分析及地理空间数据工具。

在这篇文章中，所有图片均由作者创建。

第 15 天 — OSM

OpenStreetMap 是获取地图数据的首选来源之一。虽然有许多库在其基础上进行构建，但我去年最喜欢的发现是 PrettyMaps，我用它来可视化 PlanetWare 精心挑选的一些全球顶级旅游景点，例如巴黎的埃菲尔铁塔。其余的收藏可以在 这里 查看，其中包括罗马的斗兽场（Colosseum）、纽约的自由女神像（Statue of Liberty）、法国的凡尔赛宫（Château de Versailles）、纽约的中央公园（Central Park）、中国的紫禁城（Forbidden City）、捷克共和国的布拉格城堡（Prague Castle）以及悉尼歌剧院（Sydney Opera House）。

第 15 天 — OSM

第 16 天 — 大洋洲

为了可视化澳大利亚及大洋洲众多岛屿，我使用了来自Natural Earth的 10 米分辨率物理矢量数据。也就是说，我在这里想象了澳大利亚和大洋洲的海岸线。对于这个线图，我使用了我以前的一种技巧——创建这些看起来像‘光剑’、发光、受《星球大战》启发的线条——这些线条完全使用 Python 中的 Matplotlib 创建。此外，找到合适的边界框有点棘手，所以我只是参考了一些在线地图。

第 16 天 — 大洋洲

第 17 天 — 流动

当我试图搞清楚这张地图时，主题是“流动”，我有点卡住了——然后我想，好吧，来点卡路里，继续前进。

tl;dr

我决定可视化巧克力（在协调系统中其标准代码为 HS 1806）的出口流网络，使用来自Comtrade的国际贸易数据。这个数据集包含了哪些国家向哪些国家出口了什么数量、什么价值的信息。我利用这些信息创建了一个非地理但具有拓扑结构的地图，一个网络可视化。在这个网络中，每个国家是一个节点，而国家 A 与国家 B 相连，如果 A 向 B 出口巧克力产品，其中链接的大小与 2022 年交易商品的总值成比例。节点颜色对应网络社区——那些看起来比其他国家内部交易更多的国家簇，而节点大小则衡量了每个国家通过出口这些甜品所获得的总价值。未出口的国家标记为深灰色。有趣的是：我是在布达佩斯的#Flow Specialty Coffee Bar & Bistro 喝咖啡时创建了这张地图的。

第 17 天 — 流动

第 18 天 — 大气层

在这个可视化中，我创建了显示全球月度最高气温值的实际热图，使用了World Climate data：

“这是 WorldClim 2.1 版的 1970–2000 年气候数据。该版本于 2020 年 1 月发布。包括了最低、平均和最高气温、降水量、太阳辐射、风速、水汽压和总降水量的月度气候数据。此外，还有 19 个‘生物气候’变量。”

因此，在我的地图中，每一帧对应一个月，使用了 WorldClim 2.1 栅格数据集的最高气温。

第 18 天 — 大气层

第 19 天 — 5 分钟地图

没错——我们只有 5 分钟时间来制作一个地图可视化。时间不多，所以我决定使用我经常做的事情——将网络数据和地理空间数据结合起来，展示另一个道路网络，这次是关于曼哈顿及其美丽的、完全人工的方格状道路系统。与其重新运行以前的代码，我从头开始写了这个笔记本，仅用大约 5 分钟，几乎没有查找资料（除了找到正确的十六进制代码）。

第 19 天——5 分钟地图

第 20 天——户外

为了在户外冒险——同时创建地图——我决定在线查看美丽的地标，让我感觉就像身临其境——当我坐在 Python 终端前时。因此，在我的户外地图中，我结合了可视化气候变化和宏大的ESA 的 Sentinel 数据的新更新，制作了这段乌普萨拉冰川，阿根廷的动画。这段影像基于查询了 2016 年到今天每半年最少云层的真实色彩图像。由于文件大小问题，这里提供了一个静态快照，您可以在这里找到动画。

第 20 天——户外

第 21 天——栅格

在前一天的卫星影像之后，我回到了栅格遥感数据。具体来说，我重新发布了我以前的一个作品，其中使用了旧的 Sentinel API 来计算苏黎世市的 NDVI 指数，展示了卫星图像波段的一个基本但非常实用的转换。根据#Wikipedia上的理论：

“𝘕𝘰𝘳𝘮𝘢𝘭𝘪𝘻𝘦𝘥 𝘥𝘪𝘧𝘧𝘦𝘳𝘦𝘯𝘤𝘦 𝘷𝘦𝘨𝘦𝘵𝘢𝘵𝘪𝘰𝘯 𝘪𝘯𝘥𝘦𝘹：𝘛𝘩𝘦 𝘯𝘰𝘳𝘮𝘢𝘭𝘪𝘻𝘦𝘥 𝘥𝘪𝘧𝘧𝘦𝘳𝘦𝘯𝘤𝘦 𝘷𝘦𝘨𝘦𝘵𝘢𝘵𝘪𝘰𝘯 𝘪𝘯𝘥𝘦𝘹 𝘪𝘴 𝘢 𝘴𝘪𝘮𝘱𝘭𝘦 𝘨𝘳𝘢𝘱𝘩𝘪𝘤𝘢𝘭 𝘪𝘯𝘥𝘪𝘤𝘢𝘵𝘰𝘳 𝘵𝘩𝘢𝘵 𝘤𝘢𝘯 𝘣𝘦 𝘶𝘴𝘦𝘥 𝘵𝘰 𝘢𝘯𝘢𝘭𝘺𝘻𝘦 𝘳𝘦𝘮𝘰𝘵𝘦 𝘴𝘦𝘯𝘴𝘪𝘯𝘨 𝘮𝘦𝘢𝘴𝘶𝘳𝘦𝘮𝘦𝘯𝘵𝘴, 𝘰𝘧𝘵𝘦𝘯 𝘧𝘳𝘰𝘮 𝘢 𝘴𝘱𝘢𝘤𝘦 𝘱𝘭𝘢𝘵𝘧𝘰𝘳𝘮, 𝘢𝘴𝘴𝘦𝘴𝘴𝘪𝘯𝘨 𝘸𝘩𝘦𝘵𝘩𝘦𝘳 𝘰𝘳 𝘯𝘰𝘁 𝘵𝘩𝘦 𝘵𝘢𝘳𝘨𝘦𝘵 𝘣𝘦𝘪𝘯𝘨 𝘰𝘵𝘩𝘦𝘳 𝘵𝘩𝘢𝘯 𝘢𝘯𝘺 𝘵𝘢𝘳𝘨𝘦𝘵 𝘣𝘦𝘪𝘯𝘨 𝘰𝘣𝘴𝘦𝘳𝘷𝘦𝘥 𝘤𝘰𝘯𝘵𝘢𝘪𝘯𝘴 𝘭𝘪𝘷𝘦 𝘨𝘳𝘦𝘦𝘯 𝘷𝘦𝘨𝘦𝘵𝘢𝘵𝘪𝘰𝘯。”

实际上，NDVI 是一个+1.0 到-1.0 的浮动值变量，我们可以通过结合卫星图像的红色和近红外波段来计算。由于绿色表面的光吸收特性，这样的结果会使得绿色的像素值更高。

由于我使用了免费的#sentinel卫星图像，我的像素大小为 10x10 米。在可视化中，我使用了绿色-黄色-红色的颜色图来展示从+1.0 到 0 的 NDVI 值。负值通常表示水体，接近零的正值是建筑区域的良好参考，而值在 0.1 到 0.5 之间表示稀疏，以上则表示密集的绿色植被。

第 21 天——栅格

这是我第三周参与#30DayMapChallenge 的总结——还有 10 天就结束了，所以准备好迎接更多地图吧！

第二周的概况请点击这里查看！

第一周的概况请点击这里查看！

我的（非常）个人数据仓库

原文：towardsdatascience.com/my-very-personal-data-warehouse-fitbit-activity-analysis-with-duckdb-8d1193046133?source=collection_archive---------3-----------------------#2023-05-31

使用 DuckDB 分析 Fitbit 活动数据

Simon Aubury

·

关注 发表在 Towards Data Science · 13 分钟阅读 · 2023 年 5 月 31 日

–

图片由Jake Hills提供，发布在Unsplash

可穿戴健身追踪器已成为我们生活中不可或缺的一部分，收集并跟踪我们日常活动、睡眠模式、位置、心率等数据。我已经使用 Fitbit 设备 6 年来监测我的健康。然而，我一直觉得数据分析功能不足 — 尤其是在我想跟踪长期健身目标的进展时。我的个人健身活动数据存档中埋藏了哪些见解？要开始探索，我需要一种有效的方法来对成千上万的记录不完善的 JSON 和 CSV 文件进行数据分析……额外加分的是分析过程中不需要将我的数据从笔记本电脑上移走。

进入 DuckDB — 一个轻量级、免费但强大的分析数据库，旨在简化数据分析工作流 — 它在本地运行。在这篇博客中，我想使用 DuckDB 来探索我的 Fitbit 数据，并分享使用 Seaborn 数据可视化的各种数据格式分析方法以及绘制我的健康和健身目标的途径。

导出 Fitbit 数据存档

首先，我需要获取我所有的历史健身数据。通过遵循 导出您的账户存档 的说明，Fitbit 使得导出您账户生命周期中的 Fitbit 数据变得相当简单。

使用导出 Fitbit 数据存档的说明 — 作者截图。

您需要确认您的请求……并保持耐心。我的存档创建了超过三天 — 但我最终收到了含有下载 ZIP 文件说明的电子邮件。该文件应包含由我的 Fitbit 或相关服务记录的所有个人健身活动。解压存档后会显示出大量的文件 — 例如，我在解压 79MB 文件后总共有 7,921 个文件。

数以千计的嵌套文件中的一小部分 — 作者截图。

让我们开始查看存档中可用的数据种类。

为什么选择 DuckDB？

有许多优秀的博客 (1,2,3) 描述了 DuckDB — TL;DR 摘要是 DuckDB 是一个开源的内存 OLAP 数据库，专为分析查询而构建。它本地运行，支持广泛的 SQL，并能直接在 Pandas 数据、Parquet、JSON 数据上运行查询。额外加分的是它与 Python 和 R 的无缝集成。它的极速处理能力和大部分内存处理使其成为构建个人数据仓库的好选择。

Fitbit 活动数据

我查看的第一个文件集合是活动数据。物理活动和广泛的锻炼信息似乎存储在编号的文件中，例如 Physical Activity/exercise-1700.json。

我无法弄清楚文件编号实际意味着什么，我猜这些编号只是用于一组锻炼文件的递增整数。在我的数据导出中，最早的文件从 0 开始，经过 6 年的时间到达文件编号 1700。里面是一个记录数组，每个记录都有一个活动的描述。记录似乎会根据活动的不同而变化——这是一个“步行”的示例。

"activityName" : "Walk",
  "averageHeartRate" : 79,
  "calories" : 122,
  "duration" : 1280000,
  "steps" : 1548,
  "startTime" : "01/06/23 01:08:57",
  "elevationGain" : 67.056,
  "hasGps" : false,
  : : : :
  "activityLevel" : [
    { "minutes" : 1, "name" : "sedentary"},
    { "minutes" : 2, "name" : "lightly"},
    { "minutes" : 6, "name" : "fairly"},
    { "minutes" : 6, "name" : "very"
  }]

这些物理活动数据是我笔记本电脑上 7,921 个文件中的一个。幸运的是，DuckDB 可以使用 read_json 函数从 JSON 文件中读取（并自动检测模式），让我可以通过一个 SQL 语句将所有锻炼文件加载到 physical_activity 表中。值得注意的是，我需要指定日期格式掩码，因为 Fitbit 导出的日期格式非常 美国风格 😕。

CREATE OR REPLACE TABLE physical_activity
as
SELECT 
  startTime + INTERVAL 11 hours as activityTime
, activityName
, activityLevel
, averageHeartRate
, calories
, duration / 60000 as duration_minutes
, steps
, distance
, distanceUnit
, tcxLink
, source
FROM read_json('./Physical Activity/exercise-*.json'
, format='array'
, timestampformat='%m/%d/%y %H:%M:%S');

这个 SQL 命令从磁盘读取物理活动数据，将活动和持续时间转换并加载到内存中的 DuckDB 表中。

将物理活动数据加载到数据框中

我想了解我每个月的时间花费情况。由于活动数据存储在非常细粒度的级别上，我使用了 DuckDB SQL time_bucket 函数将 activityTime 时间戳截断为每月的桶。将分组的物理活动数据加载到数据框中可以通过这个汇总 SQL 完成，查询结果可以通过 << 运算符导入 Pandas 数据框中。

activity_df <<
  select time_bucket(interval '1 month', activityTime) as activity_day
  , activityName
  , sum(duration_minutes) as duration
  from physical_activity
  where activityTime between '2022-09-01' and '2023-05-01'
  group by 1, 2
  order by 1;

这个 SQL 查询将我的活动数据（骑车、步行、跑步等）分组到每月的桶中，让我可以真实地反映我投入到物理活动中的时间。

绘制每月活动分钟数

我现在想要通过视觉方式探索我的活动数据——所以我们来获取 Fitbit 数据并制作一些统计图形。我将使用 Python 的 Seaborn 数据可视化库直接从 activity_df 数据框中创建一个每月活动分钟数的条形图。

import matplotlib.pyplot as plt
import seaborn as sns
from matplotlib.dates import DateFormatter
plt.figure(figsize=(15, 6))
plt.xticks(rotation=45)

myplot =sns.barplot(data=activity_df, x="activity_day", y="duration", hue="activityName")
myplot.set(xlabel='Month of', ylabel='Duration (min)', title='Monthly Activity Minutes')
plt.legend(loc="upper right", title='Activity') 
plt.show()

执行此操作会生成这个条形图。

锻炼活动细分——作者截图。

看起来我的主要活动仍然是步行，而且我在 2023 年的新年决心是更频繁地跑步，但实际上并没有发生（还？）。

睡眠

关于 三分之一的成年人睡眠不足，所以我想探索我的长期睡眠模式。在我的 Fitbit 档案中，睡眠数据似乎被记录在以日期命名的文件中，例如Sleep/sleep-2022-12-28.json。每个文件包含一个月的数据，但混淆的是，文件的日期为事件发生前的月份。例如，文件sleep-2022-12-28.json似乎包含了 2023 年 1 月 2 日至 2023 年 1 月 27 日的数据。不管怎样 — 文件命名的奇怪之处暂且不提，我们可以探讨文件的内容。在记录中有一个扩展的“levels”块，详细描述了睡眠类型（清醒、浅睡、快速眼动、深睡）。

"logId" : 39958970367,
  "startTime" : "2023-01-26T22:47:30.000",
  "duration" : 26040000,
  :: :: ::
  "levels": 
    "summary" : {
      {
      "light": { "count": 30, "minutes": 275},
      "rem": { "count": 4, "minutes": 48 },
      "wake": { "count" : 29, "minutes" : 42 },
      "deep" : { "count" : 12, "minutes" : 75}
      }
    }

如果查看一些较旧的文件（可能是用我以前的 Fitbit Surge 设备创建的），会发现睡眠类型的分类有所不同（躁动、不清醒、睡眠）。

"logId" : 18841054316,
  "startTime" : "2018-07-12T22:42:00.000",
  "duration" : 25440000,
  :: :: ::
  "levels" : {
    "summary" : {
      "restless" : {"count" : 9, "minutes" : 20 },
      "awake" : { "count" : 2, "minutes" : 5 },
      "asleep" : { "count" : 0,   "minutes" : 399}
    }
  }

无论数据模式如何，我们都可以使用 DuckDB JSON 读取器将记录读入一个表格中。

CREATE OR REPLACE TABLE sleep_log
as
select dateOfSleep 
, levels
from read_json('./Sleep/sleep*.json'
, columns={dateOfSleep: 'DATE', levels: 'JSON'}
, format='array') ;

睡眠数据的模式变化

我想处理我所有的睡眠数据，并处理记录睡眠的模式变化（很可能是因为我更换了 Fitbit 设备的型号）。一些记录的时间标记在$.awake上，这与$.wake类似（但不完全相同）。

我使用了 SQL 中的 coalesce 函数 — 它返回第一个计算结果为非 NULL 值的表达式，以结合类似类型的睡眠阶段。

sleep_log_df <<
  select dateOfSleep
  , cast(coalesce(json_extract(levels, '$.summary.awake.minutes'), json_extract(levels, '$.summary.wake.minutes')) as int) as min_wake
  , cast(coalesce(json_extract(levels, '$.summary.deep.minutes'), json_extract(levels, '$.summary.asleep.minutes')) as int) as min_deep
  , cast(coalesce(json_extract(levels, '$.summary.light.minutes'), json_extract(levels, '$.summary.restless.minutes')) as int) as min_light
  , cast(coalesce(json_extract(levels, '$.summary.rem.minutes'), 0) as int) as min_rem
  from sleep_log
  where dateOfSleep between '2023-04-01' and '2023-04-30'
  order by 1;

使用 DuckDB，我可以通过 json_extract 提取嵌套 JSON 中的时长阶段，以生成一个 sleep_log_df 数据框，将所有历史睡眠阶段进行分组。

绘制睡眠活动图

我们可以将每日睡眠日志制作成堆叠条形图，显示每晚清醒、浅睡、深睡和 REM 睡眠的分类。

import matplotlib.pyplot as plt
import seaborn as sns
import matplotlib.dates as mdates

#create stacked bar chart
fig, axes = plt.subplots(figsize=(15,6))
myplot = sleep_log_df.set_index('dateOfSleep').plot(ax=axes, kind='bar', stacked=True, color=['chocolate', 'palegreen', 'green', 'darkblue'])
myplot.set(xlabel='Date', ylabel='Duration (min)', title='Sleep')
axes.xaxis.set_major_locator(mdates.DayLocator(interval=7))
plt.legend(loc="upper right", labels = ['Awake', 'Deep', 'Light', 'REM']) 
plt.xticks(rotation=45)
plt.show()

加载一个月的睡眠数据让我能够进行更广泛的睡眠时长分析。

每晚的睡眠周期时长 — 作者截图。

将多晚的睡眠数据绘制在一个图表上，使我能够开始理解星期几和周期性事件如何影响我的睡眠时长和质量。

心率

心率数据被非常频繁地捕捉（每10–15 秒一次），存储在名为Physical Activity/heart_rate-2023-01-26.json的每日文件中。这些文件非常大 — 每天约有 70,000 行 — 所有数据都包装在一个数组中。

[{{"dateTime": "01/25/25 13:00:07", "value": {"bpm": 54, "confidence": 2}},
  {"dateTime": "01/25/25 13:00:22", "value": {"bpm": 54, "confidence": 2}},
  {"dateTime": "01/25/25 13:00:37", "value": {"bpm": 55, "confidence": 2}},
  : : : : : :
  {"dateTime": "01/26/26 12:59:57", "value": {"bpm": 55, "confidence": 3}
}]

我的理论是文件名表示用户的时区。例如，在我的时区（GMT+11），命名为heart_rate-2023-01-26.json的数据覆盖了 26 日 00:00（AEST）至 23:59（AEST） - 如果文件中的日期为 GMT，则逻辑上是合理的。

转换 JSON 文件

到目前为止，我已经成功处理了包含 DuckDB 函数的 Fitbit 数据。然而，在处理这些巨大的心率文件时，我遇到了问题。当尝试处理 JSON 文件中的大数组记录时，DuckDB 给出了这个错误。

(duckdb.InvalidInputException) “INTERNAL Error: Unexpected yyjson tag in ValTypeToString”

我认为这个错误信息是一个突如其来的提醒，告诉我期望一个 JSON 数组有这么多元素是不合理的。解决办法是预处理文件，使其不是 JSON 记录数组，而是转换为换行符分隔的 JSON 或 ndjson。

{"dateTime": "01/25/23 13:00:07", "value": {"bpm": 54, "confidence": 2}
{"dateTime": "01/25/23 13:00:22", "value": {"bpm": 54, "confidence": 2}
{"dateTime": "01/25/23 13:00:37", "value": {"bpm": 55, "confidence": 2}
  : : : : : :
{"dateTime": "01/26/23 12:59:57", "value": {"bpm": 55, "confidence": 3}

为了将心率数组记录转换为换行符分隔的 JSON，我使用了一点小巧的 Python 代码来转换每个文件。

import glob
import json
import ndjson
import re

for json_src_file in sorted(glob.glob('./Physical Activity/steps-*.json')):
  json_dst_file = re.sub('\.[a-z]*$', '.ndjson', json_src_file)
  print(f'{json_src_file} -->  {json_dst_file}')
  with open(json_src_file) as f_json_src_file:
    json_dict =json.load(f_json_src_file) 
    with open(json_dst_file, 'w') as outfile:
      ndjson.dump(json_dict, outfile)

这将查找每个 .json 文件，读取内容并将其转换为换行符分隔的 JSON，并用 .ndjson 文件扩展名创建新文件。这将一个包含 70,000 条记录的数组转换为一个包含 70,000 行的文件——每条 JSON 记录现在存储在新的一行上。

将心率数据加载到表中

使用新转换的 ndjson 文件，我现在准备将心率数据加载到 DuckDB 表中。注意使用 timestampformat='%m/%d/%y %H:%M:%S'); 来描述日期中的前导月份（例如 “01/25/23 13:00:07”)

CREATE OR REPLACE TABLE heart_rate
as
SELECT dateTime + INTERVAL 11 hours as hr_date_time
, cast(value->'$.bpm' as integer) as bpm
FROM read_json('./Physical Activity/*.ndjson'
, columns={dateTime: 'TIMESTAMP', value: 'JSON'}
, format='newline_delimited'
, timestampformat='%m/%d/%y %H:%M:%S');

我们可以通过将格式设置为 ’newline_delimited’ 来加载所有 .ndjson 文件。注意我们可以通过 JSON 提取来提取 BPM（每分钟心跳次数）并将其转换为整数。

DuckDB 在处理 JSON 时非常快速 — 作者截图。

值得在这里强调 DuckDB 的速度是多么惊人——加载 1200 万条记录仅用了 2.8 秒！

将心率加载到数据框中

在加载了 1200 万条心率测量记录后，我们将 5 月 21 日的一天数据加载到数据框中。

hr_df << 
  SELECT time_bucket(interval '1 minutes', hr_date_time) as created_day
  ,  min(bpm) as bpm_min
  ,  avg(bpm) as bpm_avg
  ,  max(bpm) as bpm_max
  FROM heart_rate
  where hr_date_time between '2023-05-21 00:00' and '2023-05-21 23:59'
  group by 1;

这个 DuckDB 查询将心率的变异性聚合到 1 分钟的时间块中；在每个周期内进行最小值、平均值和最大值的分类。

绘制心率图

我可以使用这样的图来绘制心率（顺便炫耀一下，我确实在早上 6 点去跑步了）

import matplotlib.pyplot as plt
from matplotlib.dates import DateFormatter
plt.figure(figsize=(15, 6))
plt.xticks(rotation=45)

myplot = sns.lineplot(data=hr_df, x="created_day", y="bpm_min")
myplot = sns.lineplot(data=hr_df, x="created_day", y="bpm_avg")
myplot = sns.lineplot(data=hr_df, x="created_day", y="bpm_max")
myFmt = DateFormatter("%H:%M")
myplot.xaxis.set_major_formatter(myFmt)
myplot.set(xlabel='Time of day', ylabel='Heart BPM', title='Heart rate')
plt.show()

一天中的心率 — 作者截图。

以细粒度探索心率使我能够跟踪我的健身目标——特别是如果我坚持我的常规跑步计划。

步骤

步骤记录在名为 Physical Activity/steps-2023-02-26.json 的每日文件中。这似乎是对一天中每个时间段块（每 5 到 10 分钟）内的步骤的详细计数。

[{
  "dateTime" : "02/25/23 13:17:00",
  "value" : "0"
},{
  "dateTime" : "02/25/23 13:52:00",
  "value" : "5"
},{
  "dateTime" : "02/25/23 14:00:00",
  "value" : "0"
},{
:: :: ::
},{
  "dateTime" : "03/24/23 08:45:00",
  "value" : "15"
}]

为了将步骤聚合为每日统计，我需要将 GMT 转换为我的本地时区（GMT+11）。

steps_df <<
select cast(time_bucket(interval '1 day', dateTime + INTERVAL 11 hours	) as DATE) as activity_day
, sum(value) as steps
from read_json('./Physical Activity/steps-2023-04-27.ndjson'
, auto_detect=True
, format='newline_delimited'
, timestampformat='%m/%d/%y %H:%M:%S') 
group by 1;

将每日步骤聚合到 steps_df 数据框中，使我能够探索长期的活动趋势，因为我努力超越 10,000 步以实现健康益处的提升。

绘制每日步骤

现在我们可以将数据框绘制为每日步数

import matplotlib.pyplot as plt
from matplotlib.dates import DateFormatter
plt.figure(figsize=(15, 6))
plt.xticks(rotation=45)

myplot = sns.barplot(data=steps_df, x="activity_day", y="steps")
myplot.set(xlabel='Day', ylabel='Steps', title='Daily steps')
plt.show()

每日步数统计——作者截图。

这显示我还需要努力达到我的每日步数目标——这对我的新年健身决心来说又是一次打击。

GPS 映射

Fitbit 将 GPS 记录的活动存储为TCX (训练中心 XML) 文件。这些 XML 文件在下载的 ZIP 文件中没有，但我们在身体活动文件中有其位置的参考，我可以像这样查询。

select tcxLink 
from physical_activity
where tcxLink is not null;

tcxLink 字段是对身体活动文件中位置的 URL 参考。

每个 TCX 文件的 URL——作者截图。

我们可以直接在浏览器中使用这个 URL（登录 Fitbit 网站后）来下载 GPS XML 文件。查看 TCX 文件内部，我们会发现每隔几秒钟就有低级别的 GPS 位置数据。

TCX GPS XML 文件样本内容——作者截图。

好消息是数据中有一些明显的字段，如纬度、经度和时间。不太好的消息是数据是 XML 格式的，因此我们需要在加载到 DuckDB 之前预处理这些文件，因为目前 XML 格式不被文件读取器支持。我们可以通过另一段 Python 代码将 XML 文件转换为 JSON 文件，循环遍历每个 .tcx 文件。

这里有一些复杂的 XML 嵌套，位置数据位于 TrainingCenterDatabase/Activities/Activity/Lap 下。

import glob
import json
import ndjson
import xmltodict
import re

for xml_src_file in sorted(glob.glob('MyFitbitData/tcx/*.tcx')):
    json_dst_file = re.sub('\.[a-z]*$', '.ndjson', xml_src_file)
    print(f'{xml_src_file} -->  {json_dst_file}')
    with open(xml_src_file) as f_xml_src_file:
        # erase file if it exists
        open(json_dst_file, 'w') 
        data_dict = xmltodict.parse(f_xml_src_file.read())
        # Loop over the "laps" in the file; roughly every 1km
        for lap in data_dict['TrainingCenterDatabase']['Activities']['Activity']['Lap']:
            data_dict_inner = lap['Track']['Trackpoint']
            # append file
            with open(json_dst_file, 'a') as outfile:
                ndjson.dump(data_dict_inner, outfile)
                outfile.write('\n')

加载 GPS 地理空间数据

我们可以这样加载地理空间数据

route_df <<
SELECT time
, position
, cast(json_extract_string(position, '$.LatitudeDegrees') as float) as latitude
, cast(json_extract_string(position, '$.LongitudeDegrees') as float) as longitude
FROM read_json('MyFitbitData/tcx/54939192717.ndjson'
, columns={Time: 'TIMESTAMP', Position: 'JSON', AltitudeMeters: 'FLOAT', DistanceMeters: 'FLOAT', HeartRateBpm: 'JSON'}
, format='newline_delimited'
, timestampformat='%Y-%m-%dT%H:%M:%S.%f%z');

这个 DuckDB 查询扁平化了 JSON，将纬度、经度和时间转换为正确的数据类型，并加载到 route_df 数据框中。

使用 Folium 可视化 GPS 路线

拥有一个位置信息的表格并不够直观，因此我想开始在交互式地图上绘制我的跑步路线。我使用了这篇博客来帮助使用 Folium 可视化路线。修改代码帮助我绘制了自己的跑步路线，例如这是我在堪培拉度假时的一次跑步路线图。

import folium

route_map = folium.Map(
    location=[-35.275, 149.129],
    zoom_start=13,
    tiles='openstreetmap',
    width=1024,
    height=600
)

coordinates = [tuple(x) for x in route_df[['latitude', 'longitude']].to_numpy()]
folium.PolyLine(coordinates, weight=8, color='red').add_to(route_map)
display(route_map)

跑步的 Folium 地图图示——作者截图。

这生成了一个使用开放街图 瓦片的跑步图，给我提供了一个很好的交互式详细跑步地图。

数据目标和健身目标总结

我是否更接近分析我的 Fitbit 设备数据的目标——绝对是的！DuckDB 证明是一个理想的灵活轻量级分析工具，能够快速处理我的大量混乱的 Fitbit 数据。通过广泛的 SQL 支持和灵活的文件解析选项，DuckDB 能够在几秒钟内处理数百万条记录，将数据本地导入数据框，这使得 DuckDB 成为构建个人数据仓库的理想选择。

至于我的健身目标——我还有一些工作要做。我觉得我现在应该离开这个博客，因为我今天的步数目标还差一点。

Code

🛠️用于 Fitbit 活动分析的代码——github.com/saubury/duckdb-fitbit

朴素贝叶斯分类

原文：towardsdatascience.com/naive-bayes-classification-41d1fe802a1e

朴素贝叶斯分类器家族的深入解释，包括一个 Python 中的文本分类示例

Dr. Roi Yehoshua

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

–

图片由 Mediocre Studio 提供，来源于 Unsplash

朴素贝叶斯分类器是一类基于应用贝叶斯定理并假设特征之间独立性的概率分类器。

这些分类器在训练和预测中都非常快速，并且具有很高的可扩展性和可解释性。尽管它们的假设过于简化，但它们在复杂的实际问题上通常表现良好，尤其是在如垃圾邮件过滤和情感分析等文本分类任务中，其朴素假设通常成立。

朴素贝叶斯也是最早的生成模型之一（早于 ChatGPT……），它学习每个类别中的输入分布。这些模型不仅可用于预测，还可以用于生成新样本（有关生成模型与判别模型的深入讨论，请参见 这篇文章）。

在这篇文章中，我们将深入探讨朴素贝叶斯模型及其变体，然后展示如何使用 Scikit-Learn 中的实现来解决文档分类任务。

背景：贝叶斯定理

贝叶斯定理（或贝叶斯规则）是概率中的一个重要定理，它允许我们基于与事件相关的先验知识计算事件的条件概率。

从数学上讲，定理指出，对于任何事件 A 和 B：

贝叶斯规则

• P(A|B)是后验概率，即给定B的情况下事件A发生的概率。

• P(B|A)是似然性，即给定A的情况下事件B发生的概率。

• P(A)是先验概率，即A的概率，没有任何先验条件。

• P(B)是边际概率，即B的概率，没有任何先验条件。

贝叶斯定理特别适用于从结果推断原因，因为通常更容易辨别在原因存在或不存在的情况下结果的概率，而不是相反。例如，估计一个患有脑膜炎的患者会出现头痛的概率比估计一个头痛的患者是否患有脑膜炎要容易得多（因为许多其他疾病也可能导致头痛）。在这种情况下，我们可以应用贝叶斯规则如下：

示例

已知大约 25%的肺癌患者会出现胸痛。假设肺癌的发生率为每 10 万人中 50 人，全球胸痛的发生率为每 10 万人中 1,500 人。那么，一个有胸痛的患者患肺癌的概率是多少？

让我们将给定的输入以概率的形式表示。设L为患肺癌的事件，C为胸痛的事件。根据我们拥有的数据，我们知道：

• P(C|L) = 0.25

• P(L) = 50 / 100,000 = 0.0005

• P(C) = 1,500 / 100,000 = 0.015

使用贝叶斯规则，给定胸痛的情况下患肺癌的后验概率是：

即，患者患肺癌的概率只有 0.833%。

朴素贝叶斯模型

朴素贝叶斯模型是概率分类器，即它们不仅为给定样本分配一个类别标签，还提供该样本属于该类别的概率估计。例如，朴素贝叶斯模型可以预测某个电子邮件有 80%的概率是垃圾邮件，20%的概率是正常邮件。

回顾一下在监督学习问题中，我们会得到一个包含n个标记样本的训练集：D = {(x₁, y₁), (x₂, y₂), … , (xₙ, yₙ)}，其中xᵢ是一个m维向量，包含样本i的特征，而yᵢ表示该样本的标签。在分类问题中，标签可以取任何一个K类别，即，y ∈ {1, …, K}。

我们区分两种类型的分类器：

1. 确定性分类器为每个样本输出一个硬标签，而不提供类别的概率估计。这类分类器的例子包括 K-最近邻、决策树 和 SVM。

2. 概率分类器输出对 k 类的概率估计，然后根据这些概率给样本分配标签（通常是概率最高的类别的标签）。这种分类器的例子包括朴素贝叶斯分类器，逻辑回归和使用逻辑/softmax 输出层的神经网络。

给定一个样本 (x, y)，朴素贝叶斯分类器使用贝叶斯规则计算它属于类别 k 的概率（即 y = k）：

方程右侧的概率是从训练集估计得到的。

首先，类别先验概率 P(y = k) 可以通过类别 k 在训练样本中的相对频率来估计：

其中 nₖ 是属于类别 k 的样本数量，n 是训练集中的样本总数。

其次，边际概率 P(x) 可以通过对贝叶斯规则中所有类别的分子部分进行求和来计算：

由于边际概率不依赖于类别，如果我们只对分配硬标签给样本感兴趣（而不提供概率估计），则不需要计算它。

最后，我们需要估计在给定类别下特征的可能性，即 P(x|y = k)。估计这些概率的主要问题是数量太多，可能在训练集中没有足够的数据来估计所有的概率。

例如，假设 x 由 m 个二进制特征组成，例如，每个特征表示某个词是否出现在文本中。在这种情况下，为了建模 P(x|y)，我们需要从训练集中为每个类别估计 2ᵐ 个条件概率（每个 x₁, …, xₘ 的所有可能组合），因此总共有 2ᵐK 个概率。在大多数情况下，我们在训练集中没有足够的样本来估计所有这些概率，即使有，也需要耗费指数时间。

朴素贝叶斯假设

为了减少需要估计的参数数量，朴素贝叶斯模型做出了以下假设：特征在给定类别变量的情况下是相互独立的。

这个假设允许我们将概率 P(x|y = k) 写作每个个体特征在给定类别下的条件概率的乘积：

朴素贝叶斯假设

例如，在垃圾邮件过滤任务中，朴素贝叶斯假设意味着诸如“rich”和“prince”这样的词独立地对预测邮件是否是垃圾邮件作出贡献，无论这些词之间是否存在任何可能的关联。

朴素假设在应用领域如文本分类和推荐系统中大致成立，在这些领域中，特征通常彼此独立。

朴素贝叶斯假设显著减少了需要从数据中估计的参数数量。例如，在输入x包含m个二元特征的情况下，它将模型的参数数量从 2ᵐ减少到每个类别m。

MAP（最大后验）

基于朴素贝叶斯假设，我们现在可以将类别后验概率写作如下：

如果我们只关心将类别标签分配给给定样本（而不关心概率），我们可以忽略分母 P(x)，并使用以下分类规则：

这称为**MAP（最大后验）**决策规则，因为它选择最大化后验概率的假设。

只要正确的类别被预测为比其他类别更可能，朴素贝叶斯将做出正确的 MAP 决策，即使概率估计不准确。这为模型提供了一些对基础朴素独立假设缺陷的鲁棒性。

请注意，如果我们假设所有先验 P(y) 是等可能的（例如，当我们没有关于哪个假设更可能的先验信息时），那么 MAP 决策规则等同于 MLE（最大似然估计）决策规则，它选择最大化给定模型 P(x|y) 的数据似然性的模型。（你可以在这篇文章中了解更多关于最大似然的信息。）

参数估计

我们现在剩下的任务是估计每个特征 j 和每个类别 k 的条件概率 P(xⱼ|y = k)。这个估计依赖于特征的类型（例如，离散或连续）以及我们假设其具有的概率分布。

对特征分布的假设称为事件模型。每个事件模型导致不同类型的朴素贝叶斯分类器。在接下来的部分，我们将讨论不同的事件模型以及如何在每种模型中估计模型参数。

伯努利朴素贝叶斯

在伯努利事件模型中，特征xⱼ被建模为具有伯努利分布的独立二元变量，即每个特征xⱼ在给定样本中出现的概率为pⱼ，不出现的概率为 1 − pⱼ。

例如，在文本分类任务中，每个特征xⱼ可能代表文本中词汇表中第j个单词的出现或缺失。

在伯努利事件模型中，概率P(xⱼ|y = k)是通过特征j在类别k的样本中出现的频率来估计的：

其中nⱼₖ是类别k中包含特征xⱼ的样本数量，nₖ是类别k中样本的总数。

分类朴素贝叶斯

分类事件模型是伯努利事件模型对V个类别（而不是仅两个类别）的扩展。在这个模型中，我们假设每个特征是一个类别（离散）变量，可以取V个可能的类别中的一个，其概率为pᵢ，其中所有概率的总和为 1。

在这个事件模型中，我们需要估计每个特征xⱼ和每个类别v的概率P(xⱼ = v|y = k)。与之前的模型类似，我们通过特征j在类别k的样本中取值v的频率来估计这个概率：

其中nⱼᵥₖ是类别k中特征xⱼ取值v的样本数量，nₖ是类别k中样本的总数。

示例：顾客购买预测

假设我们有一个包含商店顾客过去购买数据的表格：

训练集

表中的每一行包含顾客的年龄、是否是学生、收入水平、信用评级以及是否购买了产品。

一个具有以下属性的新顾客到达商店：

<Age = Young, Student = Yes, Income = Low, Credit = Excellent>

你需要预测这个顾客是否会购买产品。

我们首先通过计算 Buys = Yes（10 行中的 6 行）和 Buys = No（10 行中的 4 行）的行数来计算类别先验概率。

然后，我们计算每个类别中特征的可能性：

因此，类别后验概率为：

α是归一化因子（α = 1 / P(x））。

由于P(Buys = Yes|x) > P(Buys = No|x)，我们的预测是顾客会购买产品。

如果我们想得到顾客购买产品的实际概率，我们可以首先使用两个后验概率之和为 1 的事实来找到归一化因子：

然后，我们可以将其代入“Buy = Yes”的后验概率中：

客户购买该产品的概率是 54.21%。

多项式朴素贝叶斯

在多项式事件模型中，我们假设数据集只有一个类别特征x，它可以取m个类别中的一个，每个特征向量（x₁, …, xₘ）是一个直方图，其中xⱼ计算了x在特定实例中取值j的次数。

这个事件模型在处理文本文件时特别有用，其中m是词汇表中的词数，每个特征xⱼ表示词汇表中第j个词在文档中出现的次数。这种表示方法称为词袋模型：

词袋模型

在这个事件模型中，我们估计概率P(x = v|y = k)为特征x在类别k的样本中取得值v的频率：

其中nᵥₖ是类别k中x取值v的样本数量，nₖ是类别k中样本的总数量。

示例：垃圾邮件过滤器

我们希望基于一个包含 100 封邮件的训练集构建一个垃圾邮件过滤器：其中 80 封为正常邮件，20 封为垃圾邮件。每种类型邮件中的词语计数在以下表格中给出：

正常邮件和垃圾邮件中的词语计数

一封新的邮件到达，内容为“rich friend need money”。这是正常邮件还是垃圾邮件？

首先计算类别的先验概率：

接下来，我们估计每种类型邮件中词语的可能性。正常邮件中的词语总数为：80 + 5 + 5 + 20 + 50 + 40 = 200，而垃圾邮件中的词语总数为：20 + 20 + 30 + 25 + 25 = 120。因此，词语的可能性为：

因此，类别的后验概率是：

因此，我们的预测是该邮件是垃圾邮件。

拉普拉斯平滑

如果某个特征和某个类别在训练集中从未一起出现，那么它的可能性估计将为零。由于特征的可能性是相乘的，这将抹去我们从其他特征中获得的所有信息。

例如，由于垃圾邮件中缺少“movie”这个词，如果我们收到一封包含“movie”的邮件，它将被自动认为是正常邮件，即使邮件中的其他词语非常“垃圾”。

作为一个具体的例子，考虑一封内容为“movie rich rich prince”的邮件。这封邮件将被分类为正常邮件，尽管“rich”和“prince”这两个词与垃圾邮件高度相关，因为垃圾邮件的后验概率为零：

为了处理这个问题，我们在所有概率估计中加入一个小的样本修正（称为伪计数），以确保概率不会被设定为零。

在多项式朴素贝叶斯中，修正如下应用：

其中α是平滑参数，n是训练集中的样本总数。设置α = 1 称为拉普拉斯平滑（最常见的），而α < 1 称为利德斯通平滑。

同样，在分类朴素贝叶斯中，修正如下应用：

其中nⱼ是特征j的可能类别数。

重新审视我们的垃圾邮件过滤器示例，让我们通过为所有词汇添加伪计数 1 来应用拉普拉斯平滑：

对词频应用拉普拉斯平滑

这次，包含文本“movie rich rich prince”的电子邮件将被分类为垃圾邮件，因为：

高斯朴素贝叶斯

直到现在我们假设所有特征都是离散的。朴素贝叶斯模型如何处理连续特征？

处理连续特征的主要方法有三种：

1. 将特征值离散化，并获得伯努利或分类分布特征。

2. 假设特征按照某个已知概率分布（通常是正态分布）分布，并从训练集中估计该分布的参数（例如，正态分布的均值和方差）。

3. 使用KDE（核密度估计）来估计特征的概率密度函数，利用给定的样本点作为核。

在高斯朴素贝叶斯中，我们采用第二种方法，假设特征的似然性是高斯分布的：

其中μⱼₖ是类别k中所有样本的xⱼ值的均值，σⱼₖ是这些值的标准差（这些是分布真实参数的最大似然估计）。

上述事件模型也可以组合使用，以处理异质数据集，即包含不同类型特征的数据集（例如，既有分类特征也有连续特征）。

Scikit-Learn 中的朴素贝叶斯分类器

模块 sklearn.naive_bayes 提供了上述提到的四种朴素贝叶斯分类器的实现：

1. BernoulliNB实现了伯努利朴素贝叶斯模型。

2. CategoricalNB实现了分类朴素贝叶斯模型。

3. MultinomialNB 实现了多项式朴素贝叶斯模型。

4. GaussianNB 实现了高斯朴素贝叶斯模型。

前三种类别接受一个名为alpha的参数，该参数定义了平滑参数（默认为 1.0）。

文档分类示例

在接下来的演示中，我们将使用MultinomialNB来解决文档分类任务。我们将使用的数据集是20 newsgroups dataset，该数据集包含 18,846 篇新闻组帖子，几乎均匀地划分为 20 个不同主题。这个数据集在机器学习中的文本应用研究中被广泛使用，包括文档分类和聚类。

加载数据集

你可以使用函数fetch_20newsgroups()在 Scikit-Learn 中下载带标签的文本文档。你可以选择将所有文档作为一个组下载，或者分别下载训练集和测试集（使用subset参数）。训练集和测试集的划分基于在特定日期之前或之后发布的消息。

默认情况下，文本文档包含一些元数据，如标题（例如，帖子的日期）、页脚（签名）和其他帖子的引用。由于这些特征与文本分类任务无关，我们将通过使用remove参数将它们剥离掉：

from sklearn.datasets import fetch_20newsgroups

train_set = fetch_20newsgroups(subset='train', remove=('headers', 'footers', 'quotes'))
test_set = fetch_20newsgroups(subset='test', remove=('headers', 'footers', 'quotes'))

请注意，第一次调用此函数时，可能需要几分钟时间下载所有文档，之后它们将被缓存到本地文件夹 ~/scikit_learn_data* 中。

函数的输出是一个包含以下属性的字典：

• data — 文档集合

• target — 目标标签

• target_names — 文档类别名称

让我们将文档及其标签存储到适当的变量中：

X_train, y_train = train_set.data, train_set.target
X_test, y_test = test_set.data, test_set.target

数据探索

让我们对数据进行一些基本的探索。我们在训练集和测试集中拥有的文档数量是：

print('Documents in training set:', len(X_train))
print('Documents in test set:', len(X_test))

Documents in training set: 11314
Documents in test set: 7532

简单计算显示，60%的文档属于训练集，40%属于测试集。

让我们打印类别列表：

categories = train_set.target_names
categories

['alt.atheism',
 'comp.graphics',
 'comp.os.ms-windows.misc',
 'comp.sys.ibm.pc.hardware',
 'comp.sys.mac.hardware',
 'comp.windows.x',
 'misc.forsale',
 'rec.autos',
 'rec.motorcycles',
 'rec.sport.baseball',
 'rec.sport.hockey',
 'sci.crypt',
 'sci.electronics',
 'sci.med',
 'sci.space',
 'soc.religion.christian',
 'talk.politics.guns',
 'talk.politics.mideast',
 'talk.politics.misc',
 'talk.religion.misc']

很明显，某些类别彼此紧密相关（例如，comp.sys.mac.hardware 和 comp.sys.ibm.pc.hardware），而其他类别则高度不相关（例如，sci.electronics 和 soc.religion.christian）。

最后，让我们查看训练集中其中一份文档（例如，第一个）：

print(X_train[0])

I was wondering if anyone out there could enlighten me on this car I saw
the other day. It was a 2-door sports car, looked to be from the late 60s/
early 70s. It was called a Bricklin. The doors were really small. In addition,
the front bumper was separate from the rest of the body. This is 
all I know. If anyone can tellme a model name, engine specs, years
of production, where this car is made, history, or whatever info you
have on this funky looking car, please e-mail.

不出所料，该文档的标签是：

categories[y_train[0]]

'rec.autos'

将文本转换为向量

为了将文本文档输入到机器学习模型中，我们首先需要将它们转换为数值向量（即向量化文本）。这个过程通常涉及文本的预处理和清理，然后选择合适的数值表示来表示文本中的词。

文本预处理包括多个步骤，其中最常见的包括：

1. 清理和标准化文本。这包括去除标点符号和特殊字符，并将文本转换为小写。

2. 文本分词，即将文本拆分成单个词或术语。

3. 停用词的去除。停用词是特定语言中常用的词。例如，英语中的停用词包括“the”、“a”、“is”、“and”。这些词通常被过滤掉，因为它们不携带有用的信息。

4. 词干提取或词形还原。词干提取通过去除或替换词缀将词还原为其词汇根，而词形还原将词还原为其规范形式（词元），并考虑词的上下文（词性）。例如，词computers的词元是computer，但其词汇根是comput。

以下示例演示了这些步骤在给定句子上的应用：

文本预处理示例

清理文本后，我们需要选择如何将其向量化为数值向量。最常见的方法有：

1. 词袋（BOW）模型。在这个模型中，每个文档通过词频向量表示（类似于我们在垃圾邮件过滤器示例中使用的）。

2. TF-IDF（词频与逆文档频率的乘积）通过乘以两个指标来衡量一个词对文档的相关性：

   (a) TF（词频）——单词在文档中出现的次数。

   (b) IDF（逆文档频率）——词在整个语料库中出现的频率的倒数。

   其思想是减少在语料库中频繁出现的词的权重，同时增加稀有词的权重（从而更能指示文档的类别）。

3. 词嵌入。在这种方法中，词被映射到实值向量中，以便具有相似意义的词在向量空间中具有接近的表示。这个模型通常用于深度学习，将在未来的帖子中讨论。

Scikit-Learn 提供了以下两种变换器，支持文本预处理和向量化：

1. CountVectorizer使用词袋模型。

2. TfIdfVectorizer使用 TF-IDF 表示。

这些变换器的重要超参数包括：

• lowercase — 是否在标记化之前将所有字符转换为小写（默认为 True）。

• token_pattern — 用于定义什么是令牌的正则表达式（默认正则表达式选择两个或更多字母数字字符的令牌）。

• stop_words — 如果为‘english’，则使用内置的英语停用词列表。如果为 None（默认值），则不会使用停用词。你也可以提供自己定制的停用词列表。

• max_features — 如果不为 None，则构建一个仅包含在训练语料库中具有最高术语频率的前max_features个词汇的词汇表。否则，将使用所有特征（这是默认值）。

请注意，这些变换器不提供诸如词干提取或词形还原之类的高级预处理技术。要应用这些技术，你将需要使用其他库，如NLTK（自然语言工具包）或spaCy。

由于朴素贝叶斯模型在 TF-IDF 表示上表现更好，我们将使用 TfidfVectorizer 将训练集中的文档转换为 TF-IDF 向量：

from sklearn.feature_extraction.text import TfidfVectorizer

vectorizer = TfidfVectorizer(stop_words='english')
X_train_vec = vectorizer.fit_transform(X_train)

提取的 TF-IDF 向量的形状是：

print(X_train_vec.shape)

(11314, 101322)

也就是说，语料库的词汇表中有 101,322 个独特的令牌。我们可以通过调用向量化器的get_feature_names_out()方法来检查这些令牌：

vocab = vectorizer.get_feature_names_out()
print(vocab[50000:50010]) # pick a subset of the tokens

['innacurate' 'innappropriate' 'innards' 'innate' 'innately' 'inneficient'
 'inner' 'innermost' 'innertubes' 'innervation']

显然，在 90 年代没有自动拼写检查器 😃

TF-IDF 向量非常稀疏，平均有 67 个非零组件，超过 100,000 个：

print(X_train_vec.nnz / X_train_vec.shape[0])

66.802987449178

让我们也将测试集中的文档向量化（请注意，在测试集上我们调用transform方法而不是fit_transform）：

X_test_vec = vectorizer.transform(X_test)

构建模型

现在，让我们构建一个多项式朴素贝叶斯分类器，并将其拟合到训练集上：

from sklearn.naive_bayes import MultinomialNB

clf = MultinomialNB(alpha=0.01)
clf.fit(X_train_vec, y_train)

请注意，我们需要将平滑参数α设置为一个非常小的数字，因为 TF-IDF 值被缩放到 0 和 1 之间，因此默认的α = 1 会导致值的剧烈变化。

评估模型

接下来，让我们在训练集和测试集上评估模型。

模型在训练集上的准确率和 F1 得分为：

from sklearn.metrics import f1_score

accuracy_train = clf.score(X_train_vec, y_train)
y_train_pred = clf.predict(X_train_vec)
f1_train = f1_score(y_train, y_train_pred, average='macro')

print(f'Accuracy (train): {accuracy_train:.4f}')
print(f'F1 score (train): {f1_train:.4f}')

Accuracy (train): 0.9595
F1 score (train): 0.9622

测试集上的准确率和 F1 得分为：

accuracy_test = clf.score(X_test_vec, y_test)
y_test_pred = clf.predict(X_test_vec)
f1_test = f1_score(y_test, y_test_pred, average='macro')

print(f'Accuracy (test): {accuracy_test:.4f}')
print(f'F1 score (test): {f1_test:.4f}')

Accuracy (test): 0.7010
F1 score (test): 0.6844

与训练集相比，测试集上的得分相对较低。为了调查错误来源，让我们绘制测试文档的混淆矩阵：

from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay

cm = confusion_matrix(y_test, y_test_pred)
disp = ConfusionMatrixDisplay(confusion_matrix=cm, display_labels=clf.classes_)
fig, ax = plt.subplots(figsize=(10, 8))
disp.plot(ax=ax, cmap='Blues')

测试集上的混淆矩阵

如我们所见，大多数混淆发生在高度相关的主题之间，例如：

• 主题 0（alt.atheism）和主题 15（soc.religion.christian）之间有 74 个混淆

• 主题 18（talk.politics.misc）和主题 16（talk.politics.guns）之间有 92 个混淆

• 主题 19（talk.religion.misc）和主题 15（soc.religion.christian）之间有 89 个混淆

根据这些发现，朴素贝叶斯分类器表现得相当不错。让我们看看它与其他标准分类算法的比较。

基准测试

我们将把朴素贝叶斯模型与另外四个分类器进行基准测试：逻辑回归，KNN，随机森林和AdaBoost。

首先编写一个函数，该函数获取一组分类器，并在给定的数据集上评估它们，同时测量它们的训练时间：

import time

def benchmark(classifiers, names, X_train, y_train, X_test, y_test, verbose=True):
    evaluations = []

    for clf, name in zip(classifiers, names):
        evaluation = {}
        evaluation['classifier'] = name

        start_time = time.time()
        clf.fit(X_train, y_train)
        evaluation['training_time'] = time.time() - start_time

        evaluation['accuracy'] = clf.score(X_test, y_test)
        y_test_pred = clf.predict(X_test)
        evaluation['f1_score'] = f1_score(y_test, y_test_pred, average='macro')

        if verbose:
            print(evaluation)
        evaluations.append(evaluation)
    return evaluations

我们现在将使用这五个分类器调用这个函数：

from sklearn.linear_model import LogisticRegression
from sklearn.neighbors import KNeighborsClassifier
from sklearn.ensemble import RandomForestClassifier, AdaBoostClassifier

classifiers = [clf, LogisticRegression(), KNeighborsClassifier(), RandomForestClassifier(), AdaBoostClassifier()]
names = ['Multinomial NB', 'Logistic Regression', 'KNN', 'Random Forest', 'AdaBoost']

evaluations = benchmark(classifiers, names, X_train_vec, y_train, X_test_vec, y_test)

我们得到的输出是：

{'classifier': 'Multinomial NB', 'training_time': 0.06482672691345215, 'accuracy': 0.7010090281465746, 'f1_score': 0.6844389919212164}
{'classifier': 'Logistic Regression', 'training_time': 39.38498568534851, 'accuracy': 0.6909187466808284, 'f1_score': 0.6778246092753284}
{'classifier': 'KNN', 'training_time': 0.003989696502685547, 'accuracy': 0.08218268720127456, 'f1_score': 0.07567337211476842}
{'classifier': 'Random Forest', 'training_time': 43.847145318984985, 'accuracy': 0.6233404142326076, 'f1_score': 0.6062667217793061}
{'classifier': 'AdaBoost', 'training_time': 6.09197473526001, 'accuracy': 0.36563993627190655, 'f1_score': 0.40123307742451064}

让我们绘制分类器的准确率和 F1 分数：

df = pd.DataFrame(evaluations).set_index('classifier')

df['accuracy'].plot.barh()
plt.xlabel('Accuracy (test)')
plt.ylabel('Classifier')

测试集上的准确率

df['f1_score'].plot.barh(color='purple')
plt.xlabel('F1 score (test)')

测试集上的 F1 分数

多项式 NB 在准确率和 F1 分数上都表现最佳。注意，分类器使用的是默认参数，未经过任何调整。为了更公平的比较，应在微调其超参数后比较算法。此外，一些算法如 KNN 会遭遇维度诅咒，需要进行降维才能使其有效工作。

让我们还绘制分类器的训练时间：

df['training_time'].plot.barh(color='green')
plt.xlabel('Training time (sec)')
plt.ylabel('Classifier')

不同分类器的训练时间

多项式 NB 的训练速度如此之快，以至于我们在图中看不到它的时间！通过查看上面的函数输出，我们可以看到其训练时间仅为 0.064 秒。请注意，KNN 的训练速度也非常快（因为实际上没有构建模型），但其预测时间（未显示）非常慢。

总之，多项式 NB 在所有检查的标准中显示出优越性。

查找最具信息量的特征

朴素贝叶斯模型还允许我们获取每个类别的最具信息量的特征，即具有最高可能性的特征 P(xⱼ|y)。

MultinomialNB 类有一个名为 feature_log_prob_ 的属性，它提供了每个类别的特征的对数概率，矩阵形状为 (n_classes, n_features)。

使用这个属性，编写一个函数来找到每个类别中 10 个最有信息量的特征（词项）：

def show_top_n_features(clf, vectorizer, categories, n=10):
    feature_names = vectorizer.get_feature_names_out()

    for i, category in enumerate(categories):       
        top_n = np.argsort(clf.feature_log_prob_[i])[-n:]
        print(f"{category}: {' '.join(feature_names[top_n])}")

show_top_n_features(clf, vectorizer, categories) 

我们得到的输出是：

alt.atheism: islam atheists say just religion atheism think don people god
comp.graphics: looking format 3d know program file files thanks image graphics
comp.os.ms-windows.misc: card problem thanks driver drivers use files dos file windows
comp.sys.ibm.pc.hardware: monitor disk thanks pc ide controller bus card scsi drive
comp.sys.mac.hardware: know monitor does quadra simms thanks problem drive apple mac
comp.windows.x: using windows x11r5 use application thanks widget server motif window
misc.forsale: asking email sell price condition new shipping offer 00 sale
rec.autos: don ford new good dealer just engine like cars car
rec.motorcycles: don just helmet riding like motorcycle ride bikes dod bike
rec.sport.baseball: braves players pitching hit runs games game baseball team year
rec.sport.hockey: league year nhl games season players play hockey team game
sci.crypt: people use escrow nsa keys government chip clipper encryption key
sci.electronics: don thanks voltage used know does like circuit power use
sci.med: skepticism cadre dsl banks chastity n3jxp pitt gordon geb msg
sci.space: just lunar earth shuttle like moon launch orbit nasa space
soc.religion.christian: believe faith christian christ bible people christians church jesus god
talk.politics.guns: just law firearms government fbi don weapons people guns gun
talk.politics.mideast: said arabs arab turkish people armenians armenian jews israeli israel
talk.politics.misc: know state clinton president just think tax don government people
talk.religion.misc: think don koresh objective christians bible people christian jesus god

大多数词似乎与其对应的类别有很强的相关性。然而，也有一些像“just”和“does”这样的通用词不提供有价值的信息。这表明我们的模型可能通过更好的停用词列表来改进。实际上，Scikit-Learn 建议不要使用其默认列表，并引用其文档：“‘english’存在一些已知问题，你应考虑其他替代方案。” 😲

总结

让我们总结一下朴素贝叶斯与其他分类模型的优缺点：

优点：

• 训练和预测速度极快

• 提供类别概率估计

• 可用于二分类和多分类问题

• 需要少量的训练数据来估计其参数

• 高度可解释

• 高度可扩展（参数数量与特征数量线性相关）

• 在高维数据上表现良好

• 对噪声具有鲁棒性（噪声样本在估计条件概率时被平均处理）

• 可以处理缺失值（计算特征的似然时忽略缺失值）

• 没有超参数需要调整（除了平滑参数，通常不做更改）

缺点：

• 依赖于朴素贝叶斯假设，而该假设在许多实际领域并不成立

• 特征之间的相关性可能会降低模型的性能

• 通常被更复杂的模型超越

• 零频率问题：如果一个分类特征在训练集中未出现过，其类别将被模型赋予零概率。平滑处理可以缓解这一问题，但不能完全解决。

• 无法处理连续属性，除非进行离散化或对其分布做出假设

• 仅能用于分类任务

最终说明

这是我在 Medium 上写的最长的一篇文章。希望你阅读时的感受至少与我写作时一样愉快。如果有任何不清楚的地方，请在评论中告知我。

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

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

20 个新闻组数据集的信息：

• 引用： Mitchell, Tom (1999). Twenty Newsgroups. UCI 机器学习库。 doi.org/10.24432/C5C323.

• 许可证： Creative Commons CC BY 4.0。

从头开始的朴素贝叶斯分类器，使用 Python

原文：towardsdatascience.com/naive-bayes-classifier-from-scratch-with-python-942708211470

从理论到实践，运用贝叶斯定理

Piero Paialunga

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

–

由Joel Abraham拍摄，图片来源于Unsplash

数学和物理充满了定理、方程、原理、公理和推论。当我开始学习物理时，我记得我达到了所有课程都具有相同结构的阶段：

A. 定义基本假设

B. 使用数学构建下一个“砖块”

C. 一块块地叠加，直到所有部分汇聚成一个优雅、美丽的世界模型

让我们从我学习过的第一门物理课程开始：微积分。

1. 你从集合和数字的基本假设开始。你开始定义自然数、整数、实数和复数。

2. 从这里开始，你定义的函数不过是从空间 A（假设是 N 维实数空间）到空间 B（假设是 1 维实数空间）的映射。

3. 然后你开始研究函数。于是你开始分析它们的最小值、最大值和鞍点。你偶然（哦！）了解了**“导数”**的概念。

4. 然后你看看如何积分一个函数，那是导数的反过程。

5. 然后你将这些与微分方程结合起来。

不要误解我：这个过程是惊人的。我喜欢看到人类逻辑能带你多远。我喜欢看到非常复杂的自然事件可以通过从非常简单的概念开始逐步推导出更深的含义。

另一个宝贵的科学课程是，最初应用于**“A”**的定理，也可以应用于“B”、“C”、“D”和“E”。

令人着迷的方面在于，领域“A”和其他领域（B、C、D 和 E）不必相关。

可能最好的例子就是贝叶斯定理。

贝叶斯定理在某种程度上是一个基础且显而易见的概念。令人难以置信的酷事是，我们可以通过强调贝叶斯定理背后的理念并将其应用于机器学习来构建一些非常有趣的算法。这个工具被称为朴素贝叶斯分类器。

在这篇文章中：

1. 我们将简要介绍贝叶斯定理。我们将解释它是什么，为什么重要，以及它如何应用于机器学习。

2. 我们将看到贝叶斯定理在一个虚构的分类任务中的应用。

3. 我们将看到一个升级版的贝叶斯定理，使用所谓的高斯朴素贝叶斯分类器。

我非常兴奋。让我们开始吧！

1. 关于贝叶斯定理

有一个关于贝叶斯定理的定义一直深深地印在我的脑海里：

“贝叶斯定理是这样的定理，它证明了仅仅因为你的车是蓝色的，并不意味着全世界的车都是蓝色的但如果全世界的车都是蓝色的那么你的车必须是蓝色的。”

如果你阅读了这篇文章关于贝叶斯神经网络，你可能会认识到这与之前的定义相同，但我发誓我没有因为懒惰而重复使用同一定义！我的意思是……我确实很懒，但这不是原因。

我使用那个定义是因为它帮助我们理解事件 A 在事件 B 的条件下的发生概率与事件 B 在事件 A 的条件下的发生概率是不一样的。

让我做一个新的例子来克服我的懒惰。

假设我们有两个碗

现在，假设这两个碗里满是篮球和足球（称为足球，不是足球）球。

作者图片

现在，假设我问你这个问题：

“知道我从蓝色碗里挑了一个球，挑到一个足球球的概率是多少。”

好吧，答案很简单。概率是1，因为在蓝色碗里，我只有足球球。现在假设我随机挑选一个碗（挑选白色碗和挑选蓝色碗的概率都是 0.5）。让我问你这个问题：

“知道我挑了一个足球，选择这个球来自蓝色碗的概率是多少？”

如你所见，问题有点类似，但事件 A**（我挑了一个足球）和事件 B（我挑了蓝色碗）**的顺序正好相反。

现在，我想你可能已经猜到答案并不相同，因为我可能从白色碗中而不是蓝色碗中提取了足球。我也认为你可能猜到这个概率仍然很高，因为我在白色碗中只有一个足球。

让我们在Python中进行实验。

我可以逐行详细解释，但我觉得那样很无聊。我所做的仅仅是设置碗和球的情况，就像之前一样，然后运行概率实验N次。在这 N 次迭代结束时，我们将得到 N 个结果。然后我们需要将概率计算为频率主义概率。在这种情况下，我们将计算：

图片由作者提供

现在我们知道，如果我们在相同的实验中运行 N = 无限次，我们将收敛到分析结果，因此我们将迭代地增加 N，看看它是否收敛到任何结果。

如你所见，我们进行了 20,50,70,…, 1M 次迭代。

这个结果是什么？

图片由作者提供

嗯，看来确实存在一个分析结果，对吗？让我们来看看。

所以，贝叶斯定理告诉我们：

图片由作者提供

因此，提取到的足球来自 蓝色碗的概率等于提取一个足球从蓝色碗（注意区别！这是在你选择了蓝色碗的前提下提取一个足球的概率)的概率乘以从蓝色碗提取的概率，除以提取足球的概率。

我们知道：

图片由作者提供

因为蓝色碗中只有足球。我们还假设我们以相等的概率选择两个碗中的一个。

所以我们有：

图片由作者提供

现在，提取足球的概率是从白色碗中提取足球的概率和从蓝色碗中提取足球的概率之和。例如：

图片由作者提供

所以我们有：

图片由作者提供

因为：

图片由作者提供

所以我们有：

图片由作者提供

如果我们在之前绘制的图上绘制 5/6，我们会得到：

图片由作者提供

我们的数学计算是正确的 😄

好吧，这时我相信你们都在想：

“这与机器学习有什么关系？”

让我们来看看 😏

2. 朴素贝叶斯分类器

朴素贝叶斯分类器是朴素地将贝叶斯定理应用于机器学习分类器的过程：就是这么简单。

假设我们有一个特定的二分类问题（类别 1 和类别 2）。你有 N 个实例，每个实例都有其标签 Y。

所谓的先验概率定义如下：

图片由作者提供

现在我们真正想知道的是：给定一个特定的实例，那个实例属于类别 1 的概率是多少？属于类别 2 的概率是多少？

所以我们感兴趣的是，给定一个实例 x：

图片由作者提供

在二元数据集中，我们知道这两个数量的总和是 1，因此实际上我们只需要其中一个。现在这个数量可能看起来有些神秘，但另一个：

图片由作者提供

可以非常容易地计算（记住碗/球的例子！）。

计算 P(x)的概率也非常简单，就像我们计算 P(class 1)和 P(class 2)一样。

所以我们实际上可以计算：

图片由作者提供

这一切看起来都很好。现在我想你可能已经理解了我们为什么称之为朴素。它之所以朴素，是因为它不过是计数出现次数并使用贝叶斯逻辑来推断预测。

让我们将其应用到一个玩具数据集上。我们有这两个类别和两个特征，是使用sklearn.datasets库生成的。

太棒了。

所以：

图片由作者提供

因为我们有两个维度。

给定类别 0，我们可以轻松计算给定区域（2D 点周围的小区域）的概率。

所以：

1. 我们知道如何计算后验概率，即 P(Y|X)

2. 我们有我们的 Y 和 X 集合

3. 我们可以将规则应用于训练集，并且可以预测测试集的结果。

让我们开始吧！

1. 训练-测试分割：

2. 导入和拟合朴素贝叶斯分类器：

请注意，我们对数据（x）添加了一个偏置。这是因为我们使用的贝叶斯分类器将 x_1 特征视为分类特征，而我们将这些特征视为数值特征。一种更严格的方法是使用*LabelEncoder()*特征，但在这个特定情况下，它们是完全等价的。

3. 测试性能：

在这个非常简单的玩具示例中，性能显然是（几乎）完美的。

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

你可能已经听够了我讲的内容，我对此感到抱歉。

我只是想通过谈论高斯朴素贝叶斯来结束这篇文章。

在前面的例子中，我们将特征视为分类的。但如果它们不是呢？

好吧，我们必须假设一个特定的可能性分布：P(Y|X)。

图片由作者提供

所以它叫做高斯朴素贝叶斯，因为如上所见，似然被认为是高斯的。均值（mu_y）和方差（sigma_y 平方）是通过计算类别的均值和方差得到的。

请注意，现在我们不再需要在数据中添加偏差。

这些是结果

4. 总结

在这篇文章中，我们：

1. 承认贝叶斯定理的存在。我们对它进行了直观的介绍，并进行了一个非常简单的受控案例来查看它是如何工作的。

2. 了解贝叶斯定理如何应用于机器学习。什么是 Y，什么是 X，我们如何将它们放入贝叶斯公式中以在分类任务中进行一些预测。

3. 我们自己创建了一个数据集，并使用Categorial Naive Bayes类（来自sklearn）进行了一项玩具分类任务。它几乎完美地工作，但有一个问题，就是只能处理分类特征。

4. 我们了解了如何将分类朴素贝叶斯扩展为Gaussian Naive Bayes。正如我们所见，Gaussian Naive Bayes 使用高斯似然来工作。这种高斯似然也适用于非分类特征。

贝叶斯定理在机器学习中以多种方式应用，如贝叶斯神经网络或贝叶斯岭回归。我觉得这是一个很酷的入门示例，展示了贝叶斯定理在分类问题中的应用。我写这篇文章非常开心，希望你们喜欢它！🥰

5. 结论

如果你喜欢这篇文章，想了解更多关于机器学习的内容，或者只是想问我一些问题，你可以：

A. 在 Linkedin 上关注我，我会发布我的所有故事。

B. 订阅我的 新闻通讯。它将让你了解新故事，并给你机会通过短信与我联系，获取你可能有的所有更正或疑问。

C. 成为 会员，这样你就没有“每月最大故事数”的限制，可以阅读我（和其他数千名机器学习和数据科学顶级作者）关于最新技术的所有文章。

从零开始的朴素贝叶斯与 TensorFlow

原文：towardsdatascience.com/naive-bayes-from-scratch-with-tensorflow-6e04c5a25947

概率深度学习

路易斯·罗克

·发表于 Towards Data Science ·10 分钟阅读·2023 年 1 月 18 日

–

介绍

本文属于“概率深度学习”系列。该系列每周涵盖深度学习的概率方法。主要目标是扩展深度学习模型以量化不确定性，即了解它们不知道的内容。

在本文中，我们对使用葡萄酒样本数据集的朴素贝叶斯分类算法进行了考察。朴素贝叶斯算法是一种基于贝叶斯定理的概率机器学习技术，假设在给定目标标签的情况下特征之间是独立的。为了便于可视化类别的分离，我们将模型限制为仅使用两个特征。

我们的目标是基于选定的特征对葡萄酒样本进行分类。为实现这一目标，我们首先探索数据并选择有效区分各类别的特征。然后，我们构建类别先验分布和类别条件密度，从而能够预测具有最高概率的类别。该研究使用的数据集包含葡萄酒的各种特征，如色调、酒精、类黄酮以及一个目标类别，并且数据集来自 scikit-learn 库 [1]。

迄今为止发表的文章：

1. TensorFlow Probability 的温和介绍：分布对象

2. TensorFlow Probability 的温和介绍：可训练参数

3. 从零开始的最大似然估计，使用 TensorFlow Probability

4. 从零开始的概率线性回归，使用 TensorFlow

5. Tensorflow 中的概率性与确定性回归

6. 频率学派与贝叶斯统计在 Tensorflow 中的应用

7. 确定性与概率性深度学习

8. 从头开始使用 TensorFlow 实现朴素贝叶斯

图 1：今天的格言：在葡萄酒分类方面要保持天真？ (source)

我们使用 TensorFlow 和 TensorFlow Probability 开发我们的模型。TensorFlow Probability 是一个建立在 TensorFlow 之上的 Python 库。我们将从 TensorFlow Probability 中的基本对象开始，理解如何操作它们。接下来的几周里，我们将逐步增加复杂性，并将我们的概率模型与现代硬件（例如 GPU）上的深度学习结合起来。

如常，代码可以在我的GitHub上找到。

探索性数据

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline

import tensorflow as tf
import tensorflow_probability as tfp
tfd = tfp.distributions

from sklearn.metrics import accuracy_score
from sklearn import datasets, model_selection
from sklearn.datasets import load_wine

我们的目标是调查使用朴素贝叶斯算法根据选择的特征对葡萄酒样本进行分类。为了实现这一目标，我们首先进行数据的探索性分析。让我们开始识别两个有效区分目标变量的特征，并利用它们预测葡萄酒的类别。

dataset = load_wine(as_frame=True)
dataset = pd.concat((dataset['data'], dataset['target']), axis=1)
sns.pairplot(dataset[['alcohol','alcalinity_of_ash', 'flavanoids', 'color_intensity', 'hue', 'target']],hue='target');

图 2：分析葡萄酒数据集中的特征对。

酒精度和色调是有效区分类别的特征。因此，这些就是我们将用于构建朴素贝叶斯模型的两个特征。

sns.jointplot(x='alcohol',y='hue', hue='target', data=dataset);

图 3：按酒精度和色调分布的目标样本。

我们现在可以将数据分成训练集和测试集。

data = dataset[['alcohol', 'hue']].to_numpy()
targets = dataset[['target']].to_numpy()

label_colors = ['darkred', 'peachpuff', 'black']
x_train, x_test, y_train, y_test = model_selection.train_test_split(data, targets, test_size=0.2)

朴素贝叶斯分类器

朴素贝叶斯是一种广泛使用的概率机器学习算法，基于贝叶斯定理。它特别适用于分类任务，并以其简单性和高效性著称。尽管名字中有“朴素”，但特征之间的“朴素”独立性假设并不总是限制，并且在实践中通常能取得良好结果。在这篇文章中，我们将全面回顾朴素贝叶斯算法及其变体，并从基本原理上实现它。

我们首先简要介绍贝叶斯定理，这是朴素贝叶斯算法的基础。贝叶斯定理表明，在给定一些证据（E）的情况下，假设（H）的概率与假设的先验概率乘以证据在给定假设下的似然性成正比。朴素贝叶斯算法使用这个定理通过计算每个类别的后验概率来分类新实例，然后选择概率最高的类别。

朴素贝叶斯算法的基本原理是假设给定实例的特征在给定类别标签的情况下是条件独立的。这一假设，也称为“朴素”假设，使得算法在计算上更为高效，因为它减少了需要估计的参数数量。然而，当特征实际上并非独立时，这也可能导致准确率下降。

朴素贝叶斯算法有几种变体，每种都适用于不同类型的数据。例如，高斯朴素贝叶斯用于连续数据，而多项式朴素贝叶斯用于离散数据。伯努利朴素贝叶斯用于二元数据。在这种情况下，我们将专注于实现高斯朴素贝叶斯。

朴素贝叶斯算法已被应用于广泛的领域，包括自然语言处理、计算机视觉和生物信息学。在自然语言处理领域，它通常用于文本分类，例如垃圾邮件检测和情感分析。在计算机视觉中，它用于图像分类和目标检测。在生物信息学中，它用于蛋白质分类和基因预测。

正如我们上述所述，朴素贝叶斯分类器基于贝叶斯规则：

其中 𝑋 是输入特征，𝑌 是输出类别，𝐾 是类别的数量。更具体地说，𝑃(𝑌) 表示类别先验分布，𝑃(𝑋|𝑌) 是输入的类别条件分布，而 𝑃(𝑌|𝑋) 是给定输入特征的类别概率。

独立性假设大大简化了算法，因为我们不需要估计完整的联合分布 𝑃(𝑋|𝑌=𝑦𝑘)。相反，类别条件分布可以写作：

其中 𝑓 表示特征的数量。

先验

在朴素贝叶斯算法中，类别先验分布是一个概率分布，描述了训练数据中每个类别的概率。它是算法的一个基本组成部分，因为它用于计算给定一些证据的类别后验概率。

类别先验分布定义为给定训练数据中实例总数的类别概率。它通常表示为 𝑃(𝑌=𝑦𝑘), 其中 𝑘 是类别标签。类别先验分布通过训练数据中每个类别的相对频率来估计。例如，如果训练数据中有 100 个实例，其中 60 个属于类别 A，那么类别 A 的先验概率估计为 P(Y=A) = 0.6。

类别先验分布在朴素贝叶斯算法中起着至关重要的作用，因为它用于计算在给定一些证据的情况下某一类别的后验概率。后验概率的计算是将类别先验和在给定类别下证据的似然度相乘，并通过证据的边际似然度进行归一化。换句话说，类别先验分布作为一个权重因子，调整似然函数的相对重要性。

然而，如果类别先验分布是从有偏的训练数据中估计的，它可能会导致算法性能不佳，特别是当测试数据来自不同的分布时。这被称为类别不平衡问题，可以通过使用过采样、欠采样或合成数据生成等技术来缓解。

类别先验分布是属于类别 𝑘 的数据示例的比例。我们可以将其写成以下形式：

其中，𝑛 表示第 𝑛 个数据集示例，𝑁 是数据集中示例的总数，𝛿 是克罗内克δ函数（当类别匹配时返回 1，否则返回 0）。它返回一个与 𝑃(𝑌=𝑦𝑘) 相对应的分类分布。

def prior_fn(y):
    n_classes = np.unique(y).shape[0]
    counts = np.zeros(n_classes)
    for c_k in range(n_classes):
        counts[c_k] = np.sum(np.where(y==c_k, 1, 0))
        priors = counts/np.sum(counts)
    dist = tfd.Categorical(probs=priors)
    return dist

prior = prior_fn(y_train)
prior

<tfp.distributions.Categorical 'Categorical' batch_shape=[] event_shape=[] dtype=int32>

让我们绘制我们的先验分布。

plt.bar([0, 1, 2], prior.probs.numpy(), color=label_colors)
plt.xlabel("Class")
plt.ylabel("Prior probability")
plt.title("Class prior distribution")
plt.xticks([0, 1, 2],)
plt.show()

图 4：每个目标类别的先验概率。

似然性

在朴素贝叶斯算法中，类别条件密度是描述给定类别标签下每个特征的似然性的概率分布。它们用于计算在给定一些证据的情况下某一类别的后验概率，是算法的一个基本组成部分。类别条件密度定义为给定类别标签下每个特征的概率密度函数（pdf）。它们通常表示为 𝑃(𝑋𝑖|𝑌=𝑦𝑘), 其中 𝑋𝑖 是一个特征，𝑘 是类别标签。类别条件密度是通过各种技术从训练数据中估计出来的，具体取决于数据的类型。例如，对于连续数据，类别条件密度可以使用高斯分布估计，而对于离散数据，可以使用多项式分布或伯努利分布进行估计。正如我们之前所述，在我们的案例中，我们有连续特征，因此我们将探索高斯方法。

类别条件密度在朴素贝叶斯算法中起着关键作用，因为它们用于计算给定类别标签下证据的似然性。这一似然性是通过评估证据的每个特征的类别条件密度来计算的，然后将它们相乘。类别条件密度作为一个权重因子，调整每个特征在分类任务中的相对重要性。

现在是定义𝑃(𝑋|𝑌)——输入的类别条件分布的时候了。在这种情况下，我们使用单变量高斯分布（请记住独立性假设）：

其中𝜇𝑖𝑘和𝜎𝑖𝑘是需要估计的参数。使用最大似然估计，估计值就是每个类别样本数据点的均值和方差：

def class_conditionals_fn(x, y):
    n_classes = np.unique(y).shape[0]
    n_features = x.shape[1]
    counts = np.zeros(n_classes)
    mean_feature_given_class = []
    std_feature_given_class = []
    for c_k in range(n_classes):
        mean_feature_given_class.append(np.mean(x[np.squeeze(y==c_k)], axis=0))
        std_feature_given_class.append(np.std(x[np.squeeze(y==c_k)], axis=0))

    class_cond = tfd.MultivariateNormalDiag(loc = np.asarray(mean_feature_given_class).reshape(n_classes, n_features),
                             scale_diag=np.asarray(std_feature_given_class).reshape(n_classes, n_features))

    return class_cond

class_conditionals = class_conditionals_fn(x_train, y_train)
class_conditionals

<tfp.distributions.MultivariateNormalDiag 'MultivariateNormalDiag' batch_shape=[3] event_shape=[2] dtype=float64>

下方的等高线图展示了类别条件密度。请注意每个分布的等高线如何对应于具有对角协方差矩阵的高斯分布，因为模型假设在给定类别的情况下每个特征是独立的。

def contour_plot(x0_range, x1_range, prob_fn, batch_shape, colors, levels=None, num_points=100):
    x0 = np.linspace(x0_range[0], x0_range[1], num_points)
    x1 = np.linspace(x1_range[0], x1_range[1], num_points)
    X0, X1= np.meshgrid(x0, x1)
    Z = prob_fn(np.expand_dims(np.array([X0.ravel(), X1.ravel()]).T, 1))
    Z = np.array(Z).T.reshape(batch_shape, *X0.shape)
    for batch in np.arange(batch_shape):
        if levels:
            plt.contourf(X0, X1, Z[batch], alpha=0.2, colors=colors, levels=levels)
        else:
            plt.contour(X0, X1, Z[batch], colors=colors[batch], alpha=0.3)

plt.figure(figsize=(10, 6))
plot_data(x_train, y_train, alpha=0.3)
x0_min, x0_max = x_train[:, 0].min()-0.2, x_train[:, 0].max()+0.2
x1_min, x1_max = x_train[:, 1].min()-0.2, x_train[:, 1].max()+0.2
contour_plot((x0_min, x0_max), (x1_min, x1_max), class_conditionals.prob, 3, label_colors)
plt.title("Training set with class-conditional density contours")
plt.show()

图 5：带有类别条件密度等高线的训练集。

在执行上述计算后，算法的最后一步是预测新的数据输入𝑋̃ :=(𝑋̃ 1,…,𝑋̃ 𝑓)的类别𝑌̂。可以通过以下方式完成：

def predict_class(prior, class_conditionals, x):
    log_prob_list = []
    for sample in x:
        cond_probs = class_conditionals.log_prob(sample)
        joint_likelihood = tf.add(prior.probs.numpy(), cond_probs)
        norm_factor = tf.math.reduce_logsumexp(joint_likelihood, axis=-1, keepdims=True)
        log_prob = joint_likelihood - norm_factor
        log_prob_list.append(log_prob)
    return np.argmax(np.asarray(log_prob_list), axis=-1)

predictions = predict_class(prior, class_conditionals, x_test)

结果

在这篇文章中，我们应用了朴素贝叶斯算法来根据选定的特征对葡萄酒样本进行分类。具体来说，我们使用了两个特征：色调和酒精，来预测葡萄酒的类别。我们的结果表明，该模型在这项任务中的准确率超过了 91%。

accuracy = accuracy_score(y_test, predictions)
print("Test accuracy: {:.4f}".format(accuracy))

Test accuracy: 0.9167

为了进一步分析模型的性能，我们还绘制了模型的决策区域，即分隔不同类别的边界。决策区域有助于可视化算法执行的类别分离。如我们所见，模型能够相当有效地分隔数据集中的三个类别。

值得注意的是，朴素贝叶斯算法假设特征之间是独立的，这在实际场景中可能不成立。特征之间的相关性可以帮助提高模型的准确性。因此，考虑在模型特征之间引入相关性可能有助于提升性能。此外，还可以考虑允许特征之间相关性的其他算法来改善结果。

plt.figure(figsize=(10, 6))
plot_data(x_train, y_train)
x0_min, x0_max = x_train[:, 0].min()-0.2, x_train[:, 0].max()+0.2
x1_min, x1_max = x_train[:, 1].min()-0.2, x_train[:, 1].max()+0.2
contour_plot((x0_min, x0_max), (x1_min, x1_max), 
             lambda x: predict_class(prior, class_conditionals, x), 
             1, label_colors, levels=[-0.5, 0.5, 1.5, 2.5, 3.5],
             num_points=200)
plt.title("Training set with decision regions")
plt.show()

图 6：训练集决策区域。

结论

在这篇文章中，我们使用 TensorFlow Probability 从零开始实现了朴素贝叶斯算法。我们将其应用于使用葡萄酒样本数据集的分类任务。我们选择了两个特征，色调和酒精，来预测葡萄酒的类别，并且达到了超过 91%的准确率。我们还可视化了模型的决策区域，这有助于理解算法执行的类别分离。

这个简单的例子展示了朴素贝叶斯算法在分类任务中的简单性和有效性。然而，朴素贝叶斯算法假设特征之间是独立的，这在实际场景中可能并不成立。

参考文献和资料

[1] — 葡萄酒数据集

[2] — Coursera: 深度学习专业课程

[3] — Coursera: TensorFlow 2 深度学习 专业课程

[4] — TensorFlow 概率指南与教程

[5] — TensorFlow 博客中的 TensorFlow 概率帖子

命名实体与新闻

原文：towardsdatascience.com/named-entities-and-the-news-ac9009601095?source=collection_archive---------5-----------------------#2023-07-04

探索命名实体在荷兰新闻数据集中的应用

Felix van Deelen

·

关注 发表在 Towards Data Science ·10 min read·Jul 4, 2023

–

一个新闻推荐系统的例子，确实可以从命名实体识别（NER）中受益。来源：文章由NOS提供，照片由Rick L在Unsplash（左），作者创建的图像（中），文章由NOS，照片由Cristina Anne Costello在Unsplash（右）。

在 NOS——荷兰公共广播基金会——我们的编辑团队每天撰写数百篇新闻文章。这些文章向荷兰公民提供新闻，同时从自然语言处理的角度来看，也形成了一个有趣且高质量的数据集。在这篇博客中，我作为 NOS 的数据科学家，报告了通过将命名实体识别（NER）应用于我们荷兰新闻文章的数据集所进行的几个实验，并提出了在新闻背景下应用 NER 的几个想法。

什么是命名实体？

命名实体（NE）是一种特殊的词，指代具有专有名称的现实世界对象，例如人物、地点或组织。存在自动识别这些类型词语的模型，这些模型称为命名实体识别（NER）模型。右侧图示展示了一个应用于我们文章摘录的 NER 模型，其中 NE 被突出显示并标注了 NE 类型。

在荷兰，有一些预训练模型可用，如spaCy [1]、Flair [2] 或NTLK [3]。我们对这三种模型进行了定性评估，通过将它们应用于我们文章的随机样本并手动检查结果。由此我们决定在剩余实验中使用 spaCy。此模型可能识别的所有 NE 类型的概述见下图 1 左侧。

图 1：由spaCy提供的 NER 模型中的 NE 类型概述（左）。应用于从荷兰语翻译的新闻文章摘录的 NER 示例（右）。

使用来自 spaCy 的预训练模型，我们对数据集的几个子集应用了 NER。我们首先收集了一个月（2023 年 2 月）的所有文章，将数据分为新闻和体育（分别为 1,030 篇和 596 篇），然后应用 NER 以获得每种 NE 类型的总频率计数。新闻和体育的结果显示在图 2 中，立即展示了 NE 在新闻中的重要性。可以看出，仅在一个月的文章中，就提到了数万条 NE。为提供一个视角，平均每篇文章包含 404 个词，大约 10%的词是 NE。下图也显示了新闻和体育中最常提到的 NE 类型有所不同。对于新闻，大多数 NE 类型是国家，其次是组织和人物。而对于体育，最常提到的 NE 类型是人物，其次是国家和数字。这可能是因为体育涉及比分（基数）和个人运动员（人物），而新闻报道事件时，通常需要提到地点（gpe）。

图 2：新闻文章（左）和体育文章（右）中检测到的 NE 类型频率计数。

NER 为我们的数据提供了新的视角

我们对 2022 年世界杯足球赛的所有文章进行了案例研究，总共包含 482 篇文章。NER 被应用于数据集，以检测所有类型为Person的 NE。发现 2,171 个独特的 NE，其中 1,296 个只被提到过一次。在图 3A 中，我们展示了这一事件中最常提到的人物的概述。此外，对于最常提到的人物，我们创建了一个流图，以展示提及频率随时间的发展，见图 3B。这表明，如范加尔在整个比赛中频繁被提到，而其他人物则主要在特定日期被提到。这类图表可能为我们的编辑团队提供新的洞察，因为它们是 NOS 报道内容的定量反映。这些洞察是通过 NER 高效提供的。目前我们特别应用于世界杯 22，但可以想到许多不同的设置，这些图表可能会很有趣。例如，想一想在选举期间提到哪些政治家或政党，或更普遍地，某个较长时间范围内国家、城市、组织等的提及频率。

图 3A & B：通过 NER 获得的世界杯 22 中提到的人物频率计数。左侧显示总数，右侧显示随时间的发展。

关于[您的命名实体这里]

我们将案例研究扩展到使用所有 2022 年世界杯的文章，并提出了“我们能否利用 NER 为命名实体生成摘要？”的问题。我们首先开发了一个模块，该模块收集所有提到给定 NE 的文章，这可以作为用户特别感兴趣的 NE 的所有信息集合。但更有趣的是，该模块收集了所有提到 NE 的句子，从而生成了该集合的摘要。例如，我们将该模块应用于安德里斯·诺普特，荷兰国家队的守门员。从图 3 中可以看出，诺普特在赛事中被提及的频率相当高。应用该模块后，生成了一个很好的概述，展示了我们守门员的非凡故事，如下所示，翻译自荷兰语。

-------------------------------------------------- -------------------------------------------------- --------------------
2022-11-11
   - Noppert joining as a penalty killer?
-------------------------------------------------- -------------------------------------------------- --------------------
2022-11-16
   - sc Heerenveen goalkeeper Andries Noppert is the nineteenth premier league player in Qatar.
-------------------------------------------------- -------------------------------------------------- --------------------
2022-11-20
   - 'Don't worry about Qatar and Ecuador' and 'Failure on goal is a gamble' Analysts Leonne Stentler and Pierre van Hooijdonk agree.
   - Van Gaal does not say anything about Noppert's base place, but hints at Gakpo 'at 10' According to various media, 28-year-old Andries Noppert, who plays for sc Heerenveen, would make his debut for the Orange squad against Senegal on Monday.
-------------------------------------------------- -------------------------------------------------- --------------------
2022-11-21
   - Is Noppert the base goalkeeper now?
   - Noppert: 'This is what you dream of as a boy' Goalkeeper Andries Noppert turned out not to suffer from stage fright against Senegal.
   - Will Noppert succeed first World Cup debutant Schoenaker?
   - Goalkeeper Andries Noppert makes his debut in Orange and can look back on a successful first international match.
-------------------------------------------------- -------------------------------------------------- --------------------
2022-11-22
   - 'Disarming' Noppert takes the stage: 'In the Netherlands we are all whining' The 28-year-old goalkeeper of sc Heerenveen made his debut on Monday in the World Cup match against Senegal in the Dutch national team.
-------------------------------------------------- -------------------------------------------------- --------------------
2022-11-23
   - Noppert?
-------------------------------------------------- -------------------------------------------------- --------------------
2022-11-24
   - The Foggia episode of Orange keeper Noppert: 'He smoked like a chimney' Andries Noppert is suddenly a well-known Dutchman after the World Cup match of the Netherlands against Senegal.
-------------------------------------------------- -------------------------------------------------- --------------------
2022-11-25
   - Jurriën Timber, Virgil van Dijk and Nathan Aké had their defenses well organized and Andries Noppert once again proved to be a reliable goalkeeper.
-------------------------------------------------- -------------------------------------------------- --------------------
2022-12-03
   - View the reactions of Virgil van Dijk and Andries Noppert here: In that team, one of the important players is just back in his familiar spot in the attack.
   - Andries Noppert made a good save with his left leg.
-------------------------------------------------- -------------------------------------------------- --------------------
2022-12-07
   - Noppert lives soberly towards Argentina: 'Messi can also miss penalties, can't he?'
-------------------------------------------------- -------------------------------------------------- --------------------
2022-12-09
   - So yes..." Noppert's fairy tale ended It could have been so beautiful for sc Heerenveen goalkeeper Andries Noppert, but the keeper on the other side, Emiliano Martinez, became the great hero.
   - The Argentinian wingback Molina ran away from the back of his Dutch colleague Blind, Virgil van Dijk was just too late to correct and Molina passed Andries Noppert.
-------------------------------------------------- -------------------------------------------------- --------------------
2022-12-16
   - Six striking World Cup facts: Amrabat conquers, Modric dribbles, Noppert saves Remarkable statistics everywhere during the World Cup in Qatar.
-------------------------------------------------- -------------------------------------------------- --------------------
2022-12-18
   - Andries Noppert (Netherlands) Vermeulen: "The same goes for Noppert, of course.

一个 NE 感知的推荐系统

到目前为止，我们已经看到 NE 在新闻文章中非常普遍，应用 NER 可以提供一些有趣的见解。我们认为还有一个实验很值得在此博客中分享，涉及研究问题“我们能否利用 NER 改进基于内容的推荐系统？”。之前我们开发了一个基于内容的推荐系统，该系统最近已集成到我们的新闻应用中。通过在线和离线测试，我们比较了各种模型和优化，现在我们观察到应用中的点击率有所增加。这一切都是好消息，但我们一直在寻找进一步改进推荐系统的方法。我们收到编辑团队的反馈，称推荐系统对包含人名或地点名的文章感到困惑，这些名字在荷兰语中也是常见词汇。在接下来的部分中，我们报告了一个使用 NER 解决这种歧义的实验。

实验

我们当前的推荐系统基于使用 TF-IDF 进行文本向量化的余弦相似度。这基本上意味着它在识别相似文章时严重依赖词语重叠，但对稀有词语赋予更高的相关性。可以想象，当词语具有多重含义时，这种方法的效果可能不佳，这在 NE 的情况下尤为如此。举个例子，考虑一篇关于高尔夫球手泰格·伍兹的文章：一个基本的推荐系统可能会找到提到动物老虎或关于树林的文章。这些显然不会是有用的推荐。我们假设通过在推荐系统中引入 NE 感知，即通过按类型标注文本中的 NE，可以解决这个问题。在这种情况下，令牌将不再重叠，如图 5 所示。

图 5：当前推荐系统（基础）与 NE 感知系统的示例对比。当前系统因两篇文章中都提到了“tiger”一词而将这两篇文章关联起来，而 NE 感知系统则解决了这一歧义。来源：文章由NOS提供，照片由Rick L在Unsplash拍摄（左），文章由NOS提供，照片由Cristina Anne Costello在Unsplash拍摄（右）。

我们使用 NE 类型person、location、organisation及其组合来实现 NE 感知。我们使用由编辑团队手动标注的测试集来评估各种模型，该测试集包含了哪些文章相关的信息。该测试集包含 14,541 篇独特的文章，平均每篇文章链接到约 2 篇其他文章。作为评估指标，我们计算了在排序推荐中的策划链接文章的平均排名。

图 6 显示了我们的基础模型与各种 NE 感知模型的结果。可以看出，实际上我们的基础模型优于所有类型的 NE 感知模型。理论上，引入 NE 感知应能改善推荐系统，但实际上我们发现它引入了比解决的歧义更多的歧义。我们详细检查了各种模型的输出，发现我们受限于 NER 模型的性能。spaCy NER 模型在其自身测试集上的 F-score 为 0.77，但在应用于其他数据集时，该得分可能更低，因此模型偶尔可能不准确。从对一些 NE 感知模型输出的手动检查中，我们看到，结合 TF-IDF 时，错误检测到的 NE 的影响相当强。对于错误检测到的 NE 的文章，输出推荐经常包含相同的错误分类 NE。例如，我们看到一篇包含词汇hindsight的文章被错误分类为 Person 类型的 NE，导致推荐中包含相同错误分类的 NEhindsight。虽然在这种情况下 NER 并未启用，但推荐是有意义的，因为 TF-IDF 会给像hindsight_Person这样在语料库中非常稀有的标记分配更高的相关性。我们的结论是，目前预训练的荷兰 NER 模型不够准确，无法融入我们的推荐系统。

图 6：我们当前推荐系统（基础版）与 NE 感知系统的性能比较。

我们未来可能会从自己微调预训练模型中受益。目前，我们探索了另一种解决 NE 模糊性的方法，使用诸如类别和关键词等元数据作为无噪音但关联性较弱的 NE 测量，这显著改善了我们的推荐系统。

结论

在本博客中，我们探讨了将命名实体识别应用于荷兰新闻数据集时可以做些什么。我们发现它在获取数据集的一般见解方面效果良好，例如构建 NE 频率图和流图。然而，当应用于我们的推荐系统时，我们发现模型的准确性不够。虽然引入 NE 感知解决了一些 NE 模糊性问题，但同时也引入了 NE 检测错误的新模糊性。未来我们可能会尝试微调预训练模型或从头开始训练自己的模型，或者如果你有任何建议，请在评论中告诉我们！

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

参考文献 [1] spaCy NER 模型：spacy.io/models/nl#nl_core_news_lg

[2] Flair NER 模型：huggingface.co/flair/ner-dutch-large

[3] NLTK NER 模型：www.nltk.org/book/ch07.html

关于 NOS NOS 是荷兰的一个独立公共媒体组织，通过电视、广播、网站和移动应用等平台报道新闻和体育。我们有专门的专业团队为多个品牌创建数字服务。本博客中描述的研究是在 NOS 数据团队的成员身份下进行的，该团队负责探索新闻背景下新兴数据科学和 AI 技术的应用。

Python 标准库中的 NaN 值

原文：towardsdatascience.com/nan-values-in-the-python-standard-library-798d9ed946c0

PYTHON 编程

NaN 意味着 Not-a-Number。你可以在数值库中使用它——但也可以在 Python 标准库中使用。

Marcin Kozak

·发布在 Towards Data Science ·11 分钟阅读·2023 年 10 月 7 日

–

照片由 cyrus gomez 提供，来源于 Unsplash

NaN 代表 Not-a-Number。因此，NaN 对象表示它的名字所传达的意思——某物不是一个数字。它可以是一个缺失值，也可以是数值变量中的非数值值。由于我们不应该在纯数值容器中使用非数值值，因此我们将这样的值标记为 not-a-number，即 NaN。换句话说，我们可以说 NaN 代表一个 缺失的数值。

在本文中，我们将讨论 Python 标准库中可用的 NaN 对象。

NaN 值在数值数据中经常出现。如果你对这个值的细节感兴趣，你可以在这里找到它们：

[## NaN - 维基百科

在计算中，NaN（代表 Not a Number）是数字数据类型的一种特定值（通常是浮点数）……

en.wikipedia.org](https://en.wikipedia.org/wiki/NaN?source=post_page-----798d9ed946c0--------------------------------)

在本文中，我们不会讨论 NaN 值的所有细节。¹ 相反，我们将讨论一些如何在 Python 中处理 NaN 值的示例。

每种编程语言对 NaN 值有自己独特的方法。在以计算为重点的编程语言中，NaN 值是基本的。例如，在 R 语言中，你有 NULL（相当于 Python 的 None）、NA（表示 not available）和 NaN（表示 not-a-number）：

来自 R 会话的截图。图片由作者提供。

在 Python 中，你有None和多个表示NaN的对象。值得注意的是，Pandas 区分NaN和NaT，后者代表缺失的时间。本文将讨论标准库中的NaN值；主流数值 Python 框架中的NaN（以及NaT）——例如 NumPy 和 Pandas——将在未来的文章中讨论。

如果你没有在 Python 中处理过数值数据，可能根本没有遇到过NaN。然而，NaN值在 Python 编程中无处不在，因此了解如何处理它们非常重要。

Python 中NaN的介绍

当你处理一个list对象时，可以使用数字值和非数字值。因此，它可以是[1, 2, "three"]或[1, 2, None]等。你甚至可以这样做：

>>> [1, "two", list, list(), list.sort]
[1, 'two', <class 'list'>, [], <method 'sort' of 'list' objects>]

所以，列表可以接受任何对象。如果你想对这些列表进行数字计算，可以做到，但需要调整代码：

>>> x = [1, 2, "three"]
>>> sum(x)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: unsupported operand type(s) for +: 'int' and 'str'
>>> sum(xi for xi in x if isinstance(xi, (float, int)))
3

如你所见，非数字对象在这里可以保持原样，你仍然可以进行数字计算——这不是简单或简洁的代码，但它有效。

NaN——非数字——代表一个缺失的数值。

然而，对于只接受特定数字类型对象的容器类型来说情况并非如此。这包括array.array、numpy.array或pandas系列（以及pandas数据框中的数值列）。如果定义为数字，它们不接受非数字值，唯一的例外是：NaN值。

NaN表示一个非数字值，但正如你将很快看到的，它的类型是数字类型——准确地说是float。为什么？原因很简单，因为这样你可以在数字容器中使用NaN。

为什么不去掉它们呢？

为什么不干脆去掉所有这些值呢？干嘛要费这个劲呢？

NaN值的一个常见用例是数据分析和可视化。例如，考虑一个数据集，其中几列有某些行的缺失值。你不能从数据框中删除一个单元格，因此你可以选择保留所有这些行并以某种方式处理缺失值，或者删除所有包含一个或多个NaN值的行。删除缺失值的行是一个常见的做法，但它有代价：它会删除某些列中的非NaN值，通常不明智地丢弃我们已经拥有的信息。

NaN 值的另一个用途是在错误处理。例如，如果一个函数期望一个数值输入，但接收到一个字符串或其他非数值的值，它可能会返回一个 NaN 值以表示输入无效。我们很快会看到一个例子。这使得调用代码可以优雅地处理错误，而不是引发异常或返回意外结果。通过使用 NaN 值来表示错误或缺失的数据，可以对可能包含无效或缺失值的数据集进行计算和处理。你可以返回 None，但在 Python 中，None 可能表示各种不同的含义，而 NaN 传达了更具体的信息，直接与值的数值特性相关——这不是一个数字。

当你处理数值和工具时，你应该知道如何使用 NaN。然而，当你的应用是通用性质的，因此不需要数值框架（如 NumPy 或 Pandas）时，你会发现 NaN 常常可以被忽略，或者用 None 表示。如果这样可以使代码更简单而不牺牲任何功能，考虑这样做。

示例

NaN 值可以意味着各种事情：

• 一个普通的缺失值——它没有提供，未传输，等等。在你的笔记本中，你会将其标记为“NA”或“N/A”：不适用。不适用，即你无法应用它。它是缺失的，因此我们需要将其标记为缺失。你可以使用 NaN。

• 一个函数的结果得到了数值类型参数的错误值。函数返回 NaN，而不是抛出错误。

• 一个错误。这可能是输入错误；信不信由你，输入错误比我们大多数人想象的要频繁得多，它们可能会对后续分析产生很大影响。例如，许多人仍然认为菠菜是铁的良好来源。实际上并不是，那为什么这么多人认为如此呢？这源于一个输入错误——一个错误的十进制点。你可以使用 NaN 来表示有错误的数据元素——除非你确定可以纠正这个错误。

• 一个注释。这是一个被错误地输入到数值变量中的字符串值。当输入数据的人想解释某个特定值缺失的原因时，比如“读取不清”或“我睡过头了”。虽然这些仍然是缺失的数据，但它们提供的信息比单纯的空白值要多。有时这些信息很重要，但有时则不然。然而，对于数值计算来说，这样的注释值通常是微不足道的。因此，如果你需要使用数值容器来表示这个变量，你可以使用 NaN 来表示注释。

这些是四个例子，但也有其他情况。虽然每种情况略有不同，但从数值计算的角度来看，它们都是相同的：这个值不是一个数字。我们需要对它做点什么，使用 NaN 是一个常见的选择。

标准库中的 NaN

Python 提供了几种类型的 NaN 值，我们将在下文中讨论它们。本文重点关注标准库，但请注意，如果你使用数值框架，它很可能有自己对 NaN 的实现（或说是表示）以及与之配套的函数/方法。²

尽管 Python 标准库不是进行数值计算的最佳工具，但它确实提供了数值容器和专用工具。一个数值容器的例子是 array 模块及其 array.array 容器类型。虽然它不是直接使用的最佳工具，但它使你可以高效地使用 Cython，而无需使用像 NumPy 这样的非标准库工具。一个标准库中的专用数值工具的例子是 math 模块：

[## math — 数学函数 — Python 3.10.8 文档

该模块提供了 C 标准定义的数学函数的访问。这些函数不能被使用…

docs.python.org

在 Python 标准库中使用 NaN 值有两种方式：float("nan") 和 math.nan。我读过许多关于 Python 的书，但我不记得看到过这两个值的提及。我的记忆并不完美，但我怀疑即使这些值在一些书中被提到，也没有得到足够的关注。因此，我认为许多数据科学家，甚至是数据科学领域之外的 Python 开发者，可能对 float("nan") 和 math.nan 并不知晓，尽管他们可能熟悉 np.nan，这是在 NumPy 数组和 pandas DataFrame 中表示 NaN 值的标准方式（见下文）。一个可能的原因是这两个值不像 np.nan 那样广泛使用。

这两个 NaN 对象都是值，都是 float 类型：

float(“nan”) 和 math.nan 的类型都是 float。图像由作者提供。

顺便提一下，现在你不应该感到惊讶，了解到 np.nan 的类型也是 float。

记住两个 NaN 值的比较是很重要的：

>>> float("nan") is float("nan")
False
>>> float("nan") == float("nan")
False
>>> math.nan is math.nan
False
>>> math.nan == math.nan
False

这是因为我们只知道 NaN 不是一个数字，但我们不知道它是什么类型的值。在一种情况下，它可以是一个字符串；在另一种情况下，它可以是一个不同的字符串；还有一种，它可以是一个长字典；再一种，它可以是一个缺失值，因为 NaN 经常用于表示 NA。因此，我们不能假设两个 NaN 值彼此相等。这在处理数值向量和矩阵时可能会产生相当大的差异：

>>> [1, 2, 3] == [1, 2, 3]
True
>>> [1, 2, float("nan")] == [1, 2, float("nan")]
False

然而，如果我们创建一个新的NaN对象，我们将看到：

>>> NaN = float("nan")
>>> NaN is NaN
True
>>> NaN == NaN
False
>>> NaN = math.nan
>>> NaNmath = math.nan
>>> NaNmath is NaNmath
True
>>> NaNmath == NaNmath
False

你有没有注意到，即使is比较返回True，==比较却返回False？所以，这个对象是它自己，但它却不等于自己……

记住使用新定义的哨兵值如NaN或NaNmath时的这种行为。我知道这很诱人，我自己也做过不止一次。因此，只有当这种行为是你想要实现的时，才这样做。

让我们回到这个例子：

>>> x = [1, 2, "three"]
>>> sum(xi for xi in x if isinstance(xi, (float, int)))
3

让我们看看NaN值的实际应用。我们可以用NaN值替换掉"three"，而不是调整sum()函数。为此，我们可以使用以下函数：³

from collections.abc import Sequence
from typing import Any

def use_nan(__x: Sequence[float | Any]) -> Sequence[float]:
    """Replace non-numerical values with float("nan").

    >>> NaN = float("nan")
    >>> use_nan([1, 2, 3])
    [1, 2, 3]
    >>> use_nan([1., 2., 3.])
    [1.0, 2.0, 3.0]
    >>> use_nan([1, 2., 3.])
    [1, 2.0, 3.0]
    >>> use_nan([1, 2, "str"])
    [1, 2, nan]
    >>> use_nan((1, 2, str))
    (1, 2, nan)
    >>> use_nan(1., 2, Any, str, (1, 2,)))
    (1.0, 2, nan, nan, nan)
    """
    return type(__x)([xi
                      if isinstance(xi, (float, int))
                      else float("nan")
                      for xi in __x])

现在，让我们在使用sum()函数之前使用该函数，如上所述，它不接受非数值型的值：

>>> x = use_nan([1, 2, float("nan")])
>>> sum(x)
nan
>>> import math
>>> sum([1, 2, math.nan])
nan

哈？发生了什么？我们使用了NaN值来使sum()工作，确实它不像之前那样抛出错误。但它只是返回了nan……

从数学角度来看，这完全合理：将一个数字加到不是数字的东西上不会得到一个数字，对吧？这就是我们上面得到nan的原因。但这真的是我们想要实现的结果吗？

这要视情况而定。通常，我们可以选择如何处理NaN值。最典型的方法是删除它们。这通过从数据框中删除整行或整列，或从变量中删除单元格来完成。另一种方法——在统计学中常用——是用其他值填充缺失值，这称为插补。

本文并不打算详细讨论这些方法。你可以在许多统计学书籍中阅读这些方法，也可以在各种文章中找到它们；下面的两篇文章描述了在 Python 中使用这些方法：

[## 3 种处理 Python 中缺失值的终极方法]

在插补和数据删除之间明智地选择

towardsdatascience.com](/3-ultimate-ways-to-deal-with-missing-values-in-python-ac5a17c53787?source=post_page-----798d9ed946c0--------------------------------) [## 如何处理 Python 中的缺失数据？ [5 个简单步骤说明]]

学习如何处理 Python 中的缺失数据。识别、评估并处理缺失数据，以便你能最大限度地发挥……

www.analyticsvidhya.com](https://www.analyticsvidhya.com/blog/2021/05/dealing-with-missing-values-in-python-a-complete-guide/?source=post_page-----798d9ed946c0--------------------------------)

如上所述，标准库中的方法确实能处理 NaN 值，但它们会简单地返回 nan，这是 float("nan") 的 repr 和 str 表示。因此，我们需要手动从容器中移除不是数字的值。不幸的是，鉴于 NaN 值的比较方式，以下方法将无效

>>> x = use_nan([1, 2, "three"])
>>> sum(xi for xi in x if xi is not float("nan"))
nan
>>> sum(xi for xi in x if xi != float("nan"))
nan

那么，nan 又来了。这是怎么回事？

我们已经知道发生了什么：NaN 在与其他 NaN 值进行比较时返回 False，这发生在 if xi is not NaN 和 if xi != NaN 中。因此，我们需要一个专门的函数来检查 NaN 值。标准库提供了这样的一个函数，在 math 模块中：

>>> sum(xi for xi in x if not math.isnan(xi))
3

结论

我们讨论了在 Python 标准库中使用 NaN 值。这些知识应该足够你在标准库工具中处理 NaN 值。然而，数值框架可以实现自己的 NaN 值。例如 NumPy 的 np.nan。

正如之前提到的，我不认为作为 Python 程序员，你会频繁使用标准库中的这两个 NaN 哨兵。然而，你应该知道如何处理它们，因为有时你可能需要使用它们，即使你在使用像 NumPy 这样的数值框架。此外，我认为仅为了使用 np.nan 而安装 NumPy 并不明智。我希望这篇文章能帮助你处理这样的情况。

脚注

¹ 值得注意的是，像 None 一样，NaN 值也是哨兵值：

[## 哨兵值 — 维基百科

在计算机编程中，哨兵值（也称为标志值、三值、恶意值、信号值或…

en.wikipedia.org

² 一个例子是 np.nan 和处理 NaN 值数据的 Numpy 函数，如 np.nansum()、np.nanmean()、np.nanmax()、np.nanmin() 和 np.nanstd()。

³ 函数的文档字符串包含多个 doctest；你可以在以下文章中阅读关于这个出色的文档测试工具的信息，它也可以用于单元测试：

## 使用 doctest 进行 Python 文档测试：简便方法

doctest 允许进行文档、单元和集成测试，以及测试驱动开发。

[towardsdatascience.com

感谢阅读。如果你喜欢这篇文章，你可能也会喜欢我写的其他文章；你可以在 这里 查看。如果你想加入 Medium，请使用我下面的推荐链接：

[## 使用我的推荐链接加入 Medium - Marcin Kozak

作为 Medium 会员，你的一部分会员费会转给你阅读的作者，并且你可以完全访问所有故事……

medium.com

学习 Transformers 代码优先：第一部分 — 设置

原文：towardsdatascience.com/nanogpt-learning-transformers-code-first-part-1-f2044cf5bca0

使用 nanoGPT 作为起点的 4 部分 Transformers 探索

Lily Hughes-Robinson

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

–

图片由 Josh Riemer 提供，来源于 Unsplash

不知道你是否和我一样，有时查看代码比阅读论文更简单。当我在开发 AdventureGPT 时，我首先阅读了 BabyAGI 的源代码，这是一种用大约 600 行 Python 实现的 ReAct 论文。

最近，我通过优秀的 Cognitive Revolution Podcast 的 第 33 集 了解到了一篇名为 TinyStories 的最新论文。TinyStories 尝试展示经过数百万（而不是数十亿）参数训练的模型，在足够高质量的数据下也能有效。在论文中的微软研究人员的案例中，他们使用了从 GPT-3.5 和 GPT-4 生成的合成数据，这些数据生成的零售成本大约为 $10k。数据集和模型可以从作者的 HuggingFace repo 获取。

听到一个模型可以在 30M 及更少参数上进行训练，我感到很着迷。作为参考，我在搭载 GTX 1660 Ti 的 Lenovo Legion 5 笔记本上进行所有模型训练和推理。即使仅仅是推理，大多数拥有超过 3B 参数的模型也太大，无法在我的机器上运行。我知道有付费的云计算资源，但我是在空闲时间学习这些内容，实际上只能负担得起通过 API 调用产生的适度 OpenAI 账单。因此，能够在我的普通硬件上训练模型的想法让我立刻振奋起来。

我开始阅读 TinyStories 论文，很快意识到他们在模型训练中使用了已经停用的GPT Neo模型。我开始深入研究代码，看看是否能理解它，并意识到我需要从更小的东西开始。为了提供背景，我主要是一名后端软件工程师，具有足够的机器学习经验以便在听到别人谈论神经网络时不会完全迷失。我离一个真正的 ML 工程师还很远，这使我在搜索引擎中输入了“gpt from scratch”以寻找更温和的入门介绍。我找到下面的视频，一切都发生了变化。

这是我一直在寻找的。除了视频中链接的基本代码库，还有一个名为nanoGPT的精简版本，目前仍在积极开发中。更重要的是，训练代码和模型代码每个大约有 300 行 python。对我来说，这比视频更令人兴奋。我关闭了视频，开始仔细研究源代码。nanoGPT 利用了我从未使用过的 PyTorch。它还涉及了足够的数学和机器学习术语，让我这位新手感到紧张。这将是一个比我预期更大的工程。

理解某件事的最佳方法之一是写下来。因此，我计划深入研究 nanoGPT 代码库中的代码，阅读著名的“Attention is All You Need”论文，并以自下而上的实践方式学习变换器。无论我在这个过程中学到什么，我都希望在这个系列中写下来。如果你想跟随学习，可以将 nanoGPT 代码库克隆到你的机器上（模型甚至可以在 CPU 上训练，所以没有硬件借口），并进行跟随。

克隆代码库后，我做的第一件事是按照 README 的指示训练最简单的模型，即使用tiny_shakespeare 数据集的字符级生成模型。这里有一个脚本用于准备训练数据，一个脚本用于实际训练，还有一个采样脚本用于输出生成的文本。通过几个终端命令和一个多小时的训练，我得到一个能够输出类似莎士比亚风格文本的简单模型。

遵循说明是好的，但我直到将其修改为适合自己的用例后才真正理解某些东西。我的目标是使用 TinyStories 数据集训练一个类似的字符级模型。这需要创建我自己的数据准备脚本，以使数据集准备好进行训练。让我们深入探讨一下这个问题。

nanoGPT 有两种类型的数据准备脚本：一种用于 GPT-2 风格模型，一种用于字符级模型。我从 HuggingFace 仓库下载的 GPT-2 模型中提取了一些代码，其余的都来自 tiny_shakespeare 字符级脚本。这里一个重要的点是，tiny_shakespeare 仅有 1MB 多一点，仅包含 40k 行莎士比亚的作品。而 TinyStories 压缩后超过 3GB，包含 39.7M 个故事。对 tiny_shakespeare 进行分词和切片的方法不能直接转移，至少在我的笔记本电脑上（拥有 32GB RAM）是这样。我在尝试 pythonic、易读的数据准备方法时多次崩溃了我的机器。最终脚本使用了一些技巧，我将在下面详细说明。

首先，我处理数据列表的首选方案是 列表推导式，这是一种从现有列表生成新列表并进行修改的语法。在这种情况下，列表推导式的问题是，3GB 的压缩文本在 RAM 中变得接近 10GB。现在，列表推导式需要在 RAM 中多次复制列表。对于小数据不是问题，但对 TinyStories 来说不可行。

数据准备脚本的输出是一个压缩的 NumPy 数组，包含训练和验证数据的字符级编码，以及一个包含唯一字符完整列表和编码/解码映射的元数据 pickle，以将这些字符转换为数字。以此为参考，一旦找到并映射到数字，我们不需要其他任何东西，除了最终的编码数字数组。最有效的内存使用方法是通过简单的 for 循环迭代数据，同时分段构建这些输出。为此，在循环前初始化一个变量，然后在每次交互中更新。这样可以防止在 RAM 中保存数据集的多个版本，只输出我们需要的内容。最终的词汇生成代码如下：

chars_dataset = set([])
len_dataset = 0

# get all the unique characters that occur in this text as well as total length for training data
desc = "Enumerate characters in training set"
for story in tqdm(dataset['train']['text'], desc):
    chars = list(set(story))

    for char in chars:
        chars_dataset.add(char)

    len_dataset += len(story)

也就是说，将 30.7M 个故事（超过 40 亿字符）编码为数字的数组仍然占用大量 RAM，因为 Python 是动态存储整数的。这里引入 NumPy，它具有更高效的数组存储，可以指定整数的确切大小。除了高效的存储，NumPy 还有内存高效的数组连接，可以用于逐步构建最终的编码数组，而不是一次性完成。

我对脚本的最后修饰是使用 tqdm 为每一步添加一个进度条，最后我准备好运行脚本。所以，我让它过夜运行，早上回来时，脚本仍在运行，估计剩余计算时间超过 100 小时。

这时我真正意识到：30.7M 的故事对于语言模型来说虽然小巧，但绝对不是一个可以在单线程上处理的玩具数据集。是时候引入大招：并行化了。并行化带来了许多复杂性和开销，但性能提升是值得的。幸运的是，有许多方法可以并行化 Python 代码。许多解决方案需要对串行执行的脚本进行重大重写或复杂的抽象。经过一番挖掘，我找到了一种方法，让我可以保持大部分脚本不变，但仍然运行多个进程以利用所有线程。

Ray 是一个可以轻松地将 Python 方法并行化的库，可以在本地或集群中运行。它处理任务队列中的任务并启动工作进程来处理队列。如果这引起了你的兴趣，下面有一个很好的 ray 指南。

## 现代并行和分布式 Python：Ray 的快速教程

Ray 是一个用于并行和分布式 Python 的开源项目。

towardsdatascience.com

在选择并行化内容时，encode 函数似乎是一个不错的候选者。它有明确的输入和输出，对这些输入没有副作用，而且是计算时间中最大的一部分之一。将现有代码适配 ray 变得非常简单：通过装饰器使函数对 ray 可用，功能调用稍微更改以添加远程属性，并有一个函数来启动所有数据的执行。以下是最初在我的代码库中呈现的样子：

import ray

ray.init()

…

# given all the unique characters within a dataset, 
# create a unique mapping of characters to ints
stoi = { ch:i for i,ch in enumerate(chars_dataset) }

@ray.remote
def encode(s):
    return [stoi[c] for c in s]

…

encoded_stories = []
for story in dataset[‘train’][‘text’]:
    encoded_stories.append(encode.remote(story))

ray.get(encoded_stories)

…

拥有了所有 CPU 的力量，我继续前进，却立即导致了我的笔记本电脑崩溃。由于 ray 使用的本地分布式调用栈，整个数据集在内存中存在了几次。仅仅将整个数据集放入队列就导致了内存不足错误。我感到恼火，借此机会买了更多的 RAM（64GB，来了！），但在 RAM 发货期间继续调整代码。

下一个合乎逻辑的步骤是将 ray 处理的请求批量化成适合合理内存的大小。添加批量逻辑相对简单，并且在最终的代码库中有体现，链接将在文章末尾提供。实际上，令人感兴趣的是实验批量大小。一开始，我选择了一个随机的批量大小（5000），效果不错，但我很快意识到在每个批次中，单线程代码占用了相当多的时间。

实际上，在观察我首选的系统监视器时，我看到一个核心被占用数分钟，最后我所有笔记本的核心都亮起了几秒钟，然后又回到只有一个核心被使用。这使我尝试调整批量大小，希望能更快地喂养那些饥饿的 CPU 核心，并让它们保持更长时间的活动。减少批量大小没有帮助，因为每批中的同步代码用于从完整数据集中切分和准备一个批次。这些代码无法并行化，因此每个批次生成数据块的启动成本时间较大。这使我尝试了相反的方法，即增加数据块大小，以使核心保持更长时间的活跃。这有效，因为数据块生成的时间不论数据块大小如何都相同，但每个数据块处理了更多的数据。结合将编码后处理移到 ray 函数中，我能够在短短几个小时内处理 30%的训练数据集，全部在一台笔记本电脑上完成。

最终，在再过几个小时后，我拥有了一个完全准备好的自定义数据集，以供字符级模型使用。我很高兴我不需要求助于昂贵的云计算来处理训练集，如果 RAM 增加不起作用，这将是我的下一个步骤。更重要的是，我深入了解了为字符级模型创建/处理数据集的含义。

在本系列的下一篇文章中，我将检查实际的模型代码，尽我所能地进行解释，并链接大量外部资源以提供额外的信息，以弥补我的知识不足。一旦文章写好，我将回来在这里提供一个链接。与此同时，我已经在下面链接了我数据集准备脚本的最终版本，您可以跟随并了解在有限计算平台上处理较大数据集所需的内容。

[## nanoGPT/data/tinystories_char/prepare.py at master · oaguy1/nanoGPT

最简单、最快速的中型 GPT 训练/微调的代码库。 - nanoGPT/data/tinystories_char/prepare.py…

github.com](https://github.com/oaguy1/nanoGPT/blob/master/data/tinystories_char/prepare.py?source=post_page-----f2044cf5bca0--------------------------------)

此外，系列的第二部分已经上线！点击下面阅读！

## 学习变换器代码第一部分 — GPT 近距离了解

深入研究通过 nanoGPT 生成预训练变换器

[towardsdatascience.com