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

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

enter image description here

enter image description here

论文

L. K. Gallos, C. Song 和 H. A. Makse,度相关性的缩放及其对无标度网络扩散的影响,Phys. Rev. Lett. 100, 248701 (2008)。

您可以通过此链接阅读文章或阅读相关Google图书

问题

enter image description here

我的问题是如何使用Python计算网络的Eb(k)?我的问题是我无法复现作者的结果。我使用Condense Matter数据进行测试,Eb(k)的结果显示在上面的图中。 您可以看到我的图中一个问题是Eb(k)远大于1!!!我还尝试了互联网(As level数据)和WWW数据,但问题仍然存在。毫无疑问,我的算法或代码出了严重的问题。您可以复现我的结果,并与作者进行比较。非常感谢您提供的解决方案或建议。我将在下面介绍我的算法和Python脚本。

我按照以下步骤进行:

  1. 对于每条边,找到k=k的边以及k' > 3k的边。这些边的概率表示为P(k, k')
  2. 对于节点,获取度数大于b*k的节点比例,表示为p(k'),因此我们也可以得到k'*p(k')
  3. 获取分子P1:p1 = \sum P(k, k')/k'*P(k')
  4. 获取分母P2:P2 = \sum P(k')
  5. 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数据进行测试,您可以从此网址下载: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

1
一个区别是,似乎发表的论文使用了随着k变大而增长的箱子。这就是为什么发表的论文在x轴上具有对数刻度的等距分布符号,而你的则越来越密集。大于1的问题是另一回事。我稍后会尝试查看,但希望到那时有人已经解决了它。 - Joel
他们可以在www、互联网和蛋白质数据中使用对数分组。 - Frank Wang
1
请注意,网络首先使用[10] C. Song,L.K. Gallos,S. Havlin和H.A. Makse的盒覆盖方法进行“重新标准化”,J. Stat。Mech.(2007)P03006。 - Aric
1
在他们的图表标题中,他们说:“为了显示不变性,数据已经垂直移动。” 你关于你的图表的评论:“你可以看到我的图表中一个问题是Eb(k)远大于1!”我并不清楚这个表达式不能远大于1。 - Joel
只有互联网数据已经被重新调整以显示不变性。他们通过从0.01开始而不是0来垂直移动y轴。 - Frank Wang
3个回答

3
根据论文,Eb(k)的目的是获得相关指数epsilon:“[我们]引入一个标度不变量Eb(k),以简化epsilon的估计”(第二页,第一列底部)。
我没有找到使Eb(k)<1的方法,但我发现了一个可以正确计算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,而是它在k时给出约为0.5的对数斜率,而在论文中,该斜率约为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

Cond-mat网络的Eb(k)图

斜率是正确的(只有一位小数点后,这是该论文中提供的所有信息),因此可以正确计算ε。

关于Gamma

我从将斜率1.2和ε值1.6相加得到了γ=2.8的值(这遵循论文中方程4的结果)。我还使用了powerlaw Python模块进行了快速的合理性检查,以确定这个γ是否适合。

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

这个输出

2.84571139756

因此,对于伽玛值的正确值应为2.8(四舍五入后)。

使用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

WWW data

原论文中的斜率是0.15。通过查看论文中的图3(gamma-epsilon图)得到gamma值为2.6。

总结

我不确定为什么论文中的Eb(k)要比1小这么多。我相信论文中有一些未明确说明的重新缩放。然而,我能够使用Eb(k)恢复出正确的epsilon值。只要您能正确计算epsilon值,就不需要过于担心。


这是WWW数据的链接:http://www3.nd.edu/~networks/resources/www/www.dat.gz - Frank Wang
还要记得对数据进行对数分箱,可以使用下面的对数分箱函数。https://dev59.com/wloT5IYBdhLWcg3wlgTn#38408979 - Frank Wang
我怀疑作者只是将这一行更改为:pk = float(degree_dict[k])/node_number,变成了:pk = float(degree_dict[k])。 - Frank Wang
感谢您的反馈。我将测试WWW数据,并研究您的建议。 - bpachev
1
@FrankWang 作者可能已经用pk = float(degree_dict[k])/node_number替换了pk = float(degree_dict[k])。然而,这只会将所有内容缩小,无法解决得到错误斜率的问题。 - bpachev
显示剩余3条评论

0

看起来你正在使用离散分布计算条件概率,因此会得到许多零值,这会导致问题。

在论文中(第二页第二栏顶部),他们似乎正在使用幂律拟合数据,用一个漂亮的平滑函数替换嘈杂的离散值。我猜这也是为什么他们用积分而不是求和来写E_b的原因。

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


这并没有回答关于如何进行计算的提问。 - Joel
@pat 这是相互的 :) - Joel

0

如果考虑使用对数分组的数据,可以采用以下函数。

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)是更好的选择。

1
你应该把这个放到你的问题中。 - EvilTak
你的日志分箱函数在我尝试使用它处理我的数据时抛出了零除错误。我正在努力找出问题所在。你能用它正常工作吗? - bpachev

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接