# Generating random samples from a power law and testing them with R igraph

I haven't looked into power laws before, so the question may very well have an obvious, embarrassing answer. In some ways it dovetails with this post.

The problem may be in the method, or in the interpretation of the power.law.fit() function in the igraph r package.

I used the following method from Wolfram MathWorld and also reflected here:

$$f(y)= [(x_1^{n+1}-x_0^{n+1})\,y+x_0^{n+1}]^{1/(n+1)}$$

where $x_0$ and $x_1$ come about from using inverse transform sampling to generate a power-law sample from a uniform distribution $f_X(x)=Cx^n, x\in [x_o,x_1].$ These values arise during the calculation of $C$ to makes sure that the pdf integrates to 1.

I'm plotting the results because the result does suggest a power law distribution on visual inspection:

x1 = 5           # Maximum value
x0 = 0.1         # Min value can't be zero; otherwise X^0^(neg) is 1/0.
alpha = -2.5     # It has to be negative.
y = runif(1e5)   # Number of samples
x = ((x1^(alpha+1) - x0^(alpha+1))*y + x0^(alpha+1))^(1/(alpha+1))
hist(x, prob = T, breaks=40, ylim=c(0,10), xlim=c(0,1.2), border=F,
col="yellowgreen", main="Power law density")
lines(density(x), col="chocolate", lwd=1)
h = hist(x, prob=T, breaks=40, plot=F)
plot(h$count, log="xy", type='l', lwd=1, lend=2, xlab="", ylab="", main="Density in logarithmic scale")  Here is the summary of the data: > summary(x) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.1000 0.1208 0.1584 0.2590 0.2511 4.9388  The problem is in interpreting the results of applying power.law.fit() to the generated data in x. Aside from the fact that each time I run this function on x it takes from 5 to 10 minutes to return results, these return the minimum value,$0.1,$and the alpha value,$-2.5\$ without a glitch, yet they seem to indicate that the vector does not come from a power law distribution.

I have run the test for set.seed(0) through set.seed(4) with the following results:

1. set.seed(0); logLik 106703.4; KS.stat 0.005310278; KS.p 0.007109891;
2. set.seed(1); logLik 106737.5; KS.stat 0.005025666; KS.p 0.01280314;
3. set.seed(2); logLik 105858.9; KS.stat 0.005552832; KS.p 0.004218761;
4. set.seed(3); logLik 107102.7; KS.stat 0.005953606; KS.p 0.001674007;
5. set.seed(4); logLik 107004.3; KS.stat 0.00607452;  KS.p 0.001247162.


I presume that the key to interpret these results is according to here:

logLik Numeric scalar, the log-likelihood of the fitted parameters.

KS.stat Numeric scalar, the test statistic of a Kolmogorov-Smirnov test that compares the fitted distribution with the input vector. Smaller scores denote better fit.

KS.p Numeric scalar, the p-value of the Kolmogorov-Smirnov test. Small p-values (less than 0.05) indicate that the test rejected the hypothesis that the original data could have been drawn from the fitted power-law distribution.

THE QUESTION: Is there something incorrect, not power-law distributed, in the generated values, or am I misusing, or misinterpreting the testing with igraph?

• With that sample size is it reasonable to expect it to fit any distribution you choose? Oct 29, 2017 at 15:18
• @mdewey shouldn't this be the more reason to get a high p value, suggesting that the sample could indeed come from a power law distribution? Oct 29, 2017 at 15:35
• I would have thought that a power law distribution (or at least a Pareto distribution) would not have a maximum value parameter. I might generate data with something like xmin=0.1; alpha=2.5; y=runif(1e5); z=xmin/(1-y)^(1/alpha) Oct 29, 2017 at 18:18
• @Henry That may very well be the problem. If you have time, take a look at this Wolfram MathWorld post. Oct 29, 2017 at 18:21