-
Notifications
You must be signed in to change notification settings - Fork 0
/
power_law.R
49 lines (45 loc) · 1.42 KB
/
power_law.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
# Check direct variant to solve xmin problem
FitPowerLaw <- function(x) {
suppressMessages(require(VGAM))
x <- as.integer(x)
# Range of scaling parameters
vec <- seq(1.5, 3.5, 0.01)
zvec <- zeta(vec)
xmins <- sort(unique(x))
# limit <- c()
xmins <- xmins[-length(xmins)]
xmax <- max(x)
dat <- matrix(0, nrow=length(xmins), ncol=2)
z <- x
for (i in 1:length(xmins)) {
xmin <- xmins[i]
z <- z[z >= xmin]
n <- length(z)
# Use maximization of likelihood function to estimate alpha
if (xmin == 1) {
zdiff <- rep(1, length(vec))
} else {
zdiff <- apply(X=rep(t(1:(xmin-1)), length(vec))^-
t(kronecker(t(array(1, xmin - 1)), vec)), MARGIN=2, FUN=sum)
}
L <- -vec * sum(log(z)) - n * log(zvec - zdiff)
I <- which.max(L)
# Compute KS statistic
fit <- cumsum((((xmin:xmax)^-vec[I])) /
(zvec[I] - sum((1:(xmin-1))^-vec[I])))
# Dirty, dirty, ...
cdi <- cumsum(hist(z, c(min(z) - 1, (xmin + 0.5):xmax, max(z) + 1),
plot=FALSE)$counts / n)
dat[i, ] <- c(max(abs(fit - cdi)), vec[I])
}
D <- min(dat[, 1])
I <- which.min(dat[, 1])
xmin <- xmins[I]
n <- sum(x >= xmin)
alpha <- dat[I, 2]
# Correction for finite sample size
alpha <- alpha * (n - 1) / n + 1 / n
pval <- 1 - .C("pkolmogorov2x", p = as.double(D), as.integer(n), PACKAGE = "stats")$p
rval <- list(statistic=D, p.value=pval, xmin=xmin, n=n, alpha=alpha)
return(rval)
}