rlog {copula}

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

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\approx 1.

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).

Value

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

Note

In the copula package, the Log(p) distribution is needed only for generating Frank copula observations, namely in copFrank@V0(), where p = 1 - exp(-\theta), i.e., p = -expm1(-theta) and Ip = exp(-theta).

For large \theta it would be desirable to pass -theta to rlog() instead of p. This has not yet been implemented.

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

## 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)