幂律度分布指数计算：powerlaw库遇0值的替代求解方案

尝试使用powerlaw库拟合度分布曲线，但数据中存在0值导致拟合失败，请问是否有替代方法来计算该幂律度分布的指数？

以下是绘制度分布的代码：

import matplotlib.pyplot as plt

degree_dist = [0.0, 0.0, 0.087, 0.105, 0.088, 0.081, 0.072, 0.058, 0.04, 0.037, 0.035, 0.035, 0.021, 0.018, 0.015, 0.024,
               0.015, 0.02, 0.018, 0.01, 0.01, 0.016, 0.005, 0.005, 0.008, 0.005, 0.008, 0.008, 0.009, 0.005, 0.002,
               0.003, 0.003, 0.005, 0.003, 0.004, 0.007, 0.006, 0.002, 0.006, 0.007, 0.003, 0.004, 0.001, 0.001, 0.004,
               0.003, 0.006, 0.002, 0.006, 0.002, 0.001, 0.003, 0.003, 0.001, 0.003, 0.003, 0.002, 0.004, 0.001, 0.003,
               0.003, 0.005, 0.0, 0.001, 0.002, 0.001, 0.0, 0.002, 0.001, 0.001, 0.001, 0.0, 0.0, 0.0, 0.001, 0.001,
               0.0, 0.0, 0.001, 0.001, 0.002, 0.001, 0.001, 0.0, 0.001, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.001, 0.001, 0.0,
               0.0, 0.002, 0.0, 0.0, 0.002, 0.002, 0.002, 0.001, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
               0.0, 0.001]
degrees = range(len(degree_dist))
plt.scatter(degrees, degree_dist)
plt.xscale('log')
plt.yscale('log')
plt.xlabel('Degrees')
plt.ylabel('Distribution')
plt.title('Degree distribution')
plt.show()

该度分布呈现典型的长尾特征。

方法1：用原始节点度数据拟合（推荐） #

powerlaw库的设计逻辑是直接输入所有节点的度数列表，而非预计算的概率分布。预计算分布中的0值会干扰拟合，而原始度数据可以直接过滤无效值后使用。

代码示例 #

场景1：有原始节点度数据

import powerlaw

# 假设node_degrees是从网络中提取的所有节点度数列表
filtered_degrees = [d for d in node_degrees if d > 0]  # 过滤度数为0的孤立节点

# 拟合幂律分布
fit = powerlaw.Fit(filtered_degrees)
alpha = fit.power_law.alpha
print(f"幂律指数α: {alpha:.2f}")

# 绘制拟合结果与原始分布
fit.plot_pdf(color='b', linewidth=2, label='原始分布')
fit.power_law.plot_pdf(color='r', linestyle='--', label=f'幂律拟合 α={alpha:.2f}')
plt.xlabel('Degree')
plt.ylabel('Probability Density')
plt.legend()
plt.show()

场景2：只有预计算的分布数据

可以反向生成模拟的原始度数据：

import powerlaw
import numpy as np

# 假设总节点数为1000，可根据实际情况调整
total_nodes = 1000
node_degrees = []
for degree, prob in enumerate(degree_dist):
    if prob > 0:
        # 按概率生成对应度数的节点列表
        node_degrees.extend([degree] * int(prob * total_nodes))

# 过滤无效值后拟合
filtered_degrees = [d for d in node_degrees if d > 0]
fit = powerlaw.Fit(filtered_degrees)
print(f"幂律指数α: {fit.power_law.alpha:.2f}")

方法2：双对数坐标下的线性回归 #

幂律分布满足 p(k) ∝ k^(-α)，取双对数后转化为线性关系 log(p(k)) = -α·log(k) + C，通过线性回归拟合直线，斜率的绝对值即为幂律指数。

代码示例 #

import numpy as np
from sklearn.linear_model import LinearRegression
import matplotlib.pyplot as plt

# 过滤分布为0和度数为0的无效点
valid_k = []
valid_p = []
for k, p in zip(degrees, degree_dist):
    if p > 0 and k >= 1:
        valid_k.append(np.log10(k))
        valid_p.append(np.log10(p))

# 转换为模型所需的数组格式
X = np.array(valid_k).reshape(-1, 1)
y = np.array(valid_p)

# 线性回归拟合
model = LinearRegression()
model.fit(X, y)

alpha = -model.coef_[0]
print(f"幂律指数α: {alpha:.2f}")

# 绘制拟合结果
plt.scatter(valid_k, valid_p, label='原始数据')
plt.plot(X, model.predict(X), color='r', linestyle='--', label=f'拟合线 α={alpha:.2f}')
plt.xlabel('log₁₀(Degree)')
plt.ylabel('log₁₀(Distribution)')
plt.legend()
plt.show()

注意：此方法适合快速估算，精度略低于powerlaw库的专业拟合，但实现简单。

方法3：使用SciPy的统计函数拟合 #

SciPy的powerlaw模块也可用于拟合，需注意其参数定义与powerlaw库略有差异。

代码示例 #

from scipy.stats import powerlaw
import numpy as np

# 同方法1的过滤后原始度数据
filtered_degrees = [d for d in node_degrees if d > 0]

# 拟合分布，SciPy的形状参数a对应幂律指数α-1
params = powerlaw.fit(filtered_degrees, floc=0)
alpha = params[0] + 1
print(f"幂律指数α: {alpha:.2f}")

关键注意事项 #

1. 优先使用原始节点度数据，预计算分布容易引入离散化误差和无效0值。
2. 必须过滤度数为0的节点，幂律分布仅描述k≥k_min的区间，k=0不在拟合范围内。
3. 可通过powerlaw库的fit.xmin获取最佳拟合的最小度数阈值，仅对k≥xmin的部分拟合，结果更准确。

内容的提问来源于stack exchange，提问作者Dinh Duc Vu