Sampling transformation - rexp vs rweibull

Question

I am working with different sampling functions, and I am wondering why these two formulations do not give the same result

n=2
set.seed(1)
rweibull(n,shape = 1,scale = 1)
# [1] 1.3261078 0.9885284
set.seed(1)
rexp(n,rate = 1)
# [1] 0.7551818 1.1816428

when this is equivalent:

x <- c(0, rlnorm(50))
all.equal(dweibull(x, shape = 1), dexp(x))

Is it a problem of reverse transformation sampling ?

If yes, why ?

Thanks,

Answer 1

First of all, d* denotes the density function, which is deterministic and that's why you can achieve the same results in the following code

x <- c(0, rlnorm(50))
all.equal(dweibull(x, shape = 1), dexp(x)

---

However, r* gives random samples of a given distribution, which depends on how the randomness is triggered.

• The source code of rweibull

#include "nmath.h"

double rweibull(double shape, double scale)
{
    if (!R_FINITE(shape) || !R_FINITE(scale) || shape <= 0. || scale <= 0.) {
    if(scale == 0.) return 0.;
    /* else */
    ML_ERR_return_NAN;
    }

    return scale * pow(-log(unif_rand()), 1.0 / shape);
}

• The source code of rexp

#include "nmath.h"

double rexp(double scale)
{
    if (!R_FINITE(scale) || scale <= 0.0) {
    if(scale == 0.) return 0.;
    /* else */
    ML_ERR_return_NAN;
    }
    return scale * exp_rand(); // --> in ./sexp.c
}

where exp_rand() is used to generate exponential random variable, which is NOT simply the same as pow(-log(unif_rand()), 1.0 / shape) shown in rweibull.c, but more complicated like below

#include "nmath.h"

double exp_rand(void)
{
    /* q[k-1] = sum(log(2)^k / k!)  k=1,..,n, */
    /* The highest n (here 16) is determined by q[n-1] = 1.0 */
    /* within standard precision */
    const static double q[] =
    {
    0.6931471805599453,
    0.9333736875190459,
    0.9888777961838675,
    0.9984959252914960,
    0.9998292811061389,
    0.9999833164100727,
    0.9999985691438767,
    0.9999998906925558,
    0.9999999924734159,
    0.9999999995283275,
    0.9999999999728814,
    0.9999999999985598,
    0.9999999999999289,
    0.9999999999999968,
    0.9999999999999999,
    1.0000000000000000
    };

    double a = 0.;
    double u = unif_rand();    /* precaution if u = 0 is ever returned */
    while(u <= 0. || u >= 1.) u = unif_rand();
    for (;;) {
    u += u;
    if (u > 1.)
        break;
    a += q[0];
    }
    u -= 1.;

    if (u <= q[0])
    return a + u;

    int i = 0;
    double ustar = unif_rand(), umin = ustar;
    do {
    ustar = unif_rand();
    if (umin > ustar)
        umin = ustar;
    i++;
    } while (u > q[i]);
    return a + umin * q[0];
}