生物医学网络科学Biomedical Network Science
Example03 / 10

用最大似然估计度指数，并与最小二乘比较

所属专题
小世界与无标度网络模型
下载
Notebook · .py

对 BA 网络（理论 γ=3）的度序列，用离散数据常用的最大似然近似 γ = 1 + n / Σ ln(k_i / (k_min - 0.5)) （只用 k_i >= k_min 的 n 个节点） 依次取 k_min = 2, 4, 6, 8；再与“直方图 + 双对数最小二乘”（不推荐的做法）比较。 输出：k_min 增大时 MLE 逐渐接近 3（BA 的度分布只在大 k 处才趋于 k^-3），但所用节点变少、方差变大； 最小二乘结果明显偏离 3，因为它被尾部大量“只出现一次”的大度值拉偏。 Clauset 等（2009）选择使 KS 距离最小的 k_min 来权衡这种偏差与方差。

Output输出图

用最大似然估计度指数，并与最小二乘比较：输出图 1

Output文字输出

k_min=2: 使用 3000 个节点, MLE γ = 2.40
k_min=4: 使用  904 个节点, MLE γ = 2.67
k_min=6: 使用  426 个节点, MLE γ = 2.75
k_min=8: 使用  250 个节点, MLE γ = 2.82
直方图最小二乘估计 γ = 1.92

Code代码

# 对 BA 网络（理论 γ=3）的度序列，用离散数据常用的最大似然近似
#     γ = 1 + n / Σ ln(k_i / (k_min - 0.5))      （只用 k_i >= k_min 的 n 个节点）
# 依次取 k_min = 2, 4, 6, 8；再与“直方图 + 双对数最小二乘”（不推荐的做法）比较。
# 输出：k_min 增大时 MLE 逐渐接近 3（BA 的度分布只在大 k 处才趋于 k^-3），但所用节点变少、方差变大；
#       最小二乘结果明显偏离 3，因为它被尾部大量“只出现一次”的大度值拉偏。
# Clauset 等（2009）选择使 KS 距离最小的 k_min 来权衡这种偏差与方差。
import networkx as nx
import numpy as np

G = nx.barabasi_albert_graph(3000, 2, seed=1)
k = np.array([d for _, d in G.degree()], dtype=float)


def gamma_mle(k, k_min):
    kk = k[k >= k_min]
    return 1 + len(kk) / np.sum(np.log(kk / (k_min - 0.5)))


for k_min in [2, 4, 6, 8]:
    print("k_min=%d: 使用 %4d 个节点, MLE γ = %.2f" % (k_min, (k >= k_min).sum(), gamma_mle(k, k_min)))

vals, counts = np.unique(k, return_counts=True)
pk = counts / counts.sum()
slope, intercept = np.polyfit(np.log(vals), np.log(pk), 1)
print("直方图最小二乘估计 γ = %.2f" % (-slope))

# 作图：左为直方图 + 最小二乘直线；右为 CCDF 与不同 k_min 的 MLE 拟合
import matplotlib.pyplot as plt

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
ax1.loglog(vals, pk, "o", ms=4, color="#4393c3", label="histogram P(k)")
ax1.loglog(vals, np.exp(intercept) * vals ** slope, "-", color="#d6604d", label="least squares: gamma = %.2f" % -slope)
ax1.set_xlabel("k"); ax1.set_ylabel("P(k)"); ax1.legend(fontsize=8)
ax1.set_title("Tail of single counts biases least squares", fontsize=10)
ks = np.sort(np.unique(k))
ccdf = np.array([(k >= x).mean() for x in ks])
ax2.loglog(ks, ccdf, "o", ms=4, color="0.4", label="CCDF P(K >= k)")
for k_min, ls in [(2, "--"), (8, "-")]:
    g = gamma_mle(k, k_min)
    xx = ks[ks >= k_min]
    ax2.loglog(xx, (k >= k_min).mean() * (xx / k_min) ** (1 - g), ls, lw=2, label="MLE k_min=%d: gamma = %.2f" % (k_min, g))
ax2.set_xlabel("k"); ax2.set_ylabel("P(K >= k)"); ax2.legend(fontsize=8)
ax2.set_title("BA network (theory gamma = 3)", fontsize=10)
plt.tight_layout(); plt.show()

Run在浏览器中运行

首次加载 Python 环境约需 20–40 秒；可以修改参数后重新运行。

用最大似然估计度指数，并与最小二乘比较