欧洲的 M&M 巧克力真的比美国的 M&M 巧克力更好吃吗？
一次过度热情地应用科学和数据可视化来解答我们所有人都曾问过的问题
https://erinhwilson.medium.com/?source=post_page---byline--248e747a4c44--------------------------------https://towardsdatascience.com/?source=post_page---byline--248e747a4c44-------------------------------- Erin Wilson
·发表于Towards Data Science ·16 分钟阅读·2024 年 1 月 23 日
–
一张特别甜美的箱型图。图像来源：作者。
(哦，只有我在问这个问题吗…？嗯。如果你有一分钟时间，请享受这篇探索性的数据分析——涵盖实验设计、统计学和互动可视化——虽然有点过于认真，但仍然应用于解决一场国际争论。)
1. 引言
1.1 背景与动机
巧克力在世界各地都受到喜爱。从古代在亚马逊盆地采摘有机可可的传统，到瑞士山脉中的巧克力工匠雕刻食用艺术，再到美国宾夕法尼亚州赫尔希市的巨大工厂每天生产 7000 万个巧克力“亲吻”，巧克力的多样化形式和风味已经融入了许多文化和习俗中。尽管巧克力产品的质量差异很大，但一种广为人知、常见且易于分享的巧克力形式便是 M&M 巧克力。它们常出现在便利店的收银台和酒店的自动售货机中，五颜六色的小圆粒是大家喜爱的零食，其包装设计几乎能适应任何美国商业节日的需求。
2022 年，我在丹麦生活时听到一个令人担忧的说法：欧洲生产的 M&M 口味与美国生产的 M&M 不同，甚至可以说是“更好”。虽然我承认欧洲的精美巧克力确实非常美味，并且通常优于美国巧克力，但我不确定这个说法是否同样适用于 M&M。我得知许多欧洲人认为美国巧克力有一种“不愉快”或“酸味”的味道，这主要归因于丁酸，这种化合物是由于牛奶在加入牛奶巧克力之前的处理方式不同所产生的。
但说实话，这种差异对 M&M 来说能有多大影响呢？M&M!？我想象中，不论在哪里生产，M&M 都会保留相对加工/大规模生产/廉价糖果的味道。作为唯一一位访问一个由国际科学家组成的多元化实验室，进行前沿生物可持续性研究的美国人，我决定拿出我的数据科学工具箱，调查这个 M&M 口味现象。
1.2 以往的研究
引用一位欧洲女性的话，她在纽约旅行时品尝了一颗美国的 M&M 后说道：
“它们尝起来真恶心，像呕吐一样。我不明白人们怎么能吃这个。我把袋子里的剩下部分都扔掉了。”
呕吐？真的吗？根据我的经验，在美国长大的孩子们吃 M&M 时完全没有任何顾虑。小时候，我家里总会在高流量区域放置装满 M&M 的碗，以便随时提供糖分。显然，美国的 M&M 是可以吃的。那么，它们真的与欧洲版的 M&M 有显著不同，或者说更差吗？
为了回应那位匿名欧洲女性的尖锐评论，我和另外两位在丹麦的美国游客一起品尝了在 Lyngby Storcenter Føtex 购买的 M&M。我们希望能体验到隐藏在我们童年中的 M&M 口味的巨大改善。但奇怪的是，我们没有察觉到明显的口味改善。
不幸的是，之前的两项研究都未能进行严格控制和随机抽样的 M&M 口味对比试验。因此，我们转向了科学。
1.3 研究目标
本研究旨在弥补以往研究的不足，并探讨以下问题：
-
是否存在全球共识认为欧洲的 M&M 确实比美国的 M&M 更好？
-
欧洲人是否真的能察觉在不知道自己吃的是哪种 M&M 的情况下，来自美国和欧洲的 M&M 之间的区别？还是这是欧洲人之间的一个大规模协调的谎言，目的是让美国人感到尴尬？
-
美国人真的对美国和欧洲 M&M 的口味感到“盲目”吗？或者他们能品尝出区别，但只是不把这种区别描述为“口味的改善”？
-
这些所谓的口味差异是否能被其他大陆的公民察觉？如果能，他们是否认为某种口味明显更好？
2. 方法
2.1 实验设计与数据收集
参与者通过诱导——呃，邀请他们参加一个社交聚会（并承诺提供免费食物），这个聚会恰好与测试地点位于同一地点。参与者同意暂停社交并加入研究后，会被安排到测试站，由一位受过训练的实验员引导他们完成以下步骤：
-
参与者坐在桌前，面前有两个杯子：一个空的，另一个装满了水。每只手一个杯子，参与者被要求闭上眼睛，并在实验期间保持闭眼状态。
-
实验员随机用勺子取出一颗 M&M，将其送到参与者的空杯中，参与者被要求吃下这颗 M&M（眼睛仍然闭着）。
-
在每次吃下 M&M 后，实验员通过询问参与者是否觉得 M&M 的味道是：特别好，特别差，还是正常，来收集味觉反应。
-
每位参与者总共接受了 10 颗 M&M（5 颗欧洲的，5 颗美国的），每次一颗，顺序由 random.org 随机确定。
-
在吃每颗 M&M 之间，参与者被要求喝一口水以帮助“清洁味蕾”。
-
收集的数据：对于每位参与者，实验员记录了参与者的原始洲（如果不明确，参与者被要求列出他们童年时最有记忆的吃糖果的洲）。对于每颗 M&M，实验员记录了M&M 的来源（“丹麦”或“美国”）、M&M 的颜色和参与者的味觉反应。实验员还被鼓励记录参与者在测试过程中说出的有趣句子，这些记录被归类为备注（数据可通过此处获取）。
2.2 材料采购与参与者招募
本研究购买了两袋 M&M 巧克力。美国来源的 M&M（“美国 M&M”）是在旧金山国际机场购买的，并由作者的父母带到丹麦探访她时带来。欧洲来源的 M&M（“丹麦 M&M”）则是在位于哥本哈根北部的 Lyngby 市的一家 Føtex 超市购买的。
实验在两个主要时间点进行。前 14 位参与者是在 2022 年 8 月在丹麦的 Lyngby 市进行测试的。测试对象主要是作者在丹麦技术大学（DTU）诺和诺德基金会生物可持续性中心遇到的朋友和室友，他们参加了一个“告别派对”，而实验程序被巧妙地融入其中。在旅行期间，一些其他的朋友和家人也在丹麦进行测试（例如，在火车上）。
剩下的 37 名参与者于 2022 年 10 月在美国华盛顿州西雅图进行测试，主要是在华盛顿大学计算机科学博士项目的研究生举办的“TGIF 欢乐时光”期间。这批参与者大多是来自保罗·G·艾伦计算机科学与工程学院（UW CSE）的学生和工作人员，他们响应了每周五的召集，来到艾伦中心的中庭领取免费的零食和饮料。
图 1. 参与者分布图。在首次抽样事件中，Lyngby 的参与者主要来自北美和欧洲，少数来自亚洲、南美或澳大利亚。我们在西雅图的第二次抽样事件大大增加了参与者数量，主要来自北美和亚洲，欧洲的参与者也有所增加。这两个事件都没有招募到来自非洲的参与者。图表由 Altair 制作。
尽管本研究旨在分析全球趋势，但不幸的是，数据仅从 51 名参与者中收集，作者能够将他们吸引到研究地点，而这些数据并不均衡，也不能代表地球六大有人居住的大陆（图 1）。我们希望在未来的工作中改进招募策略。目前，我们使用此数据集的分析能力仅限于来自北美、欧洲和亚洲的个体的反应趋势，这些数据高度偏向于作者在 2022 年末接触到的子社区。
2.3 风险
尽管我们未获得正式的人体实验批准，但此实验仍存在一些小的风险：参与者被告知，由于参与本研究，他们可能会摄入较高的糖分，并可能会经历“令人不悦的味道”。没有预见到其他风险。
然而，在实验后，我们不幸地观察到，当参与者得知他们的味觉反应偏向他们未曾预期的 M&M 类型时，常常会感到自豪感受挫。尤其是欧洲参与者，当他们得知自己或未婚夫/妻的偏好倾向于美国 M&M 时，这种自豪感的受挫似乎最为严重，尽管这并未通过定量测量，且只能通过轶事证据确认。
3. 结果与讨论
3.1 “美国 M&M”与“丹麦 M&M”总体反应
3.1.1 分类反应分析 — 整体数据集
在我们的第一次分析中，我们统计了“差”、“普通”和“好”口味反应的总数，并报告了每种反应在每种 M&M 类型中所占的百分比。来自丹麦的 M&M 比美国的 M&M 更频繁地获得“好”反应，但也更频繁地获得“差”反应。来自美国的 M&M 最常被报告为“普通”味道（图 2）。这可能是因为参与者中来自北美的人数较多，而美国 M&M 是默认的，因此更多被认为是“普通”，而丹麦 M&M 则更常被认为比基准更好或更差。
图 2. 整个数据集的定性味觉反应分布。计算每种 M&M 类型的“差”、“普通”或“好”反应的百分比。图表使用 Altair 制作。
现在让我们来分析一些统计数据，比如进行卡方（X2）检验，比较我们观察到的分类味觉反应分布。使用 scipy.stats chi2_contingency函数，我们构建了每种 M&M 类型的“好”、“普通”和“差”反应的观察计数列联表。利用 X2 检验来评估原假设，即两种 M&M 之间没有差异，我们得到了检验统计量的p-值为 0.0185，这在常见的p-值截断标准 0.05 下是显著的，但在 0.01 下则不显著。所以这只是一个坚实的“也许”，具体取决于你是否希望这个结果显著。
3.1.2 定量反应分析 — 整个数据集
X2 检验有助于评估分类反应是否存在差异，但接下来，我们希望确定两种 M&M 类型之间的相对味觉排名。为此，我们将味觉反应转换为定量分布，并计算了味觉评分。简而言之，“差”=1，“普通”=2，“好”=3。对于每个参与者，我们对他们品尝的每种类型的 5 颗 M&M 的味觉评分进行了平均，保持每种 M&M 类型的独立味觉评分。
图 3. 整个数据集的定量味觉评分分布。计算每个参与者对每种 M&M 类型的平均味觉评分的核密度估计。图表使用 Seaborn 制作。
拿到每种 M&M 类型的平均味觉评分后，我们使用 scipy.stats ttest_ind（“T 检验”）来评估美国和丹麦 M&M 的味觉评分均值是否不同（原假设为均值相同）。如果均值有显著差异，那么这将证明某种 M&M 被认为比另一种更美味。
我们发现美国 M&M 和丹麦 M&M 的平均味觉评分非常接近（见图 3），且没有显著差异（T 检验：p = 0.721）。因此，在所有参与者中，我们没有观察到两种 M&M 类型在感知味觉上的差异（或者，如果你喜欢解析三重否定：“我们不能拒绝原假设，即两者之间没有差异”）。
但如果我们按照参与者的出生大陆来区分，这个结果会有所不同吗？
3.2 “美国 M&M”和“丹麦 M&M”在各大洲的反应
我们在将参与者按其大陆来源分组后重复了上述 X2 和 T 检验分析。由于澳大利亚和南美组的样本量较小，我们将其合并为一个组，以尽量保护数据隐私。由于即使是合并后的澳大利亚/南美组（n=3）样本量仍然较小，我们将避免对该组趋势进行分析，但为了完整性和参与者的享受，我们将该组的数据包含在几个图表中。
3.2.1 类别反应分析 — 按大陆
在图 4 中，我们展示了每个大陆组的口味反应计数（上面面板，注意交互式图例）和反应百分比（下面面板）。北美和亚洲的趋势与整个数据集类似：参与者报告丹麦 M&M 时，“好”的比例高于美国 M&M，但也报告丹麦 M&M 时，“不好”的比例更高。美国 M&M 最常被报告为“正常”（图 4）。
相反，欧洲参与者报告美国 M&M 时，“不好”的比例接近 50%，而“好”的比例仅为 18%，这是最负面和最不积极的反应模式（排除了样本量较少的澳大利亚/南美组）。
图 4. 按大陆的定性口味反应分布。上面面板：口味反应的计数 — 点击图例进行交互式筛选！下面面板：每种 M&M 类型的口味反应百分比。图形由 Altair 制作。
这在条形图中表现得很突出，然而，只有北美在评估每个大陆的两种 M&M 类型之间的口味反应差异时，具有显著的 X2 p-值（p = 0.0058）。欧洲的 p-值或许在某些圈子中被认为是“接近显著”，但我们即将进行更多的假设检验，并且应当注意多重假设检验（表 1）。在这里出现假阳性结果将是灾难性的。
在比较两大洲之间相同 M&M 类型的口味反应模式时，有几点有趣的发现。首先，我们发现，在评估丹麦 M&M 时，不同大陆之间的口味差异并不显著——全球对来自欧洲的 M&M 的感觉普遍一致（右栏 X2 p-值，表 2）。为了更容易地可视化这一比较，我们在图 4 中重新组织了条形图，将其按 M&M 类型分组（图 5）。
图 5. 按 M&M 类型的定性口味反应分布，结果以百分比形式报告。（与图 4 相同的数据，但重新排列）。图形由 Altair 制作。
然而，当比较各大洲对美国 M&M 的反应时，我们发现了更大的差异。我们发现一个配对的差异是显著的：欧洲和北美的参与者对美国 M&M 的评价截然不同（p = 0.000007）（表 2）。看起来这种差异不太可能是随机产生的（表 2 左列）。
3.2.2 定量反应分析 — 按大洲划分
我们再次将分类特征转化为定量分布，以评估各大洲对 M&M 类型的相对偏好。对于北美，我们看到两种 M&M 类型的口味评分均值其实相当接近，但美国 M&M 的“正常”评分周围密度较高（图 6A）。欧洲的分布在均值上保持了一定的分离（尽管这种差异不算显著），美国 M&M 的评分较低（图 6B）。亚洲参与者的口味评分分布最为相似（图 6C）。
在重新调整以比较各大洲对相同 M&M 类型的口味评分的定量方法时，只有北美和欧洲参与者对美国 M&M 的评分在 T 检验基础上有显著差异（p = 0.001）（图 6D），尽管我们现在真的面临多重假设检验的风险！如果你打算认真对待这个分析，请保持谨慎。
图 6. 按大洲分的定量口味评分分布。每种 M&M 类型的平均口味评分的核密度估计。A. 北美对每种 M&M 的反应比较。B. 欧洲对每种 M&M 的反应比较。C. 亚洲对每种 M&M 的反应比较。D. 各大洲对美国 M&M 的比较。E. 各大洲对丹麦 M&M 的比较。图由 Seaborn 绘制。
此时，我开始觉得或许欧洲人并不是在胡说八道。我并不是说他们所说的那么夸张，但也许确实存在某种差异……在某种程度上，北美参与者也感知到了这种差异，但他们对欧洲来源的 M&M 的评价并不总是正面或负面一致。
3.3 M&M 口味对齐图
在我们迄今为止的分析中，我们没有考虑到参与者在 M&M 欣赏度上的基线差异。例如，假设人物 1 将所有丹麦 M&M 评分为“好”，所有美国 M&M 评分为“正常”，而人物 2 则将所有丹麦 M&M 评分为“正常”，所有美国 M&M 评分为“差”。他们对丹麦 M&M 和美国 M&M 的相对偏好是相同的，但人物 2 可能并不像人物 1 那样喜欢 M&M，且原始评分的平均化使得相对偏好信号变得模糊。
受桌面角色扮演游戏（如《龙与地下城©™》）中使用的合法/混乱与善良/邪恶对齐图的启发，在图 7 中，我们建立了一个 M&M 对齐图，帮助确定参与者在 M&M 享受类别中的分布。
图 7. M&M 享受对齐图。x 轴表示参与者对美国 M&M 巧克力的平均口味得分；y 轴表示参与者对丹麦 M&M 巧克力的平均口味得分。图表由 Altair 制作。
值得注意的是，右上象限，即两种 M&M 巧克力均被视为“好”到“正常”的区域，主要由北美参与者和少数亚洲参与者占据。所有欧洲参与者都位于图表的左半部分，其中美国 M&M 巧克力被视为“正常”到“差”，但欧洲人对于丹麦 M&M 巧克力的感知有些分裂，分布在上下半部分，其中对丹麦 M&M 巧克力的感知从“好”到“差”不等。
图 7 的交互版本已提供，供读者探索不同 M&M 对齐区域的参与者计数。
图 7（交互版）：点击并用鼠标刷过散点图，查看不同 M&M 享受区域中各大洲的计数。图表由 Altair 制作。
3.4 参与者口味反应比率
接下来，为了去除 M&M 享受的基线影响，并专注于参与者之间对两种 M&M 巧克力类型的相对偏好，我们取了每个参与者的美国 M&M 口味得分平均值与其丹麦 M&M 口味得分平均值的对数比值。
方程式 1：计算每个参与者总体 M&M 巧克力偏好比率的方程。
因此，正分数表示偏好美国 M&M 巧克力，而负分数表示偏好丹麦 M&M 巧克力。
平均而言，欧洲参与者对丹麦 M&M 巧克力的偏好最强，亚洲人也表现出对丹麦 M&M 巧克力的轻微偏好（图 8）。对于那两位在得知自己对美国 M&M 巧克力有轻微偏好时感到自豪感下降的欧洲人，别担心：你们并不是认为美国 M&M 巧克力是“好”的，而是将其评为比丹麦 M&M 巧克力稍微“不那么差”（参见交互版图 7 中的参与者 ID 4 和 17）。如果你坚持认为 M&M 是一个不值得复制的糟糕美国发明，转而回归食用手工制作的欧洲巧克力，你的荣誉可能会恢复。
图 8. 按大陆分布的参与者 M&M 偏好比率。偏好比率按方程 1 计算。正数表示相对偏好美国 M&M 巧克力，负数表示相对偏好丹麦 M&M 巧克力。图表由 Seaborn 制作。
北美参与者在偏好比率上差异较大：有些人偏好中立，接近 0，有些人则强烈偏好熟悉的美国 M&M，而少数人则中等程度地偏好丹麦 M&M。据 anecdotal 经验，北美人中那些偏好欧洲 M&M 的人，似乎表现出某种夸大的自豪感，仿佛他们的结果象征着一种高雅的品味。
总体而言，通过 T 检验比较 M&M 偏好比率的分布，显示出欧洲和北美参与者之间可能存在显著的均值差异（p = 0.049），但拜托，这是我报告的第 20 个 p 值——这个可能太接近了，难以确定。
3.5 味觉不一致性与“完美分类者”
对每位参与者，我们通过计算他们对每种 M&M 类型的反应标准差的平均值来评估他们的味觉评分一致性，并将其与他们的偏好比率进行对比绘图（见图 9）。
图 9. 参与者的味觉一致性与偏好比率。横轴表示参与者相对的 M&M 偏好比率。纵轴表示他们对美国 M&M 和丹麦 M&M 评分的标准差的平均值。纵轴上的 0 值表示反应的一致性完美，而更高的值表示反应不一致。图表由 Altair 制作。
大多数参与者的评分在一定程度上存在不一致，针对同一类型的 M&M 在 5 个样本中的评分不同。如果欧洲产和美国产的 M&M 味觉差异并不十分明显，那么这种不一致是可以预期的。最不一致的是那些对同一类型的 M&M 给出“好”、“正常”和*“坏”*不同评分的参与者（例如，纵轴上的高点，味觉评分的标准差较大），这表明他们的味觉感知能力较差。
有趣的是，四位参与者——每个大洲组的一位——表现得完全一致：他们对每种 M&M 类型的 5 颗 M&M 的味觉反应相同，导致平均标准差为 0.0（见图 9 底部）。排除掉其中一位仅将所有 10 颗 M&M 评为“正常”的参与者，其余三位看起来是“完美分类者”——他们要么将一种类型的所有 M&M 评为“好”，另一种类型评为“正常”，要么将一种类型的所有 M&M 评为“正常”，另一种类型评为“坏”。也许这些人是“超级味觉者”。
3.6 M&M 颜色
另一种可能的解释是个体味觉反应的不一致性，可能是由于 M&M 颜色所带来的可感知的味觉差异。从视觉上看，美国生产的 M&M 明显比丹麦生产的 M&M 更加光滑和鲜艳，而丹麦 M&M 则显得有些“斑驳”（见图 10A）。在实验过程中记录了 M&M 的颜色，尽管平衡取样并未正式纳入实验设计，但颜色似乎大致均匀取样，唯一的例外是蓝色美国 M&M，其样本量过大（见图 10B）。
图 10. M&M 颜色。A. 每种类型的 M&M 颜色照片。也许在我那没有专业灯光的照片中很难看出，但肉眼观察，美国的 M&M 似乎更亮且颜色更均匀，而丹麦的 M&M 则呈现出较暗且带有斑驳色彩。是我一个人感觉这样，还是你们也能听到欧洲人说“它们更亮是因为你们在食物中加了那些我们在这里禁用的额外化学物质！”B. 在实验过程中每种颜色的 M&M 分布。蓝色的美国 M&M 并不是故意被过度抽样的——它们一定特别亮眼，令实验者难以抗拒。图表由 Altair 制作。
我们简要地根据颜色可视化了可能的口味反应差异（见图 11），然而我们认为数据不足以支持明确的结论。毕竟，平均来说，每个参与者可能只会品尝到 6 种 M&M 中的 5 种颜色一次，而 1 种颜色则根本没有品尝到。我们将进一步的 M&M 颜色研究留待未来的工作。
图 11. 各种颜色和类型的 M&M 口味反应概况。反应结果以“差”、“正常”和“好”回应的百分比表示，尽管并不是所有的 M&M 都是平均抽样的。图表由 Altair 制作。
3.7 多彩的评论
我们向每位参与者保证，在本实验中没有“正确的答案”，所有感受都是有效的。尽管有些参与者非常认真，有时会花超过一分钟时间细细品味每颗 M&M 并像品酒师一样评价它，许多参与者似乎将实验视为一场竞争（这有时会导致骄傲的情绪膨胀或萎缩）。实验者在记录 M&M 反应的同时也做了笔记和摘录，其中有些内容有些“生动”。我们提供了一个匆忙制作的每种 M&M 类型的词云，供娱乐之用（见图 12），但我们提醒大家，不要在没有认真情感分析的情况下过度解读它们。
图 11. 从每种 M&M 类型的笔记栏中生成的简单词云。提醒一下——这些词云还没有经过情感分析，并且记录了一些不太恰当的语言。图表由 WordCloud 制作。
4. 结论
总体来看，并没有出现“全球共识”认为欧洲的 M&M 比美国的 M&M 更好。然而，欧洲参与者似乎更强烈地对美国的 M&M 表达负面反应，而北美的参与者似乎对是否更喜欢来自美国还是欧洲的 M&M 观点较为分歧。亚洲参与者的偏好趋势通常介于北美和欧洲之间。
因此，我承认，欧洲人并没有在关于 M&M 的事情上进行大规模的联合欺骗。大多数欧洲参与者偏向丹麦 M&M 这一结果令人信服，尤其是考虑到我是亲自收集大部分味觉反应数据的实验者。如果他们找到了欺骗的方式，做得足够巧妙，超出了我的被动察觉，以至于我没有注意到。然而，基于这项研究，似乎强烈的“呕吐味”并不是普遍被感知的，而且在同时品尝两种 M&M 类型时，不会出现在非欧洲人中。
我们希望这项研究能带来启发！我们期待在未来扩展这项工作，改进参与者抽样，增加来自其他大洲的不同类型 M&M，并深入探讨颜色可能带来的味觉差异。
感谢所有参与并为了科学而吃 M&M 的人！
图表和分析可以在 github 上找到：github.com/erinhwilson/mnm-taste-test
本文由 Erin H. Wilson 博士[1,2,3]撰写，她决定在答辩和开始下一份工作之间的这段时间，最好用来进行这项非常有价值的分析。希望大家能看出这篇文章是幽默的——我并没有对不喜欢美式 M&M 的欧洲人怀有任何负面情绪，只是很享受这个机会，带着点调侃的心态，玩笑式地讽刺我们在充满热情的数据分析中所展开的激烈辩论。
感谢 Matt、Galen、Ameya 和 Gian-Marco 在数据收集方面的帮助！
[1] 曾在华盛顿大学保罗·G·艾伦计算机科学与工程学院担任博士生
[2] 曾在丹麦技术大学诺和诺德基金会生物可持续发展中心担任访问博士生
[3] LanzaTech 的未来数据科学家
机器学习模型是否存储受保护的内容？
概念验证
https://medium.com/@nathanReitinger?source=post_page---byline--abec357c6b70--------------------------------https://towardsdatascience.com/?source=post_page---byline--abec357c6b70-------------------------------- Nathan Reitinger
·发表于Towards Data Science ·5 分钟阅读·2024 年 5 月 6 日
–
从 chatGPT 到 Stable Diffusion，人工智能（AI）正经历一个类似于1970 年代的夏天，AI 的盛况可与那个时代的辉煌相提并论。然而，这一欢庆并非没有遭遇反对。从好莱坞到卢浮宫，人工智能似乎唤醒了一个沉睡的巨人——一个渴望保护曾经看似专属于人类的世界：创造力。
对于那些渴望保护创造力的人来说，人工智能似乎有一个致命弱点：训练数据。事实上，所有最佳模型都需要一个高质量、涵盖全球的数据源——但这意味着什么呢？
首先，高质量意味着人为创造的。尽管非人工创造的数据自从计算机自我对弈的概念被战争游戏推广以来取得了许多进展，计算机科学文献却表明，如果完全去除人的因素（即模型腐化或模型崩塌），模型质量随着时间的推移会下降。简单来说：人类数据是这些模型的命脉。
第二，全球性意味着全球性。如果你把它放到网上，你应该假设模型已经在训练中使用了它：那个你原本希望只有你和 Tom 记得的 Myspace 帖子（已被吸收），那个你高兴地忘记的图片封存记忆，直到 PimEyes 迫使你重新记起它（已被吸收），以及那些你希望只是梦境的深夜 Reddit 争论（已被吸收）。
像 LLaMa、BERT、Stable Diffusion、Claude 和 chatGPT 这样的模型都是在大量由人类创作的数据上进行训练的。而一些、许多或大多数人类创作的表达方式——尤其是那些恰好固定在计算机可以访问并学习的有形介质上的表达——具有版权保护的资格。
Anderson v. Stability AI；Concord Music Group, Inc. v. Anthropic PBC；Doe v. GitHub, Inc.；Getty Images v. Stability AI；{Tremblay, Silverman, Chabon} v. OpenAI；纽约时报诉微软
虽然可能是偶然的，这些模型无法生存的数据正是大多数受到版权保护的数据。这也催生了我们今天看到的巨大的版权斗争。
在这些诉讼中产生的许多问题中，最紧迫的一个问题是模型本身是否存储受保护的内容。这个问题似乎相当明显，因为我们怎么能说模型——仅仅是由数字（即权重）和架构组成的集合——“存储”了什么？正如 Murray 教授所说：
当前关于视觉生成型 AI 系统的辩论中的许多参与者抓住了这样一个观点：生成型 AI 系统已在包含实际版权保护的图像文件（如.jpg、.gif、.png 等）数据集和基础模型上进行训练，这些文件是从互联网上抓取的，数据集或基础模型一定已经制作并存储了这些作品的副本，并且生成型 AI 系统以某种方式进一步选择并复制了这些数据集中的个别图像，并以某种方式将这些图像的重大可版权部分复制并纳入到最终生成的图像中，供最终用户使用。这是一种魔法般的思维。
Michael D. Murray, 26 SMU 科技与法律评论 259, 281 (2023)
然而，模型本身似乎在某些情况下，确实会记住训练数据。
以下示例来自HuggingFace 上的 Gradio Space，该平台允许用户选择一个模型，查看输出，并从该模型的训练数据中检查生成的图像与其训练数据中任何图像的相似度。由于 MNIST 数字易于机器解析、易于人类从相似性角度理解，并且具有易于分类的优点——这使得相似性搜索只考虑相同数字的图像（提高效率），因此使用了 MNIST 数字。
让我们看看它是如何工作的！
以下图像的相似度得分为 0.00039。RMSE 代表均方根误差，是评估两张图像相似度的一种方式。事实上，还有许多其他相似性评估方法，但 RMSE 能很好地判断一张图像是否为副本（即，我们这里并不是在寻找法律定义的相似性）。举个例子，RMSE 值小于 0.006 时，图像已接近“复制”范围，而 RMSE 值小于 0.0009 时，则进入完美复制的领域（肉眼无法分辨）。
🤗 一个生成几乎完全相同训练数据副本的模型（RMSE 为 0.0003）🤗
要使用Gradio 空间，请按照以下三个步骤操作（如果空间处于休眠状态，可以选择构建该空间）：
-
步骤 1：选择要使用的预训练模型类型
-
步骤 2：点击“提交”按钮，模型将为您生成一张图像（28x28 的灰度图像）
-
步骤 3：Gradio 应用程序会在该模型的训练数据中搜索，识别与生成图像最相似的图像（从 60K 个示例中筛选）
如上所示，左侧生成的图像（AI 创作）几乎与右侧的训练数据完全相同，当使用“FASHION-diffusion-oneImage”模型时，结果正是如此。这是有道理的。该模型仅对FASHION 数据集中的一张图像进行了训练。同样的情况也适用于“MNIST-diffusion-oneImage”模型。
尽管如此，即使是训练了更多图像（例如 300 张、3000 张或 60000 张图像）的模型，也能产生非常相似的输出。这个示例来自一个生成对抗网络（GAN），它在完整的 60K 图像数据集（仅限训练）上进行了训练，数据集包括MNIST 手写数字。作为背景，生成对抗网络（GAN）通常生成的图像比扩散模型记忆性差：
RMSE 为 0.008
这是另一个使用扩散模型，并在 60K MNIST 数据集上训练的图像（即，支持稳定扩散的模型类型）：
RMSE 为 0.004
随时可以自己尝试使用Gradio 空间，探索模型，或者如果有问题可以联系我！
总结： 这个小型示例的重点是，机器学习模型没有什么神秘或绝对的版权豁免。机器学习模型确实可以并且会生成与其训练数据相同的图像——换句话说，模型确实会存储受保护的内容，因此可能会遇到版权问题。当然，也有许多反驳的论点（我正在进行的工作！）；这个演示应该仅作为存储的轶事性证据，可能是开发者在这一领域工作的“金丝雀”。
输入到模型中的内容和从模型中得到的结果同样重要，尤其对于某些执行特定任务的模型来说更是如此。我们需要小心并关注我们的“黑箱”，因为这个类比往往并不成立。你无法自己解读模型所持有的权重集合，并不意味着你可以摆脱所有形式的责任或审查。
— @nathanReitinger，敬请关注该领域的进一步工作！
Unless otherwise noted, all images are by the author
不要过度思考“离群值”，改用学生 t 分布
使用 R 和 Brms 的贝叶斯方法
https://medium.com/@damanrique?source=post_page---byline--b6c584b91d5c--------------------------------https://towardsdatascience.com/?source=post_page---byline--b6c584b91d5c-------------------------------- Daniel Manrique-Castano
·发表于Towards Data Science ·15 分钟阅读·2024 年 3 月 30 日
–
对许多研究人员而言，离群值是可以极大改变分析过程的“异类波”，或者“混淆”一些预期效果。我更倾向于使用“极端观测值”这一术语，并将离群值留给那些并非研究群体真正一部分的观测值。例如，在我的研究领域（脑缺血研究），离群值指的是没有缺血的动物（尽管应该有），而极端观测值则是那些小范围或大范围缺血，且与其他动物有显著差异的个体。
传统的（频率学派）统计模型是建立在高斯分布的坚实基础上的。这存在一个显著的局限性：一个固有的假设，认为所有数据点会围绕一个中央均值按照可预测的模式聚集（基于中心极限定理）。在柏拉图的理念世界中，这或许是正确的，但我们这些生物医学领域的科学家深知，考虑到有限的样本（动物数量），我们很难依赖这一假设来进行观察。
高斯分布对极端观测值非常敏感，其使用使得科学家们认为，去除极端观测值是获得“更清晰”或“更干净”结果的最佳方式（无论这是什么意思）。正如我曾在一篇文章中作为审稿人 2 所评论的，“问题不在于那些可能‘隐藏’你效果的极端观测值，而在于你使用了一个（我认为）不适合你目的的统计模型”。
应该注意的是，没有哪种统计模型是“正确”或“合适”的，但我们可以估计，在给定数据的情况下，有些统计模型比其他模型更有可能生成观察到的数据（生成模型）。
幸运的是，没人强迫我们受限于高斯模型的假设，对吧？我们还有其他选择，比如学生 t 分布（1）。我认为它是一个更具适应性的工具，用来在真实世界的生物医学数据的动荡海洋中航行。学生 t 分布提供了一种强大的替代方案，承认我们的数据可能包含极端观测值，而这些值是正常的生物学反应，任何情境下都可以预见到。可能会有一些患者或动物对治疗没有反应或反应过度，而我们的建模方法能够识别这些反应作为数据的一部分，这非常重要。因此，本教程将通过**brms**包在 R 语言中的应用（2）来探讨使用学生 t 分布的建模策略——这是贝叶斯建模的强大助手。
学生 t 分布背后是什么？
学生 t 分布不过是具有更重尾部的高斯分布。换句话说，我们可以说高斯分布是学生 t 分布的一个特例。高斯分布由均值（μ）和标准差（σ）定义。而学生 t 分布则增加了一个额外的参数，即自由度（df），它控制分布的“厚度”。这个参数赋予远离均值的事件更高的概率。这一特性对于小样本量，尤其是生物医学领域尤其有用，因为在这些领域，正态分布的假设是值得怀疑的。注意，当自由度增大时，学生 t 分布趋近于高斯分布。我们可以使用密度图来可视化这一点：
# Load necessary libraries
library(ggplot2)
# Set seed for reproducibility
set.seed(123)
# Define the distributions
x <- seq(-4, 4, length.out = 200)
y_gaussian <- dnorm(x)
y_t3 <- dt(x, df = 3)
y_t10 <- dt(x, df = 10)
y_t30 <- dt(x, df = 30)
# Create a data frame for plotting
df <- data.frame(x, y_gaussian, y_t3, y_t10, y_t30)
# Plot the distributions
ggplot(df, aes(x)) +
geom_line(aes(y = y_gaussian, color = "Gaussian")) +
geom_line(aes(y = y_t3, color = "t, df=3")) +
geom_line(aes(y = y_t10, color = "t, df=10")) +
geom_line(aes(y = y_t30, color = "t, df=30")) +
labs(title = "Comparison of Gaussian and Student t-Distributions",
x = "Value",
y = "Density") +
scale_color_manual(values = c("Gaussian" = "blue", "t, df=3" = "red", "t, df=10" = "green", "t, df=30" = "purple")) +
theme_classic()
图 1：具有不同自由度的高斯分布和学生 t 分布的比较。
请注意在图 1 中，随着自由度的降低，均值周围的山峰变得越来越小，因为概率质量转移到尾部，而尾部则更厚。这个特性使得学生 t 分布对离群值的敏感度降低。关于这一点的更多细节，你可以查看这个博客。
加载所需的包
我们加载所需的库：
library(ggplot2)
library(brms)
library(ggdist)
library(easystats)
library(dplyr)
library(tibble)
library(ghibli)
探索性数据可视化
所以，让我们跳过数据模拟，直接进入正题。我将使用我从进行转棒测试的小鼠中获得的真实数据。
首先，我们将数据集加载到环境中，并设置相应的因素水平。数据集包含动物的 ID、一个分组变量（基因型）、两个不同测试日的指示符（天数），以及同一天的不同实验。对于本文，我们仅建模其中一个实验（实验 3）。其他实验将保留在未来的关于建模变异性的文章中。
如数据处理所示，我们的建模策略将基于基因型和天数作为Trial3分布的分类预测变量。
在生物医学科学中，分类预测变量或分组因素比连续预测变量更为常见。该领域的科学家喜欢将样本分为不同的组或条件，并应用不同的处理方法。
data <- read.csv("Data/Rotarod.csv")
data$Day <- factor(data$Day, levels = c("1", "2"))
data$Genotype <- factor(data$Genotype, levels = c("WT", "KO"))
head(data)
数据框
让我们使用Raincloud 图来初步查看数据，正如吉列尔梅·A·弗朗奇博士在这篇精彩的博客文章中所展示的。
edv <- ggplot(data, aes(x = Day, y = Trial3, fill=Genotype)) +
scale_fill_ghibli_d("SpiritedMedium", direction = -1) +
geom_boxplot(width = 0.1,
outlier.color = "red") +
xlab('Day') +
ylab('Time (s)') +
ggtitle("Rorarod performance") +
theme_classic(base_size=18, base_family="serif")+
theme(text = element_text(size=18),
axis.text.x = element_text(angle=0, hjust=.1, vjust = 0.5, color = "black"),
axis.text.y = element_text(color = "black"),
plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
legend.position="bottom")+
scale_y_continuous(breaks = seq(0, 100, by=20),
limits=c(0,100)) +
# Line below adds dot plots from {ggdist} package
stat_dots(side = "left",
justification = 1.12,
binwidth = 1.9) +
# Line below adds half-violin from {ggdist} package
stat_halfeye(adjust = .5,
width = .6,
justification = -.2,
.width = 0,
point_colour = NA)
edv
图 2：探索性数据可视化。
图 2 与吉列尔梅·A·弗朗奇博士原始图形有所不同，因为我们绘制了两个因素而不是一个。然而，图形的性质是相同的。请注意红点，这些点可以被视为极端观察值，它们会将集中趋势的测量（尤其是均值）拉向一个方向。我们还观察到方差不同，因此对sigma的建模也能提供更好的估计。我们现在的任务是使用brms包来建模输出。
使用 brms 拟合统计模型
在这里，我们使用Day和Genotype作为相互作用的分类预测变量，来拟合Trial 3的分布。让我们首先拟合一个典型的高斯模型，它类似于频率学派框架中的普通最小二乘法（OLS）模型，因为我们使用的是默认的平坦brms 先验。先验超出了本文的讨论范围，但我保证我们将在未来的博客中讨论它们。
一旦我们获得了高斯模型的结果，我们可以将它们与学生 t 模型的大量结果进行比较。然后，我们将sigma添加到方程中，以考虑数据方差的差异。
在高斯分布下拟合一个“典型的”（频率学派）模型
我们的高斯模型是在典型（且常常不正确）的同方差假设下构建的（3）。换句话说，我们假设所有组的方差相同（或非常相似）。作为研究人员，我不记得曾看到过这种假设。
Gaussian_Fit1 <- brm(Trial3 ~ Day * Genotype,
data = data,
family = gaussian(),
# seed for reproducibility purposes
seed = 8807,
control = list(adapt_delta = 0.99),
# this is to save the model in my laptop
file = "Models/20240222_OutliersStudent-t/Gaussian_Fit1.rds",
file_refit = "never")
# Add loo for model comparison
Gaussian_Fit1 <-
add_criterion(Gaussian_Fit1, c("loo", "waic", "bayes_R2"))
模型诊断
在继续之前，进行一些简单的模型诊断是个好主意，以便将实际观察与我们模型的预测进行比较。我们可以通过几种方式做到这一点，但最常见的方法是绘制完整的密度图。我们可以使用brms中的pp_check函数来实现这一点。
set.seed(8807)
pp_check(Gaussian_Fit1, ndraws = 100) +
labs(title = "Gaussian model") +
theme_classic()
图 3：高斯模型的诊断
图 3 表明我们的观察值（深蓝色）与模型预测没有显著差异。下面，我为您提供了额外的代码，以检查其他pp_check的替代方案及其各自的图表。
set.seed(88071)
pp_check(Gaussian_Fit1, group = "Genotype", type = "dens_overlay_grouped", ndraws = 100) +
labs(title = "Density by Genotype") +
theme_classic()
pp_check(Gaussian_Fit1, type = "stat_grouped", group = "Genotype", stat = "var", binwidth = 3) +
coord_cartesian(xlim = c(0, 300)) +
ggtitle("Grouped variance") +
theme_classic()
pp_check(Gaussian_Fit1, type = "stat", stat = "var", binwidth = 3) +
coord_cartesian(xlim = c(0, 600)) +
ggtitle("How well we captured the variace") +
theme_classic()
pp_check(Gaussian_Fit1, type = "stat", stat = "mean", binwidth = 2) +
coord_cartesian(xlim = c(0, 50)) +
ggtitle("How well we captured the mean") +
theme_classic()
检查高斯分布的结果
现在，我们使用bayestestR包中的describe_posterior函数（4）来查看结果：
describe_posterior(Gaussian_Fit1,
centrality = "mean",
dispersion = TRUE,
ci_method = "HDI",
test = "rope",
)
这里我们专注于“截距”，即 1 DPI 时 WT 的值，以及“GenotypeKO”，即同一时间点 KO 动物的估计差异。我们看到 WT 动物在转棒上的时间约为 37 秒，而它们的 KO 同类则少于一秒（0.54）更多。作为该领域的研究人员，我可以说，这个差异是没有意义的，基因型对转棒表现没有影响。即使是“天数”的效应，2.9，在这个模型下对我来说似乎也没有意义。我们可以使用brms中的神奇函数conditional_effects轻松地可视化这些估计。
# We create the graph for convex hull
Gaussian_CondEffects <-
conditional_effects(Gaussian_Fit1)
Gaussian_CondEffects <- plot(Gaussian_CondEffects,
plot = FALSE)[[3]]
Gaussian_CondEffects +
geom_point(data=data, aes(x = Day, y = Trial3, color = Genotype), inherit.aes=FALSE) +
Plot_theme +
theme(legend.position = "bottom", legend.direction = "horizontal")
图 8：高斯模型的条件效应
在图 8 中，我们可以看到交互项的估计值和不确定性。我已经自定义了这个图，并添加了一些 ggplot 元素，您可以在原始的Quarto Notebook中查看。请注意，尽管第一天的离散度比第二天大，但两个时间点的不确定性是相似的。我们将在文章结尾的一个小片段中解决这一点。
现在让我们看看，当我们使用学生-t 分布对相同的数据建模时，我们的理解发生了多大的变化。
拟合我们的假设：使用学生 t 分布的模型
现在是时候在我们的brms模型中使用学生 t 分布了。
Student_Fit <- brm(Trial3 ~ Day * Genotype,
data = data,
family = student,
# seed for reproducibility purposes
seed = 8807,
control = list(adapt_delta = 0.99),
# this is to save the model in my laptop
file = "Models/20240222_OutliersStudent-t/Student_Fit.rds",
file_refit = "never")
# Add loo for model comparison
Student_Fit <-
add_criterion(Student_Fit, c("loo", "waic", "bayes_R2"))
模型诊断
我们像之前一样绘制模型诊断：
图 9：学生 t 分布的模型诊断
图 9 显示了观察值和预测值的平均形状及峰值相匹配。需要注意的是，我们的模型似乎预测出了低于 0 的值。这是一个重要的研究问题，我们暂时跳过。不过，这确实暗示了使用信息性先验或设定下界为 0 的分布族，如 log_normal、hurdle_lognormal 或 zero_inflated_poisson，具体取决于情况。Andrew Heiss (5) 在这方面提供了一个 很好的例子。
检查学生 t 分布的结果
让我们来看看后验分布：
describe_posterior(Student_Fit,
centrality = "mean",
dispersion = TRUE,
ci_method = "HDI",
test = "rope",
)
在这个模型下，我们可以看到我们的估计值有所变化，我会说变化是适度的。我们的截距估计（1 天时的 WT）减少了 7 秒。那为什么会这样呢？因为我们在开始时发现的极端值对数据的集中趋势度量影响较小。因此，这是对第 1 天“典型”WT 动物的更准确度量。我们还观察到天数效应的显著增加，比我们最初的高斯估计多了将近 10 秒。重要的是，我们的 KO 基因型效应似乎更加显著，从我们高斯模型中的 0.52 增加到我们学生 t 模型中的 5.5，约增加了 10 倍。从我的角度来看，鉴于这些数据的背景，两个模型之间的差异是显著的。
让我们使用 conditional_effects 以图形方式查看：
Student_CondEffects <-
conditional_effects(Student_Fit)
Student_CondEffects <- plot(Student_CondEffects,
plot = FALSE)[[3]]
Student_CondEffects +
geom_point(data=data, aes(x = Day, y = Trial3, color = Genotype), inherit.aes=FALSE) +
Plot_theme +
theme(legend.position = "bottom", legend.direction = "horizontal")
图 10：学生 t 模型的条件效应
我们能得到更好的估计吗？对于这个具体的例子，我认为我们可以。从一开始就很容易注意到数据的方差差异，尤其是在我们比较第一天和第二天的图形时。我们通过使用学生 t 分布改善了估计，进一步的改进可以通过开发一个异方差性模型来预测 sigma（残差方差）。
这样，模型并不假设你的残差方差在分组变量间是相等的，而是将其作为可以由预测变量建模的响应。
这是我们留到最后的小点。
使用学生 t 分布预测 sigma
我们使用 brms 中的 bf 函数将 sigma 作为响应变量。在这种情况下，我们将使用相同的预测变量 Day 和 Genotype 来建模这个参数。
Student_Mdl2 <- bf (Trial3 ~ Day * Genotype,
sigma ~ Day * Genotype)
Student_Fit2 <- brm(
formula = Student_Mdl2,
data = data,
family = student,
# seed for reproducibility purposes
seed = 8807,
control = list(adapt_delta = 0.99),
# this is to save the model in my laptop
file = "Models/20240222_OutliersStudent-t/Student_Fit2.rds",
file_refit = "never")
# Add loo for model comparison
Student_Fit2 <-
add_criterion(Student_Fit2, c("loo", "waic", "bayes_R2"))
模型诊断
图 11：带有 sigma 的学生 t 分布模型诊断
图 11 看起来很好，除了 0 以下的不舒服预测值。对于这个情况，我判断这不会强烈偏倚估计及其不确定性。然而，这是我在进行实际研究时会考虑的一个方面。
检查带有预测 sigma 的学生 t 分布的结果
现在，让我们来看看后验分布。
与另外两个拟合模型相比，我们看到更多的参数，因为现在模型中将 sigma 的反应作为一个主效应包含在内。在这种方案下，我们看到截距更接近高斯模型的截距，并且基因型（GenotypeKO）的效应减少了一半。
然而，有一点需要注意。在我们的第一个 Student-t 模型中，截距的不确定性为 24.1–37.4。另一方面，在最后一个模型中，不确定性增加到 24.3–46.1。这意味着当我们考虑不同的方差时，我们对这个（和其他）参数的信心减少了。例如，天数的情况也是如此，从 1.2–18.9 变化为-5.6–18.1。此时，我们对第二天与转棒上花费时间增加之间的关系的信心减弱了。
不用担心，统计建模的目的是提供对测量中不确定性的最佳量化，这正是我们现在所做的。当然，当样本中有极端值并且这些极端值也属于我们的总体时，我们的不确定性会增加。
在这个例子中，我们看到，考虑数据中不同的方差给了我们一个完全不同的结果理解。
最后，我们可以看到，绘制在对数尺度上的 sigma 在不同天数和基因型之间有显著变化：
Student_CondEffects2 <-
conditional_effects(Student_Fit2)
Student_CondEffects2 <- plot(Student_CondEffects2,
plot = FALSE)[[3]]
Student_CondEffects2 +
geom_point(data=data, aes(x = Day, y = Trial3, color = Genotype), inherit.aes=FALSE) +
Plot_theme +
theme(legend.position = "bottom", legend.direction = "horizontal")
Student_CondEffects3 <-
conditional_effects(Student_Fit2, dpar = "sigma")
Student_CondEffects3 <- plot(Student_CondEffects3,
plot = FALSE)[[3]]
Student_CondEffects3 +
Plot_theme +
theme(legend.position = "bottom", legend.direction = "horizontal")
图 12：带有 sigma 的 Student-t 模型的条件效应
图 13：sigma 的条件效应
我们在第二个图中看到的是 sigma，它有效地解释了这个参数在不同天数和基因型之间的方差。我们看到在第一天不确定性较高，特别是对于野生型小鼠，而在第二天这个参数是类似的。
我们可以通过比较三个模型的样本外预测来总结本文。
模型比较
我们使用 WAIC 标准进行模型比较（6），用于估计样本外预测误差。通过同时考虑观察数据的对数似然和有效参数数量，它在模型拟合和复杂性之间提供了平衡。与其他一些标准不同，WAIC 本质上考虑了参数的后验分布，而不是依赖于点估计，因此特别适合贝叶斯分析。
给定一个数据集和一个贝叶斯模型，WAIC 的计算公式为：
WAIC=−2×(LLPD−pWAIC)
其中：LLPD 是对数逐点预测密度，计算为每个观察数据点在后验样本中的对数似然的平均值。WAIC 是有效参数数量，通过对数似然的平均值和后验样本中的平均对数似然值之间的差异来计算。
我们使用 performance 包中的 compare_performance 函数，该包是 easystats 环境的一部分（4，7，8）。
Fit_Comp <-
compare_performance(
Gaussian_Fit1,
Student_Fit,
Student_Fit2,
metrics = "all")
Fit_Comp
输出结果显示，我们预测 sigma 的学生 t 模型在样本外预测中受到的惩罚最小（WAIC = 497）。请注意，该模型中没有 sigma 的估计值，因为它作为响应变量包含在内。此表还显示，学生 t 模型的残差方差（sigma）比高斯模型小，这意味着方差被预测变量更好地解释。我们可以将相同的结果可视化为图表：
Fit_Comp_W <-
loo_compare(
Gaussian_Fit1,
Student_Fit,
Student_Fit2,
criterion = "waic")
# Generate WAIC graph
Fit_Comp_WAIC <-
Fit_Comp_W[, 7:8] %>%
data.frame() %>%
rownames_to_column(var = "model_name") %>%
ggplot(
aes(x = model_name,
y = waic,
ymin = waic - se_waic,
ymax = waic + se_waic)
) +
geom_pointrange(shape = 21) +
scale_x_discrete(
breaks=c("Gaussian_Fit1",
"Student_Fit",
"Student_Fit2"),
labels=c("Gaussian_Fit1",
"Student_Fit",
"Student_Fit2")
) +
coord_flip() +
labs(x = "",
y = "WAIC (score)",
title = "") +
Plot_theme
Fit_Comp_WAIC
图 14：通过 WAIC 进行的模型比较
图 14 显示，我们的最后一个模型在样本外预测中受到的惩罚最小。
你可以在我的 GitHub 网站 找到此帖的更新版本。如果这段经历对你有帮助，或者你有任何建设性的意见，请告诉我。
除非另有说明，所有图像均由作者使用 R 代码生成。
参考文献
1.M. Ahsanullah, B. M. G. Kibria, M. Shakil, 正态分布与学生 t 分布及其应用（Atlantis Press，2014；dx.doi.org/10.2991/978-94-6239-061-4）。
2. P.-C. Bürkner, Brms: 一款用于贝叶斯多层次模型的 R 包，使用 Stan。80 (2017)，doi：10.18637/jss.v080.i01。
3. K. Yang, J. Tu, T. Chen, 同方差性：线性回归中的一个被忽视的关键假设。General Psychiatry。32，e100148 (2019)。
4. D. Makowski, M. S. Ben-Shachar, D. Lüdecke, bayestestR：描述贝叶斯框架中效应及其不确定性、存在性和显著性。4，1541 (2019)。
5. A. Heiss, 使用贝叶斯 Beta 回归和零膨胀 Beta 回归模型建模比例的指南 (2021)，(可在 dx.doi.org/10.59350/7p1a4-0tw75 获得)。
6. A. Gelman, J. Hwang, A. Vehtari, 理解贝叶斯模型的预测信息标准。Statistics and Computing。24，997–1016 (2013)。
7. D. Lüdecke, M. S. Ben-Shachar, I. Patil, P. Waggoner, D. Makowski, Performance: 用于评估、比较和测试统计模型的 R 包。6，3139 (2021)。
8. D. Makowski, M. Ben-Shachar, D. Lüdecke, bayestestR：描述贝叶斯框架中效应及其不确定性、存在性和显著性。Journal of Open Source Software。4，1541 (2019)。
在某些情况下，不要将过滤条件放在“WHERE”子句中
如果我们在 LEFT JOIN ON 之后放置条件，会发生什么？
https://christophertao.medium.com/?source=post_page---byline--744829da6248--------------------------------https://towardsdatascience.com/?source=post_page---byline--744829da6248-------------------------------- Christopher Tao
·发表在Towards Data Science ·7 分钟阅读·2024 年 2 月 26 日
–
只要您点击进入这篇文章，我相信您应该了解 SQL。您还必须明白，在 SELECT 查询中，我们应该将条件放在 WHERE 子句中。然而，让我问你一个问题，看看你是否能立即回答。
如果我们在 LEFT JOIN … ON …子句中放置过滤条件会发生什么？
SELECT *
FROM Employee e LEFT JOIN Department d
ON e.dept_id = d.id
AND e.name = 'Chris'
如果您对上述查询的行为不确定，或者认为它等同于以下查询，请阅读我的文章，我会告诉您它们之间的区别。
SELECT *
FROM Employee e LEFT JOIN Department d
ON e.dept_id = d.id
WHERE e.name = 'Chris'
1. 验证结果
为了演示目的，我创建了两个带有简单虚拟数据的表如下。
员工表
我们真的需要深度学习来进行海岸监测吗？
深入探讨机器学习与传统海岸侵蚀监测方法的比较
https://conorosullyds.medium.com/?source=post_page---byline--f706b1cdd47d--------------------------------https://towardsdatascience.com/?source=post_page---byline--f706b1cdd47d-------------------------------- Conor O’Sullivan
·发表于 Towards Data Science ·14 分钟阅读·2024 年 9 月 10 日
–
照片由 thiago japyassu 提供，来自 Unsplash
深度学习（DL）是解决这个问题的唯一途径。这是我阅读的许多研究中隐含的假设。我总是倾向于同意这种看法。但，可能是因为如果没有它，我的博士研究将毫无意义。
幸运的是，我越是阅读，越能意识到遥感充满了机器学习可以提供帮助的问题。这些问题包括监测空气质量、估算土壤湿度、评估作物健康状况和追踪自然灾害。我的研究领域——海岸侵蚀监测，也是如此。
海岸线很长！这意味着我们需要自动化某些任务，才能有效地监测整个海岸线。同时，由于土地开发、云层覆盖以及海上风浪等因素引起的噪音，使得传统的确定性方法可能会失败。正是在处理这些变化时，机器学习才能发挥其优势。
深度学习，作为机器学习的一个子领域，已成为遥感中的一项宝贵工具，为前所未有的挑战提供了解决方案，并在遥感应用中创造了新的机会，
成为数据科学家需要学位吗？
不需要，但它肯定会有所帮助。
https://medium.com/@egorhowell?source=post_page---byline--7fb8ba5ba5f9--------------------------------https://towardsdatascience.com/?source=post_page---byline--7fb8ba5ba5f9-------------------------------- Egor Howell
·发布于Towards Data Science ·阅读时间 8 分钟·2024 年 5 月 29 日
–
我经常被问到：“我需要学位才能进入数据科学领域吗？”
简短的回答是：不需要，但它肯定会有所帮助。
让我在本文中详细解释我的意思。希望它能为这个问题提供一些启示，或者给你一些有数据和研究支持的例子。
作为免责声明，本文中的一切仅代表我的个人观点，并非专业建议。
数据科学家的背景
“大学或学院是否还值得上？”这个问题已经被回答并讨论得不堪重负。
总是给出相同的回答：“除非你从事需要学位的工作，比如医生或律师，否则它不值得。”
我同意这个说法，因为它在技术上是正确的，但这个问题远比这复杂。
根据 Coursera 的这篇文章，以及 Zippia 的研究，51%的数据科学家拥有学士学位，34%拥有……
你真的了解 Python 中的 *args 吗？
图片来源：Miguel Á. Padriñán 来自 Pixabay
*args 的全面指南与实用示例
https://christophertao.medium.com/?source=post_page---byline--1e3402c77190--------------------------------https://towardsdatascience.com/?source=post_page---byline--1e3402c77190-------------------------------- Christopher Tao
·发表于 Towards Data Science ·阅读时间 9 分钟·2024 年 1 月 29 日
–
作为 Python 中最独特的语法之一，*args 在编程过程中能为我们提供极大的灵活性和便利。我会说，它们体现了所谓的“Pythonic”风格和 Python 的禅意。
然而，我发现学习者往往难以理解它们。在本文中，我将尽我所能解释这个 Python 中标志性的概念，并根据我的知识提供实际应用案例。希望能帮助你更好地理解它。
1. “*args” 到底是什么？
*args 代表“参数”。它允许我们将任意数量的 位置参数（稍后解释）传递给一个函数。在函数内部，我们可以将所有位置参数存储为元组。因此，我们可以在函数中对这些参数元组进行任意操作。
这是一个简单的 *args 示例。
def add_up(*numbers):
result = 0
for num in numbers:
result += num
return result
print(add_up(1, 2, 3))
医生利用多模态数据；医疗 AI 也应该如此
整合多模态数据使新一代医疗 AI 系统能够更好地捕捉医生的思维和决策过程
https://medium.com/@fimafurman?source=post_page---byline--6475c247b243--------------------------------https://towardsdatascience.com/?source=post_page---byline--6475c247b243-------------------------------- Fima Furman
·发布于 Towards Data Science ·阅读时间 9 分钟·2024 年 9 月 25 日
–
多模态医疗 AI 的数据类型和应用示意图。图片由作者提供。
多模态 AI 模型利用来自各种格式的数据，如文本、图像和音频，为用户提供更全面的医疗情况理解。这些模型由于能够处理和整合多种数据类型，能够呈现出比任何单一数据类型更为全面的健康图景，因此正在迅速发展。随着变换器架构和大型语言模型（LLMs）的崛起，这些模型广泛适用于多种数据模式，开发人员获得了新的工具来综合这些数据格式。谷歌的Gemini 多模态 AI及其他前沿的生成型 AI 模型能够无缝理解并综合文本、视频、图像、音频和代码（遗传或计算）等数据格式。尽管在过去几年中，医疗 AI 取得了令人振奋的进展，但其应用进展缓慢，现有应用程序通常只针对非常具体和狭窄的使用案例。医疗 AI 的未来在于多模态应用，因为它们反映了医生的临床过程，医生在做出评估时必须考虑许多因素和数据来源。能够在这个充满巨大潜力的领域执行的开发者和公司，将在 AI 辅助医疗的未来中占据重要角色。
利用多模态数据的好处
医学数据本质上是多模态的，AI 系统应当反映这一现实。在评估患者时，医生会利用多种数据来源，如病历记录、医学影像、音频录音和基因序列。传统上，AI 应用程序被设计来处理这些单一数据类型中的特定、狭义的任务。例如，某个 AI 系统可能擅长识别 CT 扫描中的肺结节，但它无法将这些数据与患者的报告症状、家族史和基因信息结合起来，帮助医生诊断肺癌。相比之下，多模态 AI 应用可以整合多种数据类型，将 LLMs 的灵活性与专业 AI 系统的专长相结合。研究表明，这些系统在传统 AI 任务中的表现也优于单模态 AI 系统，准确度提升了 6–33%。
多模态 AI 模型还致力于打破医学专业之间的壁垒。由于医学领域的专业化不断增加，研究和数据的不断扩展，医学已经发展成一个碎片化的格局，放射学、内科和肿瘤学等不同领域可能会独立运作。处理复杂疾病的患者通常需要跨多个专家团队的协作，然而，由于沟通不畅，关键的见解可能会丧失。多模态 AI 模型通过跨专业领域获取知识，弥合这些鸿沟，确保患者能够从所有相关领域的最新医学进展中受益。
不同医学数据模态概述
医学数据占全球所有数据的 30%以上，且有多种形式。以下是其中一些最突出的数据形式（非详尽列表）：
医学影像
医学影像在医疗诊断和治疗规划中扮演着至关重要的角色，以至于它已经成为一个独立的专业（放射学）。CT 扫描和 X 光通常用于可视化骨骼结构并检测骨折或肿瘤，而超声波在监测胎儿发育和检查软组织方面至关重要。医生使用病理切片图像分析组织样本，以检测癌症等疾病。像卷积神经网络（CNN）这样的 AI 算法，通过处理大量带标签的图像，学习识别这些图像中的模式和异常。这些工具帮助放射科医生和其他医生更快速、更准确地解读图像。
组学
近年来，由于测序成本的下降，组学数据，包括基因组学、转录组学和蛋白质组学，迅速增长。它通过提供对疾病分子机制的洞察，彻底改变了个性化医学。在一个多模态医学 AI 系统中，组学数据可以帮助更好地理解患者对某些疾病的易感性以及对治疗选项的潜在反应。例如，BRCA 基因中的特定突变表明，患者更有可能发展为某些类型的癌症。
患者与电子健康记录（EHR）笔记
传统上，由于缺乏结构，患者笔记（临床观察、治疗计划等）一直是分析的难点。然而，大型语言模型（LLMs）可以利用这些笔记提取见解、识别模式，并支持大规模数据分析，这在过去是无法实现的。例如，LLMs 可以阅读临床试验中潜在患者的笔记，并识别符合资格要求的患者——这是以前需要大量人工操作的任务。
可穿戴设备数据
健康监测传感器，如可穿戴健身追踪器，可以实时测量心率、血压、睡眠模式和血糖水平等生命体征。人工智能应用可以分析这些时间序列数据，发现趋势并预测健康事件。此类应用可以通过提供个性化的健康建议，帮助患者，并帮助医生在医院之外监控患者的健康状况。
音频记录
音频记录，如心脏和肺部听诊，通常用于诊断某些形式的疾病。医生通过心脏听诊来标记心脏杂音的范围和强度，而肺部听诊可以帮助识别肺炎等疾病。人工智能系统可以分析这些音频记录，检测异常，并协助更快且成本更低的诊断。
病理学
病理数据来源于组织样本和显微图像，在诊断癌症等疾病中起着关键作用。人工智能算法可以分析这些数据源，识别异常的细胞结构、分类组织类型，并检测出提示疾病的模式。通过处理大量病理数据，人工智能可以帮助病理学家做出更准确的诊断，标出潜在的关注区域，甚至预测疾病进展。事实上，哈佛医学院和麻省理工学院的研究人员最近推出了一个多模态生成式 AI 助手用于人体病理学，帮助病理学家处理常见的医疗任务。
集成标注图像与文本数据的示例。图像由作者提供。
多模态人工智能模型的应用
多模态算法有潜力开启人工智能驱动医疗应用的新范式。多模态人工智能的一个有前景的应用是个性化医学，在这种应用中，系统通过整合患者的病情、病史、生活方式和基因组等数据，预测最有效的治疗方案。例如，考虑一个旨在为肺癌患者识别最有效治疗方案的应用程序。该应用可以考虑患者的基因组特征、病理（组织样本）图像和记录、放射学图像（肺部 CT 扫描）和记录，以及病史临床记录（收集吸烟史和环境影响等因素）。通过整合所有这些数据源，应用程序可以推荐最适合患者独特体质的治疗方案。这种方法已经在黄等人研究中显示出有前景的结果，研究人员能够根据患者的基因表达谱预测他们对常规化疗药物的反应，准确率超过 80%。这种方法将有助于最大化治疗效果，并减少通常与寻找合适药物或干预措施相关的试错过程。
另一个关键的应用场景是提高诊断和预后速度与准确性。通过整合医学影像、实验室结果和病历记录等数据源，多模态医学人工智能系统可以帮助医生获得全面的洞察。例如，Tempus Next 利用超声心动图和心电图的波形数据、电子健康记录文本数据以及腹部放射影像（CT 扫描、超声波检查）来帮助心脏病专家诊断和预测患者的心脏问题风险，如腹主动脉瘤和房颤。Optellum 的虚拟结节诊所采用类似的方法，通过使用 CT 扫描和临床记录帮助诊断肺癌。这类应用不仅提高了诊断的准确性，还节省了医生的时间，从而帮助解决持续的医生短缺问题，并降低了医疗成本。
多模态 AI 还将通过集成来自可穿戴设备、家庭监测系统和患者自我报告的记录的数据，推动远程患者监测和远程医疗的重大进展，为患者健康状态提供持续的实时洞察。这一能力对慢性病管理尤其重要，因为持续的监测可以检测到恶化的早期迹象，并及时采取干预措施。例如，一个 AI 系统可以监测来自Eight Sleep Pod的患者睡眠数据以及来自Levels（连续血糖监测）的血糖数据，以识别糖尿病前期患者的恶化情况。医生可以利用这一早期预警提出积极建议，帮助患者避免进一步的健康下降。这项技术有助于减少医院再入院率，改善慢性病的整体管理，使医疗服务更具可及性，并减轻整个医疗系统的负担。
构建多模态 AI 模型的方法
构建多模态系统的不同方法的示意图。图片由作者提供。
研究人员目前正在尝试不同的方法来构建多模态医疗 AI 系统，研究仍处于初步阶段。谷歌团队正在探索的三种主要系统开发方法是：
-
工具使用 — 在这种方法中，一个主控 LLM 将不同数据源的分析任务外包给专门的软件子系统，这些子系统是针对特定数据形式训练的。例如，LLM 可能会将胸部 X 光片转交给放射学 AI 系统，心电图分析交给专门的波形分析系统，然后将这些响应与患者记录整合，以评估心脏健康。这种方法允许子系统之间具有灵活性和独立性，使得每个特定任务都能使用最佳的工具。
-
模型嫁接 — 这种方法涉及将每个相关领域的专门神经网络进行调整，并直接集成到 LLM 中。例如，一个经过训练用于解读医学图像的神经网络可以通过将其输出直接映射到 LLM 的输入空间，将其嫁接到 LLM 中。这种方法利用了现有的优化模型，允许模块化开发，尽管它需要为每个特定的模型和领域创建适配器。
-
通用系统 — 最具雄心的方法是构建一个能够本地处理所有数据模态的单一集成系统。这种方法使用统一模型，例如Med-PaLM M，它将语言模型与视觉编码器结合，处理多种数据类型。尽管这种方法最大化了灵活性和信息传递，但也伴随着更高的计算成本，并可能在领域专门化和系统调试性方面面临挑战。
实现多模态 AI 模型的挑战
尽管构建多模态人工智能模型充满了巨大潜力，但在实现可行系统时仍面临诸多挑战。以下是一些挑战：
-
数据标注 — 为了启用监督学习，机器学习算法需要由专家人工标注数据，并正确识别特征。跨领域识别专家以标注不同类型的数据模态可能具有挑战性。模型构建者应考虑与具有跨模态专业知识的数据标注服务提供商合作，例如Centaur Labs。
-
避免偏见 — 部署 AI 系统在医疗环境中的最大风险之一是它们可能加剧现有的偏见和医疗不平等。多模态系统可能进一步加深这种偏见，因为代表性不足的群体在系统为其构建的一个或多个模态中更可能缺失数据。为了避免偏见，模型构建者应考虑减少 AI 应用中偏见的技术。
-
监管 — 数据隐私法规如 HIPAA 对患者数据的共享和使用施加了严格的控制，这使得开发人员在整合和关联跨不同模态的数据时面临挑战。这需要额外的开发工作以确保合规性。
-
采纳与信任 — 许多传统的 AI 系统发现影响最大的障碍是推动医疗用户社区的采纳和信任。医生们担心 AI 输出的准确性和一致性，不希望在使用这些系统来为患者护理提供决策前，因过度信任而危及患者健康。多模态 AI 模型将面临类似的采纳障碍。开发人员必须与此类系统的最终用户密切合作，建立信任，确保系统能够融入现有的临床工作流程。
-
数据格式共享标准缺乏 — 对于许多数据格式（例如，组织图像），不同提供者之间没有标准化的数据共享协议。这种互操作性的缺失可能会阻碍整合开发强大 AI 模型所需的数据源。为了加快在（目前）尚未标准化的医疗数据领域中运行 AI 系统的开发和应用，研究与开发社区应制定通用的数据共享标准/框架，并确保各机构遵守这些标准。
结论
多模态 AI 代表了医疗应用的未来，通过整合和全面利用数据，提供了革命化医疗保健的潜力，提升应用的灵活性、准确性和能力。如果这些应用能够有效开发和部署，它们有望降低医疗成本，扩大可及性，并提供更高质量的患者护理和结果。
知识和技术的最重大进展通常来自于跨领域的见解融合。以列奥纳多·达·芬奇为例，他将自己在绘画和流体力学方面的知识应用于心脏和生理学的研究。医疗人工智能也不例外。通过将计算机科学的发现融入医学，开发者开启了一波初步的突破。现在，融合多种数据模态的前景将带来第二波创新，这一波创新将由日益智能化的人工智能系统推动。
文档提取是 GenAI 的杀手级应用
未来已经到来，你不会看到致命的机器人。你会看到用于繁琐办公室工作的优秀自动化。
https://urimerhav.medium.com/?source=post_page---byline--9e9c816e6caf--------------------------------https://towardsdatascience.com/?source=post_page---byline--9e9c816e6caf-------------------------------- Uri Merhav
·发表于Towards Data Science ·阅读时间 9 分钟·2024 年 8 月 13 日
–
几乎十年前，我在 LinkedIn 的著名数据标准化团队担任机器学习工程师。从我加入的那一天到离开的那一天，我们仍然无法自动读取一个人的个人资料，并可靠地理解某人在所有语言和地区中的职级和职位。
乍一看，这看起来很简单。“软件工程师”已经足够明确了，对吧？那如果有人只写“助理”呢？如果他在沃尔玛，那可能是一个低级零售工人；如果他在律师事务所，那可能是一个高级律师。但你大概已经知道了——你知道什么是Java Fresher吗？什么是Freiwilliges Soziales Jahr？这不仅仅是了解德语——它翻译为“志愿社会年”。但是，什么是代表这一角色的标准职称呢？如果你有一大堆已知的职位名称，你会把它映射到哪里？
我加入了 LinkedIn，我离开了 LinkedIn。我们取得了一些进展，但即使是最简单的常规文本——一个人的简历，依然难以理解。
非常困难变得简单
你可能不会感到惊讶地发现，这个问题对于像 GPT-4 这样的 LLM 来说是微不足道的
对 GPT 来说轻松简单（来源：我和 GPT）
等等，我们是公司，不是一个聊天终端上的人，我们需要结构化的输出。
（来源：GPT）
啊，这样好多了。你可以对那些最复杂且具有文化特定性的提问进行重复练习。更棒的是，当你得到一个完整的个人资料时，你可以重复这个练习，这样你就能获得更多的背景信息；使用代码时，你可以在商业环境中稳定地使用这些结果，而不仅仅是进行一次性的聊天。通过更多的工作，你可以将结果强制转换为一个标准的可接受职位标题分类法，这样它就能被索引。毫不夸张地说，如果你复制并粘贴某个人的全部简历并正确提示 GPT，你将超越十年前一些相当聪明的人的最佳成果，这些人花了多年时间在这方面工作。
高价值办公室工作 == 理解文档
标准化简历的具体例子很有趣，但它仍然局限于技术一直在努力的领域——一个自然应用人工智能工具的技术网站。我认为这里有更深层的机会。全球 GDP 的大部分来自办公室工作，这些工作归结为将专家级人类智能应用于从文档中反复提取洞见，且需要考虑背景信息。以下是一些复杂度逐渐增加的例子：
-
费用管理就是读取发票并将其转化为标准化视图，显示支付了什么、何时支付、使用了哪种货币以及属于哪个费用类别。可能这个决策是基于关于业务、报销人等背景信息做出的。
-
医疗保险理赔裁定过程就是阅读一堆杂乱的发票和临床记录，并判断：“好吧，总共有一次胸部 X 光检查，包含了一些重复项目，总费用为 800 美元，且它对应健康保险政策中的 1-C 类别”。
-
贷款核准员可能会查看一堆申请人的银行账单并回答一系列问题。同样，这之所以复杂，仅仅因为输入信息杂乱无章。实际的决策过程像是：“现金的平均流入和流出是多少，多少用于贷款偿还，其中有多少是一时性的支出，多少是实际的经常性收入”。
关于文本推理是 LLM 的强项。
到现在为止，大型语言模型（LLMs）因易发生幻觉（即胡乱编造）而臭名昭著。现实情况更为复杂：幻觉实际上在某些情境下是可预测的结果，而在其他情况下则几乎可以保证不会发生。
产生幻觉的地方是当你让它回答事实性问题，并期望模型仅凭其对世界的固有知识“知道”答案时。大语言模型（LLM）在自我反思自己对世界的知识方面表现得很差——它们能做到这一点更多是一个非常偶然的结果。它们并没有专门为此任务进行训练。它们的训练目标是生成可预测的文本序列补全。当 LLM 被绑定到输入文本并需要回答关于该文本内容的问题时，它不会产生幻觉。如果你将这篇博客文章复制并粘贴到 chatGPT 中，问它是否教你做美国苹果派，你会 100%得到正确答案。对于 LLM 来说，这是一个非常可预测的任务，它看到一段文本，并尝试预测一个有能力的数据分析师如何用预定义的字段和预定义的结果来填写，这些结果之一是{"is cooking discussed": false}。
以前作为 AI 顾问，我们多次解决涉及从文档中提取信息的项目。事实证明，在保险、金融等领域这方面有很多用途。客户之间对 LLM 的恐惧（“LLM 会产生幻觉”）与实际上摧毁我们的原因（我们没有正确提取表格，所有错误都源于此）之间存在很大的差距。LLM 确实失败了——当我们没有以清晰且不含歧义的方式呈现输入文本时，它们会失败。构建能够推理文档的自动化管道需要两个必要的成分：
-
完美的文本提取，将输入文档转换为干净、易懂的纯文本。这意味着需要处理表格、勾选框、手写注释、可变文档布局等。现实世界表单的复杂性需要转化为 LLM 能够理解的清晰纯文本。
-
健壮的架构，明确规定你希望从给定文档类型中获得哪些输出，如何处理边缘案例，使用什么数据格式等等。
文本提取比看起来要复杂得多
这就是导致 LLM 崩溃和产生荒谬输出的原因：
-
输入中有复杂的格式，例如双列布局，你从例如 PDF 中从左到右复制并粘贴文本，完全把句子从上下文中剥离开来。
-
输入中有复选框、勾选标记、手写注释，而你在转换为文本时完全忽略了这些
-
更糟糕的是：你认为可以绕过转换为文本的步骤，指望只是粘贴一张文档的图片，让 GPT 自行推理。这会让你进入幻觉之城。只需让 GPT 转录一张带有空白单元格的表格图片，你就会看到它高兴地胡乱发挥，乱编一些东西。
时常记住现实世界文档的混乱程度是很有帮助的。这里有一份随意的税表：
当然，真实的税表上这些字段通常是填写完整的，且常常是手写的
或者这是我的简历
来源：我的简历
或者一个公开的示例实验报告（这是 Google 的首页结果）
来源：research gate，公共领域图像
顺便说一句，最糟糕的事情就是让 GPT 的多模态能力来转录一个表格。如果你敢试试——一开始看起来没问题，但它会为某些表格单元格随意编造内容，完全脱离上下文等。
如果这个世界有问题，那就建立一个 SaaS 公司来解决它。
当我们需要理解这些类型的文档时，我和我的联合创始人Nitai Dean感到困惑，因为没有现成的解决方案可以帮助我们理解这些文本。
有些人声称能够解决这个问题，比如 AWS Textract。但在我们测试的任何复杂文档中，它们都会犯很多错误。接着就是那些必须处理的小细节，比如识别勾选框、单选按钮、删除线文本、表单上的手写涂鸦等。
所以，我们建立了Docupanda.io——它首先生成任何你输入的页面的清晰文本表示。左边是原始文档，右边是文本输出。
来源：docupanda.io
表格也类似地处理。在背后，我们只是将表格转换为人类和 LLM 可读的 Markdown 格式：
来源：docupanda.io
使用 LLM 来理解数据的最后一块拼图是生成并遵循严格的输出格式。虽然我们可以让 AI 将输出格式化为 json，但为了对数据应用规则、推理、查询等——我们需要让它以规律的方式运作。数据需要符合预定义的插槽集，我们将用内容填充这些插槽。在数据领域，我们称之为模式。
构建模式是一个试错过程……这是一个 LLM 可以做的事情。
我们需要模式的原因是，数据没有规律性是毫无意义的。如果我们在处理患者记录，并且它们映射到“male”、“Male”、“m”和“M”——那我们做得很糟糕。
那么，如何构建一个模式呢？在教科书中，你可能通过长时间地坐在那里，盯着墙壁，来定义你想要提取的内容。你坐在那里，思考你的医疗数据操作，并说“我想提取患者姓名、日期、性别和他们的医生姓名。哦，性别必须是 M/F/Other。”
在现实中，为了定义从文档中提取什么内容，你得盯着文档看……很多。你一开始可能像上面那样，但是当你看文档时，你会发现其中一个文档列出了多个医生的名单，而不是一个。还有些文档也列出了医生的地址。一些地址包含单元号和楼栋号，所以你可能需要为此预留一个字段。事情就是这么一件接一件地发生。
我们意识到，能够准确定义你想要提取的所有内容，是既不简单、又困难，但使用 AI 是完全可以解决的。
这是 DocuPanda 的关键部分。我们不是仅仅要求一个 LLM 为每个文档即兴生成输出，而是建立了一个机制，让你可以：
-
使用自由语言指定你需要从文档中提取的内容。
-
让我们的 AI 映射 多个 文档，并找出一个能回答所有问题、并适应实际文档中观察到的漏洞和不规则之处的模式。
-
根据反馈调整模式，以适应你的业务需求。
最终你得到的是一个强大的 JSON 模式 —— 一个模板，明确指出你想从每个文档中提取的内容，并且能够处理数十万份文档，从中提取所有答案，同时遵循像始终以相同格式提取日期、尊重一组预定义类别等规则。
来源：docupanda.io
还有更多！
就像任何兔子洞一样，总是有比最初看到的更多的东西。随着时间的推移，我们发现需要更多的东西：
-
组织经常需要处理大量匿名文档，因此我们会自动对它们进行分类，并决定应用哪种模式。
-
文档有时是多个文档的拼接，你需要一个智能的解决方案将一篇非常长的文档拆分成其原子化的独立组件。
-
使用生成的结果查询正确的文档非常有用。
如果从这篇文章中能得到一个启示，那就是你应该研究如何利用 LLM 来以规范的方式理解文档。如果有两个启示，那就是你也应该试试 Docupanda.io。我之所以构建它，是因为我相信它。也许这就足够成为尝试它的理由？
未来的办公室工作人员（来源：unsplash.com）
使用大语言模型进行文档解析——附带代码
你将不再考虑使用正则表达式
https://zoumanakeita.medium.com/?source=post_page---byline--9229fda09cdf--------------------------------https://towardsdatascience.com/?source=post_page---byline--9229fda09cdf-------------------------------- Zoumana Keita
·发布于Towards Data Science ·14 分钟阅读·2024 年 7 月 25 日
–
动机
多年来，正则表达式一直是我解析文档的首选工具，我相信对于许多其他技术人员和行业来说也是如此。
尽管正则表达式在某些情况下非常强大且成功，但它们在面对现实世界文档的复杂性和多变性时常常显得力不从心。
另一方面，大语言模型提供了一种更强大且灵活的方法来处理多种类型的文档结构和内容类型。
系统的一般工作流程
理解正在构建的系统的主要组件总是很重要的。为了简单起见，我们聚焦于一个科研论文处理的场景。
使用 LLM 进行文档解析工作流程（作者：Zoumana Keita）
-
工作流程总体上包括三个主要组件：输入、处理和输出。
-
首先，文档（在此案例中为 PDF 格式的科研论文）提交进行处理。
-
处理组件的第一个模块从每个 PDF 中提取原始数据，并将其与包含大语言模型指令的提示结合起来……
使用 MkDocs 记录 Python 项目
使用 Markdown 快速为你的项目创建一个漂亮的文档页面
https://gustavorsantos.medium.com/?source=post_page---byline--60e26b64380e--------------------------------https://towardsdatascience.com/?source=post_page---byline--60e26b64380e-------------------------------- Gustavo R Santos
·发布于 Towards Data Science ·阅读时间 9 分钟·2024 年 11 月 22 日
–
图片由 OpenAI 的 DALL·E 生成。openai.com. 使用 MkDocs 进行 Python 文档编写。
介绍
项目文档是必需的。非常必要，我必须强调这一点。
在我职业生涯的开始，我通过痛苦的经历学到了一个重要的教训，那就是一个项目必须有文档。
让我们回到过去——2000 年代——那时我在大型美国公司担任客户代表。我是一个团队的一员，我和我的同事们大约在同一个月加入了公司。因此，一段时间内，我们不需要担心，因为刚刚开始新工作的人不可能在几周或几个月后就请假。
然而，经过一段时间，这种情况不可避免地会发生。而且我们每个人都被分配去互相备份。正是在那时，文档开始在我的职业生涯中扮演重要角色。
那天，第一个人请了几天假，我慌了！我开始工作时，不知道该做什么，甚至不知道从哪里开始。任务源源不断地涌来，而我则在努力弄清楚如何处理它们。
最终，一切都顺利解决了。我终于弄明白了并继续前进。但从那天起，我知道无论是休假还是团队成员的变动，如晋升等，文档都必须到位。
数据驱动的故事讲述是否需要客观？
在数据驱动故事中找到效率与吸引力之间的平衡
https://medium.com/@michalszudejko?source=post_page---byline--7053a29f3d47--------------------------------https://towardsdatascience.com/?source=post_page---byline--7053a29f3d47-------------------------------- Michal Szudejko
·发表于Towards Data Science ·阅读时长：11 分钟·2024 年 7 月 29 日
–
来源：图片由作者提供。
我一直想学会解魔方。 我在 10 到 12 岁的时候有一个魔方。我从未成功过。即使成功了，那也只是通过“调整”方块之间的贴纸。所以，那并不算真正的成功。
问题是什么？我曾以为解魔方是一种智力和逻辑的练习。我错了。多亏了 YouTube 上的一些教程，30 年后，我终于学会了如何解魔方。事实证明，你可以通过遵循一系列算法来解决它。没有魔法——只是记住一些规则和步骤（幸运的是，至少对我来说，还是需要一些脑力努力）。
数据驱动故事的魔力
我曾经也遇到过类似的数据驱动故事讲述问题。让我们把它放在一个更具体的情境中，就像我 10 岁时试图解决魔方一样。你可能认为数据驱动的故事讲述是用来向观众讲述一个关于数据和分析的客观故事——结论、假设、优点和缺点。观众应该知道所有的事实，才能做出充分知情的决策。这也可能包括分享那些证明是错误的分析场景。最终，观众应该…
半监督学习是否有助于训练更好的模型？
评估半监督学习如何利用未标注数据
https://medium.com/@reinhard.sellmair?source=post_page---byline--338283d1f4e9--------------------------------https://towardsdatascience.com/?source=post_page---byline--338283d1f4e9-------------------------------- Reinhard Sellmair
·发表于 Towards Data Science ·阅读时长 8 分钟·2024 年 9 月 9 日
–
作者提供的图像 — 使用必应的图像生成器创建
数据科学家面临的最常见挑战之一是缺乏足够的标注数据来训练一个可靠且准确的模型。标注数据对于监督学习任务（如分类或回归）至关重要。然而，在许多领域，获取标注数据可能成本高昂、耗时或不可行。另一方面，未标注数据通常容易收集，但它们并不能直接用于训练模型。
我们如何利用未标注数据来提高监督学习模型的性能？这正是半监督学习发挥作用的地方。半监督学习是机器学习的一个分支，它结合了标注数据和未标注数据，以训练一个比仅使用标注数据表现更好的模型。半监督学习背后的直觉是，未标注数据可以提供关于数据潜在结构、分布和多样性的有用信息，这有助于模型更好地泛化到新的和未见过的样本。
在这篇文章中，我介绍了三种可以应用于不同类型数据和任务的半监督学习方法。我还将评估它们在实际数据集上的表现，并将其与仅使用标注数据的基准进行比较。
什么是半监督学习？
半监督学习是一种机器学习方法，它同时利用已标记数据和未标记数据来训练模型。已标记数据是指那些已知输出或目标变量的样本，比如分类任务中的类别标签或回归任务中的数值。未标记数据是指那些没有已知输出或目标变量的样本。半监督学习可以利用在实际问题中通常大量存在的未标记数据，同时也利用较少的已标记数据，而后者通常更加昂贵或耗时。
使用未标记数据来训练监督学习方法的基本思想是通过监督或无监督学习方法给这些数据打上标签。尽管这些标签可能不像实际标签那样准确，但拥有大量此类数据可以提升监督学习方法的表现，相比只用已标记数据进行训练。
scikit-learn 包提供了三种半监督学习方法：
-
自训练：首先仅在已标记数据上训练分类器，以预测未标记数据的标签。在下一轮中，另一个分类器在已标记数据和从未标记数据中预测出高置信度标签的结果上进行训练。这个过程会重复进行，直到没有新的高置信度标签被预测，或者达到最大迭代次数为止。
-
标签传播：构建一个图，其中节点表示数据点，边表示它们之间的相似性。标签通过图进行迭代传播，从而使算法能够根据已标记数据与未标记数据之间的连接关系，将标签分配给未标记数据点。
-
标签传播：使用与标签传播相同的概念。区别在于标签扩散使用软分配，在这种方法中，标签是根据数据点之间的相似性进行迭代更新的。这种方法还可能“覆盖”已标记数据集的标签。
为了评估这些方法，我使用了一个糖尿病预测数据集，该数据集包含患者的年龄、BMI 等特征，以及描述患者是否患有糖尿病的标签。该数据集包含 100,000 条记录，我将其随机分为 80,000 条训练数据、10,000 条验证数据和 10,000 条测试数据。为了分析在不同标签数据量下学习方法的有效性，我将训练数据分为标签数据集和无标签数据集，其中标签大小表示有多少样本已标记。
数据集划分（图片来源：作者）
我使用验证数据评估了不同参数设置，并使用测试数据评估了各方法在参数调优后的性能。
我使用了 XGBoost 进行预测，并使用 F1 得分来评估预测性能。
基准
基准模型用于将自学习算法与不使用任何无标签数据的情况进行比较。因此，我在不同大小的标签数据集上训练了 XGB，并计算了验证数据集上的 F1 得分：
基准分数（图片来源：作者）
结果显示，对于少于 100 个样本的训练集，F1 得分相对较低，但随着样本量的增加，得分稳定提高，到达 1,000 个样本时达到了 79%的得分。更大的样本量对 F1 得分的提升几乎没有作用。
自学习
自学习是通过多次迭代来预测无标签数据的标签，这些标签将在下一次迭代中用于训练另一个模型。可以使用两种方法来选择将作为标签数据用于下一次迭代的预测：
-
阈值（默认）：选择所有置信度超过阈值的预测
-
K 最优：选择置信度最高的 k 个预测
我评估了默认参数（ST Default），并根据验证数据集调优了阈值（ST Thres Tuned）和 k 最优参数（ST KB Tuned）。这些模型的预测结果在测试数据集上进行了评估：
自学习得分（图片来源：作者）
对于小样本量（<100），默认参数（红线）的表现不如基准（蓝线）。对于更大的样本量，模型取得的 F1 得分稍微优于基准。调整阈值（绿线）带来了显著的提升。例如，在标签大小为 200 时，基准 F1 得分为 57%，而调整阈值后的算法 F1 得分为 70%。除了标签大小为 30 时的一个例外，调整 K 最优值（紫线）几乎与基准表现相同。
标签传播
标签传播有两种内置的核方法：RBF 和 KNN。RBF 核通过使用密集矩阵生成一个完全连接的图，这在大数据集上需要大量内存且耗时较长。为了考虑内存限制，我只使用了最多 3,000 个样本进行 RBF 核的训练。KNN 核则使用更节省内存的稀疏矩阵，这使我能够在最多 80,000 个样本的整个训练数据上进行训练。以下图表展示了这两种核方法的结果比较：
标签传播得分（图片由作者提供）
图表显示了不同标签传播方法在测试数据集上的 F1 得分与标签大小的关系。蓝色线表示基准线，与自我训练相同。红色线表示使用默认参数的标签传播，明显低于基准线，且所有标签大小下都表现较差。绿色线表示使用 RBF 核和调优后的参数 gamma 的标签传播。Gamma 定义了单个训练样本的影响范围。调优后的 RBF 核在小标签大小（<=100）下表现优于基准线，但在大标签大小下则表现较差。紫色线表示使用 KNN 核和调优后的参数 k 的标签传播，k 决定了使用多少个最近邻。KNN 核的表现与 RBF 核相似。
标签扩展
标签扩展是一种类似于标签传播的方法，但增加了一个控制实例应该接受多少邻居信息的参数 alpha。Alpha 的取值范围为 0 到 1，其中 0 表示实例保持其原始标签，1 表示实例完全采用邻居的标签。我还对 RBF 和 KNN 核方法进行了调优。标签扩展的结果显示在下图中：
标签扩展得分（图片由作者提供）
标签传播的结果与标签扩展非常相似，唯一显著的例外是，RBF 核方法在标签扩展中的测试得分低于基准得分，且所有标签大小下都如此，不仅仅是小标签大小。这表明，邻居标签的“覆盖”对这个数据集有相当负面的影响，该数据集可能只有少数离群点或噪声标签。另一方面，KNN 核方法不受 alpha 参数的影响。似乎这个参数只与 RBF 核方法相关。
所有方法的比较
接下来，我将所有方法的最佳参数进行对比。
最佳得分比较（图片由作者提供）
图表显示了不同半监督学习方法的测试得分与标签大小的关系。自我训练优于基线，因为它很好地利用了无标签数据。标签传播和标签扩散仅在小标签大小下超过基线，对于较大的标签大小，表现较差。
结论
结果可能因不同的数据集、分类器方法和评估指标而显著变化。半监督学习的性能依赖于许多因素，例如无标签数据的质量和数量、基本学习器的选择以及评估标准。因此，在没有适当的测试和验证的情况下，不应将这些发现推广到其他设置中。
如果你有兴趣进一步探索半监督学习，欢迎查看我的 Git 仓库并自行实验。你可以在这里找到这个项目的代码和数据。
我从这个项目中学到的一点是，参数调优对于显著提升这些方法的性能至关重要。通过优化参数，自我训练在任何标签大小下都优于基线，并且 F1 分数提高了最多 13%！标签传播和标签扩散仅在非常小的样本量下有所改进，但用户必须非常小心，以免得到比不使用任何半监督学习方法更差的结果。
在招聘过程中使用大型语言模型（LLM）会让你成为候选人中的骗子吗？
雇主们，放弃 AI 检测工具，改问一个重要的问题。
https://christineegan42.medium.com/?source=post_page---byline--99a05678536b--------------------------------https://towardsdatascience.com/?source=post_page---byline--99a05678536b-------------------------------- Christine Egan
·发表于 Towards Data Science ·12 分钟阅读·2024 年 1 月 6 日
–
我在 LinkedIn 上看到一篇来自一家咨询公司总监的帖子，描述了他如何布置一篇关于机器学习系统漂移的作文题目，以筛选潜在的候选人。
然后，基于他通过直觉建立的标准（“你可以闻到它”），他使用了四种不同的“AI 检测器”来“确认”申请者是否使用 ChatGPT 写了他们的作文回答。
[## Adam Sroka 在 LinkedIn 上发布：我自动拒绝了一个在招聘过程中使用 LLM 的候选人。我设定了一个讨论机器学习系统和漂移的任务。…|…
我自动拒绝了一个在招聘过程中使用 LLM 的候选人。我设定了一个讨论机器学习系统和漂移的任务。这个…
“涉嫌”机器人生成的作文的标准包括：
-
奇怪的句子结构
-
奇怪的类比
-
重复
-
从一个英语方言切换到另一个方言（在同一申请中的不同写作样本中）
有一个标准显然缺失：准确性。
这么做的理由是，使用 AI 工具是在试图颠覆候选人筛选过程。不用说，评论区简直乱套（而且非常典型 LinkedIn 风格）。
你的公司有数据战略吗？
这个复杂度矩阵可以显示你需要前往的方向
https://medium.com/@kminoguem?source=post_page---byline--e8944dbec4ae--------------------------------https://towardsdatascience.com/?source=post_page---byline--e8944dbec4ae-------------------------------- 凯特·米诺格
·发布于Towards Data Science ·7 分钟阅读·2024 年 4 月 8 日
–
如今每家公司似乎都沉浸于制定人工智能战略。从个体经营者到大型组织，过去 6 个月最常见的问题是——我如何在我的业务中使用人工智能？往往是在“我应该吗？”之前就已经有了这个问题。
幸运的是，随着由人工智能构建的新工具爆炸式增长，个人或公司比以往任何时候都更容易将人工智能融入到日常工作中，且我们每个人都应该不断思考，是否在用现有技术以最有效的方式开展工作。但当它成为一个重大投资或优先事项偏离时，就需要更深入的反思。
我坚信利用技术、人工智能和数据来推动现有的商业战略和优先事项。其目的是为了增强竞争优势，而不是成为优先事项本身。
当我思考人工智能及其成功潜力时，我会从三个方面来看——是的，当然有那些似乎得到所有关注和热议的机器学习和模型，但同样重要（如果不是更重要的话）的是那些构建或使用它的人，最后，还有支撑它的数据和基础设施。
所以，如果你是一个商业领袖，正在考虑 2024 年你应该制定什么样的人工智能战略，那么我建议你首先考虑，你当前的数据战略是什么？
狗狗排便指南针
狗狗行为的贝叶斯分析
https://datawondering.com/?source=post_page---byline--f95a4b9f2bf9--------------------------------https://towardsdatascience.com/?source=post_page---byline--f95a4b9f2bf9-------------------------------- Dima Sergeev
·发表于Towards Data Science ·阅读时间：22 分钟·2024 年 11 月 25 日
–
tl;dr
狗狗是朝北和朝南排便吗？结果证明它们确实是！想学会如何使用指南针应用程序、贝叶斯统计和一只狗（狗狗不包括在内）在家里测量这个现象吗？那就跟我一起看看吧！
引言
这是我的狗。它的名字叫 Auri，是一只 5 岁的凯文犬（Cavalier King Charles Spaniel）。
Auri（图由作者提供）
和许多其他狗主人一样，在我们散步时，我注意到 Auri 在需要去厕所时有一个非常独特的仪式。当他找到一个合适的地点时，他会开始围绕某个东西转圈，就像一个指南针。
起初，我只是觉得这种行为很有趣。毕竟，谁知道狗狗在想些什么呢？但过了一段时间，我记得曾读过一篇 2013 年的研究论文，标题是“Dogs are sensitive to small variations of the Earth’s magnetic field”。这项研究在相对较大的狗狗样本中进行，证实了“狗狗倾向于在北南轴对齐的情况下排便”。
这似乎是个有趣的研究课题！ 我心想。多么幸运，我正好有一个完美的实验对象——我自己心爱的狗狗。我决定复现这些发现，并通过 Auri 这个 N=1 的意外研究参与者来验证（或推翻！）这个假设。
就这样，我开始了长达几个月的数据收集之旅，记录了超过 150 次的“对齐”行为，如果你明白我的意思的话。
数据收集
对于我的研究，我需要记录每次 Auri 排便时的指南针数据。得益于现代科技的进步，我们不仅有iPad 的计算器应用，还可以在手机上使用相当精确的指南针。于是我决定使用这个。
方法非常简单。每次我的狗安静下来准备享受私人时光时，我就打开指南针应用程序，把手机与 Auri 的身体对齐，并截图。在原始论文中，作者优雅地将这种对齐称为胸椎（肩胛骨之间）朝向头部的指南针方向。非常科学。其实就是意味着指南针的箭头应指向与狗头相同的方向。
狗对地球磁场微小变化非常敏感，Vlastimil Hart 等人
总之，这就是我在几个月的时间里总共做了大约 150 次的事情。每当我似乎在拍摄我家狗狗的排泄行为时，我几乎能感觉到路人们混合着困惑和好奇的目光。但这值得吗？让我们来看看吧！
分析
我将在这里简要讨论数据提取和预处理，然后直接进入圆形分布和假设检验的部分。
和往常一样，所有代码都可以在我的 GitHub 上找到：Data Wondering。
如何处理应用截图？
数据收集后，我得到了若干张指南针应用截图：
指南针应用截图（作者提供的图片）
因为我懒得一张张图片地查看并耐心地写下指南针的度数，所以我决定将所有这些图片发送到我的笔记本，并自动化这个过程。
任务很简单——我只需要从这些图片中获取屏幕底部的大数字。幸运的是，有很多小型的预训练网络可以做基本的光学字符识别（OCR）。我选择了一个名为[easyocr](https://github.com/JaidedAI/EasyOCR)的包，虽然它稍慢一点，但免费且容易使用。
我将为你展示一个简单的例子，演示如何使用easyocr与[opencv](https://opencv.org/)结合，从单张截图中提取数字。
首先，我们加载图像并显示它：
import cv2
import os
image_dir = '../data/raw/'
img = cv2.imread(image_dir + 'IMG_6828.png')
plt.imshow(img)
plt.show()
图片由作者提供
接下来，我将图像转换为灰度，以去除任何颜色信息并减少噪声：
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
plt.imshow(gray, cmap='gray')
plt.show()
图片由作者提供
接着，我将放大感兴趣的区域：
图片来源：作者
最后，我使用easyocr从图片中提取数字。它提取了数字和置信度评分，但我只对数字感兴趣。
import easyocr
reader = easyocr.Reader(['en'])
result = reader.readtext(gray[1850:2100, 200:580])
for bbox, text, prob in result:
print(f"Detected text: {text} with confidence {prob}")
>> Detected text: 340 with confidence 0.999995182215476
就是这样！我写了一个简单的 for 循环来遍历所有截图，并将结果保存到 CSV 文件中。
这是完整预处理笔记本的链接：数据预处理。
图片来源：作者
转动轮盘：圆形分布
我通常不处理圆形分布，因此我必须做一些阅读。与我们常见的常规数据不同，圆形数据有一个特殊的属性：分布的“端点”是连接的。
例如，如果你考虑一天中的小时分布，你会发现 23:00 到 00:00 之间的距离与 00:00 到 01:00 之间的距离是相同的。或者，在指南针角度的情况下，359°和 0°之间的距离与 0°和 1°之间的距离是相同的。
即使计算样本均值也不是直截了当的。360°和 0°之间的标准算术均值将是 180°，尽管 360°和 0°指向完全相同的方向。
在我的情况下，计算算术均值和正确的均值时，我得到了几乎完全相反的估算。我使用这个不错的库中的辅助函数将角度转换为弧度：pingouin，并使用circ_mean函数计算均值。
from pingouin import circ_mean
arithmetic_mean = data['radians'].mean()
circular_mean = circ_mean(data['radians'])
print(f"Arithmetic mean: {arithmetic_mean:.3f}; Circular mean: {circular_mean:.3f}")
>> Arithmetic mean: 0.082; Circular mean: 2.989
接下来，我想可视化指南针分布。我使用冯·米塞斯分布来建模圆形数据，并使用matplotlib绘制极坐标图。
冯·米塞斯分布是圆形分布的正态分布类比。它由两个参数定义：均值位置μ和集中度κ。集中度参数控制分布的扩展，类似于方差的倒数。当κ为 0 时，分布是均匀的，随着κ的增加，分布会围绕均值收缩。
让我们导入必要的库并定义辅助函数：
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import vonmises
from pingouin import convert_angles
from typing import Tuple, List
def vonmises_kde(series: np.ndarray, kappa: float, n_bins: int = 100) -> Tuple[np.ndarray, np.ndarray]:
"""
Estimate a von Mises kernel density estimate (KDE) over circular data using scipy.
Parameters:
series: np.ndarray
The input data in radians, expected to be a 1D array.
kappa: float
The concentration parameter for the von Mises distribution.
n_bins: int
The number of bins for the KDE estimate (default is 100).
Returns:
bins: np.ndarray
The bin edges (x-values) used for the KDE.
kde: np.ndarray
The estimated density values (y-values) for each bin.
"""
bins = np.linspace(-np.pi, np.pi, n_bins)
kde = np.zeros(n_bins)
for angle in series:
kde += vonmises.pdf(bins, kappa, loc=angle)
kde = kde / len(series)
return bins, kde
def plot_circular_distribution(
data: pd.DataFrame,
plot_type: str = 'kde',
bins: int = 30,
figsize: tuple = (4, 4),
**kwargs
) -> None:
"""
Plot a compass rose with either KDE or histogram for circular data.
Parameters:
-----------
data: pd.DataFrame
DataFrame containing 'degrees'and 'radians' columns with circular data
plot_type: str
Type of plot to create: 'kde' or 'histogram'
bins: int
Number of bins for histogram or smoothing parameter for KDE
figsize: tuple
Figure size as (width, height)
**kwargs: dict
Additional styling arguments for histogram (color, edgecolor, etc.)
"""
plt.figure(figsize=figsize)
ax = plt.subplot(111, projection='polar')
ax.set_theta_zero_location('N')
ax.set_theta_direction(-1)
# add cardinal directions
directions = ['N', 'E', 'S', 'W']
angles = [0, np.pi / 2, np.pi, 3 * np.pi / 2]
for direction, angle in zip(directions, angles):
ax.text(
angle, 0.45, direction,
horizontalalignment='center',
verticalalignment='center',
fontsize=12,
weight='bold'
)
if plot_type.lower() == 'kde':
x, kde = vonmises_kde(data['radians'].values, bins)
ax.plot(x, kde, color=kwargs.get('color', 'red'), lw=2)
elif plot_type.lower() == 'histogram':
hist_kwargs = {
'color': 'teal',
'edgecolor': 'black',
'alpha': 0.7
}
hist_kwargs.update(kwargs)
angles_rad = np.deg2rad(data['degrees'].values)
counts, bin_edges = np.histogram(
angles_rad,
bins=bins,
range=(0, 2*np.pi),
density=True
)
widths = np.diff(bin_edges)
ax.bar(
bin_edges[:-1],
counts,
width=widths,
align='edge',
**hist_kwargs
)
else:
raise ValueError("plot_type must be either 'kde' or 'histogram'")
ax.xaxis.grid(True, linestyle='--', alpha=0.5)
ax.yaxis.grid(True, linestyle='--', alpha=0.5)
ax.set_yticklabels([])
plt.show()
现在，让我们加载数据并绘制图表：
data = pd.read_csv('../data/processed/compass_degrees.csv', index_col=0)
data['radians'] = convert_angles(data['degrees'], low=0, high=360)
plot_circular_distribution(data, plot_type='histogram', figsize=(6, 6))
plot_circular_distribution(data, plot_type='kde', figsize=(5, 5))
狗狗排便的圆形直方图（图片来源：作者）
从直方图中可以明显看出，Auri 在选择缓解方向时有自己的偏好。北方方向有明显的峰值，南方方向则有一个低谷。很棒！
狗狗排便的圆形 KDE（图片来源：作者）
通过 KDE 图，我们可以获得分布的更平滑表示。好消息是，它离均匀圆形非常远。
该是进行统计验证的时候了！
统计学显著的排便
就像圆形数据需要特别的可视化和分布处理一样，它也需要特别的统计检验。
我将使用之前提到的 pingouin 库中的几个测试。我将使用的第一个测试是 Rayleigh 检验，它是一个用于检验圆形数据均匀性的测试。原假设认为数据在圆周上均匀分布，备择假设认为数据不均匀。
from pingouin import circ_rayleigh
z, pval = circ_rayleigh(data['radians'])
print(f"Z-statistics: {z:.3f}; p-value: {pval:.6f}")
>> Z-statistics: 3.893; p-value: 0.020128
好消息，大家！p 值小于 0.05，我们拒绝原假设。Auri 的战略性排便位置并非随机的！
唯一的缺点是该测试假设分布只有一个模态，并且数据是从冯·米塞斯分布中采样的。唉，那我们换个方法试试吧。Auri 的数据显然有多个模态。
接下来是 V 检验。该测试检查数据是否具有特定的均值方向且非均匀。从文档中我们得知：
V 检验比雷 leigh 检验有更高的效能，并且如果有理由相信某个特定的均值方向，推荐使用 V 检验。
完美！让我们试试。
从分布来看，很明显 Auri 更喜欢南方方向。我将均值方向设置为 π 弧度（南方），并运行测试。
from pingouin import circ_vtest
v, pval = circ_vtest(data['radians'], dir=np.pi)
print(f"V-statistics: {v:.3f}; p-value: {pval:.6f}")
>> V-statistics: 24.127; p-value: 0.002904
现在我们开始有进展了！p 值接近零，我们拒绝原假设。Auri 是一个统计学上显著的南向排便者！
贝叶斯排便：数学部分
现在换个完全不同的方式。让我们尝试一种贝叶斯方法。
一开始，我决定查看均值方向的估计值如何随着样本量的增加而变化。这个想法很简单：我会从一个圆形均匀的先验分布开始，并随着每一个新的数据点更新它。
我们需要定义一些概念，所以让我们开始进入数学部分。如果你不喜欢方程式，可以跳到下一部分，那里有很酷的可视化！
1. 冯·米塞斯分布
冯·米塞斯分布 p(θ∣μ,κ) 的概率密度函数为：
其中：
-
μ 是均值方向
-
κ 是集中参数（类似于正态分布中方差的倒数）
-
I₀(κ) 是第一类修正贝塞尔函数，确保分布是标准化的。
2. 先验和似然性
假设我们有：
- 先验 分布：
- 似然性 对于一个新的观察值 θₙ：
我们想使用贝叶斯定理更新我们的先验：
其中
3. 在 von Mises 形式中乘以先验和似然
两个 von Mises 分布，参数分别为 (μ1, κ1) 和 (μ2, κ2)，它们的乘积会得到另一个 von Mises 分布，具有更新后的参数。我们一步步来看：
给定：
并且
后验与乘积成正比：
使用余弦和公式的三角恒等式：
这变为：
4. 转换为极坐标形式以获得后验
最后阶段！上面的表达式是一个伪装的 von Mises 分布。我们可以将其重新写成极坐标形式，从而估计更新后的平均方向和集中度参数。
令：
现在后验表达式简化为：
让我们暂停一下，仔细看看这个简化的表达式。
- 请注意，C cos(θ)+S sin(θ) 是两个向量 (C,S) 和 (cos(θ),sin(θ)) 的点积，我们可以将其表示为：
其中 ϕ 是向量 (C,S) 和 (cos(θ),sin(θ)) 之间的角度。
2. 向量的大小是：
3. 向量 (cos(θ),sin(θ)) 和正 x 轴之间的角度就是 θ，而 (C,S) 和正 x 轴之间的角度，由定义为：
4. 那么两个向量之间的角度是：
将我们的发现代入简化后的后验表达式：
或者
其中
- kappa_post 是后验的集中度参数：
- mu_post 是后验的平均方向
哇，我们做到了！后验也是一个 von Mises 分布，具有更新后的参数 (mu_post, kappa_post)。现在，我们可以通过每个新的观测值更新先验，并观察平均方向的变化。
贝叶斯小插曲：有趣的部分
欢迎回到那些跳过数学部分的朋友们，恭喜那些完成了的朋友！现在让我们编写贝叶斯更新代码，并可视化结果。
首先，让我们定义一些辅助函数，用于可视化后验分布。稍后我们将用它来创建一个漂亮的动画。
import imageio
from io import BytesIO
def get_posterior_distribution_image_array(
mu_grid: np.ndarray,
posterior_pdf: np.ndarray,
current_samples: List[float],
idx: int,
fig_size: Tuple[int, int],
dpi: int,
r_max_posterior: float
) -> np.ndarray:
"""
Creates the posterior distribution and observed samples histogram on a polar plot,
converts it to an image array, and returns it for GIF processing.
Parameters:
-----------
mu_grid (np.ndarray):
Grid of mean direction values for plotting the posterior PDF.
posterior_pdf (np.ndarray):
Posterior probability density function values for the given `mu_grid`.
current_samples (List[float]):
List of observed angle samples in radians.
idx (int):
The current step index, used for labeling the plot.
fig_size (Tuple[int, int]):
Size of the plot figure (width, height).
dpi (int):
Dots per inch (resolution) for the plot.
r_max_posterior (float):
Maximum radius for the posterior PDF plot, used to set plot limits.
Returns:
np.ndarray: Image array of the plot in RGB format, suitable for GIF processing.
"""
fig = plt.figure(figsize=fig_size, dpi=dpi)
ax = plt.subplot(1, 1, 1, projection='polar')
ax.set_theta_zero_location('N')
ax.set_theta_direction(-1)
ax.plot(mu_grid, posterior_pdf, color='red', linewidth=2, label='Posterior PDF')
# observed samples histogram
n_bins = 48
hist_bins = np.linspace(-np.pi, np.pi, n_bins + 1)
hist_counts, _ = np.histogram(current_samples, bins=hist_bins)
# normalize the histogram counts
if np.max(hist_counts) > 0:
hist_counts_normalized = hist_counts / np.max(hist_counts)
else:
hist_counts_normalized = hist_counts
bin_centers = (hist_bins[:-1] + hist_bins[1:]) / 2
bin_width = hist_bins[1] - hist_bins[0]
# set the maximum radius to accommodate both the posterior pdf and histogram bars
r_histogram_height = r_max_posterior * 0.9
r_max = r_max_posterior + r_histogram_height
ax.set_ylim(0, r_max)
# plot the histogram bars outside the circle
for i in range(len(hist_counts_normalized)):
theta = bin_centers[i]
width = bin_width
hist_height = hist_counts_normalized[i] * r_histogram_height
if hist_counts_normalized[i] > 0:
ax.bar(
theta, hist_height, width=width, bottom=r_max_posterior,
color='teal', edgecolor='black', alpha=0.5
)
ax.text(
0.5, 1.1, f'Posterior Distribution (Step {idx + 1})',
transform=ax.transAxes, ha='center', va='bottom', fontsize=18
)
ax.set_yticklabels([])
ax.grid(linestyle='--')
ax.yaxis.set_visible(False)
ax.spines['polar'].set_visible(False)
plt.subplots_adjust(top=0.85, bottom=0.05, left=0.05, right=0.95)
# saving to buffer for gif processing
buf = BytesIO()
plt.savefig(buf, format='png', bbox_inches=None, pad_inches=0)
buf.seek(0)
img_array = plt.imread(buf)
img_array = (img_array * 255).astype(np.uint8)
plt.close(fig)
return img_array
现在我们准备写更新循环了。记住我们需要设置先验分布。我将从一个圆形均匀分布开始，它等价于一个冯·米塞斯分布，集中度参数为 0。对于 kappa_likelihood，我设置了一个固定的中等集中度参数 2。这将使得后验更新更加明显。
# initial prior parameters
mu_prior = 0.0 # initial mean direction (any value, since kappa_prior = 0)
kappa_prior = 0.0 # uniform prior over the circle
# fixed concentration parameter for the likelihood
kappa_likelihood = 2.0
posterior_mus = []
posterior_kappas = []
mu_grid = np.linspace(-np.pi, np.pi, 200)
# vizualisation parameters
fig_size = (10, 10)
dpi = 100
current_samples = []
frames = []
for idx, theta_n in enumerate(data['radians']):
# compute posterior parameters
C = kappa_prior * np.cos(mu_prior) + kappa_likelihood * np.cos(theta_n)
S = kappa_prior * np.sin(mu_prior) + kappa_likelihood * np.sin(theta_n)
kappa_post = np.sqrt(C**2 + S**2)
mu_post = np.arctan2(S, C)
# posterior distribution
posterior_pdf = np.exp(kappa_post * np.cos(mu_grid - mu_post)) / (2 * np.pi * i0(kappa_post))
# store posterior parameters and observed samples
posterior_mus.append(mu_post)
posterior_kappas.append(kappa_post)
current_samples.append(theta_n)
# plot posterior distribution
r_max_posterior = max(posterior_pdf) * 1.1
img_array = get_posterior_distribution_image_array(
mu_grid,
posterior_pdf,
current_samples,
idx,
fig_size,
dpi,
r_max_posterior
)
frames.append(img_array)
# updating priors for next iteration
mu_prior = mu_post
kappa_prior = kappa_post
# Create GIF
fps = 10
frames.extend([img_array]*fps*3) # repeat last frame a few times to make a "pause" at the end of the GIF
imageio.mimsave('../images/posterior_updates.gif', frames, fps=fps)
就是这样！该代码将生成一个 GIF，展示每次新观测后，后验分布的更新。这里是辉煌的结果：
后验分布更新（图像由作者提供）
随着每次新的观测，后验分布变得越来越集中于真实的均值方向。如果我能用 Auri 的轮廓替代红线，那就完美了！
我们可以进一步可视化后验均值方向和集中度参数的历史。让我们绘制它们：
# Convert posterior_mus to degrees
posterior_mus_deg = np.rad2deg(posterior_mus) % 360
n_samples = data.shape[0]
true_mu = data['degrees'].mean()
# Plot evolution of posterior mean direction
fig, ax1 = plt.subplots(figsize=(12, 6))
color = 'tab:blue'
ax1.set_xlabel('Number of Observations')
ax1.set_ylabel('Posterior Mean Direction (Degrees)', color=color)
ax1.plot(range(1, n_samples + 1), posterior_mus_deg, marker='o', color=color)
ax1.tick_params(axis='y', labelcolor=color)
ax1.axhline(true_mu, color='red', linestyle='--', label='Sample Distribution Mean Direction')
ax1.legend(loc='upper left')
ax1.grid(True)
ax2 = ax1.twinx() # instantiate a second axes that shares the same x-axis
color = 'tab:orange'
ax2.set_ylabel('Posterior Concentration Parameter (kappa)', color=color) # we already handled the x-label with ax1
ax2.plot(range(1, n_samples + 1), posterior_kappas, marker='o', color=color)
ax2.tick_params(axis='y', labelcolor=color)
fig.tight_layout() # otherwise the right y-label is slightly clipped
sns.despine()
plt.title('Evolution of Posterior Mean Direction and Concentration Over Time')
plt.show()
后验均值、kappa 演变（图像由作者提供）
图表展示了后验均值方向和集中度参数随着每次新观测如何变化。均值方向最终收敛到样本值，而集中度参数随着估计的确定性增加而上升。
贝叶斯因子：PyMC 用于狗狗指南针
我最后想尝试的是使用贝叶斯因子方法进行假设检验。贝叶斯因子背后的思想非常简单：它是两个竞争假设/模型的边际似然比。
通常，贝叶斯因子定义为：
其中：
-
p(D∣Mi) 和 p(D∣Mj) 是在 i 和 j 假设下的数据边际似然
-
p(Mi∣D) 和 p(Mj∣D) 是给定数据后的模型后验概率
-
p(Mi) 和 p(Mj) 是模型的先验概率
结果是一个数值，告诉我们一个假设比另一个假设更有可能。解释贝叶斯因子的方式有很多种，其中一种常见的方法是使用哈罗德·杰弗里斯的杰弗里斯尺度：
你可能会问，模型是什么？很简单！它们是具有不同参数的分布。我将使用 PyMC 来定义模型，并从中采样后验分布。
首先，我们重新引入零假设。我仍然假设它是一个圆形均匀的冯·米塞斯分布，且kappa=0，但这次我们需要计算在此假设下数据的似然性。为了简化后续的计算，我们将计算对数似然。
# Calculate log likelihood for H0
log_likelihood_h0 = vonmises.logpdf(data['radians'], kappa=0, loc=0).sum()
接下来，是时候构建备择模型了。首先从一个简单的场景开始：单峰南方方向，在这个场景下，我假设分布集中在 180°或π弧度处。
单峰南方
让我们在 PyMC 中定义模型。我们将使用冯·米塞斯分布，并设置固定位置参数μ=π，同时为非负浓度参数κ设置半正态先验。这使得模型能够从数据中学习浓度参数，并检查南方方向是否更受偏好。
import pymc as pm
import arviz as az
import arviz.data.inference_data as InferenceData
from scipy.stats import halfnorm, gaussian_kde
with pm.Model() as model_uni:
# Prior for kappa
kappa = pm.HalfNormal('kappa', sigma=10)
# Likelihood
likelihood_h1 = pm.VonMises('angles', mu=np.pi, kappa=kappa, observed=data['radians'])
# Sample from posterior
trace_uni = pm.sample(
10000, tune=3000, chains=4,
return_inferencedata=True,
idata_kwargs={'log_likelihood': True})
这给我们提供了一个简单的模型，我们也可以将其可视化：
# Model graph
pm.model_to_graphviz(model_uni)
PyMC 模型图（作者提供的图片）
这是浓度参数κ的后验分布：
az.plot_posterior(trace_uni, var_names=['kappa'])
plt.show()
后验 kappa 分布（作者提供的图片）
剩下的就是计算备择模型的对数似然值和贝叶斯因子。
# Posterior samples for kappa
kappa_samples = trace_uni.posterior.kappa.values.flatten()
# Log likelihood for each sample
log_likes = []
for k in kappa_samples:
# Von Mises log likelihood
log_like = vonmises.logpdf(data['radians'], k, loc=np.pi).sum()
log_likes.append(log_like)
# Log-mean-exp trick for numerical stability
log_likelihood_h1 = np.max(log_likes) +\
np.log(np.mean(np.exp(log_likes - np.max(log_likes))))
BF = np.exp(log_likelihood_h1 - log_likelihood_h0)
print(f"Bayes Factor: {BF:.4f}")
print(f"Probability kappa > 0.5: {np.mean(kappa_samples > 0.5):.4f}")
>> Bayes Factor: 32.4645
>> Probability kappa > 0.5: 0.0649
因为我们是将备择模型的似然除以原假设模型的似然，所以贝叶斯因子表明数据在备择假设下的可能性增加了多少。在这种情况下，我们得到了 32.46，这是非常强的证据，表明数据不是均匀分布在圆周上，而是偏向南方方向。
然而，我们还计算了浓度参数kappa大于 0.5 的概率。这是一种简单的方法，用来检查分布是否显著不同于均匀分布。在单峰南方模型下，这个概率只有 0.0649，意味着分布仍然相当分散。
让我们尝试另一个模型：双峰南北混合模型。
双峰南北混合模型
这次我假设分布是双峰的，峰值分别位于 0°和 180°，正如我们在罗盘玫瑰图上看到的那样。
为了实现这一点，我需要使用两个具有不同固定均值方向和共享浓度参数的冯·米塞斯分布的混合。
首先，让我们定义一些辅助函数：
# Type aliases
ArrayLike = Union[np.ndarray, pd.Series]
ResultDict = Dict[str, Union[float, InferenceData.InferenceData]]
def compute_mixture_vonmises_logpdf(
series: ArrayLike,
kappa: float,
weights: npt.NDArray[np.float64],
mus: List[float]
) -> float:
"""
Compute log PDF for a mixture of von Mises distributions
Parameters:
-----------
series: ArrayLike
Array of observed angles in radians
kappa: float
Concentration parameter
weights: npt.NDArray[np.float64],
Array of mixture weights
mus: List[float]
Array of means for each component
Returns:
--------
float: Sum of log probabilities for all data points
"""
mixture_pdf = np.zeros_like(series)
for w, mu in zip(weights, mus):
mixture_pdf += w * vonmises.pdf(series, kappa, loc=mu)
return np.log(np.maximum(mixture_pdf, 1e-300)).sum()
def compute_log_likelihoods(
trace: az.InferenceData,
series: ArrayLike,
mus: List[float]
) -> np.ndarray:
"""
Compute log likelihoods for each sample in the trace
Parameters:
-----------
trace: az.InferenceData
The trace from the PyMC3 model sampling.
series: ArrayLike
Array of observed angles in radians
"""
kappa_samples = trace.posterior.kappa.values.flatten()
weights_samples = trace.posterior.weights.values.reshape(-1, 2)
# Calculate log likelihood for each posterior sample
log_likes = []
for k, w in zip(kappa_samples, weights_samples):
log_like = compute_mixture_vonmises_logpdf(
series,
kappa=k,
weights=w,
mus=mus
)
log_likes.append(log_like)
# Calculate marginal likelihood using log-sum-exp trick
log_likelihood_h1 = np.max(log_likes) + np.log(np.mean(np.exp(log_likes - np.max(log_likes))))
return log_likelihood_h1
def posterior_report(
log_likelihood_h0: float,
log_likelihood_h1: float,
kappa_samples: ArrayLike,
kappa_threshold: float = 0.5
) -> str:
"""
Generate a report with Bayes Factor and probability kappa > threshold
Parameters:
-----------
log_likelihood_h0: float
Log likelihood for the null hypothesis
log_likelihood_h1: float
Log likelihood for the alternative hypothesis
kappa_samples: ArrayLike
Flattened posterior samples of the concentration parameter
kappa_threshold: float
Threshold for computing the probability that kappa > threshold
Returns:
--------
summary: str
A formatted string containing the summary statistics.
"""
BF = np.exp(log_likelihood_h1 - log_likelihood_h0)
summary = (
f"Bayes Factor: {BF:.4f}\n"
f"Probability kappa > {kappa_threshold}: {np.mean(kappa_samples > kappa_threshold):.4f}"
)
return summary
现在回到模型：
mu1 = 0 # 0 degrees
mu2 = np.pi # 180 degrees
with pm.Model() as model_mixture_bimodal_NS:
# Priors for concentration parameters
kappa = pm.HalfNormal('kappa', sigma=10)
# Priors for component weights
weights = pm.Dirichlet('weights', a=np.ones(2))
# Define the von Mises components
vm1 = pm.VonMises.dist(mu=mu1, kappa=kappa)
vm2 = pm.VonMises.dist(mu=mu2, kappa=kappa)
# Mixture distribution
likelihood = pm.Mixture(
'angles',
w=weights,
comp_dists=[vm1, vm2],
observed=data['radians']
)
# Sample from the posterior
trace_mixture_bimodal_NS = pm.sample(
10000, tune=3000, chains=4, return_inferencedata=True, idata_kwargs={'log_likelihood': True})
# Get kappa samples
kappa_samples = trace_mixture_bimodal_NS.posterior.kappa.values.flatten()
再次，让我们可视化模型图和浓度参数κ的后验分布：
# Model graph
pm.model_to_graphviz(model_mixture_bimodal_NS)
PyMC 模型图（作者提供的图片）
# Posterior Analysis
az.plot_posterior(trace_mixture_bimodal_NS, var_names=['kappa'])
plt.show()
后验 kappa 分布（作者提供的图片）
最后，让我们计算贝叶斯因子和浓度参数κ大于 0.5 的概率：
log_likelihood_h1 = compute_log_likelihoods(trace_mixture_bimodal_NS, data['radians'], [mu1, mu2])
print(posterior_report(log_likelihood_h0, log_likelihood_h1, kappa_samples))
>> Bayes Factor: 214.2333
>> Probability kappa > 0.5: 0.9110
太棒了！ 我们的两个指标都表明这个模型更适合数据。贝叶斯因子表明有决定性证据，并且大多数后验κ样本大于 0.5，均值为 0.99，正如我们在分布图上看到的那样。
在结束之前，让我们再试试其他几个模型。
双峰西南混合模型
这个模型再次假设一个双峰分布，但这次峰值位于 270°和 180°，这些方向在罗盘玫瑰图中较为常见。
mu1 = np.pi # 180 degrees
mu2 = 3 * np.pi / 2 # 270 degrees
with pm.Model() as model_mixture_bimodal_WS:
# Priors for concentration parameters
kappa = pm.HalfNormal('kappa', sigma=10)
# Priors for component weights
weights = pm.Dirichlet('weights', a=np.ones(2))
# Define the four von Mises components
vm1 = pm.VonMises.dist(mu=mu1, kappa=kappa)
vm2 = pm.VonMises.dist(mu=mu2, kappa=kappa)
# Mixture distribution
likelihood = pm.Mixture(
'angles',
w=weights,
comp_dists=[vm1, vm2],
observed=data['radians']
)
# Sample from the posterior
trace_mixture_bimodal_WS = pm.sample(
10000, tune=3000, chains=4, return_inferencedata=True, idata_kwargs={'log_likelihood': True})
# Get kappa samples
kappa_samples = trace_mixture_bimodal_WS.posterior.kappa.values.flatten()
# Posterior Analysis
az.plot_posterior(trace_mixture_bimodal_WS, var_names=['kappa'])
plt.show()
log_likelihood_h1 = compute_log_likelihoods(trace_mixture_bimodal_WS, data['radians'], [mu1, mu2])
print(posterior_report(log_likelihood_h0, log_likelihood_h1, kappa_samples))
>> Bayes Factor: 20.2361
>> Probability kappa > 0.5: 0.1329
后验 kappa 分布（作者提供的图片）
不，明显不如之前的模型好。下一个！
四态混合模型
最后一轮。也许我的狗确实喜欢与基准方向对齐？让我们尝试一个四态分布，峰值分别位于 0°、90°、180°和 270°。
mu1 = 0 # 0 degrees
mu2 = np.pi / 2 # 90 degrees
mu3 = np.pi # 180 degrees
mu4 = 3 * np.pi / 2 # 270 degrees
with pm.Model() as model_mixture_quad:
# Priors for concentration parameters
kappa = pm.HalfNormal('kappa', sigma=10)
# Priors for component weights
weights = pm.Dirichlet('weights', a=np.ones(4))
# Define the four von Mises components
vm1 = pm.VonMises.dist(mu=mu1, kappa=kappa)
vm2 = pm.VonMises.dist(mu=mu2, kappa=kappa)
vm3 = pm.VonMises.dist(mu=mu3, kappa=kappa)
vm4 = pm.VonMises.dist(mu=mu4, kappa=kappa)
# Mixture distribution
likelihood = pm.Mixture(
'angles',
w=weights,
comp_dists=[vm1, vm2, vm3, vm4],
observed=data['radians']
)
# Sample from the posterior
trace_mixture_quad = pm.sample(
10000, tune=3000, chains=4, return_inferencedata=True, idata_kwargs={'log_likelihood': True}
)
# Get kappa samples
kappa_samples = trace_mixture_quad.posterior.kappa.values.flatten()
# Posterior Analysis
az.plot_posterior(trace_mixture_quad, var_names=['kappa'])
plt.show()
log_likelihood_h1 = compute_log_likelihoods(trace_mixture_quad, data['radians'], [mu1, mu2, mu3, mu4])
print(posterior_report(log_likelihood_h0, log_likelihood_h1, kappa_samples))
>> Bayes Factor: 0.0000
>> Probability kappa > 0.5: 0.9644
后验 kappa 分布（作者提供的图片）
嗯… 其实并不是。尽管集中参数κκ大于 0.5 的概率相当高，但贝叶斯因子却是 0.0。
贝叶斯因子的优点在于它有效地惩罚过度复杂的模型，有效防止过拟合。
模型比较
让我们用信息准则总结所有模型的结果。我们将使用Widely Applicable Information Criterion（WAIC）和 Leave-One-Out Cross-Validation（LOO）来比较这些模型。
# Compute WAIC for each model
wail_uni = az.waic(trace_uni)
waic_quad = az.waic(trace_mixture_quad)
waic_bimodal_NS = az.waic(trace_mixture_bimodal_NS)
waic_bimodal_WS = az.waic(trace_mixture_bimodal_WS)
model_dict = {
'Quadrimodal Model': trace_mixture_quad,
'Bimodal Model (NS)': trace_mixture_bimodal_NS,
'Bimodal Model (WS)': trace_mixture_bimodal_WS,
'Unimodal Model': trace_uni
}
# Compare models using WAIC
waic_comparison = az.compare(model_dict, ic='waic')
waic_comparison
# Compare models using LOO
loo_comparison = az.compare(model_dict, ic='loo')
loo_comparison
# Visualize the comparison
az.plot_compare(waic_comparison)
plt.show()
WAIC 比较（作者提供的图片）
我们找到了优胜者！根据 WAIC 和 LOO，双模态南北模型是数据的最佳拟合模型。
结论
Christmas Auri（作者提供的图片）
何等的旅程！几个月前仅仅是对我狗拉屎习惯的简单观察，现在却变成了全面的贝叶斯分析。
在本文中，我展示了如何建模圆形数据，估计平均方向和集中参数，并通过新观察更新后验分布。我们还看到如何使用贝叶斯因子进行假设检验，并使用信息准则比较模型。
结果非常有趣！Auri 确实有自己的喜好，并且在南北轴上能够对齐。如果我和我的狗迷失在树林中，我知道该往哪个方向走了。只需足够大的样本量来确认！
希望您像我一样享受这段旅程。如果您有任何问题或建议，请随时联系。如果您想支持我的工作，请考虑给我买杯咖啡 ❤️
https://www.buymeacoffee.com/datawondering
我的社交媒体账号：
参考文献：
扫一扫
1380
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
