def plfit(self, nosmall=True, finite=False, quiet=False, silent=False,<EOL>usefortran=False, usecy=False, xmin=None, verbose=False,<EOL>discrete=None, discrete_approx=True, discrete_n_alpha=<NUM_LIT:1000>,<EOL>skip_consistency_check=False):
x = self.data<EOL>if any(x < <NUM_LIT:0>):<EOL><INDENT>raise ValueError("<STR_LIT>"<EOL>"<STR_LIT>"<EOL>"<STR_LIT>")<EOL><DEDENT>z = np.sort(x)<EOL>xmins,argxmins = np.unique(z,return_index=True)<EOL>self._nunique = len(xmins)<EOL>if self._nunique == len(x) and discrete is None:<EOL><INDENT>if verbose:<EOL><INDENT>print("<STR_LIT>"<EOL>"<STR_LIT>")<EOL><DEDENT>discrete = False<EOL><DEDENT>elif self._nunique < len(x) and discrete is None:<EOL><INDENT>if verbose:<EOL><INDENT>print("<STR_LIT>"<EOL>"<STR_LIT>")<EOL><DEDENT>discrete = True<EOL><DEDENT>t = time.time()<EOL>if xmin is None:<EOL><INDENT>if discrete:<EOL><INDENT>self.discrete_best_alpha(approximate=discrete_approx,<EOL>n_alpha=discrete_n_alpha,<EOL>verbose=verbose,<EOL>finite=finite)<EOL>return self._xmin,self._alpha<EOL><DEDENT>elif usefortran and fortranOK:<EOL><INDENT>kstest_values,alpha_values = fplfit.plfit(z, <NUM_LIT:0>)<EOL>if not quiet:<EOL><INDENT>print(("<STR_LIT>" % (time.time()-t)))<EOL><DEDENT><DEDENT>elif usecy and cyOK:<EOL><INDENT>kstest_values,alpha_values = cplfit.plfit_loop(z,<EOL>nosmall=False,<EOL>zunique=xmins,<EOL>argunique=argxmins)<EOL>if not quiet:<EOL><INDENT>print(("<STR_LIT>" % (time.time()-t)))<EOL><DEDENT><DEDENT>else:<EOL><INDENT>f_alpha = alpha_gen(z)<EOL>f_kstest = kstest_gen(z)<EOL>alpha_values = np.asarray(list(map(f_alpha,xmins)),<EOL>dtype='<STR_LIT:float>')<EOL>kstest_values = np.asarray(list(map(f_kstest,xmins)),<EOL>dtype='<STR_LIT:float>')<EOL>if not quiet:<EOL><INDENT>print(("<STR_LIT>" % (time.time()-t)))<EOL><DEDENT><DEDENT>if not quiet:<EOL><INDENT>if usefortran and not fortranOK:<EOL><INDENT>raise ImportError("<STR_LIT>")<EOL><DEDENT>if usecy and not cyOK:<EOL><INDENT>raise ImportError("<STR_LIT>")<EOL><DEDENT><DEDENT>sigma = (alpha_values-<NUM_LIT:1>)/np.sqrt(len(z)-argxmins)<EOL>if nosmall:<EOL><INDENT>goodvals = sigma<<NUM_LIT:0.1><EOL>nmax = argmin(goodvals)<EOL>if nmax <= <NUM_LIT:0>:<EOL><INDENT>nmax = len(xmins) - <NUM_LIT:1><EOL>if not silent:<EOL><INDENT>print("<STR_LIT>"<EOL>"<STR_LIT>"<EOL>"<STR_LIT>")<EOL><DEDENT><DEDENT><DEDENT>else:<EOL><INDENT>nmax = len(xmins)-<NUM_LIT:1><EOL><DEDENT>best_ks_index = argmin(kstest_values[:nmax])<EOL>xmin = xmins[best_ks_index]<EOL>self._alpha_values = alpha_values<EOL>self._xmin_kstest = kstest_values<EOL>if scipyOK:<EOL><INDENT>self._ks_prob_all = np.array([scipy.stats.ksone.sf(D_stat,<EOL>len(kstest_values)-ii)<EOL>for ii,D_stat in<EOL>enumerate(kstest_values)])<EOL><DEDENT>self._sigma = sigma<EOL>n = np.count_nonzero(z>=xmin)<EOL>alpha = <NUM_LIT:1.> + float(n)/sum(log(z[z>=xmin]/xmin))<EOL>try:<EOL><INDENT>if not skip_consistency_check:<EOL><INDENT>np.testing.assert_almost_equal(alpha,<EOL>alpha_values[best_ks_index],<EOL>decimal=<NUM_LIT:4>)<EOL><DEDENT><DEDENT>except AssertionError:<EOL><INDENT>raise AssertionError("<STR_LIT>"<EOL>"<STR_LIT>"<EOL>"<STR_LIT>"<EOL>"<STR_LIT>"<EOL>"<STR_LIT>"<EOL>"<STR_LIT>"<EOL>"<STR_LIT>"<EOL>.format(alpha,<EOL>alpha_values[best_ks_index]))<EOL><DEDENT><DEDENT>z = z[z>=xmin]<EOL>n = len(z)<EOL>alpha = <NUM_LIT:1.> + float(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 not finite and not silent:<EOL><INDENT>print(('<STR_LIT>' % n))<EOL><DEDENT>ks = max(abs( np.arange(n)/float(n) - (<NUM_LIT:1>-(xmin/z)**(alpha-<NUM_LIT:1>)) ))<EOL>L = n*log((alpha-<NUM_LIT:1>)/xmin) - alpha*sum(log(z/xmin))<EOL>self._likelihood = L<EOL>self._xmin = xmin<EOL>self._xmins = xmins<EOL>self._alpha= alpha<EOL>self._alphaerr = (alpha-<NUM_LIT:1>)/np.sqrt(n)<EOL>self._ks = ks<EOL>if scipyOK:<EOL><INDENT>self._ks_prob = scipy.stats.ksone.sf(ks, n)<EOL><DEDENT>self._ngtx = n<EOL>if n == <NUM_LIT:1>:<EOL><INDENT>if not silent:<EOL><INDENT>print("<STR_LIT>")<EOL><DEDENT>self._alpha = alpha = <NUM_LIT:0><EOL>self._alphaerr = <NUM_LIT:0><EOL>self._likelihood = L = <NUM_LIT:0><EOL>self._ks = <NUM_LIT:0><EOL>self._ks_prob = <NUM_LIT:0><EOL>self._xmin = xmin<EOL>return xmin,<NUM_LIT:0><EOL><DEDENT>if np.isnan(L) or np.isnan(xmin) or np.isnan(alpha):<EOL><INDENT>raise ValueError("<STR_LIT>")<EOL><DEDENT>if not quiet:<EOL><INDENT>if verbose: print("<STR_LIT>", end='<STR_LIT:U+0020>')<EOL>print("<STR_LIT>" % xmin, end='<STR_LIT:U+0020>')<EOL>if verbose: print("<STR_LIT>", end='<STR_LIT:U+0020>')<EOL>print("<STR_LIT>" % n, end='<STR_LIT:U+0020>')<EOL>if verbose: print("<STR_LIT>", end='<STR_LIT:U+0020>')<EOL>print("<STR_LIT>" % (alpha,self._alphaerr), end='<STR_LIT:U+0020>')<EOL>if verbose: print("<STR_LIT>", end='<STR_LIT:U+0020>')<EOL>print("<STR_LIT>" % L, end='<STR_LIT:U+0020>')<EOL>if verbose: print("<STR_LIT>", end='<STR_LIT:U+0020>')<EOL>print("<STR_LIT>" % (ks), end='<STR_LIT:U+0020>')<EOL>if scipyOK:<EOL><INDENT>if verbose: print("<STR_LIT>", end='<STR_LIT:U+0020>')<EOL>print("<STR_LIT>" % (self._ks_prob))<EOL><DEDENT>else:<EOL><INDENT>print()<EOL><DEDENT><DEDENT>return xmin,alpha<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, 51, 661-703 (2009). (arXiv:0706.1062)
http://arxiv.org/abs/0706.1062

There are 3 implementations of xmin estimation.  The fortran version is
fastest, the C (cython) version is ~10% slower, and the python version
is ~3x slower than the fortran version.  Also, the cython code suffers
~2% numerical error relative to the fortran and python for unknown
reasons.

There is also a discrete version implemented in python - it is
different from the continous version!

Parameters
----------
discrete : bool or None
    If *discrete* is None, the code will try to determine whether the
    data set is discrete or continous based on the uniqueness of the
    data; if your data set is continuous but you have any non-unique
    data points (e.g., flagged "bad" data), the "automatic"
    determination will fail.  If *discrete* is True or False, the
    discrete or continuous fitter will be used, respectively.
xmin : float or int
    If you specify xmin, the fitter will only determine alpha assuming
    the given xmin; the rest of the code (and most of the complexity)
    is determining an estimate for xmin and alpha.
nosmall : bool
    When on, the code rejects low s/n points.  WARNING: This option,
    which is on by default, may result in different answers than the
    original Matlab code and the "powerlaw" python package
finite : bool
    There is a 'finite-size bias' to the estimator.  The "alpha" the
    code measures is "alpha-hat" s.t. ᾶ = (nα-1)/(n-1), or α = (1 + ᾶ
    (n-1)) / n
quiet : bool
    If False, delivers messages about what fitter is used and the fit
    results
verbose : bool
    Deliver descriptive messages about the fit parameters (only if
    `quiet==False`)
silent : bool
    If True, will print NO messages
skip_consistency_check : bool
    The code will normally perform a consistency check to make sure the
    alpha value computed by the fitter matches the alpha value computed
    directly in python.  It is possible for numerical differences to
    creep in, usually at the 10^-6 or less level.  If you see an
    exception reporting this type of error, skipping the check can be
    the appropriate next step.

Returns
-------
(xmin, alpha)
The best-fit xmin and alpha values