rlog: Sampling Logarithmic Distributions

Description

Generating random variates from a Log(p) distribution with probability mass function $$p_k=\frac{p^k}{-\log(1-p)k},\ k\in\mathbf{N}, $$ where $p\in(0,1)$. The implemented algorithm is the one named “LK” in Kemp (1981).

Usage

rlog(n, p, Ip = 1 - p)

Arguments

sample size, that is, length of the resulting vector of random variates.

parameter in $(0,1)$.

$= 1 - p$, possibly more accurate, e.g, when $p\approx 1$.

Value

A vector of positive integers of length n containing the generated random variates.

Details

For documentation and didactical purposes, rlogR is a pure-R implementation of rlog. However, rlogR is not as fast as rlog (the latter being implemented in C).

References

Kemp, A. W. (1981), Efficient Generation of Logarithmically Distributed Pseudo-Random Variables, Journal of the Royal Statistical Society: Series C (Applied Statistics) 30, 3, 249--253.

Examples

Run this code

# NOT RUN {
## Sample n random variates from a Log(p) distribution and plot a
## "histogram"
n <- 1000
p <- .5
X <- rlog(n, p)
table(X) ## distribution on the integers {1, 2, ..}
## ==> The following plot is more reasonable than a  hist(X, prob=TRUE) :
plot(table(X)/n, type="h", lwd=10, col="gray70")

## case closer to numerical boundary:
lV <- log10(V <- rlog(10000, Ip = 2e-9))# Ip = exp(-theta) <==> theta ~= 20
hV <- hist(lV, plot=FALSE)
dV <- density(lV)
## Plot density and histogram on log scale with nice axis labeling & ticks:
plot(dV, xaxt="n", ylim = c(0, max(hV$density, dV$y)),
     main = "Density of [log-transformed] Log(p),  p=0.999999..")
abline(h=0, lty=3); rug(lV); lines(hV, freq=FALSE, col = "light blue"); lines(dV)
rx <- range(pretty(par("usr")[1:2]))
sx <- outer(1:9, 10^(max(0,rx[1]):rx[2]))
axis(1, at=log10(sx),     labels= FALSE, tcl = -0.3)
axis(1, at=log10(sx[1,]), labels= formatC(sx[1,]), tcl = -0.75)
# }

Run the code above in your browser using DataCamp Workspace