在题为程度相关性的缩放及其对无标度网络中扩散的影响的论文中，作者定义了 $E_b(k)$ 的数量来衡量程度相关性的程度。
纸
L. K. Gallos、C. Song 和 H. A. Makse，度相关的标度及其对无标度网络中扩散的影响，物理学。牧师莱特。 100, 248701 (2008)。
您可以阅读下面的文章this link或阅读相关 google book .
问题
我的问题是如何用 Python 计算网络的 Eb(k)？我的问题是我无法重现作者的结果。我使用 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 数据进行测试，您可以从以下网址下载: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
最佳答案
根据论文，Eb(k) 的目的是得到相关指数 epsilon:“[We] 引入一个尺度不变量 Ebk 到 简化 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，而是它为您提供了大约 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
斜率(小数点后最多 1 位，这是论文中给出的全部内容)是正确的，因此现在可以正确计算 epsilon。
关于 Gamma
我通过将 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，我就不会太担心。
关于python - 如何用 Python 计算网络的 Eb(k)？，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/38408224/