在题为 Scaling of degree correlations and its influence on diffusion in scale-free networks 的论文中，作者定义了$ E_b（k）$的数量来衡量学位相关程度 .

enter image description here

纸

L. K. Gallos，C . Song和H. A. Makse，度量相关的缩放及其对无标度网络扩散的影响，物理学 . 莱特牧师 . 100,248701（2008） .

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

问题

enter image description here

我的问题是如何使用Python计算网络的Eb（k）？我的问题是我无法重现作者的结果 . 我使用Condense Matter数据测试它 . Eb（k）的结果如上图所示 . You can see that one problem in my figure is the Eb(k) is much larger than 1! !!我也尝试了互联网（作为级别数据）和WWW数据，问题仍然存在 . 毫无疑问，我的算法或代码存在严重问题 . 您可以重现我的结果，并将其与作者进行比较 . 您的解决方案或建议非常感谢 . 我将在下面介绍我的算法和python脚本 .

我按照以下步骤操作：

对于每条边，找到k = k和k ' > 3k. The probability of these edges is denoted as P(k, k'的边
对于节点，获取度数大于b * k的节点的比例，表示为p（k '), thus we can also have 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数据测试，你可以从这个网址下载：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()

Update ：问题尚未解决 .

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

3 回答

0

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

在论文（第二列顶部，第二页）中，看起来他们正在使用适合要替换的数据的幂律具有良好平滑功能的嘈杂离散值 . 而且，我认为，这也是为什么他们用积分而不是求和来写E_b .

如果我是你，我会问论文的作者他们的代码 . 然后我会要求期刊在没有支持代码的情况下停止发表论文 .

回复于 2024-04-29T11:36:16+08:00

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

我还没有找到使Eb（k）<1的方法，但我找到了一个修正 computes epsilon correctly .

根据等式4，Eb（k）~k ^ - （epsilon-gamma）（其中度分布P（k）~k ^ -gamma，幂律） . 因此，如果我们绘制log（Eb（k））与log（k）的斜率，我们应该得到gamma - epsilon . 知道伽玛，我们就可以轻松获得epsilon .

注意，如果Eb（k）按常数缩放，则该斜率是不变的 . 因此， problem 与你的计算Eb（k） is not 它大于1，但它给你一个约为.5的对数斜率与k，而在论文中斜率约为1.2，因此你将获得 wrong epsilon .

我的算法

我开始时通过复制代码，查看代码并以相同的方式重新实现它 . 我的重新实现复制了你的结果 . 我非常有信心你正确地为E_b（k）实现了公式的离散版本 . 然而，仔细研究该论文表明作者在其代码中使用了平滑近似 .

在第二页和第二列，第一个积分的分母中的等式P（k | k ') = P(k, k'）/（k ')^(1-gamma) is stated. This is equivalent to replacing the exact probability P(k'）与度分布的平滑幂律近似（k'）^（ - gamma），以及不是平等 .

作者将这种近似描述为没有限定的平等这一事实告诉我，他们可能已经在他们的代码中使用了它 . 所以，我决定在代码中使用它们的近似值，得到下面的内容（我得到的伽玛= 2.8，对于cond-mat将在下面解释） .

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

Plot of Eb(k) for the cond-mat network

斜率（小数点后最多1位数，这是纸上给出的全部数字）是正确的，因此现在可以正确计算epsilon .

关于Gamma

从斜率1.2加到epsilon值1.6得到gamma = 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图表），我得到了2.6的伽玛值 .

结论

我不确定为什么Eb（k）在纸张图形中小于1 . 我很确定正在进行一些重新调整，这在论文中并不明确 . 但是，我能够使用Eb（k）恢复正确的epsilon值 . 只要你能够正确计算epsilon，我就不会太担心它 .

回复于 2024-04-29T11:36:16+08:00

考虑使用数据的日志分级，可以采用以下功能 .

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将是更好的选择 .