How can you draw a random sample from an RV with density $ \exp(x - \exp(x))$?

Question

How can you draw a random sample from a random variable whose density is given by $\exp(x - \exp(x))$?

I am trying to work through the blog at https://staffblogs.le.ac.uk/bayeswithstata/2015/03/27/poisson-regression-with-two-random-effects-mcmc-by-data-augmentation/

The section Simulating the unobserved data describes a way to sample from an exponential variable (using -log(runif)/mu) but then, as far as I can understand, seems to jump into slightly different Stata code for sampling from $\exp(x - \exp(x))$ ( I have no access to Stata and so cannot step through it unfortunately which may have helped my understanding).

These are my attempts to sample, in R, which appear successful but I am unsure. In the first method below, the Stata code from the link uses two calls to a random uniform generator and does a few extra calculations.

# set parameters
set.seed(1)
n = 1000000
mu = 10
Trying to translate the code from web
u = runif(n)
v = -log(u)/mu
hist(log(v) + log(mu), breaks=100, xlim=c(-10, 10), probability = TRUE, main="Method 1")
plot(function(x) exp(x - exp(x)), from=-10, to=10, add=TRUE, col="red", lwd=2)

I can also get samples that appear a good fit using the inbuilt random exponential sampler using the same approach as above.

# using in built functions
v = rexp(n, rate=mu) 
hist(log(v) + log(mu), breaks=100, xlim=c(-10, 10), probability = TRUE, main="Method 2")
plot(function(x) exp(x - exp(x)), from=-10, to=10, add=TRUE, col="red", lwd=2)

Are these correct? Why?

Answer 1

Just to formalise whuber's excellent answer in comments, taking $y = e^x$ gives $dy/dx = e^x$ so that:

$$\begin{align} f_X(x) &= \exp(x - e^x) \\[16pt] &= \exp(-e^x) \cdot e^x \\[12pt] &= e^{-y} \cdot \bigg| \frac{dy}{dx} \bigg| \\[6pt] &= \text{Exp}(y|1) \cdot \bigg| \frac{dy}{dx} \bigg|. \\[6pt] \end{align}$$

Consequently, using the rules for monotonic transformations of random variables, we can generate the random variable $X$ using the transformation:

$$X = \log Y \quad \quad \quad \quad \quad Y \sim \text{Exp}(1).$$

You can also use the relationship between the exponential and uniform distributions to generate $X$ from a continuous uniform random variable on the unit interval (which is the foundational unit of a standard PRNG):

$$X = \log Y \quad \quad \quad \quad \quad Y = -\log(U) \quad \quad \quad \quad \quad U \sim \text{U}(0,1).$$

The code you present appears to include a rate parameter, which is not part of the density function you initially describe. Your code appears to implicitly generate $Y$ from an exponential distribution with rate $\lambda = 10$ but then remove this latter part from the result when presented in the histogram --- this part is unecessary and can be removed entirely. (If you want to build the parameter into the density then you will need to generalise the form of your density function appropriately.) To implement random generation from this distribution in R you can use the following function:

rdist <- function(n) { log(-log(runif(n))) }

We can now confirm graphically that this method simulates from the stipulated density:

#Generate random variables
set.seed(1)
x <- rdist(10^6)
#Plot histogram with overlaid density
hist(x, breaks = 100, xlim = c(-10, 10), probability = TRUE,
     main = 'Histogram of randomly generated values')
plot(function(x) exp(x - exp(x)), from = -10, to = 10, add = TRUE, col = 'red', lwd = 2)