dpowbend {sads}

Puyeo's Power-bend discrete distribution

Description

Density, distribution function, quantile function and random generation for discrete version of Pueyo's power-bend distribution with parameters s and omega.

Usage

dpowbend( x, s, omega = 0.01, oM = -log(omega), log=FALSE)
ppowbend ( q, s, omega = 0.01, oM = -log(omega), lower.tail=TRUE, log.p=FALSE)
qpowbend( p, s, omega = 0.01, oM = -log(omega), lower.tail= TRUE, log.p=FALSE)
rpowbend( n, s, omega)

Arguments

x

vector of (integer x>0) quantiles. In the context of species abundance distributions, this is a vector of abundances of species in a sample.

q

vector of (integer x>0) quantiles. In the context of species abundance distributions, a vector of abundances of species in a sample.

p

vector of probabilities.

n

number of random values to return.

s

positive real s > 1; exponent of the power-bend distribution.

omega

positive real greater than 0; bending parameter of the distribution. The current implementation only accepts omega smaller than 5.

oM

real number; alternative parametrization which is used for faster fitting on the fitpowbend function. The current implementation only accepts oM greater than approximately -1.5. You can specify either 'omega' or 'oM', but not both.

log, log.p

logical; if TRUE, probabilities p are given as log(p).

lower.tail

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

Details

The power-bend density is a discrete probability distribution based on the power probability density with an added term for “bending”, defined for integer x > 0:

p(x) = \frac{x^{-s} e^{- \omega x}}{Li_s (e^\omega)}

The bending term can be seen as a finite-size correction of the power law (Pueyo 2006). Therefore, the power-bend can describe better samples taken from power-law distributions.

The function Li_s(e^\omega) is the Polylogarithm of the exponential of omega on base s, and represents the integration constant. A naive implementation of the Polylogarithm function is included in the package, and accepts values for non-integer s and omega smaller than 5, which cover the cases of biological relevance.

This distribution was proposed by S. Pueyo to describe species abundance distributions (sads) as a generalization of the logseries distribution. Fisher's logseries correspond to the power-bend with parameters s=1 and \omega = - \log(\frac{N}{N + \alpha}) where N is sample size or total number of individuals and alpha is the Fisher's alpha parameter of the logseries.

When fitted to sads, power-law distributions usually overestimates the abundance of common species, and in practice power-bend corrects this problem and usually provides a better fit to abundance distributions.

Value

dpowbend gives the (log) density of the function, ppowbend gives the (log) distribution function, qpowbend gives the quantile function.

Invalid values for parameter s or omega will result in return values NaN, with a warning. Note that integer values of s and omega values larger than 5 are currently not supported, and will also return NaN.

Author(s)

Paulo I Prado [e-mail] and Andre Chalom.

References

Pueyo, S. (2006) Diversity: Between neutrality and structure, Oikos 112: 392-405.

See Also

dpower for the power distribution; dls for the logseries distribution, which is a particular case of the power-bend, fitpowbend for maximum likelihood estimation in the context of species abundance distributions.

Examples

x <- 1:20
PDF <- dpowbend(x=x, s=2.1, omega=0.5)
CDF <- ppowbend(q=x, s=2.1, omega=0.5)
par(mfrow=c(1,2))
plot(x,CDF, ylab="Cumulative Probability", type="b",
     main="Powerbend distribution, CDF")
plot(x,PDF, ylab="Probability", type="h",
     main="Powerbend distribution, PDF")
par(mfrow=c(1,1))

## The powbend distribution is a discrete PDF, hence:
all.equal( ppowbend(10, s=2.1, omega=0.5), sum(dpowbend(1:10, s=2.1, omega=0.5)) ) # should be TRUE

## quantile is the inverse of CDF
all.equal(qpowbend(CDF, s=2.1, omega=0.5), x)

## Equivalence between power-bend and logseries
x <- 1:100
N <- 1000
alpha <- 5
X <- N/(N+alpha)
omega <- -log(X)
PDF1 <- dls(x, N, alpha)
PDF2 <- dpowbend(x, s=1, omega)
plot(PDF1,PDF2, log="xy")
abline(0,1, col="blue")