The doc is definitely unclear; noting that \alpha is the variance, and that \lambda is the rate (or inverse-scale), and hence \sigma = \frac{1}{\lambda}. Some small tweaks give me identical answers:

Thanks, @hhau ! Yes, I came to the same tweaks. And I agree that it is strange that the documentation changes how the Exponential, Normal, and DoubleExponential are parameterized.

The Double Exponential is parameterized in terms of scale on the top of the page, but on the bottom of the page it is parameterized in terms of the inverse-scale 1/\lambda.

Great idea, let’s do it !! That sure looks like the source to me. I think it should be (could cite Wikipedia or Ding and Blitzstein (2018) in the documentation or PR):

If \mu \in \mathbb{R} and \sigma \in \mathbb{R}^+, then for y \in
\mathbb{R},

Note that the double exponential distribution is parameterized in terms of the scale, in contrast to the exponential distribution (see section exponential distribution), which is parameterized in terms of inverse scale.

The double-exponential distribution can be defined as a compound exponential-normal distribution. Specifically, if \alpha \sim \mathsf{Exponential}\left( \frac{1}{2\sigma^2} \right)
and \beta \ |\ \alpha \sim \mathsf{Normal}(\mu, \sqrt{\alpha}),
then \beta \sim \mathsf{DoubleExponential}(\mu, \sigma).

This may be used to code a non-centered parameterization by taking \beta^{\text{raw}} \sim
\mathsf{Normal}(0, 1)
and defining \beta = \mu + \alpha \, \beta^{\text{raw}}.

With tweaked code:

n = 10000000
mu = 1
sigma = 2
alpha = rexp(n = n, rate = 1/(2*sigma^2))
beta_mixture = rnorm(n = n, mean = mu, sd = sqrt(alpha))
beta_marginal = rdexp(n = n, location = mu, scale = sigma)
par(mfrow = c(2,1))
plot(density(beta_marginal), xlim = c(-20,20))
plot(density(beta_mixture), xlim = c(-20,20))