用Python在现实世界数据中检测幂律

使用Python在现实世界数据中检测幂律律定律

用示例代码分解最大似然方法

这是关于幂律和重尾的系列文章中的第二篇。在之前的文章中，我给出了一个用户友好的介绍，并提出了在分析幂律时我们标准统计工具的三个问题。尽管意识可以帮助我们避免这些问题，但在实践中往往不明确某些给定数据遵循什么分布。在本文中，我将描述如何客观地检测现实世界数据中的幂律，并通过社交媒体数据共享一个具体的示例。

注意：如果您对幂律分布或重尾这类术语不熟悉，请先阅读本系列的第一篇文章。

幂律打破统计学101

在上一篇文章中，我们专注于两种类型的分布：高斯分布和幂律分布。我们发现这些分布具有完全相反的统计特性。即，幂律由罕见事件驱动，而高斯分布则不是。

这种罕见事件驱动的特性在使用许多常用统计工具（如均值、标准偏差、回归等）分析幂律时引发了一些问题。重要的是，如果数据类似于高斯分布，则可以使用常见的方法如回归和计算期望值而不必担心。然而，如果数据更像是幂律分布，这些技术可能会给出错误和误导性的结果。

• 通过自动化的工单分类提高客户支持效率
• 谷歌的BARD现在可以“观看并回答有关YouTube视频的问题”
• 这篇AI论文介绍了一种开创性的机器学习模型，用于高效预测氢燃烧：利用负设计和反应化学中的元动力学

我们还看到了一种既类似高斯分布又类似幂律分布（尽管它们具有相反的属性）的第三个（更捣蛋的）分布，称为对数正态分布。

这种不确定性对于决定如何分析给定的数据集的从业者提出了挑战。为了帮助克服这些难题，确定数据是否符合幂律、对数正态分布或其他类型的分布可能是有优势的。

对数-对数方法

拟合实际数据的幂律的一种流行方法是我称之为“对数-对数方法” [1]。这个想法源于取幂律概率密度函数的对数，如下所示。

上述推导将幂律的概率密度函数定义转化为线性方程，如下图所示。

这意味着遵循幂律的数据的直方图将遵循一条直线。在实践中，这看起来就像是为一些数据生成直方图并在对数对数图中绘制它。甚至可以进一步进行线性回归，以估计分布的α值（这里，α = -m+1）。

然而，这种方法存在着显著的局限性。这些局限性在参考文献[1]中被描述并总结如下：

• 斜率（因此α）的估计存在系统误差
• 回归误差难以估计
• 即使分布不遵循幂律，拟合看起来也可能很好
• 拟合可能不满足概率分布的基本条件，例如归一化

最大似然方法

虽然对数对数方法很容易实现，但其局限性使其不够理想。相反，我们可以转向更数学上严谨的方法，即最大似然。最大似然是一种广泛应用的统计方法，用于根据一些数据推断出模型的最佳参数。

最大似然包括两个关键步骤。步骤1：获取似然函数。步骤2：最大化似然以获得模型参数。

步骤1：编写似然函数

似然是一种特殊类型的概率。简而言之，它量化了在给定特定模型的情况下，我们的数据的概率。我们可以将其表示为我们观察到的所有数据的联合概率。对于 Pareto 分布，我们可以写成如下形式。

为了使最大化似然变得更容易，习惯上使用对数似然（它们由相同的 α 值最大化）。

步骤2：最大化似然

有了（对数）似然函数，我们现在可以将确定最佳参数选择的任务框架化为一个优化问题。为了基于我们的数据找到最优的 α 值，这归结为将 l(α) 关于 α 的导数设为零，然后解出 α。以下是其推导。

这为我们提供了所谓的最大似然估计 α 的估计值。通过这个估计值，我们可以将观察到的 x 值带入来估计 Pareto 分布的 α 值。

有了理论基础，让我们看看当应用于真实世界数据时会是什么样子（来自我的社交媒体账户）。

示例代码：将幂律拟合到社交媒体数据上

社交媒体是一个普遍存在长尾数据的领域。例如，部分创作者得到了大部分的关注，只有少数 VoAGI 博客获得了大部分的阅读量，等等。

在这里，我们将使用 Python 的 powerlaw 库来确定我各种社交媒体渠道（即 VoAGI、YouTube、LinkedIn）的数据是否真正遵循幂律分布。这些示例的数据和代码可在 GitHub 仓库 上获得。

YouTube-Blog/power-laws/2-detecting-powerlaws at main · ShawhinT/YouTube-Blog

用于补充 VoAGI 的 YouTube 视频和博文的代码。- YouTube-Blog/power-laws/2-detecting-powerlaws at main ·…

github.com

人造数据

在将基于最大似然方法应用于来自真实世界的混乱数据之前，让我们看看当我们将这种技术应用于分别从 Pareto 分布和对数正态分布生成的（真正的）人造数据时会发生什么。这将有助于我们在使用该方法处理我们不知道“真实”潜在分布类型的数据之前对期望有所了解。

首先，我们导入一些有用的库。

import numpy as npimport matplotlib.pyplot as pltimport powerlawimport pandas as pdnp.random.seed(0)

接下来，让我们从 Pareto 分布和对数正态分布生成一些数据。

# 幂律数据a = 2x_min = 1n = 1_000x = np.linspace(0, n, n+1)s_pareto = (np.random.pareto(a, len(x)) + 1) * x_min# 对数正态数据m = 10s = 1s_lognormal = np.random.lognormal(m, s, len(x)) * s * np.sqrt(2*np.pi)

为了了解这些数据的样子，绘制直方图是很有帮助的。在这里，我绘制了每个样本原始值和原始值的对数的直方图。后者的分布使得在视觉上更容易区分幂律和对数正态数据。

从上面的直方图中可以看出，对于两个分布的原始值，其分布在质的上是相似的。然而，我们可以看到对数分布中存在明显的差异。也就是说，对数幂律分布高度倾斜且不以平均值为中心，而对数对数正态分布则类似于高斯分布。

现在，我们可以使用 powerlaw 库来将每个样本拟合为幂律分布，并估计α和x_min的值。以下是我们对 Power Law 样本的拟合结果。

# 对幂律数据进行拟合results = powerlaw.Fit(s_pareto)# 打印结果print("alpha = " + str(results.power_law.alpha)) # 注意：powerlaw 库的 alpha 定义与标准定义不同，即 a_powerlawlib = a_standard + 1print("x_min = " + str(results.power_law.xmin))print('p = ' + str(compute_power_law_p_val(results)))# 计算幂律拟合的最佳最小值# alpha = 2.9331912195958676# x_min = 1.2703447024073973# p = 0.999

通过上面打印的 alpha 和 x_min 值，我们可以看到拟合对真实参数值（即 a=3，x_min=1）的估计相当准确。数值 p 表示了拟合的质量，p 值越高意味着拟合质量越好（关于该值的更多详细信息可参考第 4.1 节的参考文献 [1]）。

我们可以对对数正态分布进行类似的分析。

# 将幂律分布拟合到对数正态分布数据results = powerlaw.Fit(s_lognormal)print("alpha = " + str(results.power_law.alpha)) # 注意：powerlaw 库中的 alpha 定义与标准定义不同，即 a_powerlawlib = a_standard + 1print("x_min = " + str(results.power_law.xmin))print('p = ' + str(compute_power_law_p_val(results)))# 计算幂律拟合的最佳最小值# alpha = 2.5508694755027337# x_min = 76574.4701482522# p = 0.999

我们可以看到，对数正态样本也很好地符合幂律分布（p=0.999）。然而，请注意 x_min 值远在尾部。虽然这在某些用例中可能有所帮助，但它并没有告诉我们最适合拟合样本中所有数据的分布。

为了解决这个问题，我们可以手动将 x_min 值设为样本最小值，然后重新进行拟合。

# 修正 xmin，这样拟合必须包含所有数据results = powerlaw.Fit(s_lognormal, xmin=np.min(s_lognormal))print("alpha = " + str(results.power_law.alpha))print("x_min = " + str(results.power_law.xmin))# alpha = 1.3087955873576855# x_min = 2201.318351239509

.Fit() 方法还会自动生成对于对数正态分布的估计值。

print("mu = " + str(results.lognormal.mu))print("sigma = " + str(results.lognormal.sigma))# mu = 10.933481999687547# sigma = 0.9834599169175509

估计的对数正态参数值接近实际值（mu=10，sigma=1），因此拟合又一次完成得很好！

然而，通过修正 x_min，我们失去了质量度量 p（由于某种原因，当提供 x_min 时，该方法不会生成 p 值）。这就引出了一个问题，我应该选择哪个分布参数？幂律分布还是对数正态分布？

为了回答这个问题，我们可以通过对数似然比（R）将幂律分布与其他候选分布进行比较。正R表示幂律分布更适合，而负R表示候选分布更适合。此外，每个比较还给出了一个显著性值（p）。下面的代码块演示了这一点。

distribution_list = ['lognormal', 'exponential', 'truncated_power_law', \      'stretched_exponential', 'lognormal_positive']    for distribution in distribution_list:    R, p = results.distribution_compare('power_law', distribution)    print("power law vs " + distribution +           ": R = " + str(np.round(R,3)) +           ", p = " + str(np.round(p,3)))# power law vs lognormal: R = -776.987, p = 0.0# power law vs exponential: R = -737.24, p = 0.0# power law vs truncated_power_law: R = -419.958, p = 0.0# power law vs stretched_exponential: R = -737.289, p = 0.0# power law vs lognormal_positive: R = -776.987, p = 0.0

如上所示，在包含所有数据的对数正态样本中，每种替代分布都优于幂律分布。此外，根据似然比，对数正态分布和正值对数正态分布的拟合效果最好。

现实世界数据

现在，我们将将幂律分布库应用于我们不知道真实分布的数据上。

我们将按照与前面类似的过程，但使用真实世界的数据。这里，我们将分析以下数据：我 VoAGI 个人资料上的每月新增关注者，我所有 YouTube 视频的收益，以及我 LinkedIn 帖子的每日印象过去一年的数据。

我们将从绘制直方图开始。

从这些图中，有两个问题引起了我的注意。首先，这三个图都更像是对数正态直方图，而不是我们之前看到的幂律直方图。其次，VoAGI和YouTube的分布稀疏，这意味着它们可能没有足够的数据来得出强有力的结论。

接下来，我们将对所有三个分布应用幂律拟合，将x_min设置为每个样本中的最小值。以下是打印出的结果。

为了确定哪个分布最好，我们可以再次将幂律拟合与其他分布进行对比。以下是这些结果。

使用统计学中常用的显著性水平p<0.1，我们可以得出以下结论。VoAGI的追随者和LinkedIn的印象最适合对数正态分布，而幂律最能代表YouTube的收益。

当然，由于VoAGI的追随者和YouTube的收益数据有限（N<100），我们应该对这些数据的任何结论持保留态度。

结论

许多传统的统计工具在应用于满足幂律分布的数据时会出现故障。因此，检测实证数据中的幂律分布可以帮助从业者避免错误分析和误导性的结论。

然而，幂律是一种更一般的叫做“尾巴肥胖”的现象的特例。在本系列的下一篇文章中，我们将进一步量化任何给定数据集的尾巴肥胖程度，提供4个方便的经验法则。

帕累托、幂律和尾巴肥胖

统计学没有教给你的知识

towardsdatascience.com

资源

联系方式：我的网站 | 预约通话 | 问我任何问题

社交媒体：YouTube 🎥 | LinkedIn | Twitter

支持: 给我买杯咖啡 ☕️

数据创业者

为数据领域的创业者提供一个社区。👉加入Discord！

VoAGI.com

[1] arXiv:0706.1062 [physics.data-an]

[2] arXiv:2001.10488 [stat.OT]

[3] https://en.wikipedia.org/wiki/Likelihood_function

[4] https://en.wikipedia.org/wiki/Pareto_distribution