def plfit(x,nosmall=False,finite=False):
xmins = unique(x)<EOL>xmins = xmins[<NUM_LIT:1>:-<NUM_LIT:1>]<EOL>dat = xmins * <NUM_LIT:0> <EOL>z = sort(x)<EOL>for xm in arange(len(xmins)):<EOL><INDENT>xmin = xmins[xm]<EOL>z    = z[z>=xmin] <EOL>n    = float(len(z))<EOL>a    =  n / sum( log(z/xmin) )<EOL>if nosmall:<EOL><INDENT>if (a-<NUM_LIT:1>)/sqrt(n) > <NUM_LIT:0.1>:<EOL><INDENT>dat = dat[:xm]<EOL>xm = len(xmins)+<NUM_LIT:1><EOL>break<EOL><DEDENT><DEDENT>cx   = arange(n)/float(n)  <EOL>cf   = <NUM_LIT:1>-(xmin/z)**a  <EOL>dat[xm] = max( abs(cf-cx) )<EOL><DEDENT>D     = min(dat);<EOL>xmin  = xmins[argmin(dat)]<EOL>z     = x[x>=xmin]<EOL>n     = len(z)<EOL>alpha = <NUM_LIT:1> + n / sum( log(z/xmin) )<EOL>if finite:<EOL><INDENT>alpha = alpha*(n-<NUM_LIT:1>)/n+<NUM_LIT:1>/n<EOL><DEDENT>if n < <NUM_LIT:50> and ~finite:<EOL><INDENT>print('<STR_LIT>')<EOL><DEDENT>L = n*log((alpha-<NUM_LIT:1>)/xmin) - alpha*sum(log(z/xmin));<EOL>return xmin,alpha,L,dat<EOL>
A Python implementation of the Matlab code http://www.santafe.edu/~aaronc/powerlaws/plfit.m
from http://www.santafe.edu/~aaronc/powerlaws/

See A. Clauset, C.R. Shalizi, and M.E.J. Newman, "Power-law distributions
in empirical data" SIAM Review, to appear (2009). (arXiv:0706.1062)
http://arxiv.org/abs/0706.1062