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

Question

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?