NetworkX 2.8.8 实战：Python 构建 BA 无标度网络与幂律分析

1. 复杂网络与无标度网络基础

当我们观察互联网、社交网络或蛋白质相互作用网络时，会发现这些系统具有一个共同特征：少数节点拥有大量连接，而大多数节点只有少量连接。这种特性正是无标度网络（Scale-Free Network）的典型表现。

无标度网络由Barabási和Albert在1999年提出，其核心特征体现在度分布服从幂律分布：

P(k) ~ k^(-γ)

其中γ通常介于2到3之间。这种网络具有两个关键机制：

• 增长特性 ：网络规模随时间不断扩大
• 优先连接 ：新节点倾向于连接到已有高度节点

现实世界中许多网络都呈现无标度特性：

• 互联网拓扑结构（γ≈2.2）
• 学术引用网络（γ≈3）
• 社交网络关注关系（γ≈2.1）

2. 环境配置与NetworkX基础

2.1 安装与导入

确保已安装最新版NetworkX和科学计算套件：

pip install networkx==2.8.8 matplotlib numpy scipy

基础导入与版本检查：

import networkx as nx
import matplotlib.pyplot as plt
import numpy as np
from collections import Counter
print(f"NetworkX版本: {nx.__version__}")

2.2 网络生成基础

NetworkX提供多种网络生成方法：

# ER随机网络
G_er = nx.erdos_renyi_graph(100, 0.1)

# WS小世界网络
G_ws = nx.watts_strogatz_graph(100, 4, 0.2)

# BA无标度网络
G_ba = nx.barabasi_albert_graph(100, 2)

网络基本统计量计算：

def network_stats(G):
    stats = {
        "节点数": nx.number_of_nodes(G),
        "边数": nx.number_of_edges(G),
        "平均度": sum(dict(G.degree()).values())/nx.number_of_nodes(G),
        "集聚系数": nx.average_clustering(G),
        "平均路径长度": nx.average_shortest_path_length(G)
    }
    return stats

3. BA无标度网络构建实战

3.1 标准BA模型实现

NetworkX内置的 barabasi_albert_graph() 函数参数说明：

• n : 最终网络节点数
• m : 新节点加入时创建的边数

# 生成1000个节点的BA网络，每个新节点连接2条边
G = nx.barabasi_albert_graph(1000, 2)

# 可视化前100个节点
pos = nx.spring_layout(G.subgraph(list(G.nodes)[:100])) 
nx.draw(G.subgraph(list(G.nodes)[:100]), pos, node_size=20, alpha=0.6)
plt.title("BA网络局部可视化")
plt.show()

3.2 自定义BA模型实现

更灵活的手动实现方式：

def custom_ba_graph(n, m):
    G = nx.complete_graph(m)  # 初始为完全图
    targets = list(range(m))
    
    for i in range(m, n):
        # 优先连接：按度比例选择目标节点
        degrees = np.array([G.degree(j) for j in G.nodes()])
        probs = degrees / degrees.sum()
        targets = np.random.choice(G.nodes(), size=m, replace=False, p=probs)
        
        G.add_node(i)
        for target in targets:
            G.add_edge(i, target)
    
    return G

# 性能对比
%timeit nx.barabasi_albert_graph(1000,2)  # 约3.5ms
%timeit custom_ba_graph(1000,2)          # 约120ms

注意：自定义实现更灵活但效率较低，适合需要修改连接规则的研究场景

4. 度分布分析与幂律拟合

4.1 度分布计算与可视化

def plot_degree_distribution(G, loglog=True):
    degrees = [d for n, d in G.degree()]
    degree_counts = Counter(degrees)
    x = list(degree_counts.keys())
    y = [degree_counts[k]/len(G) for k in x]
    
    plt.scatter(x, y, alpha=0.6, label='实际分布')
    
    if loglog:
        plt.xscale('log')
        plt.yscale('log')
        plt.xlabel('度 (log)')
        plt.ylabel('概率 (log)')
        plt.title('双对数坐标下的度分布')
    else:
        plt.xlabel('度')
        plt.ylabel('概率')
        plt.title('线性坐标下的度分布')
    
    plt.grid(True, which="both", ls="-")
    return x, y

# 绘制双对数坐标图
x, y = plot_degree_distribution(G)

4.2 幂律指数γ的拟合

使用最小二乘法进行幂律拟合：

from scipy.optimize import curve_fit

def power_law(x, a, b):
    return a * np.power(x, -b)

# 过滤零值
x_fit = np.array(x)[y > 0]
y_fit = np.array(y)[y > 0]

# 拟合
params, cov = curve_fit(power_law, x_fit, y_fit)
gamma = params[1]
print(f"拟合幂律指数γ = {gamma:.3f}")

# 绘制拟合曲线
plt.plot(x_fit, power_law(x_fit, *params), 'r--', 
         label=f'拟合曲线: γ={gamma:.3f}')
plt.legend()
plt.show()

典型输出结果：

拟合幂律指数γ = 2.876

4.3 拟合质量评估

计算R²值评估拟合优度：

def r_squared(y_true, y_pred):
    ss_res = np.sum((y_true - y_pred) ** 2)
    ss_tot = np.sum((y_true - np.mean(y_true)) ** 2)
    return 1 - (ss_res / ss_tot)

y_pred = power_law(x_fit, *params)
r2 = r_squared(y_fit, y_pred)
print(f"R² = {r2:.3f}")

高质量拟合应满足：

• γ介于2-3之间
• R² > 0.9
• 在双对数坐标下呈现线性关系

5. 进阶分析与可视化技巧

5.1 累积度分布分析

累积度分布可减少尾部噪声影响：

def plot_cumulative_dist(G):
    degrees = [d for n, d in G.degree()]
    max_degree = max(degrees)
    
    # 计算累积分布
    hist = np.zeros(max_degree + 1)
    for d in degrees:
        hist[d] += 1
    cum_dist = 1 - np.cumsum(hist)/np.sum(hist)
    
    # 绘制
    plt.loglog(range(max_degree + 1), cum_dist, 'bo', alpha=0.6)
    plt.xlabel('度 (log)')
    plt.ylabel('P(K≥k) (log)')
    plt.title('累积度分布')
    plt.grid(True)
    
    # 拟合
    valid = cum_dist > 0
    x = np.arange(max_degree + 1)[valid]
    y = cum_dist[valid]
    params, _ = curve_fit(power_law, x, y)
    plt.plot(x, power_law(x, *params), 'r--', 
             label=f'γ={params[1]:.3f}')
    plt.legend()
    
plot_cumulative_dist(G)

5.2 节点重要性分析

识别网络中的关键枢纽节点：

def analyze_centrality(G):
    # 计算多种中心性指标
    metrics = {
        '度中心性': nx.degree_centrality(G),
        '介数中心性': nx.betweenness_centrality(G),
        '接近中心性': nx.closeness_centrality(G),
        '特征向量中心性': nx.eigenvector_centrality(G, max_iter=1000)
    }
    
    # 找出每种指标的前5节点
    top_nodes = {}
    for name, values in metrics.items():
        sorted_nodes = sorted(values.items(), key=lambda x: x[1], reverse=True)[:5]
        top_nodes[name] = [n[0] for n in sorted_nodes]
    
    return metrics, top_nodes

centralities, top_nodes = analyze_centrality(G)
print("关键节点识别结果:")
for metric, nodes in top_nodes.items():
    print(f"{metric}: {nodes}")

5.3 动态生长过程可视化

展示网络生长过程中的优先连接现象：

def visualize_growth(n=50, m=1):
    G = nx.Graph()
    G.add_node(0)
    
    plt.figure(figsize=(12,8))
    
    for i in range(1, n):
        # 优先连接
        degrees = np.array([G.degree(j) for j in G.nodes])
        probs = degrees / degrees.sum()
        targets = np.random.choice(G.nodes, size=m, replace=False, p=probs)
        
        G.add_node(i)
        for target in targets:
            G.add_edge(i, target)
        
        # 每10步绘制一次
        if i % 10 == 0 or i == n-1:
            plt.clf()
            pos = nx.spring_layout(G)
            node_size = [d*50 for n, d in G.degree()]
            nx.draw(G, pos, node_size=node_size, 
                   with_labels=False, alpha=0.7)
            plt.title(f"BA网络生长过程 (n={i+1})")
            plt.pause(0.5)
    
visualize_growth(n=30, m=1)

6. 实际应用与扩展

6.1 不同参数对网络的影响

比较不同m值对网络特性的影响：

m值	平均度	集聚系数	幂律指数γ
1	1.98	0.021	2.91
2	3.96	0.038	2.87
5	9.90	0.092	2.79

def compare_parameters(n=1000, m_list=[1,2,5]):
    results = []
    for m in m_list:
        G = nx.barabasi_albert_graph(n, m)
        x, y = plot_degree_distribution(G, loglog=False)
        params, _ = curve_fit(power_law, np.array(x)[y>0], np.array(y)[y>0])
        
        results.append({
            'm': m,
            'avg_degree': sum(dict(G.degree()).values())/n,
            'clustering': nx.average_clustering(G),
            'gamma': params[1]
        })
    return pd.DataFrame(results)

6.2 鲁棒性分析

测试网络对随机攻击和针对性攻击的抵抗能力：

def test_robustness(G, attack='random'):
    giant = []
    removed = []
    nodes = list(G.nodes())
    
    if attack == 'targeted':
        nodes.sort(key=lambda x: G.degree(x), reverse=True)
    
    for i in range(1, 101):
        # 移除1%的节点
        to_remove = nodes[:int(len(nodes)*0.01*i)]
        H = G.copy()
        H.remove_nodes_from(to_remove)
        
        # 计算最大连通分量比例
        if nx.number_of_nodes(H) > 0:
            largest = max(nx.connected_components(H), key=len)
            giant.append(len(largest)/nx.number_of_nodes(H))
        else:
            giant.append(0)
        
        removed.append(i)
    
    plt.plot(removed, giant, label=attack)
    plt.xlabel('移除节点比例(%)')
    plt.ylabel('最大连通分量比例')
    plt.title('网络鲁棒性测试')
    plt.legend()

plt.figure()
test_robustness(G, 'random')
test_robustness(G, 'targeted')

6.3 与其他网络模型对比

def compare_models(n=1000):
    models = {
        'ER随机网络': nx.erdos_renyi_graph(n, 4/n),
        'WS小世界': nx.watts_strogatz_graph(n, 4, 0.1),
        'BA无标度': nx.barabasi_albert_graph(n, 2)
    }
    
    fig, axes = plt.subplots(1, 3, figsize=(18,5))
    for ax, (name, G) in zip(axes, models.items()):
        x, y = plot_degree_distribution(G, loglog=False)
        ax.scatter(x, y, alpha=0.6)
        ax.set_title(name)
        ax.set_xlabel('度')
        ax.set_ylabel('概率')
    plt.tight_layout()

compare_models()