poweRlaw (version 0.80.0)

dpldis: Discrete power-law distribution

Description

Density, distribution function and random number generation for the discrete power law distribution with parameters xmin and alpha.

Usage

dpldis(x, xmin, alpha, log = FALSE)
ppldis(q, xmin, alpha, lower.tail = TRUE)
rpldis(n, xmin, alpha, discrete_max = 10000)

Value

dpldis returns the density, ppldis returns the distribution function and rpldis return random numbers.

Arguments

x, q

vector of quantiles. The discrete power-law distribution is defined for x > xmin.

xmin

The lower bound of the power-law distribution. For the continuous power-law, xmin >= 0. for the discrete distribution, xmin > 0.

alpha

The scaling parameter: alpha > 1.

log

logical (default FALSE) if TRUE, log values are returned.

lower.tail

logical; if TRUE (default), probabilities are \(P[X \le x]\), otherwise, \(P[X > x]\).

n

Number of observations. If length(n) > 1, the length is taken to be the number required.

discrete_max

The value when we switch from the discrete random numbers to a CTN approximation.

Details

The Clauset, 2009 paper provides an algorithm for generating discrete random numbers. However, if this algorithm is implemented in R, it gives terrible performance. This is because the algorithm involves "growing vectors". Another problem is when alpha is close to 1, this can result in very large random number being generated (which means we need to calculate the discrete CDF for very large values).

The algorithm provided in this package generates true discrete random numbers up to 10,000 then switches to using continuous random numbers. This switching point can altered by changing the discrete_max argument.

In order to get a efficient power-law discrete random number generator, the algorithm needs to be implemented in C.

References

Clauset, Aaron, Cosma Rohilla Shalizi, and Mark EJ Newman. "Power-law distributions in empirical data." SIAM review 51.4 (2009): 661-703.

Examples

xmin = 1; alpha = 2
x = xmin:100

plot(x, dpldis(x, xmin, alpha), type="l")
plot(x, ppldis(x, xmin, alpha), type="l", main="Distribution function")
dpldis(1, xmin, alpha)

###############################################
## Random number generation                   #
###############################################
n = 1e5
x1 = rpldis(n, xmin, alpha)
## Compare with exact (dpldis(1, xmin, alpha))
sum(x1==1)/n
## Using only the approximation
x2 = rpldis(n, xmin, alpha, 0)
sum(x2==1)/n