VGAM (version 0.9-1)

Polono: The Poisson Lognormal Distribution

Description

Density, and random generation for the Poisson lognormal distribution.

Usage

dpolono(x, meanlog = 0, sdlog = 1, bigx = 170, ...)
ppolono(q, meanlog = 0, sdlog = 1,
        isOne = 1 - sqrt( .Machine$double.eps ), ...)
rpolono(n, meanlog = 0, sdlog = 1)

Arguments

x, q

vector of quantiles.

n

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

meanlog, sdlog

the mean and standard deviation of the normal distribution (on the log scale). They match the arguments in Lognormal.

bigx

Numeric. This argument is for handling large values of x and/or when integrate fails. A first order Taylor series approximation [Equation (7) of Bulmer (1974)] is used at value

isOne

Used to test whether the cumulative probabilities have effectively reached unity.

...

Arguments passed into integrate.

Value

• dpolono gives the density, ppolono gives the distribution function, and rpolono generates random deviates.

Details

The Poisson lognormal distribution is similar to the negative binomial in that it can be motivated by a Poisson distribution whose mean parameter comes from a right skewed distribution (gamma for the negative binomial and lognormal for the Poisson lognormal distribution).

References

Bulmer, M. G. (1974) On fitting the Poisson lognormal distribution to species-abundance data. Biometrics, 30, 101--110.

See Also

lognormal, poissonff, negbinomial.

Examples

meanlog <- 0.5; sdlog <- 0.5; yy <- 0:19
sum(proby <- dpolono(yy, m = meanlog, sd = sdlog)) # Should be 1
max(abs(cumsum(proby) - ppolono(yy, m = meanlog, sd = sdlog))) # Should be 0

opar = par(no.readonly = TRUE)
par(mfrow = c(2, 2))
plot(yy, proby, type = "h", col = "blue", ylab = "P[Y=y]", log = "",
     main = paste("Poisson lognormal(m = ", meanlog,
                  ", sdl = ", sdlog, ")", sep = ""))

y <- 0:190 # More extreme values; use the approximation and plot on a log scale
(sum(proby <- dpolono(y, m = meanlog, sd = sdlog, bigx = 100))) # Should be 1
plot(y, proby, type = "h", col = "blue", ylab = "P[Y=y] (log)", log = "y",
     main = paste("Poisson lognormal(m = ", meanlog,
                  ", sdl = ", sdlog, ")", sep = "")) # Note the kink at bigx

# Random number generation
table(y <- rpolono(n = 1000, m = meanlog, sd = sdlog))
hist(y, breaks = ((-1):max(y))+0.5, prob = TRUE, border = "blue")
par(opar)