rlog | R Documentation |
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).
rlog(n, p, Ip = 1 - p)
n |
sample size, that is, length of the resulting vector of random variates. |
p |
parameter in |
Ip |
|
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 integer
s 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(-\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.
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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.