使用Scipy curve_fit拟合幂律分布结果异常的问题排查

你的curve_fit拟合结果异常，主要是以下几个操作失误导致的：

1. Bounds参数格式完全错误 #

你写的bounds = ([2, 3])不符合curve_fit的要求。这个参数需要传入两个独立的数组：第一个是所有参数的下限，第二个是所有参数的上限。因为你只拟合alpha这一个参数，正确写法应该是bounds=([2], [3])。之前的写法会被解析成“第一个参数下限2，第二个参数下限3”，但你的函数只有一个参数，导致约束逻辑混乱，最终alpha被强制卡在下限2。

2. 拟合数据的权重逻辑不匹配powerlaw库 #

powerlaw库拟合时，会自动用每个度数的出现频次作为权重（相当于把原始度数列表里的每个节点都纳入计算），但你用curve_fit时直接传了unique_degrees和counts_pdf——这相当于把每个不同度数当成权重相同的独立点，而社交网络里小度数节点的数量远多于大度数节点，这种处理会让小度数点完全主导拟合结果，最终得到偏小的alpha值。

解决办法：要么给curve_fit传入权重（用频次的平方根作为sigma，对应泊松误差），要么直接基于原始度数数据拟合。

3. 幂律函数的形式与拟合目标不匹配 #

powerlaw库默认拟合的是互补累积分布函数(CCDF)，或者会对离散度数的概率质量函数(PMF)做正确归一化。你自定义的函数是连续幂律的PDF，直接用来拟合离散度数的PMF会有系统偏差，尤其是大度数区域。更稳定的做法是拟合CCDF，因为它在对数坐标下更接近线性，拟合难度更低。

修正后的Curve_fit代码（匹配powerlaw逻辑） #

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

# 第一步：确保数据和powerlaw拟合时一致，过滤xmin=1及以上的度数
filtered_degree = [d for d in degree if d >= 1]
unique_degrees, counts = np.unique(filtered_degree, return_counts=True)
counts_pdf = counts / len(filtered_degree)  # 归一化为PMF

# 自定义离散幂律PMF函数（xmin=1时的归一化形式）
def powerlaw_pmf(x, alpha):
    # 计算离散归一化因子（近似，当x范围大时也可以用黎曼和）
    norm_factor = np.sum(np.power(np.arange(1, max(unique_degrees)+1), -alpha))
    return np.power(x, -alpha) / norm_factor

# 正确的bounds格式
bounds = ([2], [3])
p0 = [2.5]

# 传入频次平方根作为sigma，给大度数点更高的权重（小度数计数误差小，权重低）
params, _ = curve_fit(powerlaw_pmf, unique_degrees, counts_pdf, 
                      p0=p0, bounds=bounds, sigma=np.sqrt(counts)/len(filtered_degree))

print("alpha = ", params[0])

# 绘图对比
x = np.linspace(min(unique_degrees), max(unique_degrees), 1000)
y1 = powerlaw_pmf(x, *params)

plt.loglog(unique_degrees, counts_pdf, 'o', label='Degree Distribution')
plt.loglog(x, y1, c='g', label='Fitted Power Law')
plt.legend()
plt.show()

更稳定的方案：拟合CCDF #

# 计算互补累积分布函数CCDF：P(X >= x)
sorted_degrees = np.sort(filtered_degree)
ccdf = 1 - np.arange(len(sorted_degrees)) / len(sorted_degrees)
# 提取唯一度数对应的CCDF值
unique_ccdf_x = np.unique(sorted_degrees)
unique_ccdf_y = [ccdf[sorted_degrees == x][0] for x in unique_ccdf_x]

# xmin=1时的幂律CCDF形式
def powerlaw_ccdf(x, alpha):
    return np.power(x / 1, -(alpha - 1))

params, _ = curve_fit(powerlaw_ccdf, unique_ccdf_x, unique_ccdf_y, 
                      p0=p0, bounds=bounds)

print("alpha = ", params[0])

# 绘图
x = np.linspace(min(unique_ccdf_x), max(unique_ccdf_x), 1000)
y1 = powerlaw_ccdf(x, *params)

plt.loglog(unique_ccdf_x, unique_ccdf_y, 'o', label='Degree CCDF')
plt.loglog(x, y1, c='g', label='Fitted Power Law CCDF')
plt.legend()
plt.show()

内容的提问来源于stack exchange，提问作者Alexandre Bloch