alpha=2.5
xmin=1
n=1000
xr=xmin*rand(n)**(-1/(alpha-1))

figure(figsize=(6,5))
yh,bh,p=hist(xr,bins=arange(xmin,1000,.1),log=True,histtype='step')
xlabel('x')
ylabel('# occurrences in bin of size .1')
ylim(.5,300),xlim(1,200)
xscale('log');

from scipy.stats import linregress
slope, intercept, r_value, p_value, std_err = linregress(log(bh[yh>0]),log(yh[yh>0]))
print slope,intercept

figure(figsize=(6,5))
yh,bh,p=hist(xr,bins=arange(xmin,1000,.1),log=True,histtype='step')
xlabel('x')
ylabel('# occurrences in bin of size .1')
ylim(.5,300),xlim(1,200)
xscale('log')

xl=arange(1,200,.1)
plot(xl,exp(intercept)*xl**(slope),'k-',linewidth=2,label='$\\alpha={:.2f}$'.format(-slope))

legend();

figure(figsize=(6,5))
yh,bh,p=hist(xr,bins=arange(xmin,1000,.1),log=True,histtype='step')
xlabel('x')
ylabel('# occurrences in bin of size .1')
ylim(.5,300),xlim(1,200)
xscale('log')

slope, intercept = linregress(log(bh[yh>0]),log(yh[yh>0]))[:2]
xl=arange(1,200,.1)
plot(xl,exp(intercept)*xl**(slope),'k-',linewidth=2,
     label='$\\alpha={:.2f}$ (all)'.format(-slope))

for xmax in (64,32,16,8,4,2):
    retain=logical_and(yh>0,bh[:-1]<=xmax)
    slope,intercept = linregress(log(bh[retain]),log(yh[retain]))[:2]
    xl=arange(1,xmax+1,xmax-1)
    plot(xl,exp(intercept)*xl**(slope),linewidth=1,
         label='$\\alpha={:.2f}$ ({})'.format(-slope,xmax));
legend();

alpha=1+n/sum(log(xr))
sigma=(alpha-1)/sqrt(n)
print 'alpha=',round(alpha,3),', with stderr sigma=',round(sigma,3)

C=(alpha-1)*xmin**(alpha-1)

figure(figsize=(6,5))
yh,bh,p=hist(xr,bins=arange(xmin,1000,.1),log=True,histtype='step')
xlabel('x')
ylabel('# occurrences in bin of size .1')
ylim(.5,300),xlim(1,200)
xscale('log')

xl=arange(1,200,.1)
plot(xl,(.1*n)*C*xl**(-alpha),'k-',linewidth=2,label='$\\alpha={:.2f}$'.format(alpha))

legend();

alpha=2.5
xmin=1
n=1000000
xr=xmin*rand(n)**(-1/(alpha-1))

alpha=1+n/sum(log(xr))
sigma=(alpha-1)/sqrt(n)
print 'alpha=',round(alpha,4),', with stderr sigma=',round(sigma,4)
C=(alpha-1)*xmin**(alpha-1)

figure(figsize=(6,5))
yh,bh,p=hist(xr,bins=arange(xmin,20000,.1),log=True,histtype='step')
xlabel('x')
ylabel('# occurrences in bin of size .1')
ylim(.5,150000),xlim(1,20000)
xscale('log')
xl=arange(1,200,.1)
plot(xl,(.1*n)*C*xl**(-alpha),'r-',linewidth=2,label='$\\alpha={:.4f}$'.format(alpha))
title('n={:,} data points'.format(n))
legend();