度保持随机化检验——观测到的 r 是否超出度序列本身能解释的范围

# 两个“观测网络”有相同的幂律度序列：
#   网络 A：对简单配置模型再做充分的双边交换，相当于从“给定度序列的均匀随机简单图”中抽一张
#           （只含结构性异配）；
#   网络 B：在 A 的基础上做部分异配重连（p=0.5，模拟真实的“枢纽回避”机制）。
# 对每个网络生成 100 个保度、保持简单图的随机网络（双边交换，Maslov–Sneppen），比较 r。
# 解读：A 的 r 虽为负，却落在零分布之内（结构性）；B 的 r 显著低于零分布（超出度序列的额外异配）。
# 注意：直接“合并重边、删自环”得到的简单图（erased configuration model）与均匀随机简单图不是同一个
# 分布，拿它当零模型会给出有偏的 r，所以这里先用双边交换把 A 混合充分。
import numpy as np
import networkx as nx

rng = np.random.default_rng(5)
N, gamma = 1500, 2.3
deg = np.floor(2 * (1 - rng.random(N)) ** (-1 / (gamma - 1))).astype(int)
deg = np.minimum(deg, 120)
if deg.sum() % 2:
    deg[0] += 1
A = nx.Graph(nx.configuration_model(deg.tolist(), seed=1))
A.remove_edges_from(nx.selfloop_edges(A))
nx.double_edge_swap(A, nswap=20 * A.number_of_edges(), max_tries=10**7, seed=3)


def xbs_rewire(G, n_swaps, p, seed):
    rng = np.random.default_rng(seed)
    G = G.copy()
    dg = dict(G.degree())
    edges = list(G.edges())
    done = 0
    while done < n_swaps:
        i, j = rng.choice(len(edges), 2, replace=False)
        (a, b), (c, d) = edges[i], edges[j]
        nodes = [a, b, c, d]
        if len(set(nodes)) < 4:
            continue
        if rng.random() < p:
            s = sorted(nodes, key=lambda x: dg[x], reverse=True)
            new = [(s[0], s[3]), (s[1], s[2])]
        else:
            rng.shuffle(nodes)
            new = [(nodes[0], nodes[1]), (nodes[2], nodes[3])]
        if any(G.has_edge(*e) for e in new):
            continue
        G.remove_edges_from([(a, b), (c, d)])
        G.add_edges_from(new)
        edges[i], edges[j] = new
        done += 1
    return G


B = xbs_rewire(A, A.number_of_edges(), p=0.5, seed=2)


def null_r(G, n_null=100, swaps_per_edge=5, seed=0):
    vals = []
    for s in range(n_null):
        H = G.copy()
        m = H.number_of_edges()
        nx.double_edge_swap(H, nswap=swaps_per_edge * m, max_tries=100 * swaps_per_edge * m, seed=seed + s)
        vals.append(nx.degree_assortativity_coefficient(H))
    return np.array(vals)


tests = []                                                  # 保存结果用于作图
for name, G in [("网络 A（仅结构性异配）", A), ("网络 B（额外异配重连）", B)]:
    r_obs = nx.degree_assortativity_coefficient(G)
    r0 = null_r(G, n_null=100, swaps_per_edge=5, seed=100)
    tests.append((r_obs, r0))
    z = (r_obs - r0.mean()) / r0.std(ddof=1)
    p_low = (1 + np.sum(r0 <= r_obs)) / (len(r0) + 1)       # 单侧：是否比零模型更异配
    print(f"{name}: r_obs = {r_obs:+.3f}; 零模型 r = {r0.mean():+.3f} ± {r0.std(ddof=1):.3f}; "
          f"z = {z:+.1f}; 单侧经验 p = {p_low:.3f}")

# 对照：erased configuration model（不再交换）与均匀随机简单图的 r 不同
E = nx.Graph(nx.configuration_model(deg.tolist(), seed=1))
E.remove_edges_from(nx.selfloop_edges(E))
print(f"erased 配置模型 r = {nx.degree_assortativity_coefficient(E):+.3f}，"
      f"同一度序列的保度随机网络 r = {null_r(E, n_null=20, seed=500).mean():+.3f}")

# 作图：两个网络的保度零分布与观测 r
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1, 2, figsize=(10, 4), sharey=True)
for a, (r_obs, r0), ttl in zip(ax, tests, ["A: structural disassortativity only", "B: extra disassortative rewiring"]):
    a.hist(r0, bins=20, color="0.7", edgecolor="white", label="degree-preserving null (n=100)")
    a.axvline(r_obs, color="C3", lw=2, label=f"observed r = {r_obs:+.3f}")
    a.set_xlabel("degree assortativity r"); a.set_title(ttl, fontsize=10); a.legend(fontsize=8)
ax[0].set_ylabel("count")
fig.suptitle("Only network B is more disassortative than its degree sequence implies")
plt.tight_layout(); plt.show()