Search code examples
rplotpower-law

Multiple power law plots with package poweRlaw


I tried using the package poweRlaw to plot a few powerlaw fits. It seems working for a single curve. But I could not plot multiple plots on the same graph. Ref: Is there a way within this package? [P.S. I am a novice]

set.seed(1)
x1 <- ceiling(rlnorm(1000, 4))

x2 <- ceiling(rlnorm(1000, 2))

library(poweRlaw)

pl_d = pl_data$new(x1)
plot(pl_d)

#Now fit the powerlaw
m = displ$new(pl_d)
#Estimate the cut-off
#estimate_xmin(m)
aaa <- estimate_xmin(m)
aaa <- as.data.frame(aaa)
aaa <- aaa[2,1]
x_min <- min(table(x))
m$setXmin(aaa); m$mle()
#Plot the data and the PL line
#plot(m)
lines(m, col=2)

# next POWER LAW graph

#Plot the data
pl_d = pl_data$new(x2)
points(pl_d)

#Now fit the powerlaw
m = displ$new(pl_d)
#Estimate the cut-off
#estimate_xmin(m)
aaa <- estimate_xmin(m)
aaa <- as.data.frame(aaa)
aaa <- aaa[2,1]
x_min <- min(table(x))
m$setXmin(aaa); m$mle()
#Plot the data and the PL line
#points(m)
lines(m, col=3)

Solution

  • You can use the lines commands to add multiple best fitting power-laws. So, using your data:

    set.seed(1)
    x1 = ceiling(rlnorm(1000, 4))
    

    we load the package and create a discrete power-law object:

    m = displ$new(x1)
    

    Then plot the data

    plot(m)
    

    Next, we set xmin and alpha for each power-law and add it to the graph using lines:

    m$setXmin(100)
    p_100 = estimate_pars(m)
    m$setPars(p_100)
    lines(m, col=2, lwd=2)
    
    ##Line 2    
    m$setXmin(202)
    p_200 = estimate_pars(m)
    m$setPars(p_200)
    lines(m, col=3, lwd=2)
    

    enter image description here


    If you want to have multiple data sets on the same plot, the easiest way is to plot one data set, but save the data as a data frame. For example, generate some data sets:

    set.seed(1)
    x1 = ceiling(rlnorm(1000, 4))
    x2 = ceiling(rlnorm(1000, 5))
    

    Fit the power-laws

    m1 = displ$new(x1)
    m1$setXmin(estimate_xmin(m1))
    
    m2 = displ$new(x2)
    m2$setXmin(estimate_xmin(m2))
    

    Plot data set 2, but save the data as pts2

    pts2 = plot(m2)
    

    Plot and add lines as normal:

    plot(m1)
    points(pts2$x, pts2$y, col=3)
    lines(m1, col=2)
    lines(m2, col=4)
    

    to get:

    enter image description here