I'd like to see if using Stan or similar I can successfully model Laplace noise added to data through the use of a convolved Normal-Laplace distribution and MCMC sampling. In the literature I can only find this https://www.tandfonline.com/doi/full/10.1080/03610926.2015.1040510 which concerns an asymmetric Laplace distribution. I was wondering if anyone knows of a simpler expression I could use for a Normal convolved with a symmetric Laplace. And further to this, whether this is even the correct approach to modelling noise added to data (provided I know the scale of the Laplace noise added). Here is the corrected PDF in the linked article which is unfortunately behind a paywall:

$$ \begin{aligned} f(y)=& \frac{\alpha \beta}{2(\alpha+\beta)}\left[e^{\frac{1}{2} \alpha\left(-2 y+2 \mu+\alpha \sigma^{2}\right)} \operatorname{erf} c\left(\frac{\alpha \sigma}{\sqrt{2}}-\frac{y-\mu}{\sqrt{2} \sigma}\right)\right.\\ &\left.+e^{\frac{1}{2} \beta\left(2 y-2 \mu+\beta \sigma^{2}\right)} \operatorname{erf} c\left(\frac{\beta \sigma}{\sqrt{2}}+\frac{y-\mu}{\sqrt{2} \sigma}\right)\right] \end{aligned} $$

And for further reference, the original:

$$ \begin{array}{c} \phi(x)=\frac{e^{-\frac{1}{2} x^{2}}}{\sqrt{2 \pi}} \\ \left.\Phi(x)=\frac{1}{\sqrt{2 \pi}} \int_{-\infty}^{x} e^{-\frac{t^{2}}{2}} d t=\frac{1}{2}\left[1+\operatorname{erf}\left(\frac{x}{\sqrt{2}}\right)\right]\right] \\ R(z)=\frac{1-\Phi(z)}{\phi(z)}=\frac{\sqrt{2 \pi}\left[\operatorname{erfc}\left(\frac{z}{\sqrt{2}}\right)-\frac{1}{2}\right]}{e^{-\frac{z^{2}}{2}}} \\ f(y)=\frac{\alpha \beta}{\alpha+\beta} \phi\left(\frac{y-\mu}{\sigma}\right)\left[R\left(\alpha \sigma-\frac{y-\mu}{\sigma}\right)+R\left(\beta \sigma+\frac{y-\mu}{\sigma}\right)\right] \end{array}$$

Here is the PDF I get by setting $\alpha = \beta$ in the referenced article. I am presuming that the asymmetric Laplace here is centred, meaning this would only work for a $Laplace(0, \lambda)$. This results in extremely large PDF values in some cases, so I believe I have done something wrong.

$$\frac{λ}{4} \text{erfc} \left(\frac{λσ}{\sqrt{2}} - \frac{y - μ}{\sigma \sqrt{2}} \right) \left( \exp\left(\frac{\lambda}{2} \left(-2y + 2μ + λσ^2\right)\right) + \exp\left(\frac{\lambda}{2} \left(2y - 2μ + λσ^2\right)\right)\right)$$

function pdf_NL(μ, σ, λ, y)
    return (λ / 4) * erfc((λ * σ) / √2 - (y - μ) / (√2 * σ)) * (
        exp(0.5λ * (-2y + 2μ + λ * abs2(σ))) +
        exp(0.5λ * (2y - 2μ + λ * abs2(σ)))
  • 2
    $\begingroup$ Your reference explicitly answers the question in the title by providing expressions for the pdf and cdf of this convolution. Just set $\alpha=\beta.$ It's clear nothing will cancel or disappear from the formulas, so little simplification will occur. $\endgroup$ – whuber May 19 '20 at 11:45
  • $\begingroup$ @whuber Yes, this was my initial attempt too, the PDF that results through setting $\alpha = \beta = \lambda$ doesn't seem right when I pass it values though. I will edit the question to include it and some Julia code to go alongside. I am presuming here that the asymmetric Laplace is always centred, but perhaps this is where I am misunderstanding? $\endgroup$ – Harrison W. May 19 '20 at 12:15
  • $\begingroup$ Could you please rephrase your question by (a) avoiding the connection with DP which is irrelevant for the question and (b) formally and rigorously defining the convolution you are interested in? Note that the density of the convolution is NOT the sum of the two densities involved. Further, convolutions are well-suited to MCMC algorithms as they naturally invoke a latent variable. Enjoying a closed form density is thus less important. $\endgroup$ – Xi'an May 19 '20 at 12:35
  • $\begingroup$ The link is to a paper behind a paywall. $\endgroup$ – Xi'an May 19 '20 at 12:36
  • $\begingroup$ Your formula does not seem to agree with that of the paper: you are neglecting terms in the denominator. See the formula for $R.$ $\endgroup$ – whuber May 19 '20 at 12:44

Let's work it out from first principles, beginning with the hard work of computing a convolution.

As an auxiliary calculation, consider the distribution of $W=X+Y$ where $Y$ has an Exponential distribution with pdf $$f_Y(y) = e^{-y}\,\mathcal{I}(y\gt 0)$$ and $X$ has a Normal$(\mu,\sigma^2)$ distribution with pdf $f_X(x;\mu,\sigma) = \phi((x-\mu)/\sigma)/\sigma$ where $$\phi(z) = \frac{1}{\sqrt{2\pi}}\,e^{-z^2/2}$$ is the standard Normal pdf. The PDF of the sum is the convolution

$$f_W(w;\mu,\sigma) = \int_{-\infty}^\infty f_Y(y) f_X(w-y;\mu,\sigma)\,\mathrm{d}y = \int_0^\infty e^{-y} f_X(w-y;\mu,\sigma)\,\mathrm{d}y.$$

Substituting $\sigma z = w - y - \mu$ expresses this integral as

$$\eqalign{f_W(w;\mu,\sigma) &= e^{\mu-w}\,e^{\sigma^2/2}\int_{-\infty}^{(w-\mu)/\sigma} \phi(z-\sigma)\,\mathrm{d}z \\ &= e^{\mu-w+\sigma^2/2}\, \Phi\left(\frac{w-\mu}{\sigma}-\sigma\right)}\tag{1}$$

where $\Phi$ is the standard normal CDF,

$$\Phi(z) = \int_{-\infty}^z \phi(z)\,\mathrm{d}z.$$

The rest builds on this work and is relatively easy.

An asymmetric Laplace random variable $U$ is based on a mixture of a scaled exponential distribution and the negative of a scaled exponential distribution (potentially with a different scale, thereby making the mixture asymmetric). This mixture is then shifted by a specified amount. The amount of mixing is established to give the Laplace pdf a unique value at its peak--but this is unimportant.

One component of $U$ therefore can be expressed as $$U_+ = \alpha Y + \lambda$$ with a positive scale $\alpha$ and the other component as $$U_- = -\beta Y + \lambda$$ with a positive scale $\beta.$ (I apologize: I worked this out before realizing that my $\alpha$ is $1/\alpha$ in the paper and my $\beta$ is $1/\beta$ in the paper: in the end, after setting $\alpha=\beta,$ this won't matter.)

When we add $X = \sigma Z + \mu$ we obtain two components, of which the first is $$W_+ = U_+ + X = \alpha Y + \lambda + \sigma Z + \mu = \alpha\left(Y + \left[\frac{\sigma}{\alpha} Z + \frac{\lambda + \mu}{\alpha}\right]\right)$$ and the second is similarly written. To obtain its pdf, all we need to do is scale formula $(1)$ by $\alpha,$ giving

$$f_{W_+}(w;\mu,\sigma,\lambda,\alpha) = \frac{1}{\alpha}\,f_W\left(\frac{w}{\alpha};\frac{\lambda+\mu}{\alpha}, \frac{\sigma}{\alpha}\right).\tag{2}$$

Likewise, because

$$W_- = U_- + X = -\beta Y + \lambda + \sigma Z + \mu = -\beta \left(Y + \left[-\frac{\sigma}{\beta } Z + \frac{\lambda + \mu}{\beta }\right]\right)$$

and $-Z$ has the same distribution as $Z$, formula $(1)$ yields

$$\eqalign{f_{W_-}(w;\mu,\sigma,\lambda,\beta) &= \frac{1}{\beta }\,f_W\left(-\frac{w}{\beta };-\frac{\lambda+\mu}{\beta }, \frac{\sigma}{\beta }\right) \\ &= f_{W_+}(-w;-\lambda,\beta,-\mu,\sigma).}\tag{3}$$

The mixture pdf is

$$f_W(w;\mu,\sigma,\lambda,\alpha,\beta,p) = pf_{W_+}(w;\mu,\sigma,\lambda,\alpha) + (1-p) f_{W_-}(w;\mu,\sigma,\lambda,\beta).\tag{4}$$


For the Laplace-Normal distribution, use $p = \alpha / (\alpha + \beta).$ In your case $\alpha=\beta,$ which evidently "simplifies" $(4)$ a tiny bit--but a quick look at its component formulas $(2)$ and $(3)$ suggests there's not much one can do algebraically to reduce the amount of computation, so why bother?

Each of the components of the final formula $(4),$ as embodied in formulas $(1),$ $(2),$ and $(3)$ can be separately and flexibly implemented and separately tested. This makes for an easier and more reliable software implementation than attempting to combine them all into one monster combination of $\phi$ and $\Phi,$ as done in the referenced paper. As a bonus, important numerical improvements in the calculation can be implemented exactly where they are needed, making the code relatively easy to maintain. As an example, see how $f_W$ is implemented using logarithms (as f.1) in the code below.


This plot compares a histogram of one million iid draws from an asymmetric Laplace-Normal distribution with pdf $f_W(w;4,0.5,-3,2,1,2/3)$ to a calculation based directly on formulas $(1) - (4):$


The agreement is a pretty good test.


Here's the R code that generated this simulation and this plot.

n <- 1e6      # Size of simulation
mu <- 4
sigma <- 1/2
alpha <- 2
lambda <- -3
beta <- 1
# Generate data.
# set.seed(17)
X <- rnorm(n, mu, sigma)
Y <- ifelse(runif(n, 0, alpha + beta) < alpha, alpha, -beta) * rexp(n) + lambda
W <- X + Y
# Plot their histogram.
hist(W, freq=FALSE, breaks=200, cex.main=1)
# Overplot the PDF.
f.1 <- function(w, mu=0, sigma=1) {
  exp(mu - w + sigma^2/2 + pnorm((w - mu)/sigma - sigma, log=TRUE))
f.plus <- function(w, mu=0, sigma=1, lambda=0, alpha=1) {
  f.1(w / alpha, (lambda + mu) / alpha, sigma / alpha) / alpha
f.minus <- function(w, mu=0, sigma=1, lambda=0, beta=1) {
  f.plus(-w, -mu, sigma, -lambda, beta)
f <- function(w, mu=0, sigma=1, lambda=0, alpha=1, beta=1, p=1/2) {
  p * f.plus(w, mu, sigma, lambda, alpha) + (1-p) * f.minus(w, mu, sigma, lambda, beta)
f.asymmetric <- function(y, mu=0, sigma=1, lambda=0, alpha=1, beta=1) {
  f(y, mu, sigma, lambda, alpha, beta, alpha / (alpha + beta))
curve(f.asymmetric(x, mu, sigma, lambda, alpha, beta), add=TRUE, lwd=2, col="Red")

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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