Estimating the PDF
By taking the derivative, you are estimating the pdf from the cdf. A more straightforward estimate of the pdf is to slightly alter the query to make a count by reputation.
Like the cdf query you have, this data needs to be smoothed (specifically, binned) before we can estimate the distribution. A telltale sign that the data is not smoothed is that the pdf estimate in your plots are not one-to-one. That is, there are multiple scatter plot points for a given value of $\log(R)$. This means we aren't really plotting the pdf and that it doesn't make sense to fit this data to infer whether we see a power law.
What we are seeing in both of the unbinned dataset (yours above and mine below) are reputation values in that appear much more rare than they are. For instance, only one user has a reputation of 729, but that doesn't mean that a 729 rep is a 1/283,128 event (there are 283,128 users) . In fact a reputation of 729 is much closer to being a 1/280 event.
Binned Data
I took a crack at binning the data and then regressing with a linear and quadratic fit (in log-base-10 scale). A linear fit implies a power law distribution, I am not sure what a quadratic fit implies other than a pdf of the form $f(x)=b e^{-a \log^2(x)}$. For comparison, I also plotted the raw reputation vs user count in blue, which is the un-binned data that is similar to the data you are plotting.
Stack Overflow Data
I repeated this for Stack Overflow user data with this query. SO has 1.5 orders of magnitude more users, so the trend is a bit more clear.
Code
The python code to make these plots in a Jupyter Notebook is:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
%matplotlib inline
def plot_hist(df,alpha,c,label, ax=None, limits = [-1,7,-8,1], trends=False):
ax = df.plot(x='log_reputation',
y=['log_user_proportion'],
kind='scatter',
alpha=alpha,
ax=ax,
color=c,edgecolor=c,
label=label,
figsize=(15.0, 10.0))
if trends:
X = df.log_reputation
X = sm.add_constant(X)
model = sm.OLS(df.log_user_proportion, X)
results = model.fit()
df['linear_trend'] = results.fittedvalues
plt.text(df.log_reputation.min(),df.linear_trend.max(),
'{log_reputation:10.4f}x +{const:10.4f}'.format(**results.params)+
'\nr**2 = {}'.format(results.rsquared))
X = np.column_stack((df.log_reputation, df.log_reputation**2))
X = sm.add_constant(X)
model = sm.OLS(df.log_user_proportion, X)
results = model.fit()
df['quad_trend'] = results.fittedvalues
plt.text(3,-2,
'{x2:10.4f}x**2 +{x1:10.4f}x +{const:10.4f}'.format(**results.params)+
'\nr**2 = {}'.format(results.rsquared))
ax = df.plot(x='log_reputation',
y=['linear_trend','quad_trend'],
kind='line',
alpha=1,
ax=ax)
plt.axis(limits)
plt.grid(True)
return ax
def process_df(url):
df = pd.read_csv(url,names=['Reputation','UserCount'],dtype=float,skiprows=1)
df['user_proportion'] = df.UserCount/sum(df.UserCount)
df['log_user_proportion'] = np.log10(df.user_proportion)
df['log_reputation'] = np.log10(df.Reputation)
ax = plot_hist(df,.1,c='b',label='un-binned')
cut_points = np.power(10,np.arange(0,8,.1))
bins = pd.cut(df.Reputation,cut_points,right=False)
df = df.groupby(bins).user_proportion.sum().reset_index()
df['cut_points'] = cut_points[:-1]
df.dropna(inplace=True)
df['log_reputation'] = np.log10(df.cut_points)
df['log_user_proportion'] = np.log10(df.user_proportion)
plot_hist(df,1,c='m',ax=ax,label='binned',trends=True)
plt.legend()
return df
df = process_df('http://data.stackexchange.com/math/csv/715124')
plt.title('Mathematics')
df2 = process_df('https://data.stackexchange.com/stackoverflow/csv/715082')
plt.title('Stack Overflow')
