What is the PDF of a Normal convolved with a Laplace 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(σ)))
    )
end

 A: 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}$$


Comments
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.

Illustration
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.

Code
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")

