diff --git a/src/Distributions.jl b/src/Distributions.jl
index fcf699434..907480098 100644
--- a/src/Distributions.jl
+++ b/src/Distributions.jl
@@ -139,6 +139,7 @@ export
     MvNormalCanon,
     MvNormalKnownCov,
     MvTDist,
+    SymmetricMvLaplace,
     NegativeBinomial,
     NoncentralBeta,
     NoncentralChisq,
diff --git a/src/multivariate/mvlaplace.jl b/src/multivariate/mvlaplace.jl
new file mode 100644
index 000000000..9ba8a17c2
--- /dev/null
+++ b/src/multivariate/mvlaplace.jl
@@ -0,0 +1,122 @@
+"""
+
+The [symmetric multivariate Laplace distribution](https://en.wikipedia.org/wiki/Multivariate_Laplace_distribution#Symmetric_multivariate_Laplace_distribution)
+is a multidimensional generalization of the *Laplace distribution*, as described in 
+T. Eltoft, Taesu Kim and Te-Won Lee, 
+"On the multivariate Laplace distribution," 
+in IEEE Signal Processing Letters, vol. 13, no. 5, pp. 300-303, May 2006, doi: 10.1109/LSP.2006.870353.
+The symmetric multivariate Laplace distribution is defined by three parameters:
+- ``\\lambda``, which is a positive real number used to define an exponential distribution `Exp(λ)`
+- ``\\boldsymbol{\\Gamma}``, which is a k-by-k positive definite matrix with `det(Γ)=1` (as per the assumptions in the source paper)
+- ``\\boldsymbol{\\mu}``, which is a k-dimensional real-valued vector
+The expected valued of the symmetric multivariate Laplace distribution is simply ``\\boldsymbol{\\mu}``,
+whereas the covariance ``\\boldsymbol{\\Sigma} = \\lambda \\boldsymbol{\\Gamma}``.
+
+The symmetric multivariate Laplace distribution can be created by specifying either a ``\\boldsymbol{\\mu}`` and a ``\\boldsymbol{\\Sigma}`` 
+(analogously to a `MvNormal`) and the ``\\lambda`` and ``\\boldsymbol{\\Gamma}`` are calculated internally,
+or by specifying all three parameters, ``\\boldsymbol{\\mu}``, ``\\lambda`` and ``\\boldsymbol{\\Gamma}``.
+
+The probability density function of
+a k-dimensional symmetric multivariate Laplace distribution with parameters ``\\boldsymbol{\\mu}``,  ``\\lambda`` and ``\\boldsymbol{\\Gamma}`` is:
+```math
+f(\\mathbf{x}; \\boldsymbol{\\mu}, \\lambda, \\boldsymbol{\\Gamma}) = \\frac{1}{(2 \\pi)^{d/2}} \\frac{2}{\\lambda}
+\frac{\\mathrm{K}_{(d/2)-1}\\left(\\sqrt{\\frac{2}{\\lambda}q(\\mathbf{x})}\\right)}{\\left(\\sqrt{\\frac{2}{\\lambda}q(\\mathbf{x})}\\right)^{(d/2)-1}}
+```
+where ``\\mathrm{K}_m (y)`` is the Bessel function of the second kind of order ``m`` evaluated at ``y`` and
+```math
+q(\\mathbf{x} =  (\\mathbf{x} - \\boldsymbol{\\mu})^T \\Gamma^{-1} (\\mathbf{x} - \\boldsymbol{\\mu})
+```
+Note that the density function becomes undefined as ``x\\rightarrow\\boldsymbol{\\mu}`` because ``\\mathrm{K}_m (y)`` approaches
+infinity as ``y`` approaches zero. This feature is not unique to the formulation from Eltoft et al., but is also present in
+the (slightly different) version presented in https://en.wikipedia.org/wiki/Multivariate_Laplace_distribution.
+"""
+struct SymmetricMvLaplace{T<:Real,Cov<:AbstractPDMat,iCov<:AbstractMatrix,Mean<:AbstractVector} <: ContinuousMultivariateDistribution
+    μ::Mean
+    Γ::Cov
+    iΓ::Cov
+    λ::T
+end
+
+### generic methods for SymmetricMvLaplace
+length(d::SymmetricMvLaplace) = length(d.μ)
+params(d::SymmetricMvLaplace) = (d.μ, d.λ, d.Γ)
+
+insupport(d::SymmetricMvLaplace, x::AbstractVector) =
+    length(d) == length(x) && all(isfinite, x)
+
+minimum(d::SymmetricMvLaplace) = fill(eltype(d)(-Inf), length(d))
+maximum(d::SymmetricMvLaplace) = fill(eltype(d)(Inf), length(d))
+mode(d::SymmetricMvLaplace) = mean(d)
+modes(d::SymmetricMvLaplace) = [mean(d)]
+mean(d::SymmetricMvLaplace) = d.μ
+var(d::SymmetricMvLaplace) = diag(d.λ * d.Γ)
+cov(d::SymmetricMvLaplace) = d.λ * d.Γ
+invcov(d::SymmetricMvLaplace) = inv(d.λ * d.Γ)
+
+### Construction when all three parameters are specified
+function SymmetricMvLaplace(μ::AbstractVector{<:Real}, λ::Real, Γ::AbstractPDMat{<:Real})
+    size(Γ, 1) == length(μ) || throw(DimensionMismatch("The dimensions of mu and Gamma are inconsistent."))
+    isapprox(det(Γ), 1) || throw(ArgumentError("det(Gamma) is not approximately equal to 1."))
+    λ > 0 || throw(ArgumentError("λ is not a positive real number."))
+
+    R = Base.promote_eltype(μ, λ, Γ)
+    _μ, _λ, _Γ = convert(AbstractArray{R}, μ), convert(R, λ), convert(AbstractArray{R}, Γ)
+    _iΓ = inv(_Γ) 
+    SymmetricMvLaplace{R,typeof(_Γ), typeof(_iΓ), typeof(_μ)}(_μ, _Γ, _iΓ, _λ)
+end
+
+### Construction when only an overall covariance is specified
+function SymmetricMvLaplace(μ::AbstractVector{T}, Σ::AbstractPDMat{T}) where {T<:Real}
+    size(Σ, 1) == length(μ) || throw(DimensionMismatch("The dimensions of mu and Sigma are inconsistent."))
+    λ = det(Σ)^(1/size(Σ,1))
+    Γ = 1/λ * Σ
+    iΓ = inv(Γ)
+    SymmetricMvLaplace{T,typeof(Γ), typeof(iΓ), typeof(μ)}(μ, Γ, iΓ, λ)
+end
+
+function SymmetricMvLaplace(μ::AbstractVector{<:Real}, Σ::AbstractPDMat{<:Real})
+    R = Base.promote_eltype(μ, Σ)
+    SymmetricMvLaplace(convert(AbstractArray{R}, μ), convert(AbstractArray{R}, Σ))
+end
+
+# constructor with general gamma matrix and lambda
+SymmetricMvLaplace(μ::AbstractVector{<:Real}, λ::Real, Γ::AbstractMatrix{<:Real}) = SymmetricMvLaplace(μ, λ, PDMat(Γ))
+SymmetricMvLaplace(μ::AbstractVector{<:Real}, λ::Real, Γ::Diagonal{<:Real}) = SymmetricMvLaplace(μ, λ, PDiagMat(Γ.diag))
+SymmetricMvLaplace(μ::AbstractVector{<:Real}, λ::Real, Γ::Union{Symmetric{<:Real,<:Diagonal{<:Real}},Hermitian{<:Real,<:Diagonal{<:Real}}}) = SymmetricMvLaplace(μ, λ, PDiagMat(Γ.data.diag))
+SymmetricMvLaplace(μ::AbstractVector{<:Real}, λ::Real, Γ::UniformScaling{<:Real}) =
+    SymmetricMvLaplace(μ, λ, ScalMat(length(μ), Γ.λ))
+function SymmetricMvLaplace(
+    μ::AbstractVector{<:Real}, λ::Real, Γ::Diagonal{<:Real,<:FillArrays.AbstractFill{<:Real,1}}
+)
+    return SymmetricMvLaplace(μ, λ, ScalMat(size(Γ, 1), FillArrays.getindex_value(Γ.diag)))
+end
+
+# constructor with general covariance matrix
+SymmetricMvLaplace(μ::AbstractVector{<:Real}, Σ::AbstractMatrix{<:Real}) = SymmetricMvLaplace(μ, PDMat(Σ))
+SymmetricMvLaplace(μ::AbstractVector{<:Real}, Σ::Diagonal{<:Real}) = SymmetricMvLaplace(μ, PDiagMat(Σ.diag))
+SymmetricMvLaplace(μ::AbstractVector{<:Real}, Σ::Union{Symmetric{<:Real,<:Diagonal{<:Real}},Hermitian{<:Real,<:Diagonal{<:Real}}}) = SymmetricMvLaplace(μ, PDiagMat(Σ.data.diag))
+SymmetricMvLaplace(μ::AbstractVector{<:Real}, Σ::UniformScaling{<:Real}) =
+    SymmetricMvLaplace(μ, ScalMat(length(μ), Σ.λ))
+function SymmetricMvLaplace(
+    μ::AbstractVector{<:Real}, Σ::Diagonal{<:Real,<:FillArrays.AbstractFill{<:Real,1}}
+)
+    return SymmetricMvLaplace(μ, ScalMat(size(Σ, 1), FillArrays.getindex_value(Σ.diag)))
+end
+
+### generic methods for SymmetricMvLaplace
+Base.eltype(::Type{<:SymmetricMvLaplace{T}}) where {T} = T
+@inline partype(d::SymmetricMvLaplace{T}) where {T<:Real} = T
+
+function _rand!(rng::AbstractRNG, d::SymmetricMvLaplace, x::VecOrMat)
+    unwhiten!(d.Γ, randn!(rng, x))
+    x .*= sqrt.(rand(rng, Exponential(d.λ), 1, size(x,2)))
+    x .+=  d.μ
+    return x
+end
+
+function _logpdf(d::SymmetricMvLaplace, x::AbstractArray)
+    _d = length(d) / 2
+    xdif = x - d.μ
+    q = dot(xdif, d.iΓ, xdif)
+    return -_d * log(2π) + log(2/d.λ) + log(besselk(_d-1, sqrt(2/d.λ * q))) - 0.5*(_d-1)*log(0.5*d.λ*q)
+end
\ No newline at end of file
diff --git a/src/multivariates.jl b/src/multivariates.jl
index 962b90ed5..4693c8067 100644
--- a/src/multivariates.jl
+++ b/src/multivariates.jl
@@ -126,6 +126,7 @@ for fname in ["dirichlet.jl",
               "mvlognormal.jl",
               "mvtdist.jl",
               "product.jl", # deprecated
-              "vonmisesfisher.jl"]
+              "vonmisesfisher.jl",
+              "mvlaplace.jl"]
     include(joinpath("multivariate", fname))
 end
diff --git a/test/multivariate/mvlaplace.jl b/test/multivariate/mvlaplace.jl
new file mode 100644
index 000000000..5b928bbe9
--- /dev/null
+++ b/test/multivariate/mvlaplace.jl
@@ -0,0 +1,92 @@
+using Distributions
+import PDMats: ScalMat, PDiagMat, PDMat
+using LinearAlgebra, Test
+using FillArrays
+
+
+@testset "SymmetricMvLaplace tests" begin
+    mu = [1., 2., 3.]
+    mu_r = 1.:3.
+
+    C = [4. -2. -1.; -2. 5. -1.; -1. -1. 6.]
+    dv = [1.2, 3.4, 2.6]
+    for (g, μ, Σ) in [
+        (SymmetricMvLaplace(mu, C), mu, C),
+        (SymmetricMvLaplace(mu_r, C), mu_r, C),
+        (SymmetricMvLaplace(mu, Symmetric(C)), mu, Matrix(Symmetric(C))),
+        (SymmetricMvLaplace(mu_r, Symmetric(C)), mu_r, Matrix(Symmetric(C))),
+        (SymmetricMvLaplace(mu, Diagonal(dv)), mu, Matrix(Diagonal(dv))),
+        (SymmetricMvLaplace(mu, Symmetric(Diagonal(dv))), mu, Matrix(Diagonal(dv))),
+        (SymmetricMvLaplace(mu, Hermitian(Diagonal(dv))), mu, Matrix(Diagonal(dv))),
+        (SymmetricMvLaplace(mu_r, Diagonal(dv)), mu_r, Matrix(Diagonal(dv))) ]
+
+        @test mean(g)   ≈ μ
+        @test cov(g)    ≈ Σ
+        @test invcov(g) ≈ inv(Σ)
+    end
+end
+
+@testset "SymmetricMvLaplace constructors" begin
+    @testset "Providing mu and Sigma" begin
+        mu = [1., 2., 3.]
+        C = [4. -2. -1.; -2. 5. -1.; -1. -1. 6.]
+        @test typeof(SymmetricMvLaplace(mu, PDMat(Array{Float32}(C)))) == typeof(SymmetricMvLaplace(mu, PDMat(C)))
+        @test typeof(SymmetricMvLaplace(mu, Array{Float32}(C))) == typeof(SymmetricMvLaplace(mu, PDMat(C)))
+        @test SymmetricMvLaplace(mu, I) === SymmetricMvLaplace(mu, Diagonal(Ones(length(mu))))
+        @test SymmetricMvLaplace(mu, 9 * I) === SymmetricMvLaplace(mu, Diagonal(Fill(9, length(mu))))
+        @test SymmetricMvLaplace(mu, 0.25f0 * I) === SymmetricMvLaplace(mu, Diagonal(Fill(0.25f0, length(mu))))
+    end
+    @testset "Providing mu, lambda and Gamma" begin
+        mu = [1., 2., 3.]
+        C = [4. -2. -1.; -2. 5. -1.; -1. -1. 6.]
+        L = det(C)^(1/size(C,1))
+        G = 1/L * C
+        @test typeof(SymmetricMvLaplace(mu, L, PDMat(Array{Float32}(G)))) == typeof(SymmetricMvLaplace(mu, L, PDMat(G)))
+        @test typeof(SymmetricMvLaplace(mu, L, Array{Float32}(G))) == typeof(SymmetricMvLaplace(mu, L, PDMat(G)))
+        @test SymmetricMvLaplace(mu, 1, I) === SymmetricMvLaplace(mu, Diagonal(Ones(length(mu))))
+
+        L9 = det(9*I(3))^(1/3)
+        L025 = det(0.25f0*I(3))^(1/3)
+        @test SymmetricMvLaplace(mu, L9, I) === SymmetricMvLaplace(mu, L9, Diagonal(Fill(1, length(mu))))
+        @test SymmetricMvLaplace(mu, L025, I) === SymmetricMvLaplace(mu, L025, Diagonal(Fill(1, length(mu))))
+        @test SymmetricMvLaplace(mu, L9, Symmetric(Float64.(I(3)))) == SymmetricMvLaplace(mu, L9, Float64.(I(3)))
+    end
+end
+
+@testset "Symmetric MvLaplace basic functions" begin
+    d_dim = 2
+    for t in (Float32, Float64)
+        @testset "Type: $t" begin
+            mu = t.(zeros(d_dim))
+            λ = t.(1)
+            Γ = t.(I(d_dim))
+            d = SymmetricMvLaplace(mu, λ, Γ)
+    
+            @test length(d) == d_dim
+            @test params(d) == (mu, λ, Γ)
+            @test maximum(d) == repeat([Inf], d_dim)
+            @test minimum(d) == repeat([-Inf], d_dim)
+            @test mode(d) == mu
+            @test modes(d) == [mu]
+            @test var(d) == diag(λ*Γ)
+            @test eltype(d) == t
+            @test partype(d) == t
+            @test insupport(d, rand(2))
+
+            s = rand(d, 10)
+            @test size(s) == (2, 10)
+            @test eltype(s) == t
+            @test isfinite(loglikelihood(d, s))
+
+            l = Laplace(0, 0.707)
+            s_d = rand(d, 10_000_000)
+            qs = collect(0.01:0.01:0.99)
+            q_d = quantile.(Ref(s_d[1,:]), qs)
+            q_l = quantile.(Ref(l), qs)
+            @test all(isapprox.(q_d, q_l; atol=1e-2)) ### tolerance a bit high as tails are heavy and a bit stochastic
+            ### below rhs calculated by pen and paper
+            @test isapprox(logpdf(d, [1; 1]), log(1/π * besselk(0, 2))) 
+            @test isapprox(logpdf(d, [2; 2]), log(1/π * besselk(0, 4)))
+        end
+    end
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index f63a3d5e8..b52e31c6d 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -46,6 +46,7 @@ const tests = [
     "multivariate/dirichletmultinomial",
     "univariate/continuous/logitnormal",
     "multivariate/mvtdist",
+    "multivariate/mvlaplace",
     "univariate/continuous/kolmogorov",
     "edgeworth",
     "matrixreshaped", # extra file compared to /src