生物医学网络科学Biomedical Network Science
Example10 / 15

ER 网络的最大度与最小度(自然截断)

所属专题
网络结构指标与随机图
下载
Notebook · .py

对泊松度分布,令“度 >= k_max 的期望节点数”N·P(K >= k_max) ≈ 1,“度 <= k_min 的期望节点数”≈ 1, 解出 k_max 与 k_min。用 Barabási《Network Science》3.B 节的设定 N = 10^9、 = 1000 计算, 再在可模拟的规模(N = 10^5, = 10)上与实际生成的 ER 网络比较。 输出:N=10^9 时 k_max 约 1.2×10^3、k_min 约 816,都在平均度的 ±20% 以内——随机网络没有枢纽; N=10^5、=10 时泊松估计与模拟网络的最大度一致,最小度为 0(存在孤立节点,因为 < ln N)。

Output输出图

ER 网络的最大度与最小度(自然截断):输出图 1

Output文字输出

N=1e9, <k>=1000: k_max ≈ 1196, k_min ≈ 816
N=1e5, <k>=10: 泊松估计 k_max≈27, k_min≈0;模拟 k_max=27, k_min=0

Code代码

# 对泊松度分布,令“度 >= k_max 的期望节点数”N·P(K >= k_max) ≈ 1,“度 <= k_min 的期望节点数”≈ 1,
# 解出 k_max 与 k_min。用 Barabási《Network Science》3.B 节的设定 N = 10^9、<k> = 1000 计算,
# 再在可模拟的规模(N = 10^5, <k> = 10)上与实际生成的 ER 网络比较。
# 输出:N=10^9 时 k_max 约 1.2×10^3、k_min 约 816,都在平均度的 ±20% 以内——随机网络没有枢纽;
#       N=10^5、<k>=10 时泊松估计与模拟网络的最大度一致,最小度为 0(存在孤立节点,因为 <k> < ln N)。
import networkx as nx
import numpy as np
from scipy.stats import poisson


def natural_cutoffs(N, c):
    k = np.arange(0, int(c + 50 * np.sqrt(c) + 50))
    kmax = k[np.argmin(np.abs(N * poisson.sf(k - 1, c) - 1))]   # N·P(K >= k) ≈ 1
    kmin = k[np.argmin(np.abs(N * poisson.cdf(k, c) - 1))]      # N·P(K <= k) ≈ 1
    return kmax, kmin


print("N=1e9, <k>=1000: k_max ≈ %d, k_min ≈ %d" % natural_cutoffs(1e9, 1000))
N, c = 100_000, 10
G = nx.fast_gnp_random_graph(N, c / (N - 1), seed=1)
deg = np.array([d for _, d in G.degree()])
print("N=1e5, <k>=10: 泊松估计 k_max≈%d, k_min≈%d;模拟 k_max=%d, k_min=%d"
      % (*natural_cutoffs(N, c), deg.max(), deg.min()))

# 作图:模拟度分布与泊松分布(对数纵轴),竖线为自然截断的估计值
import matplotlib.pyplot as plt

kmax_est, kmin_est = natural_cutoffs(N, c)
ks, cnt = np.unique(deg, return_counts=True)
fig, ax = plt.subplots(figsize=(6.4, 4.2))
ax.semilogy(ks, cnt / N, "o", ms=5, label="ER simulation (N=1e5)")
kk = np.arange(0, ks.max() + 3)
ax.semilogy(kk, poisson.pmf(kk, c), "-", color="0.4", label="Poisson(<k>=10)")
ax.axhline(1 / N, color="0.6", ls=":", label="1/N")
for x, lab in [(kmax_est, "k_max est."), (kmin_est, "k_min est.")]:
    ax.axvline(x, color="#d6604d", ls="--")
    ax.text(x, 0.12, lab, rotation=90, ha="right", va="top", fontsize=8, color="#d6604d")
ax.set_xlabel("degree k"); ax.set_ylabel("P(k)"); ax.set_ylim(1e-6, 0.3); ax.legend(fontsize=8)
ax.set_title("No hubs: observed k_max = %d, k_min = %d" % (deg.max(), deg.min()), fontsize=10)
plt.tight_layout(); plt.show()

Run在浏览器中运行

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

ER 网络的最大度与最小度(自然截断)