如何用 Python 计算网络的 Eb(k)？

Posted 2023-02-16

【中文标题】如何用 Python 计算网络的 Eb(k)？

【英文标题】：How to calculate Eb(k) of networks with Python?

【发布时间】：2016-11-19 09:05:50

【问题描述】：

在题为程度相关性的缩放及其对无标度网络中的扩散的影响的论文中，作者定义了 $E_b(k)$ 的数量来衡量程度相关性的程度。

纸

L. K. Gallos、C. Song 和 H. A. Makse，度相关的标度及其对无标度网络中扩散的影响，物理学。牧师莱特。 100, 248701 (2008)。

您可以阅读this link后面的文章或阅读相关的google book。

问题

我的问题是我的问题是我无法重现作者的结果。我使用 Condense Matter 数据对其进行测试。 Eb(k)的结果如上图所示。 您可以看到我图中的一个问题是 Eb(k) 远大于 1！！！我也试过互联网（作为级别数据）和万维网数据，问题依然存在。毫无疑问，我的算法或代码存在严重问题。您可以复制我的结果，并将其与作者的结果进行比较。非常感谢您的解决方案或建议。下面我将介绍我的算法和python脚本。

我遵循以下步骤：

对于每条边，找到 k=k 且 k' > 3k 的边。这些边的概率记为 P(k, k') 对于节点，求度数大于b*k的节点的比例，记为p(k')，因此我们也可以有k'*p(k') 要获得分子 P1：p1 = \sum P(k, k')/k'*P(k') 求分母p2：P2 = \sum P(k') Eb(k) = p1/p2

Python 脚本

python脚本如下：

%matplotlib inline
import networkx as nx
import matplotlib.cm as cm
import matplotlib.pyplot as plt
from collections import defaultdict

def ebks(g, b):
    edge_dict = defaultdict(lambda: defaultdict(int))
    degree_dict = defaultdict(int)
    edge_degree = [sorted(g.degree(e).values()) for e in g.edges()]
    for e in edge_degree:
        edge_dict[e[0]][e[-1]] +=1
    for i in g.degree().values():
        degree_dict[i] +=1
    edge_number = g.number_of_edges()
    node_number = g.number_of_nodes()
    ebks, ks = [], []
    for k1 in edge_dict:
        p1, p2 = 0, 0
        for k2 in edge_dict[k1]:
            if k2 >= b*k1:
                pkk = float(edge_dict[k1][k2])/edge_number
                pk2 = float(degree_dict[k2])/node_number
                k2pk2 = k2*pk2
                p1 += pkk/k2pk2
        for k in degree_dict:
            if k>=b*k1:
                pk = float(degree_dict[k])/node_number
                p2 += pk
        if p2 > 0:
            ebks.append(p1/p2)
            ks.append(k1)
    return ebks, ks

我用 ca-CondMat 数据测试，你可以从这个 url 下载它：http://snap.stanford.edu/data/ca-CondMat.html

# Load the data
# Remember to change the file path to your own
ca = nx.Graph()
with open ('/path-of-your-file/ca-CondMat.txt') as f:
    for line in f:
        if line[0] != '#':
            x, y = line.strip().split('\t')
            ca.add_edge(x,y)
nx.info(ca)

#calculate ebk 
ebk, k = ebks(ca, b=3)

plt.plot(k,ebk,'r^')
plt.xlabel(r'$k$', fontsize = 16)
plt.ylabel(r'$E_b(k)$', fontsize = 16)
plt.xscale('log')
plt.yscale('log')
plt.show()

更新：问题还没有解决。

def ebkss(g, b, x):
    edge_dict = defaultdict(lambda: defaultdict(int))
    degree_dict = defaultdict(int)
    edge_degree = [sorted(g.degree(e).values()) for e in g.edges()]
    for e in edge_degree:
        edge_dict[e[0]][e[-1]] +=1
    for i in g.degree().values():
        degree_dict[i] +=1
    edge_number = g.number_of_edges()
    node_number = g.number_of_nodes()
    ebks, ks = [], []
    for k1 in edge_dict:
        p1, p2 = 0, 0
        nk2k = np.sum(edge_dict[k1].values())
        pk1 = float(degree_dict[k1])/node_number
        k1pk1 = k1*pk1
        for k2 in edge_dict[k1]:
            if k2 >= b*k1:
                pk2k = float(edge_dict[k1][k2])/nk2k
                pk2 = float(degree_dict[k2])/node_number
                k2pk2 = k2*pk2
                p1 += (pk2k*k1pk1)/k2pk2
        for k in degree_dict:
            if k>=b*k1:
                pk = float(degree_dict[k])/node_number
                p2 += pk
        if p2 > 0:
            ebks.append(p1/p2**x)
            ks.append(k1)
    return ebks, ks

【问题讨论】：

一个不同之处在于，看起来已发表的论文使用的 bin 会随着 k 的增大而增长。这就是为什么已发表的论文在 x 上用对数刻度均匀分布符号，而您的论文变得越来越密集。大于 1 的问题是另外一回事。我稍后会尝试看看，但希望到那时有人会解决它。

他们可能对 www、internet 和蛋白质数据使用 log-binning。

请注意，首先使用 [10] C. Song, L. K. Gallos, S. Havlin 和 H. A. Makse, J. Stat 中的框覆盖方法对网络进行“重新归一化”。机甲。 (2007) P03006。

在他们的图的标题中，他们说“数据已经垂直移动以显示不变性。”您对您的图的评论：“您可以看到我图中的一个问题是 Eb(k) 远大于 1！！”我不清楚这个表达式不能远大于 1。

只有互联网数据被重新规范化以显示不变性。他们从 0.01 而不是 0 开始垂直移动 y 轴。

【参考方案1】：

根据论文，Eb(k) 的目的是得到相关指数 epsilon：“[We] 引入一个尺度不变量 Ebk 到 简化 epsilon 的估计”（第二页，第一列底部）。

我还没有找到使 Eb(k) 正确计算 epsilon 的更正方法。

根据等式4，Eb(k)~k^-(epsilon-gamma)（其中度数分布P(k)~k^-gamma，一个幂律）。因此，如果我们绘制 log(Eb(k)) 对 log(k) 的斜率，我们应该得到 gamma - epsilon。知道了 gamma，我们就可以很容易地得到 epsilon。

请注意，如果 Eb(k) 按常数缩放，则该斜率是不变的。因此，您计算出的 Eb(k) 问题不是大于 1，而是它为您提供了大约 0.5 的 k 对数斜率，而在论文中，斜率约为 1.2，因此您会得到错误的 epsilon。

我的算法

我首先复制您的代码，查看它，然后以等效的方式重新实现它。我的重新实现复制了您的结果。我非常有信心您正确地实现了 E_b(k) 公式的离散版本。然而，仔细研究论文表明作者在他们的代码中使用了平滑近似。

在第二页和第二列，陈述了等式 P(k|k') = P(k, k')/ (k')^(1-gamma)。这相当于用度分布的平滑幂律近似 (k')^(-gamma) 代替第一个积分的分母中的精确概率 P(k')，并且 不是 平等。

作者将这种近似表示为无条件的等式这一事实向我表明，他们可能在他们的代码中使用了它。所以，我决定在代码中使用它们的近似值，结果如下（我得到的 cond-mat 的 gamma = 2.8 如下所述）。

def ebkss(g, b, gamma=2.8):
    edge_dict = defaultdict(lambda: defaultdict(int))
    degree_dict = defaultdict(int)
    edge_degree = [sorted(g.degree(e).values()) for e in g.edges()]
    for e in edge_degree:
        edge_dict[e[0]][e[-1]] +=1
    for i in g.degree().values():
        degree_dict[i] +=1
    edge_number = g.number_of_edges()
    node_number = g.number_of_nodes()
    ebks, ks = [], []
    for k1 in edge_dict:
        p1, p2 = 0, 0
        nk2k = np.sum(edge_dict[k1].values())
        pk1 = float(degree_dict[k1])/node_number
        k1pk1 = k1*pk1

        for k2 in edge_dict[k1]:
            if k2 >= b*k1:
                pk2k = float(edge_dict[k1][k2])/edge_number
                pk2 = float(degree_dict[k2])/node_number
                p1 += pk2k/(k2*k2**(-gamma))
        for k in degree_dict:
            if k>=b*k1:
                pk = float(degree_dict[k])/node_number
                p2 += pk
        if p2 > 0 and p1 > 0:
            ebks.append(p1/p2)
            ks.append(k1)
    return ebks, ks

结果

使用此代码：

def get_logslope(x,y):
    A = np.empty((len(x), 2))
    A[:,0] = np.log(x)
    A[:,1] = 1
    res = la.lstsq(A, np.log(y))
    return res[0]

def show_eb(ca, b, gamma):
    #calculate ebk 
    ebk, k = ebkss(ca, b=b,gamma=gamma)
    print "Slope = ", get_logslope(np.array(k), np.array(ebk) )
    plt.plot(k,ebk,'r^')
    plt.xlabel(r'$k$', fontsize = 16)
    plt.ylabel(r'$E_b(k)$', fontsize = 16)
    plt.xscale('log')
    plt.yscale('log')
    plt.show()
show_eb(ca, 3, 2.8)

我得到了这个输出：

Slope =  1.22136715547

斜率（小数点后最多一位，这是论文中给出的全部内容）是正确的，因此现在可以正确计算 epsilon。

关于伽玛

我通过将 1.2 的斜率与 1.6 的 epsilon 值相加得到 gamma = 2.8 的值（这来自论文的等式 4）。我还使用 powerlaw Python 模块进行了快速的健全性检查，以确定这个 gamma 是否合适。

import powerlaw
res = powerlaw.Fit(np.array(ca.degree().values())+1, xmin=10)
print res.alpha

这个输出

2.84571139756

因此 2.8 对 gamma 的值是正确的，直到四舍五入。

用 WWW 数据编辑

我用 WWW 数据集测试了我的方法。我最终得到了一个接近论文中的斜率，但缩放仍然关闭。 这是我的代码：

def log_binning(x, y, bin_count=50):
    max_x = np.log10(max(x))
    max_y = np.log10(max(y))
    max_base = max([max_x,max_y])
    xx = [i for i in x if i>0]
    min_x = np.log10(np.min(xx))
    bins = np.logspace(min_x,max_base,num=bin_count)
    hist = np.histogram(x,bins)[0]
    nonzero_mask = np.logical_not(hist==0)       
    hist[hist==0] = 1
    bin_means_y = (np.histogram(x,bins,weights=y)[0] / hist)
    bin_means_x = (np.histogram(x,bins,weights=x)[0] / hist)
    return bin_means_x[nonzero_mask],bin_means_y[nonzero_mask]
def single_line_read(fname):    
    g = nx.Graph()
    with open(fname, "r") as f:
        for line in f:
          a = map(int,line.strip().split(" "))
          g.add_edge(a[0], a[1])
    return g

www = single_line_read("data/www.dat")
ebk, k = ebkss(www, 3, 2.6)
lk, lebk = log_binning(np.array(k,dtype=np.float64), np.array(ebk), bin_count=70)
#print lk, lebk
print "Slope", get_logslope(lk, lebk)
plt.plot(lk,lebk/www.number_of_edges(),'r^')
plt.xlabel(r'$k$', fontsize = 16)
plt.ylabel(r'$E_b(k)$', fontsize = 16)
plt.xscale('log')
plt.yscale('log')
plt.show()

坡度 0.162453554297

原始论文的斜率为 0.15。通过查看论文中的图 3（gamma-epsilon 图表），我得到了 2.6 的 gamma 值。

结论

我不确定为什么 Eb(k) 在论文图形中比 1 小得多。我很确定正在进行一些重新调整，这在论文中没有明确说明。但是，我能够使用 Eb(k) 恢复 epsilon 的正确值。只要您能够正确计算 epsilon，我就不会太担心。

【讨论】：

这是WWW数据的链接www3.nd.edu/~networks/resources/www/www.dat.gz

还要记得对数据进行log-binning，可以使用下面的log-binning功能。 ***.com/a/38408979/912429

我怀疑作者只是将这一行：pk = float(degree_dict[k])/node_number 改为：pk = float(degree_dict[k])

感谢您的反馈。我要测试一下 WWW 数据，看看你的建议。

@FrankWang 作者可能已将 pk = float(degree_dict[k])/node_number 替换为 pk = float(degree_dict[k])。但是，这只会按比例缩小所有内容，并不能解决斜率错误的问题。

【参考方案2】：

考虑使用数据的log-binning，可以采用以下功能。

import numpy as np

def log_binning(x, y, bin_count=35):
    max_x = np.log10(max(x))
    max_y = np.log10(max(y))
    max_base = max([max_x,max_y])
    xx = [i for i in x if i>0]
    min_x = np.log10(np.min(xx))
    bins = np.logspace(min_x,max_base,num=bin_count)
    bin_means_y = (np.histogram(x,bins,weights=y)[0] / np.histogram(x,bins)[0])
    bin_means_x = (np.histogram(x,bins,weights=x)[0] / np.histogram(x,bins)[0])
    return bin_means_x,bin_means_y

如果您想对数据进行线性分箱，请使用以下函数：

def LinearBinData(x, y, number): 
    data=sorted(zip(x,y))
    rs = np.linspace(min(x),max(x),number)
    rs = np.transpose(np.vstack((rs[:-1],rs[1:])))
    ndata = []
    within = []
    for start,end in rs:
        for i,j in data:
            if i>=start and i<end:
                within.append(j)
        ndata.append([(start+end)/2.0,np.mean(np.array(within))]  )
    nx,ny = np.array(ndata).T
    return nx,ny

通常，对于缩放关系，log-binning 会是更好的选择。

【讨论】：

你应该把这个放到你的问题中。

当我尝试将它与我的数据一起使用时，您的日志分箱功能会引发零除法错误。我正在努力追查问题。它对你有用吗？

【参考方案3】：

看起来您实际上是在使用离散分布计算条件概率，因此您会得到很多零，这会产生问题。

在论文（第二列顶部，第二页）中，看起来他们正在使用适合数据的幂律来用一个很好的平滑函数替换嘈杂的离散值。我想这也是为什么他们用积分而不是求和来写 E_b 的原因。

如果我是你，我会向论文的作者询问他们的代码。然后我会要求期刊停止发表没有支持代码的论文。

【讨论】：

这并没有回答 OP 关于如何进行计算的问题。

@pat 这是相互的 :)