使用pandas和numpy来参数化堆栈溢出的用户数和声誉

Question

我注意到Stack Overflow的用户数量及其声誉遵循有趣的分布.我创建了一个pandas DF,看看我是否可以创建一个参数拟合:

import pandas as pd
import numpy as np
soDF = pd.read_excel('scores.xls')
print soDF

哪个返回:

    total_rep    users
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

如果我绘制图表,我会得到以下图表:

该分布似乎遵循幂律.因此,为了更好地可视化,我添加了以下内容:

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')

其中产生了以下内容:

有没有一种简单的方法可以使用pandas来找到最适合这些数据的方法？虽然拟合看起来是线性的,但也许多项式拟合更好,因为现在我处理的是对数尺度.

Answer 1

python,pandas以及scipy,噢,我的!

科学的python生态系统有几个免费的库.按设计,没有一个图书馆可以做任何事情 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()

这些数据代表什么？

接下来,我们需要知道我们正在合作的是什么.我们绘制的内容类似于直方图,因为它是给定信誉级别的用户数的原始计数.但是,请注意每个bin旁边的小"+"信誉表.这意味着,例如,2082个用户的信誉得分为25000 或更高.

我们的数据基本上是互补累积分布函数(CCDF)的估计,与直方图是概率分布函数(PDF)的估计相同.我们只需要通过我们样本中的用户总数对其进行标准化,以估算CCDF.在这种情况下,我们可以简单地除以第一个元素num_users.信誉永远不会小于1,因此x轴上的1对应于定义的概率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用户的总数(以及每个信誉级别的数字)将有所不同.但是,任何给定用户具有特定声誉的总概率不会发生显着变化.如果我们想要预测John Skeet的声誉(最高代表用户),当该网站达到500万注册用户时,使用概率而不是原始计数要容易得多.

天权适合幂律分布

接下来,让我们为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()

哇!他们都做得非常糟糕!首先,这是在拟合分布而不仅仅是分箱数据时使用完整系列的一个很好的理由.然而,问题的根源在于,在这种情况下,幂律分布不是很合适.乍一看,它看起来像指数分布可能更合适,但让我们留待以后.

作为一个例子,不同的幂律适合过度预测低概率观察(即具有最高代表的用户),让我们预测Jon Skeet对每个模型的声誉:

import numpy as np

# Jon Skeet's 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

Answer 2

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')

拟合二次多项式

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

print "Coefficients: ", coefficients

接下来,让我们绘制原始+拟合:

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

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