随机块模型SBM的社区检测与似然比检验

第一部分：理解问题与准备工作

明确你的任务目标：使用随机块模型识别网络中的社区结构，并运用似然比检验来评估不同社区划分方案的优劣。这适用于社交网络、生物网络等任何由节点和边构成的数据。

准备你的输入数据。你需要一个网络的邻接矩阵 $A$，其中 $A_{ij} = 1$ 表示节点 $i$ 和 $j$ 之间有连接，$A_{ij} = 0$ 则表示没有。同时，你还需要一个假定的社区划分标签向量 $z$，其中 $z_i$ 表示节点 $i$ 所属的社区编号。

确定要比较的假设。原假设 $H_0$ 是一个基准划分（例如，已知的真实划分或一个简单的划分）。备择假设 $H_1$ 是一个你认为可能更好的新划分。检验的目标是判断新划分是否显著优于基准划分。

第二部分：构建随机块模型与计算似然值

理解随机块模型的核心思想。它假设同一社区内的节点连接概率相同，不同社区间的节点连接概率也相同。对于有 $K$ 个社区的网络，模型由社区概率矩阵 $\Omega$ 定义，其中 $\omega_{rs}$ 表示属于社区 $r$ 的节点与属于社区 $s$ 的节点相连的期望概率。

计算给定划分 $z$ 和概率矩阵 $\Omega$ 下网络 $A$ 的对数似然值。这个值衡量了模型解释数据的拟合程度。计算公式如下：
$$\mathcal{L}(A, z, \Omega) = \sum_{i<j} \left[ A_{ij} \ln(\omega_{z_i, z_j}) + (1-A_{ij}) \ln(1-\omega_{z_i, z_j}) \right]$$
对于大型网络，这个求和可以高效地通过统计每个社区对之间的实际连接数和可能连接数来计算。

估计最优的概率矩阵 $\hat{\Omega}$。对于给定的划分 $z$，使似然值最大的 $\Omega$ 可以直接计算出来。令 $n_{rs}$ 为社区 $r$ 和 $s$ 之间可能的边总数（若 $r=s$，则为 $n_r(n_r-1)/2$）， $m_{rs}$ 为实际存在的边数。那么：
$$\hat{\omega}_{rs} = \frac{m_{rs}}{n_{rs}}$$

计算在最优参数下的最大对数似然值。将 $\hat{\Omega}$ 代入之前的似然函数，得到：
$$\mathcal{L}_{\max}(z) = \sum_{r\le s} \left[ m_{rs} \ln(\hat{\omega}_{rs}) + (n_{rs} - m_{rs}) \ln(1-\hat{\omega}_{rs}) \right]$$
这个 $\mathcal{L}_{\max}(z)$ 是评估划分 $z$ 优劣的关键指标。

第三部分：执行似然比检验

计算原假设划分 $z_0$ 和备择假设划分 $z_1$ 各自的最大对数似然值 $\mathcal{L}_{\max}(z_0)$ 和 $\mathcal{L}_{\max}(z_1)$。

构造似然比统计量 $\Lambda$。这个统计量衡量了两个划分在似然上的差异：
$$\Lambda = 2 \left[ \mathcal{L}_{\max}(z_1) - \mathcal{L}_{\max}(z_0) \right]$$

判断 $\Lambda$ 的显著性。在零假设（$z_0$ 为真）下，当网络规模很大时，$\Lambda$ 近似服从卡方分布。自由度 $\nu$ 等于两个模型的独立参数个数之差。对于社区划分问题，参数主要是社区概率矩阵 $\Omega$ 中的独立元素。因此：
$$\nu = \left(\frac{K_1(K_1+1)}{2}\right) - \left(\frac{K_0(K_0+1)}{2}\right)$$
其中 $K_0$ 和 $K_1$ 分别是划分 $z_0$ 和 $z_1$ 的社区数量。

执行检验。设定一个显著性水平，例如 $\alpha = 0.05$。将计算出的 $\Lambda$ 值与自由度为 $\nu$ 的卡方分布临界值进行比较。如果 $\Lambda$ 大于该临界值，则拒绝原假设，认为备择假设划分 $z_1$ 在统计上显著优于 $z_0$。

第四部分：实践指南与代码示例

选择一个合适的图分析库。在Python中，NetworkX 用于基础图操作，graph-tool 或 karateclub 等库提供了更高效的SBM实现。以下流程以 graph-tool 为例。

执行基础SBM社区检测并获取似然值。graph-tool 的 minimize_blockmodel_dl 函数可以自动寻找一个良好的划分，并返回其对数似然。

import graph_tool.all as gt

# 从边列表或邻接矩阵构建图 `g`
# ... (此处省略图构建代码)

# 执行最小描述长度优化，寻找一个SBM划分
state = gt.minimize_blockmodel_dl(g)
# 获取该划分的对数似然值
log_likelihood_h1 = state.entropy()

# 获取划分标签
z1 = state.get_blocks().a

获取基准划分的对数似然。假设你有一个基准划分的标签数组 z0。

# 使用已有的划分 `z0` 创建一个SBM状态
state0 = gt.BlockState(g, b=z0)
# 计算该固定划分下的对数似然（实际上是最小描述长度）
log_likelihood_h0 = state0.entropy()

计算似然比统计量并进行检验。

import scipy.stats as stats

# 计算似然比统计量
Lambda = 2 * (log_likelihood_h0 - log_likelihood_h1)  # 注意：graph-tool的entropy是描述长度，越小越好，等价于负的似然

# 确定自由度（需要根据你的z0和z1中的社区数K0和K1计算）
K0 = len(set(z0))
K1 = len(set(z1))
df = (K1*(K1+1)//2) - (K0*(K0+1)//2)

# 计算p值
p_value = 1 - stats.chi2.cdf(Lambda, df)

# 设定显著性水平
alpha = 0.05

# 做出判断
if p_value < alpha:
    print(f"拒绝原假设，p值={p_value:.4f}。新划分显著更优。")
else:
    print(f"无法拒绝原假设，p值={p_value:.4f}。无充分证据表明新划分更优。")

注意进行检验的前提条件。此方法适用于较大规模的网络。对于非常小的网络，卡方近似可能不准确，需要采用置换检验等更稳健的方法。