Ответ 1
Согласно статье, целью Eb (k) является получение показателя корреляции epsilon: "[We] вводим масштабно-инвариантную величину Ebk в упростить оценку epsilon "(вторая страница, нижняя часть первого столбца).
Я не нашел способ сделать Eb (k) < 1, но я нашел исправление, которое правильно вычисляет epsilon.
Согласно уравнению 4 Eb (k) ~ k ^ - (эпсилон-гамма) (где распределение степени P (k) ~ k ^ -гамма, степенной закон). Таким образом, если мы построим наклон логарифма (Eb (k)) против log (k), мы должны получить гамма-эпсилон. Зная гамму, мы можем легко получить эпсилон.
Обратите внимание, что этот наклон инвариантен, если Eb (k) масштабируется константой. Таким образом, проблема с вашим вычисленным Eb (k) не, что она больше 1, но она дает вам логарифм наклона около 0,5 с k, тогда как в документе наклон около 1,2, поэтому вы получите неправильный эпсилон.
Мой алгоритм
Я начал с копирования кода, просмотра его и повторного его реализации эквивалентным образом. Моя ре-реализация повторила ваши результаты. Я вполне уверен, что вы внедрили дискретную версию формулы для E_b (k) правильно. Однако тщательное изучение статьи предполагает, что авторы использовали гладкие аппроксимации в своем коде.
На второй странице и в столбце указано равенство P (k | k ') = P (k, k')/(k ') ^ (1-gamma). Это эквивалентно замене точной вероятности P (k ') в знаменателе первого интеграла с гладким степенным приближением (k') ^ (- гамма) распределения степени и не является равенством.
Тот факт, что авторы утверждают, что это приближение как равенство без квалификации, подсказывает мне, что они, возможно, использовали его как таковой в своем коде. Итак, я решил использовать их приближение в коде, в результате чего ниже (где я получил гамма = 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
Наклон (до десятизначной цифры после десятичной точки, который является все, что дается в документе) является правильным, и, следовательно, теперь epsilon может быть правильно рассчитан.
Об Gamma
Я получил значение gamma = 2.8 от добавления наклона 1.2 к эпсилонному значению 1,6 (это следует из уравнения 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
Наклон от оригинальной бумаги равен 0,15. Я получил гамма-значение 2,6, посмотрев на рис. 3 в статье (диаграмма гамма-эпсилон).
В заключение
Я не уверен, почему Eb (k) настолько меньше, чем 1 на графике. Я почти уверен, что происходит перемасштабирование, которое не указано в документе. Тем не менее, я смог восстановить правильное значение epsilon, используя Eb (k). До тех пор, пока вы сможете правильно вычислить epsilon, я бы не стал слишком беспокоиться об этом.