# Power law fitting package that optimizes fitting parameters based off of the KolmogorovSmirnov MC method

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"}]]:=
Module[
{$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?