Sampling Logarithmic Distributions
Generating random variates from a Log(p) distribution with probability mass function
p_k = p^k/(-log(1-p)k), k in IN,
where p in (0,1). The implemented algorithm is the one named “LK” in Kemp (1981).
rlog(n, p, Ip = 1 - p)
n |
sample size, that is, length of the resulting vector of random variates. |
p |
parameter in (0,1). |
Ip |
= 1 - p, possibly more accurate, e.g, when p ~= 1. |
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).
A vector of positive integers of length n containing the
generated random variates.
In the copula package, the Log(p) distribution is needed only
for generating Frank copula observations, namely in
copFrank@V0(), where p = 1 - exp(-θ), i.e.,
p = -expm1(-theta) and Ip = exp(-theta).
For large θ it would be desirable to pass -theta
to rlog() instead of p. This has not yet been implemented.
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.
## 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)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.