Error in stan documentation about Double Exponential?

Modeling
1

I think there may be an error in the Stan documentation about the Double Exponential ?

See the simulation below, as well as Wikipedia and Ding and Blitzstein (2018).

library(nimble)

n = 10000000
mu = 1
lambda = 2

alpha = rexp(n = n, rate = 1/lambda)
beta_mixture = rnorm(n = n, mean = mu, sd = alpha)
beta_marginal = rdexp(n = n, location = mu, scale = lambda)

par(mfrow = c(2,1))
plot(density(beta_marginal), xlim = c(-20,20))
plot(density(beta_mixture), xlim = c(-20,20))
1 Like
2

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:

library(nimble)

n = 10000000
mu = 1
lambda = 2

alpha = rexp(n = n, rate = lambda)
beta_mixture = rnorm(n = n, mean = mu, sd = sqrt(alpha))
beta_marginal = rdexp(n = n, location = mu, scale = 1 / lambda)

withr::with_par(new = list(mfrow = c(2,1)), {
  hist(beta_marginal, xlim = c(-20,20), breaks = 200, freq = F)
  hist(beta_mixture, xlim = c(-20,20), breaks = 200, freq = F)
})

Created on 2020-06-16 by the reprex package (v0.3.0)

1 Like
3

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 Exponential is parameterized in terms of inverse-scale (i.e. “rate”), except on the page about the Double Exponential, where it is parameterized in terms of the scale 1/\lambda.

The Normal is parameterized using the scale \sigma, except on the page about the Double Exponential, where it is parameterized in terms of the variance \sigma^2.

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.

3 Likes
4

I see exactly what you’re talking about now!

This is particularly unclear because it isn’t mentioned. I think this is the source for the docs that we could adjust and submit a PR?

2 Likes
5

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},

\text{DoubleExponential}(y|\mu,\sigma) = \frac{1}{2\sigma} \exp \left( - \, \frac{|y - \mu|}{\sigma} \right) .

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))
3 Likes
6

Looks good to me, the PR is here.

3 Likes