plfit/plfit_py.py
def plfit(self,nosmall=True,finite=False,quiet=False,silent=False,
              xmin=None, verbose=False):
        """
        A pure-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

        nosmall is on by default; it rejects low s/n points
        can specify xmin to skip xmin estimation

        This is only for continuous distributions; I have not implemented a
        pure-python discrete distribution fitter
        """
        x = self.data
        z = sorted(x)
        t = time.time()
        possible_xmins = sorted(set(z))
        argxmins = [z.index(i) for i in possible_xmins]
        self._nunique = len(possible_xmins)
        if xmin is None:
            av  = map(self.alpha_(z),possible_xmins)
            dat = map(self.kstest_(z),possible_xmins)
            sigma = [(a-1)/math.sqrt(len(z)-i+1) for a,i in zip(av,argxmins)]
            if nosmall:
                # test to make sure the number of data points is high enough
                # to provide a reasonable s/n on the computed alpha
                goodvals = [s<0.1 for s in sigma]
                if False in goodvals: 
                    nmax = goodvals.index(False)
                    dat = dat[:nmax]
                    possible_xmins = possible_xmins[:nmax]
                    av = av[:nmax]
                else:
                    print("Not enough data left after flagging - using all positive data.")
            if not quiet: print("PYTHON plfit executed in %f seconds" % (time.time()-t))
            self._av = av
            self._xmin_kstest = dat
            self._sigma = sigma
            # [:-1] to weed out the very last data point; it cannot be correct
            # (can't have a power law with 1 data point).
            # However, this should only be done if the ends have not previously
            # been excluded with nosmall
            if nosmall:
                xmin = possible_xmins[dat.index(min(dat))] 
            else:
                xmin = possible_xmins[dat.index(min(dat[:-1]))] 
        z     = [i for i in z if i >= xmin]
        n     = len(z)
        alpha = 1 + n / sum([math.log(a/xmin) for a in z]) 
        if finite:
            alpha = alpha*(n-1.)/n+1./n
        if n == 1 and not silent:
            print("Failure: only 1 point kept.  Probably not a power-law distribution.")
            self._alpha = 0
            self._alphaerr = 0
            self._likelihood = 0
            self._ks = 0
            self._ks_prob = 0
            self._xmin = xmin
            return xmin,0
        if n < 50 and not finite and not silent:
            print('(PLFIT) Warning: finite-size bias may be present. n=%i' % n)
        # ks = max(abs( numpy.arange(n)/float(n) - (1-(xmin/z)**(alpha-1)) ))
        ks = max( [abs( i/float(n) - (1-(xmin/b)**(alpha-1))) for i,b in zip(range(n),z)] )
        # Parallels Eqn 3.5 in Clauset et al 2009, but zeta(alpha, xmin) = (alpha-1)/xmin.  Really is Eqn B3 in paper.
        #L = n*log((alpha-1)/xmin) - alpha*sum(log(z/xmin))
        sl = sum([math.log(a/xmin) for a in z])
        L = (n*math.log((alpha-1)/xmin) - alpha*sl)
        #requires another map... Larr = arange(len(unique(x))) * log((av-1)/unique(x)) - av*sum
        self._likelihood = L
        self._xmin = xmin
        self._xmins = possible_xmins
        self._alpha= alpha
        self._alphaerr = (alpha-1)/math.sqrt(n)
        self._ks = ks  # this ks statistic may not have the same value as min(dat) because of unique()
        #if scipyOK: self._ks_prob = scipy.stats.kstwobign.sf(ks*numpy.sqrt(n))
        self._ngtx = n
        if math.isnan(L) or math.isnan(xmin) or math.isnan(alpha):
            raise ValueError("plfit failed; returned a nan")

        if not quiet:
            if verbose: print("The lowest value included in the power-law fit, ", end=' ')
            print("xmin: %g" % xmin, end=' ')
            if verbose: print("\nThe number of values above xmin, ", end=' ')
            print("n(>xmin): %i" % n, end=' ')
            if verbose: print("\nThe derived power-law alpha (p(x)~x^-alpha) with MLE-derived error, ", end=' ')
            print("alpha: %g +/- %g  " % (alpha,self._alphaerr), end=' ') 
            if verbose: print("\nThe log of the Likelihood (the maximized parameter), ", end=' ')
            print("Log-Likelihood: %g  " % L, end=' ')
            if verbose: print("\nThe KS-test statistic between the best-fit power-law and the data, ", end=' ')
            print("ks: %g" % (ks))

        return xmin,alpha
def plfit(self,nosmall=True,finite=False,quiet=False,silent=False,
              xmin=None, verbose=False):
        """
        A pure-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

        nosmall is on by default; it rejects low s/n points
        can specify xmin to skip xmin estimation

        This is only for continuous distributions; I have not implemented a
        pure-python discrete distribution fitter
        """
        x = self.data
        z = sorted(x)
        t = time.time()
        possible_xmins = sorted(set(z))
        argxmins = [z.index(i) for i in possible_xmins]
        self._nunique = len(possible_xmins)
        if xmin is None:
            av  = map(self.alpha_(z),possible_xmins)
            dat = map(self.kstest_(z),possible_xmins)
            sigma = [(a-1)/math.sqrt(len(z)-i+1) for a,i in zip(av,argxmins)]
            if nosmall:
                # test to make sure the number of data points is high enough
                # to provide a reasonable s/n on the computed alpha
                goodvals = [s<0.1 for s in sigma]
                if False in goodvals: 
                    nmax = goodvals.index(False)
                    dat = dat[:nmax]
                    possible_xmins = possible_xmins[:nmax]
                    av = av[:nmax]
                else:
                    print("Not enough data left after flagging - using all positive data.")
            if not quiet: print("PYTHON plfit executed in %f seconds" % (time.time()-t))
            self._av = av
            self._xmin_kstest = dat
            self._sigma = sigma
            # [:-1] to weed out the very last data point; it cannot be correct
            # (can't have a power law with 1 data point).
            # However, this should only be done if the ends have not previously
            # been excluded with nosmall
            if nosmall:
                xmin = possible_xmins[dat.index(min(dat))] 
            else:
                xmin = possible_xmins[dat.index(min(dat[:-1]))] 
        z     = [i for i in z if i >= xmin]
        n     = len(z)
        alpha = 1 + n / sum([math.log(a/xmin) for a in z]) 
        if finite:
            alpha = alpha*(n-1.)/n+1./n
        if n == 1 and not silent:
            print("Failure: only 1 point kept.  Probably not a power-law distribution.")
            self._alpha = 0
            self._alphaerr = 0
            self._likelihood = 0
            self._ks = 0
            self._ks_prob = 0
            self._xmin = xmin
            return xmin,0
        if n < 50 and not finite and not silent:
            print('(PLFIT) Warning: finite-size bias may be present. n=%i' % n)
        # ks = max(abs( numpy.arange(n)/float(n) - (1-(xmin/z)**(alpha-1)) ))
        ks = max( [abs( i/float(n) - (1-(xmin/b)**(alpha-1))) for i,b in zip(range(n),z)] )
        # Parallels Eqn 3.5 in Clauset et al 2009, but zeta(alpha, xmin) = (alpha-1)/xmin.  Really is Eqn B3 in paper.
        #L = n*log((alpha-1)/xmin) - alpha*sum(log(z/xmin))
        sl = sum([math.log(a/xmin) for a in z])
        L = (n*math.log((alpha-1)/xmin) - alpha*sl)
        #requires another map... Larr = arange(len(unique(x))) * log((av-1)/unique(x)) - av*sum
        self._likelihood = L
        self._xmin = xmin
        self._xmins = possible_xmins
        self._alpha= alpha
        self._alphaerr = (alpha-1)/math.sqrt(n)
        self._ks = ks  # this ks statistic may not have the same value as min(dat) because of unique()
        #if scipyOK: self._ks_prob = scipy.stats.kstwobign.sf(ks*numpy.sqrt(n))
        self._ngtx = n
        if math.isnan(L) or math.isnan(xmin) or math.isnan(alpha):
            raise ValueError("plfit failed; returned a nan")

        if not quiet:
            if verbose: print("The lowest value included in the power-law fit, ", end=' ')
            print("xmin: %g" % xmin, end=' ')
            if verbose: print("\nThe number of values above xmin, ", end=' ')
            print("n(>xmin): %i" % n, end=' ')
            if verbose: print("\nThe derived power-law alpha (p(x)~x^-alpha) with MLE-derived error, ", end=' ')
            print("alpha: %g +/- %g  " % (alpha,self._alphaerr), end=' ') 
            if verbose: print("\nThe log of the Likelihood (the maximized parameter), ", end=' ')
            print("Log-Likelihood: %g  " % L, end=' ')
            if verbose: print("\nThe KS-test statistic between the best-fit power-law and the data, ", end=' ')
            print("ks: %g" % (ks))

        return xmin,alpha
def plfit ( self , nosmall = True , finite = False , quiet = False , silent = False , xmin = None , verbose = False ) : x = self . data z = sorted ( x ) t = time . time ( ) possible_xmins = sorted ( set ( z ) ) argxmins = [ z . index ( i ) for i in possible_xmins ] self . _nunique = len ( possible_xmins ) if xmin is None : av = map ( self . alpha_ ( z ) , possible_xmins ) dat = map ( self . kstest_ ( z ) , possible_xmins ) sigma = [ ( a - 1 ) / math . sqrt ( len ( z ) - i + 1 ) for a , i in zip ( av , argxmins ) ] if nosmall : # test to make sure the number of data points is high enough # to provide a reasonable s/n on the computed alpha goodvals = [ s < 0.1 for s in sigma ] if False in goodvals : nmax = goodvals . index ( False ) dat = dat [ : nmax ] possible_xmins = possible_xmins [ : nmax ] av = av [ : nmax ] else : print ( "Not enough data left after flagging - using all positive data." ) if not quiet : print ( "PYTHON plfit executed in %f seconds" % ( time . time ( ) - t ) ) self . _av = av self . _xmin_kstest = dat self . _sigma = sigma # [:-1] to weed out the very last data point; it cannot be correct # (can't have a power law with 1 data point). # However, this should only be done if the ends have not previously # been excluded with nosmall if nosmall : xmin = possible_xmins [ dat . index ( min ( dat ) ) ] else : xmin = possible_xmins [ dat . index ( min ( dat [ : - 1 ] ) ) ] z = [ i for i in z if i >= xmin ] n = len ( z ) alpha = 1 + n / sum ( [ math . log ( a / xmin ) for a in z ] ) if finite : alpha = alpha * ( n - 1. ) / n + 1. / n if n == 1 and not silent : print ( "Failure: only 1 point kept.  Probably not a power-law distribution." ) self . _alpha = 0 self . _alphaerr = 0 self . _likelihood = 0 self . _ks = 0 self . _ks_prob = 0 self . _xmin = xmin return xmin , 0 if n < 50 and not finite and not silent : print ( '(PLFIT) Warning: finite-size bias may be present. n=%i' % n ) # ks = max(abs( numpy.arange(n)/float(n) - (1-(xmin/z)**(alpha-1)) )) ks = max ( [ abs ( i / float ( n ) - ( 1 - ( xmin / b ) ** ( alpha - 1 ) ) ) for i , b in zip ( range ( n ) , z ) ] ) # Parallels Eqn 3.5 in Clauset et al 2009, but zeta(alpha, xmin) = (alpha-1)/xmin.  Really is Eqn B3 in paper. #L = n*log((alpha-1)/xmin) - alpha*sum(log(z/xmin)) sl = sum ( [ math . log ( a / xmin ) for a in z ] ) L = ( n * math . log ( ( alpha - 1 ) / xmin ) - alpha * sl ) #requires another map... Larr = arange(len(unique(x))) * log((av-1)/unique(x)) - av*sum self . _likelihood = L self . _xmin = xmin self . _xmins = possible_xmins self . _alpha = alpha self . _alphaerr = ( alpha - 1 ) / math . sqrt ( n ) self . _ks = ks # this ks statistic may not have the same value as min(dat) because of unique() #if scipyOK: self._ks_prob = scipy.stats.kstwobign.sf(ks*numpy.sqrt(n)) self . _ngtx = n if math . isnan ( L ) or math . isnan ( xmin ) or math . isnan ( alpha ) : raise ValueError ( "plfit failed; returned a nan" ) if not quiet : if verbose : print ( "The lowest value included in the power-law fit, " , end = ' ' ) print ( "xmin: %g" % xmin , end = ' ' ) if verbose : print ( "\nThe number of values above xmin, " , end = ' ' ) print ( "n(>xmin): %i" % n , end = ' ' ) if verbose : print ( "\nThe derived power-law alpha (p(x)~x^-alpha) with MLE-derived error, " , end = ' ' ) print ( "alpha: %g +/- %g  " % ( alpha , self . _alphaerr ) , end = ' ' ) if verbose : print ( "\nThe log of the Likelihood (the maximized parameter), " , end = ' ' ) print ( "Log-Likelihood: %g  " % L , end = ' ' ) if verbose : print ( "\nThe KS-test statistic between the best-fit power-law and the data, " , end = ' ' ) print ( "ks: %g" % ( ks ) ) return xmin , alpha
A pure-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

        nosmall is on by default; it rejects low s/n points
        can specify xmin to skip xmin estimation

        This is only for continuous distributions; I have not implemented a
        pure-python discrete distribution fitter
A pure - Python implementation of the Matlab code http : // www . santafe . edu / ~aaronc / powerlaws / plfit . m from http : // www . santafe . edu / ~aaronc / powerlaws /