demo/moby.R

In poweRlaw: Analysis of Heavy Tailed Distributions

library("poweRlaw")
## Load in the full Moby dick data set
data(moby)

##Create a distribution object
m = displ$new(moby)

##Set parameters explicitly
m$setXmin(1);m$setPars(2.2)

##Plot the data with fitted lines
plot(m)
lines(m, col=2)

##Now calculate the mle estimate and add to the plot
mle_pars = estimate_pars(m)
est = estimate_pars(m, pars=seq(1.5, 3.5, 0.001))
m$setPars(est)

lines(m, col=3)

##Estimate xmin
(est = estimate_xmin(m))
m$setXmin(est)
lines(m, col=4)

##Uncertainity in xmin
bs1 = bootstrap(m, no_of_sims=100, threads=2)
bs2 = bootstrap_p(m, no_of_sims=100, threads=2)

##This can take a while
bs1 = bootstrap(m, no_of_sims=5000, threads=4, pars = seq(1.6, 2.5, 0.01))
bs2 = bootstrap(m, no_of_sims=5000, threads=4)

##Plot uncertainity in xmin and alpha
par(mfrow=c(1, 2))
hist(bs1$bootstraps[,2], breaks="fd")
hist(bs1$bootstraps[,3], breaks="fd")

hist(bs2$bootstraps[,2], breaks="fd")
hist(bs2$bootstraps[,3], breaks="fd")

apply(bs1[[2]], 2, sd)
apply(bs2[[3]], 2, sd)

bootstrap(m, no_of_sims = 4, threads=4)