powerlaw Python Toolbox

When comparing two distributions where one is a specific case of the other (e.g., power_law and truncated_power_law), the p-value calculation must be adjusted. This is known as a nested comparison.

How to use:

• Automatic detection: By default, powerlaw attempts to decide if distributions are nested based on their names. If the name of one distribution is contained within the other, it may trigger nested logic.
• Explicit control: Use the nested=True keyword argument in Fit.distribution_compare or powerlaw.loglikelihood_ratio to force the nested calculation method.
• Manual comparison: You can compare distributions without a Fit object by using powerlaw.loglikelihood_ratio directly on the loglikelihoods of individual distribution instances.

# Using the Fit object with explicit nested flag
fit.distribution_compare('power_law', 'truncated_power_law', nested=True)

# Manual comparison without a Fit object
data = powerlaw.load_test_dataset('fires')

loglikelihoods_pl = powerlaw.Power_Law(data=data).loglikelihoods()
loglikelihoods_exp = powerlaw.Exponential(data=data).loglikelihoods()

# Use nested=False if you know they are not nested
ratio = powerlaw.loglikelihood_ratio(loglikelihoods_pl, loglikelihoods_exp, nested=False)

The powerlaw package uses specific mathematical derivations for the exponential distribution to perform analytics and generate random samples. The distribution is defined by the parameter $\lambda$ over a domain $[x_{min}, x_{max}]$.

Probability Distribution Function (PDF)

The PDF follows the form: $p(x) = C e^{-\lambda x}$

Where the normalization constant $C$ is: $C = - \lambda [ e^{-\lambda x_{max}} - e^{-\lambda x_{min}} ]^{-1}$

For a semi-infinite domain ($x_{max} \rightarrow \infty$), it simplifies to: $C = \lambda e^{\lambda x_{min}}$

Note: Storing the normalization constant $C$ separately from unnormalized PDF values can lead to overflow or underflow issues.

Cumulative Distribution Function (CDF)

The CDF $c(x)$ is defined as: $c(x) = \frac{ e^{-\lambda x} - e^{-\lambda x_{min}} }{ e^{-\lambda x_{max}} - e^{-\lambda x_{min}} }$

For a semi-infinite domain ($x_{max} \rightarrow \infty$), it simplifies to: $c(x) = 1 - e^{-\lambda (x - x_{min})}$

To impose relational rules between parameters (e.g., Lambda < alpha), use the parameter_constraints keyword argument in a Distribution object or the Fit class.

Constraints are defined using a dictionary with the following keys:

• type: Either 'eq' (equality) or 'ineq' (inequality).
  • 'eq': The constraint function must evaluate to zero.
  • 'ineq': The constraint function must evaluate to a non-negative (positive or zero) value.
• fun: A function that implements the constraint. Crucially, this function must take a single argument: the distribution object itself.
• dists (optional): A list of distribution names (e.g., ['power_law']) to which the constraint applies. If omitted, the constraint applies to all distributions.

Important Implementation Detail: Avoid using boolean expressions (like return x1 + x2 == 5) inside your constraint function. Instead, return the mathematical difference (like return 5 - x2 - x1) so the optimizer can determine how close it is to satisfying the constraint.

# Example: Inequality constraint requiring at least N points in the data
def constraint(dist):
    N = 100
    # dist.data is already cropped from xmin to xmax
    return len(dist.data) - N 

constraint_dict = {
    "type": 'ineq',
    "fun": constraint,
    "dists": ['power_law']
}

fit = powerlaw.Fit(data, parameter_constraints=constraint_dict)

Unlike general-purpose fitting libraries (e.g., scipy, lmfit) that minimize a cost function (like Mean Squared Error) to fit a mathematical function $f(x, \vec p)$ to specific data points, powerlaw is specialized for statistical distributions.

Key differences include:

1. PDF-based Fitting: Instead of fitting a function to raw data values, powerlaw fits the Probability Density Function (PDF) to the probabilities of finding those samples. The cost function compares the PDF of the original data to the PDF of the candidate distribution.
2. Parameter Constraints: Because powerlaw assumes the model is a proper probability distribution, it automatically handles normalization. For example, in a power law, you only need to fit the exponent ($\alpha$), as the normalization constant is mathematically constrained and does not need to be an arbitrary parameter.
3. Analytical vs. Numerical: Many distributions in powerlaw use analytical expressions for Maximum Likelihood Estimation (MLE). This allows for $\mathcal{O}(1)$ parameter calculation without the need for expensive numerical optimization (like scipy.optimize.minimize).
4. Statistical Comparison: Because the library understands the underlying distributions, it enables comparisons using likelihood measures and distribution distances (e.g., Kolmogorov-Smirnov distance) rather than just simple residual/error analysis.

Fitting a discrete power law requires numerical optimization, which can be slow when iterating over different xmin values. To speed this up, you can use the estimate_discrete option:

• Automatic: When discrete=True, estimate_discrete is True by default if xmin >= 6 and no xmax is provided. This uses an analytic estimate of $\alpha$ that is accurate to within ~1%.
• Explicitly True: Forces the use of the fast analytic estimation technique. Use this with caution, as it may produce incorrect fits if your data violates the underlying assumptions.
• Explicitly False: Forces the use of the slower, exact numerical optimization method.

Note: This option can be set during powerlaw.Fit initialization or when calling Distribution.generate_random to control the speed/accuracy tradeoff of simulated data generation.

# Uses the default discrete estimation technique (if xmin >= 6)
fit = powerlaw.Fit(data, discrete=True)

# Forces the fast estimation technique (use with caution!)
fit = powerlaw.Fit(data, discrete=True, estimate_discrete=True)

If you do not specify an xmin when creating a powerlaw.Fit object, the library automatically executes an optimization process to find the best xmin value. This process involves fitting the data starting from each unique value in the dataset and selecting the one that results in the minimal distance between the data and the fit.

Distance Metrics

You can choose the evaluation metric using the xmin_distance keyword argument. Supported schemes include:

• Kolmogorov-Smirnov distance (D)
• Kuiper distance (V)
• Anderson-Darling distance (A^2)

Performance and Progress

Finding xmin is computationally expensive. To monitor progress, set verbose=2 in the Fit constructor to display a progress bar.

Restricting the Search Range

You can restrict the range of values the algorithm searches through for the optimal xmin by passing a tuple (min_bound, max_bound) to the xmin argument.

# Use find_xmin() to get the best value for xmin automatically
fit = powerlaw.Fit(data)

# Use find_xmin() with a progress bar
fit = powerlaw.Fit(data, verbose=2)

# Restrict the search for xmin to the range between 50 and 100
fit = powerlaw.Fit(data, xmin=(50, 100))

The package generates random numbers following a power law distribution using inverse transform sampling. This involves passing uniform random numbers $r \in [0, 1)$ through the inverse CDF $c^{\dagger}(r)$.

Inverse CDF Formulas

Depending on the parameters, the package uses different expressions for $c^{\dagger}(r)$:

1. For $\alpha < 1$ (requires finite $x_{max}$): $c^{\dagger}(r) = [ r x_{max}^{1 - \alpha} + (1 - r) x_{min}^{1-\alpha} ]^{\frac{1}{1 - \alpha}}$

2. For $\alpha > 1$ and semi-infinite domain ($x_{max} \rightarrow \infty$): $c^{\dagger}(r) = x_{min} (1 - r)^{\frac{1}{1 - \alpha}}$

Sampling Strategies

• Unbounded/Semi-infinite sampling: Generate $r \in [0, 1)$ and use the simplified expression for $\alpha > 1$ to get samples in $[x_{min}, \infty)$.
• Bounded sampling: To get samples in $[x_{min}, a)$ where $a < x_{max}$, you can restrict the uniform random numbers to the range $[0, c(a))$. This allows using the same simplified inverse CDF expression for both bounded and unbounded samples when $\alpha > 1$.
• Mandatory finite $x_{max}$: If $\alpha < 1$, you must use the full expression with a finite $x_{max}$ and can use uniform random numbers in the range $[0, 1)$.