Используя pandas и numpy для параметризации стека переполнения количества пользователей и репутации

Ответ 1

NumPy имеет множество функций для подгонки. Для полиномиальных подстановок мы используем numpy.polyfit (документация).

Инициализировать свой набор данных:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

data = [k.split() for k in '''0           1  4364226
1         200   269110
2         500   158824
3        1000    90368
4        2000    48609
5        3000    32604
6        5000    18921
7       10000     8618
8       25000     2802
9       50000     1000
10     100000      334'''.split('\n')]

soDF = pd.DataFrame(data, columns=('index', 'total_rep', 'users'))

soDF['total_rep'] = pd.to_numeric(soDF['total_rep'])
soDF['users'] = pd.to_numeric(soDF['users'])

soDF['log_total_rep'] = soDF['total_rep'].apply(np.log10)
soDF['log_users'] = soDF['users'].apply(np.log10)
soDF.plot(x='log_total_rep', y='log_users')

Установите многочлен 2-й степени

coefficients = np.polyfit(soDF['log_total_rep'] , soDF['log_users'], 2)

print "Coefficients: ", coefficients

Далее, давайте нарисуем оригинал + fit:

polynomial = np.poly1d(coefficients)
xp = np.linspace(-2, 6, 100)

plt.plot(soDF['log_total_rep'], soDF['log_users'], '.', xp, polynomial(xp), '-')

Ответ 2

python, pandas и scipy, о мой!

В научной экосистеме питона имеется несколько бесплатных библиотек. Ни одна библиотека не делает все по дизайну. pandas предоставляет инструменты для обработки табличных данных и таймсерий. Однако он намеренно не включает тип функциональности, который вы ищете.

Для установки статистических распределений обычно используется другой пакет, например scipy.stats.

Однако в этом случае у нас нет "сырых" данных (т.е. длинной последовательности оценок репутации). Вместо этого мы имеем нечто похожее на гистограмму. Поэтому нам нужно будет соответствовать вещам на чуть более низком уровне, чем scipy.stats.powerlaw.fit.

Автономный пример

На данный момент оставьте pandas полностью. Здесь нет никаких преимуществ в использовании, и мы в любом случае быстро завершим преобразование данных в другие структуры данных. pandas отлично, он просто переполняет эту ситуацию.

Как быстрый автономный пример для воспроизведения вашего сюжета:

import matplotlib.pyplot as plt

total_rep = [1, 200, 500, 1000, 2000, 3000, 5000, 10000,
             25000, 50000, 100000]
num_users = [4364226, 269110, 158824, 90368, 48609, 32604, 
             18921, 8618, 2802, 1000, 334]

fig, ax = plt.subplots()
ax.loglog(total_rep, num_users)
ax.set(xlabel='Total Reputation', ylabel='Number of Users',
       title='Log-Log Plot of Stackoverflow Reputation')
plt.show()

Что представляют эти данные?

Затем нам нужно знать, с чем мы работаем. То, что мы построили, похоже на гистограмму, так как это необработанные подсчеты количества пользователей на заданном уровне репутации. Однако обратите внимание на маленькую "+" рядом с каждой ячейкой таблицы репутации. Это означает, что, например, 2082 пользователей имеют рейтинг репутации 25000 или больше.

Наши данные в основном представляют собой оценку функции кумулятивного распределения (CCDF), в том же смысле, что гистограмма представляет собой оценку функции распределения вероятности (PDF). Нам просто нужно нормализовать его на общее количество пользователей в нашем примере, чтобы получить оценку CCDF. В этом случае мы можем просто делить на первый элемент num_users. Репутация никогда не может быть меньше 1, поэтому 1 по оси x соответствует вероятности 1 по определению. (В других случаях нам нужно оценить это число.) В качестве примера:

import numpy as np
import matplotlib.pyplot as plt

total_rep = np.array([1, 200, 500, 1000, 2000, 3000, 5000, 10000,
                      25000, 50000, 100000])
num_users = np.array([4364226, 269110, 158824, 90368, 48609, 32604, 18921,
                      8618, 2802, 1000, 334])

ccdf = num_users.astype(float) / num_users.max()

fig, ax = plt.subplots()
ax.loglog(total_rep, ccdf, color='lightblue', lw=2, marker='o',
          clip_on=False, zorder=10)
ax.set(xlabel='Reputation', title='CCDF of Stackoverflow Reputation',
       ylabel='Probability that Reputation is Greater than X')
plt.show()

Вы можете удивиться, почему мы конвертируем вещи в "нормализованную" версию. Самый простой ответ заключается в том, что он более полезен. Это позволяет нам говорить то, что напрямую не связано с нашим размером выборки. Завтра общее количество пользователей Stackoverflow (и номеров на каждом уровне репутации) будет отличаться. Однако общая вероятность того, что любой данный пользователь имеет определенную репутацию, не будет существенно изменяться. Если мы хотим предсказать репутацию Джона Скита (наивысший представитель пользователя), когда сайт попадает на 5 миллионов зарегистрированных пользователей, гораздо проще использовать вероятности вместо необработанных счетчиков.

Наивное соответствие степенного распределения

Далее, давайте подпишем распределение мощности по закону CCDF. Опять же, если бы у нас были "сырые" данные в виде длинного списка показателей репутации, было бы лучше использовать статистический пакет для этого. В частности, scipy.stats.powerlaw.fit.

Однако у нас нет необработанных данных. CCDF степенного распределения принимает вид ccdf = x**(-a + 1). Поэтому мы поместим строку в лог-пространстве, и мы можем получить параметр a для распределения из a = 1 - slope.

На данный момент используйте np.polyfit для соответствия строке. Нам нужно обработать преобразование назад и вперед из лог-пространства самостоятельно:

import numpy as np
import matplotlib.pyplot as plt

total_rep = np.array([1, 200, 500, 1000, 2000, 3000, 5000, 10000,
                      25000, 50000, 100000])
num_users = np.array([4364226, 269110, 158824, 90368, 48609, 32604, 18921,
                      8618, 2802, 1000, 334])

ccdf = num_users.astype(float) / num_users.max()

# Fit a line in log-space
logx = np.log(total_rep)
logy = np.log(ccdf)
params = np.polyfit(logx, logy, 1)
est = np.exp(np.polyval(params, logx))

fig, ax = plt.subplots()
ax.loglog(total_rep, ccdf, color='lightblue', ls='', marker='o',
          clip_on=False, zorder=10, label='Observations')

ax.plot(total_rep, est, color='salmon', label='Fit', ls='--')

ax.set(xlabel='Reputation', title='CCDF of Stackoverflow Reputation',
       ylabel='Probability that Reputation is Greater than X')

plt.show()

У нас возникла непосредственная проблема с этим подходом. Согласно нашей оценке, вероятность того, что у пользователей будет репутация 1, будет более 1, что невозможно.

Проблема состоит в том, что мы можем выбрать polyfit наилучший y-перехват для нашей линии. Если мы посмотрим на params в нашем коде выше, это второе число:

In [11]: params
Out[11]: array([-0.81938338,  1.15955974])

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

Фиксация y-перехвата в подгонке

Прежде всего, обратите внимание, что log(1) --> 0. Поэтому мы фактически хотим заставить y-перехват в лог-пространстве равным 0 вместо 1.

Проще всего это сделать, используя np.linalg.lstsq для решения задач вместо np.polyfit. Во всяком случае, вы бы сделали что-то похожее на:

import numpy as np
import matplotlib.pyplot as plt

total_rep = np.array([1, 200, 500, 1000, 2000, 3000, 5000, 10000,
                      25000, 50000, 100000])
num_users = np.array([4364226, 269110, 158824, 90368, 48609, 32604, 18921,
                      8618, 2802, 1000, 334])

ccdf = num_users.astype(float) / num_users.max()

# Fit a line with a y-intercept of 1 in log-space
logx = np.log(total_rep)
logy = np.log(ccdf)
slope, _, _, _ = np.linalg.lstsq(logx[:,np.newaxis], logy)

params = [slope, 0]
est = np.exp(np.polyval(params, logx))

fig, ax = plt.subplots()
ax.loglog(total_rep, ccdf, color='lightblue', ls='', marker='o',
          clip_on=False, zorder=10, label='Observations')

ax.plot(total_rep, est, color='salmon', label='Fit', ls='--')

ax.set(xlabel='Reputation', title='CCDF of Stackoverflow Reputation',
       ylabel='Probability that Reputation is Greater than X')

plt.show()

Хммм... Теперь у нас есть новая проблема. Наша новая линия не очень подходит нашим данным. Это общая проблема с степенными распределениями.

Используйте только "хвосты" в подгонке

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

Тот факт, что только хвост распределения следует степенному закону, объясняет, почему мы не смогли хорошо подобрать наши данные, когда мы зафиксировали y-перехват.

Существует множество различных модифицированных моделей степенного права для того, что происходит вблизи вероятности 1, но все они следуют степенному правилу справа от некоторого значения отсечки. Основываясь на наших наблюдаемых данных, похоже, что мы могли бы соответствовать двум строкам: один справа от репутации ~ 1000 и один слева.

Имея это в виду, позвольте забыть о левой стороне вещей и сосредоточьтесь на "длинном хвосте" справа. Мы будем использовать np.polyfit, но исключаем левые три точки из подгонки.

import numpy as np
import matplotlib.pyplot as plt

total_rep = np.array([1, 200, 500, 1000, 2000, 3000, 5000, 10000,
                      25000, 50000, 100000])
num_users = np.array([4364226, 269110, 158824, 90368, 48609, 32604, 18921,
                      8618, 2802, 1000, 334])

ccdf = num_users.astype(float) / num_users.max()

# Fit a line in log-space, excluding reputation <= 1000
logx = np.log(total_rep[total_rep > 1000])
logy = np.log(ccdf[total_rep > 1000])

params = np.polyfit(logx, logy, 1)
est = np.exp(np.polyval(params, logx))

fig, ax = plt.subplots()
ax.loglog(total_rep, ccdf, color='lightblue', ls='', marker='o',
          clip_on=False, zorder=10, label='Observations')

ax.plot(total_rep[total_rep > 1000], est, color='salmon', label='Fit', ls='--')

ax.set(xlabel='Reputation', title='CCDF of Stackoverflow Reputation',
       ylabel='Probability that Reputation is Greater than X')

plt.show()

Протестируйте различные настройки

В этом случае у нас есть некоторые дополнительные данные. Давайте посмотрим, насколько хорошо каждая различная посадка предсказывает рейтинг репутации топ-5:

import numpy as np
import matplotlib.pyplot as plt

total_rep = np.array([1, 200, 500, 1000, 2000, 3000, 5000, 10000,
                      25000, 50000, 100000])
num_users = np.array([4364226, 269110, 158824, 90368, 48609, 32604, 18921,
                      8618, 2802, 1000, 334])

top_5_rep = [832131, 632105, 618926, 596889, 576697]
top_5_ccdf = np.array([1, 2, 3, 4, 5], dtype=float) / num_users.max()

ccdf = num_users.astype(float) / num_users.max()

# Previous fits
naive_params = [-0.81938338,  1.15955974]
fixed_intercept_params = [-0.68845134, 0]
long_tail_params = [-1.26172528, 5.24883471]

fits = [naive_params, fixed_intercept_params, long_tail_params]
fit_names = ['Naive Fit', 'Fixed Intercept Fit', 'Long Tail Fit']


fig, ax = plt.subplots()
ax.loglog(total_rep, ccdf, color='lightblue', ls='', marker='o',
          clip_on=False, zorder=10, label='Observations')

# Plot reputation of top 5 users
ax.loglog(top_5_rep, top_5_ccdf, ls='', marker='o', color='darkred',
          zorder=10, label='Top 5 Users')

# Plot different fits
for params, name in zip(fits, fit_names):
    x = [1, 1e7]
    est = np.exp(np.polyval(params, np.log(x)))
    ax.loglog(x, est, label=name, ls='--')

ax.set(xlabel='Reputation', title='CCDF of Stackoverflow Reputation',
       ylabel='Probability that Reputation is Greater than X',
       ylim=[1e-7, 1])
ax.legend()

plt.show()

Ничего себе! Они все делают довольно ужасную работу! Во-первых, это хорошая причина, чтобы использовать полную серию при установке дистрибутива, а не только биндинговых данных. Однако корень проблемы заключается в том, что распределение степенного закона в этом случае не очень хорошо подходит. На первый взгляд, похоже, что экспоненциальное распределение может быть лучше, но оставьте это позже.

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

import numpy as np

# Jon Skeet actual reputation
skeet_prob = 1.0 / 4364226
true_rep = 832131

# Previous fits
naive_params = [-0.81938338,  1.15955974]
fixed_intercept_params = [-0.68845134, 0]
long_tail_params = [-1.26172528, 5.24883471]

fits = [naive_params, fixed_intercept_params, long_tail_params]
fit_names = ['Naive Fit', 'Fixed Intercept Fit', 'Long Tail Fit']

for params, name in zip(fits, fit_names):
    inv_params = [1 / params[0], -params[1]/params[0]]
    est = np.exp(np.polyval(inv_params, np.log(skeet_prob)))

    print '{}:'.format(name)
    print '    Pred. Rep.: {}'.format(est)
    print ''

print 'True Reputation: {}'.format(true_rep)

Это дает:

Naive Fit:
    Pred. Rep.: 522562573.099

Fixed Intercept Fit:
    Pred. Rep.: 4412664023.88

Long Tail Fit:
    Pred. Rep.: 11728612.2783

True Reputation: 832131

Ответ 3

Прочитав отличные объяснения Джо Кингтона и Джоса Полфлита, я решил, что по моим данным, 5 точек данных из хвоста дистрибутива (включая главного пользователя), чтобы узнать, могу ли я найти один, хороший достаточно поместиться, используя только полиномиальное соответствие.

Оказывается, что полином 6-й степени отлично работает в центре и в хвосте дистрибутива с меньшими шагами.

В приведенной ниже таблице показаны данные и полиномиальное соответствие, которое кажется почти идеальным:

Добавление конечных точек данных

Это мой df с некоторыми дополнительными точками данных из хвоста дистрибутива:

0           1  4364226
1         200   269110
2         500   158824
3        1000    90368
4        2000    48609
5        3000    32604
6        5000    18921
7       10000     8618
8       25000     2802
9       50000     1000
10     100000      334
11     193000      100
12     261000       50
13     441000       10
14     578000        5
15     833000        1

Это мой код:

soDF['log_total_rep'] = soDF['total_rep'].apply(np.log10)
soDF['log_users']     = soDF['users'].apply(np.log10)
coefficients = np.polyfit(soDF['log_total_rep'] , soDF['log_users'], 6)
polynomial = np.poly1d(coefficients)
print polynomial

Что возвращает это:

          6           5          4          3          2
-0.00258 x + 0.04187 x - 0.2541 x + 0.6774 x - 0.7697 x - 0.2513 x + 6.64

График выполняется с помощью этого кода:

xp = np.linspace(0, 6, 100)
plt.figure(figsize=(18,6))
plt.title('Stackoverflow Reputation', fontsize =15)
plt.xlabel('Log reputation', fontsize =15)
plt.ylabel('Log probability that reputation is greater than X', fontsize = 15)
plt.plot(soDF['log_total_rep'], soDF['log_users'],'o', label ='Data')
plt.plot(xp, polynomial(xp), color='red', label='Fit', ls='--')
plt.legend(loc='upper right', fontsize = 15)

Тестирование параметрического соответствия

Чтобы проверить соответствие в центре и в хвостах, я выбираю следующие профили для пользователей с рангом 150, 25 и 5:

Это мой код:

total_users = 4407194
def predicted_rank(total_rep):
    parametric_rank_position   = 10**polynomial(np.log10(total_rep))
    parametric_rank_percentile = parametric_rank_position/total_users
    print "Position is " + str(int(parametric_rank_position)) + ", and rank is top " +  "{:.4%}".format(parametric_rank_percentile)

Итак, для Иоахима Зауера это результат:

predicted_rank(165671)
Position is 133, and rank is top 0.0030%

Отключено на 17 позиций. Для Эрика Липперта:

predicted_rank(374507)
Position is 18, and rank is top 0.0004%

Выкл. на 7 позиций. Для Марка Гравелла:

predicted_rank(579042)
Position is 4, and rank is top 0.0001%

Отключено на 1 позицию. Чтобы проверить центр распространения, я тестирую свой собственный:

predicted_rank(1242)
Position is 75961, and rank is top 1.7236%

который приближается к реальному ранга 75630.