2
$\begingroup$

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?

enter image description here

$\endgroup$
5
  • 2
    $\begingroup$ Write $y=e^x$ so that $\exp\left(x-e^x\right)\mathrm dx = e^{-y}\mathrm dy.$ Draw a random variate $Y$ from the latter distribution and set $X=\log Y.$ See stats.stackexchange.com/questions/591943 for details of the built-in R implementation and the code at stats.stackexchange.com/a/492176/919 for a "hand-rolled" solution (in one short line). If you're unsure about your code, just do a simple goodness of fit test. A chi-squared test will work fine and a probability plot will be simple and helpful. $\endgroup$
    – whuber
    May 10, 2023 at 19:54
  • $\begingroup$ Thanks @whuber. Do I understand that your first sentence is on the algebra for using a change of variables to show that the log of the standard exponential ($exp(1)$) is distributed $exp(𝑥−𝑒^𝑥)$. ( I also think i see now why I was adding on log(mu), as this same value was subtracted in the log(v) step) $\endgroup$ May 10, 2023 at 20:54
  • $\begingroup$ Yes. It looks like your solution is the same as the one I posted in the second link above, with the calculations broken between the preliminary drawing of u and its recentering within the call to hist. $\endgroup$
    – whuber
    May 10, 2023 at 20:56
  • 1
    $\begingroup$ Thanks again @whuber . The links were useful in getting me over the line. $\endgroup$ May 10, 2023 at 21:01
  • 1
    $\begingroup$ I just verified that ks.test will work here. It reverts to an asymptotic calculation of the p-value for large datasets. (You really don't need such large datasets for this purpose unless you are concerned about possible errors far out in their tails.) This is the ideal application for the KS test because you have definite hypotheses about both distributions! $\endgroup$
    – whuber
    May 10, 2023 at 21:33

1 Answer 1

4
$\begingroup$

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)

enter image description here

$\endgroup$
1
  • $\begingroup$ Thanks Ben, I appreciate you clarifying and writing this up. (re the adding in and removing of the rate; this was due to the link sampling conditional on rate / counts, but this was of course not what i asked here) $\endgroup$ May 11, 2023 at 13:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.