Ответ 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, я бы не стал слишком беспокоиться об этом.

Ответ 2

Учитывая использование лог-бинания данных, можно воспользоваться следующей функцией.

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

Обычно для отношения масштабирования логический бининг будет лучшим выбором.

Ответ 3

Похоже, вы на самом деле вычисляете условную вероятность с использованием дискретных распределений, поэтому вы получаете много нулей, что создает проблемы.

В документе (вверху второго столбца, вторая страница) похоже, что они используют закон мощности, соответствующий данным, чтобы заменить шумные дискретные значения приятной гладкой функцией. И это также, я полагаю, почему они пишут E_b в терминах интегралов, а не суммирования.

Если бы я был вами, я бы попросил авторов статьи для их кода. И тогда я попрошу журнал прекратить публикацию статей без поддержки кода.