I'm exploring power law fitting and would like to get to a point of making a power law package like what was discussed in this post, where OP references the R poweRlaw package. As per the comments/answer in the post, MMa can find a distribution (e.g. ZipfDistribution) however it's not clear to me how to optimize the fit. I mean specifically the package's utility in choosing a starting value for x, i.e. $x_{min}$.

Here is a little background to show where I am at, and so others trying to do the same can find utility from this post:

The probability distribution of a power law follows:

$$p(x) \propto x^{-\alpha}$$

For discrete data, e.g. in the case of graph vertex degrees, the Maximum Likelihood Estimator (MLE) can be used to approximate, $\hat\alpha$ with a standard error, $\sigma$:

$$\hat\alpha \space \backsimeq 1+N[\sum_{i=1}^{N}\ln\space \frac{k_{i}}{k_{min}-0.5}]^{-1}$$ $$\sigma=\frac{\hat{\alpha}-1}{N^{0.5}}$$

There's a great 2013 MMa journal article written by Todd Silvestri that shows how you can do the power law fitting exercise based off of a random walk of the internet (i.e. emulating Barabasi's world wide web work). I am working off of Todd's code:

PLExponentEstimated[g_?DirectedGraphQ, kmin_Integer /; kmin >= 6,distType_String /; 
StringMatchQ[distType, {"InDegree","OutDegree"}]]:=
{$FunctionName = "PLExponentEstimated", vertexDegrees, ki, M, exponentMLE, expStandardError},

Switch[distType,(*Picks whether to choose In-Degree)
  "InDegree", vertexDegrees = VertexInDegree[g]
  "OutDegree", vertexDegrees = VertexOutDegree[g];
ki = Select[vertexDegrees, # >= kmin &];(*defined the values for ki*)
M = Length[ki];(*values in the graph, N*)

exponentMLE = 1 + M*Total[Map[N[Log[#/(kmin - 1/2)]] &, ki]]^-1;(*MLE estimator*)

expStandardError = (exponentMLE - 1)/Sqrt[M];(*standard error*)
{exponentMLE, expStandardError} (*output*)

The inputs to the function are a graph and a starting value ($x_{min}$), and the MLE for $\alpha$ & $\sigma$ are the outputs.

This is what I want to accomplish:

  1. In the article Todd chooses $x_{min} \geq$ 6... however like the poweRlaw package, I want to instead estimate $x_{min}$ by comparing the distance between the data and the fitted model CDF (Clauset et al 09): $$D= max_{x\geq x_{min}}|S(x)-P(x)|$$ Where S(x) and P(x) are the CDFs of the data and the model, resp (for $x\geq x_{min}$).
  2. Then I want to plot the fitted power law distribution over the complementary cumulative distribution function (cCDF). Two fits on the same plot: one for all values $x_{0},x_{max}$, and the other for the optimized $x_{min}$ to $x_{max}$ (this isn't too challenging, but I'm adding it in the spirit of creating a complete analysis tool).

My general approach:

I'm not sure how to do this.. but so far (1) I can estimate an $\alpha$ for a graph using Todd's code, and (2) Create a complementary CDF also using the code with some minor changes of mine:

g = RandomGraph[BarabasiAlbertGraphDistribution[100, 2]] //DirectedGraph
g[g_, k_] := Sum[Count[VertexInDegree[g], j]/VertexCount[g], {j, k,Max[VertexInDegree[g]]}]
DeleteCases[#, {x_, y_} /; x y == 0] &@Table[{k, g[sp, k]}, {k, 0, Max[VertexInDegree[sp]]}]

I don't know how to write the iterative loop that estimates the optimal $x_{min}$. However, I do know the R package uses an equivalent to the KolmogorovSmirnovTest with the built in MonteCarlo method (Sjoerd explains this well here) to find the optimum $x_{min}$.

My Question:

  1. The R poweRlaw package uses a for loop to iterate different values of $\alpha$ and $x_{min}$ and then outputs the respective values when minimum $D$ is found. Is there another, perhaps most efficient way of doing this in MMa? Or should it be a loop, and how to do it?
  2. Will MMa's ZipfDistribution work for $P(x)$ or should I define the function I wish to work with?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.