csgillespie / poweRlaw

This package implements both the discrete and continuous maximum likelihood estimators for fitting the power-law distribution to data. Additionally, a goodness-of-fit based approach is used to estimate the lower cutoff for the scaling region.
109 stars 24 forks source link

Output of 'lines' incorrect? #102

Open Brennen-Fagan opened 2 months ago

Brennen-Fagan commented 2 months ago

Hi Colin,

Some collaborators and I noticed something odd with the plots of the fitted discrete power-law distribution. I've also attached the results of using Aaron Clauset's plplot for contrast, which seems to replicate what I've done in blue instead of what lines does in red. Any idea why there's such a significant difference?

Data from Correlates of War.

library(poweRlaw)
dat <- read.csv(text = paste0(
  "d\n",
  "1000\n130000\n19283\n7527\n6000\n2600\n1300\n264200\n2000\n22500\n10000\n",
  "1000\n1000\n20000\n1000\n4481\n310000\n1000\n44100\n204313\n4000\n285000\n",
  "13868\n10079\n12100\n1000\n15000\n2000\n3685\n3003\n4000\n151831\n1000\n",
  "1000\n10000\n20000\n82000\n60500\n8578031\n11750\n13246\n100000\n11000\n",
  "50000\n40000\n1000\n3200\n60000\n92661\n2100\n20000\n1000000\n1726\n",
  "28000\n16634907\n151798\n1400\n3500\n8000\n910084\n2370\n3221\n2426\n",
  "1122\n1800\n1853\n1021442\n7061\n19600\n13866\n5368\n1900\n6525\n11223\n",
  "14439\n1500\n2700\n10500\n8000\n3000\n21000\n1250000\n1001\n1655\n8000\n",
  "4000\n41466\n5240\n14000\n1500\n120000\n5002\n1172\n4002\n7173"
))

# Fit a discrete power law.
dat_pl <- displ$new(dat$d)
dat_pl$setXmin(estimate_xmin(dat_pl, xmax = Inf))

# Plot.
plot(dat_pl) # Plot Data.
lines(dat_pl, col = "red") # Plot power law fit.
lines( # Plot adjusted complementary cumulative density:
  dat_pl$dat,
  ( # Complementary Density.
    1 - ppldis(dat_pl$dat,
      xmin = dat_pl$getXmin(),
      alpha = dat_pl$getPars()
    )
  ) * ( # Scale by data in tail only.
    1 - sum(dat_pl$dat < dat_pl$getXmin()) / length(dat_pl$dat)
  ),
  col = "blue"
)

Created on 2024-09-16 with reprex v2.1.1

warsizes_cow_MATLAB